Loading observation_sim/instruments/chip/effects.py +53 −55 Original line number Original line Diff line number Diff line Loading @@ -559,22 +559,22 @@ def CTEModelColRow(img, trail_direction='up', direction='column', threshold=27): # ---------- For Cosmic-Ray Simulation ------------ # ---------- For Cosmic-Ray Simulation ------------ # ---------- Zhang Xin ---------------------------- # ---------- Zhang Xin ---------------------------- def getYValue(collection, x): def getYValue(collection, x): index = 0; index = 0 if (collection.shape[1] == 2): if (collection.shape[1] == 2): while(x > collection[index, 0] and index < collection.shape[0]): while(x > collection[index, 0] and index < collection.shape[0]): index= index + 1; index = index + 1 if (index == collection.shape[0] or index == 0): if (index == collection.shape[0] or index == 0): return 0; return 0 deltX = collection[index, 0] - collection[index-1, 0]; deltX = collection[index, 0] - collection[index-1, 0] deltY = collection[index, 1] - collection[index-1, 1]; deltY = collection[index, 1] - collection[index-1, 1] if deltX == 0: if deltX == 0: return (collection[index, 1] + collection[index-1, 1])/2.0 return (collection[index, 1] + collection[index-1, 1])/2.0 else: else: a = deltY/deltX; a = deltY/deltX return a * (x - collection[index-1, 0]) + collection[index-1, 1]; return a * (x - collection[index-1, 0]) + collection[index-1, 1] return 0; return 0 def selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size): def selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size): Loading @@ -582,43 +582,43 @@ def selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_s normalRay = 0.90 normalRay = 0.90 nnormalRay = 1-normalRay nnormalRay = 1-normalRay max_nrayLen = 100 max_nrayLen = 100 pixelNum = int(xLen * yLen * cr_pixelRatio * normalRay ); pixelNum = int(xLen * yLen * cr_pixelRatio * normalRay) pixelNum_n = int(xLen * yLen * cr_pixelRatio * nnormalRay) pixelNum_n = int(xLen * yLen * cr_pixelRatio * nnormalRay) CRPixelNum = 0; CRPixelNum = 0 maxValue = max(attachedSizes[:, 1]) maxValue = max(attachedSizes[:, 1]) maxValue += 0.1; maxValue += 0.1 cr_event_num = 0; cr_event_num = 0 CRs = np.zeros(pixelNum); CRs = np.zeros(pixelNum) while (CRPixelNum < pixelNum): while (CRPixelNum < pixelNum): x = CR_max_size * np.random.random(); x = CR_max_size * np.random.random() y = maxValue * np.random.random(); y = maxValue * np.random.random() if (y <= getYValue(attachedSizes, x)): if (y <= getYValue(attachedSizes, x)): CRs[cr_event_num] = np.ceil(x); CRs[cr_event_num] = np.ceil(x) cr_event_num = cr_event_num + 1; cr_event_num = cr_event_num + 1 CRPixelNum = CRPixelNum + round(x); CRPixelNum = CRPixelNum + round(x) while (CRPixelNum < pixelNum + pixelNum_n): while (CRPixelNum < pixelNum + pixelNum_n): nx = np.random.random()*(max_nrayLen-CR_max_size)+CR_max_size nx = np.random.random()*(max_nrayLen-CR_max_size)+CR_max_size CRs[cr_event_num] = np.ceil(nx); CRs[cr_event_num] = np.ceil(nx) cr_event_num = cr_event_num + 1; cr_event_num = cr_event_num + 1 CRPixelNum = CRPixelNum + np.ceil(nx); CRPixelNum = CRPixelNum + np.ceil(nx) return CRs[0:cr_event_num]; return CRs[0:cr_event_num] def defineEnergyForCR(cr_event_size, seed=12345): def defineEnergyForCR(cr_event_size, seed=12345): import random import random sigma = 0.6 / 2.355; sigma = 0.6 / 2.355 mean = 3.3; mean = 3.3 random.seed(seed) random.seed(seed) energys = np.zeros(cr_event_size); energys = np.zeros(cr_event_size) for i in np.arange(cr_event_size): for i in np.arange(cr_event_size): energy_index = random.normalvariate(mean, sigma); energy_index = random.normalvariate(mean, sigma); energys[i] = pow(10, energy_index); energys[i] = pow(10, energy_index) return energys; return energys def convCR(CRmap=None, addPSF=None, sp_n=4): def convCR(CRmap=None, addPSF=None, sp_n=4): sh = CRmap.shape sh = CRmap.shape Loading Loading @@ -668,13 +668,13 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 # Return: an 2-D numpy array # Return: an 2-D numpy array # attachedSizes = np.loadtxt('./wfc-cr-attachpixel.dat'); # attachedSizes = np.loadtxt('./wfc-cr-attachpixel.dat'); np.random.seed(seed) np.random.seed(seed) CR_max_size = 20.0; CR_max_size = 20.0 cr_size = selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size); cr_size = selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size) cr_event_size = cr_size.shape[0]; cr_event_size = cr_size.shape[0] cr_energys = defineEnergyForCR(cr_event_size,seed = seed); cr_energys = defineEnergyForCR(cr_event_size, seed=seed) CRmap = np.zeros([yLen, xLen]); CRmap = np.zeros([yLen, xLen]) # produce conv kernel # produce conv kernel from astropy.modeling.models import Gaussian2D from astropy.modeling.models import Gaussian2D Loading @@ -690,27 +690,25 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 convKernel = addPSF/addPSF.sum() convKernel = addPSF/addPSF.sum() # --------------------------------- # --------------------------------- for i in np.arange(cr_event_size): for i in np.arange(cr_event_size): xPos = round((xLen - 1)* np.random.random()); xPos = round((xLen - 1)* np.random.random()); yPos = round((yLen - 1)* np.random.random()); yPos = round((yLen - 1)* np.random.random()); cr_lens = int(cr_size[i]); cr_lens = int(cr_size[i]); if cr_lens == 0: if cr_lens == 0: continue; continue pix_energy = cr_energys[i]/gain/cr_lens; pix_energy = cr_energys[i]/gain/cr_lens pos_angle = 1/2*math.pi*np.random.random(); pos_angle = 1/2*math.pi*np.random.random() crMatrix = np.zeros([cr_lens+1, cr_lens + 1]) crMatrix = np.zeros([cr_lens+1, cr_lens + 1]) for j in np.arange(cr_lens): for j in np.arange(cr_lens): x_n = int(np.cos(pos_angle)*j - np.sin(pos_angle)*0); x_n = int(np.cos(pos_angle)*j - np.sin(pos_angle)*0) if x_n < 0: if x_n < 0: x_n = x_n + cr_lens+1 x_n = x_n + cr_lens+1 y_n = int(np.sin(pos_angle)*j + np.cos(pos_angle)*0); y_n = int(np.sin(pos_angle)*j + np.cos(pos_angle)*0) if x_n < 0 or x_n > cr_lens or y_n < 0 or y_n > cr_lens: if x_n < 0 or x_n > cr_lens or y_n < 0 or y_n > cr_lens: continue; continue crMatrix[y_n, x_n] = pix_energy; crMatrix[y_n, x_n] = pix_energy crMatrix_n = convCR(crMatrix, convKernel, sp_n) crMatrix_n = convCR(crMatrix, convKernel, sp_n) # crMatrix_n = crMatrix # crMatrix_n = crMatrix Loading Loading
observation_sim/instruments/chip/effects.py +53 −55 Original line number Original line Diff line number Diff line Loading @@ -559,22 +559,22 @@ def CTEModelColRow(img, trail_direction='up', direction='column', threshold=27): # ---------- For Cosmic-Ray Simulation ------------ # ---------- For Cosmic-Ray Simulation ------------ # ---------- Zhang Xin ---------------------------- # ---------- Zhang Xin ---------------------------- def getYValue(collection, x): def getYValue(collection, x): index = 0; index = 0 if (collection.shape[1] == 2): if (collection.shape[1] == 2): while(x > collection[index, 0] and index < collection.shape[0]): while(x > collection[index, 0] and index < collection.shape[0]): index= index + 1; index = index + 1 if (index == collection.shape[0] or index == 0): if (index == collection.shape[0] or index == 0): return 0; return 0 deltX = collection[index, 0] - collection[index-1, 0]; deltX = collection[index, 0] - collection[index-1, 0] deltY = collection[index, 1] - collection[index-1, 1]; deltY = collection[index, 1] - collection[index-1, 1] if deltX == 0: if deltX == 0: return (collection[index, 1] + collection[index-1, 1])/2.0 return (collection[index, 1] + collection[index-1, 1])/2.0 else: else: a = deltY/deltX; a = deltY/deltX return a * (x - collection[index-1, 0]) + collection[index-1, 1]; return a * (x - collection[index-1, 0]) + collection[index-1, 1] return 0; return 0 def selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size): def selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size): Loading @@ -582,43 +582,43 @@ def selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_s normalRay = 0.90 normalRay = 0.90 nnormalRay = 1-normalRay nnormalRay = 1-normalRay max_nrayLen = 100 max_nrayLen = 100 pixelNum = int(xLen * yLen * cr_pixelRatio * normalRay ); pixelNum = int(xLen * yLen * cr_pixelRatio * normalRay) pixelNum_n = int(xLen * yLen * cr_pixelRatio * nnormalRay) pixelNum_n = int(xLen * yLen * cr_pixelRatio * nnormalRay) CRPixelNum = 0; CRPixelNum = 0 maxValue = max(attachedSizes[:, 1]) maxValue = max(attachedSizes[:, 1]) maxValue += 0.1; maxValue += 0.1 cr_event_num = 0; cr_event_num = 0 CRs = np.zeros(pixelNum); CRs = np.zeros(pixelNum) while (CRPixelNum < pixelNum): while (CRPixelNum < pixelNum): x = CR_max_size * np.random.random(); x = CR_max_size * np.random.random() y = maxValue * np.random.random(); y = maxValue * np.random.random() if (y <= getYValue(attachedSizes, x)): if (y <= getYValue(attachedSizes, x)): CRs[cr_event_num] = np.ceil(x); CRs[cr_event_num] = np.ceil(x) cr_event_num = cr_event_num + 1; cr_event_num = cr_event_num + 1 CRPixelNum = CRPixelNum + round(x); CRPixelNum = CRPixelNum + round(x) while (CRPixelNum < pixelNum + pixelNum_n): while (CRPixelNum < pixelNum + pixelNum_n): nx = np.random.random()*(max_nrayLen-CR_max_size)+CR_max_size nx = np.random.random()*(max_nrayLen-CR_max_size)+CR_max_size CRs[cr_event_num] = np.ceil(nx); CRs[cr_event_num] = np.ceil(nx) cr_event_num = cr_event_num + 1; cr_event_num = cr_event_num + 1 CRPixelNum = CRPixelNum + np.ceil(nx); CRPixelNum = CRPixelNum + np.ceil(nx) return CRs[0:cr_event_num]; return CRs[0:cr_event_num] def defineEnergyForCR(cr_event_size, seed=12345): def defineEnergyForCR(cr_event_size, seed=12345): import random import random sigma = 0.6 / 2.355; sigma = 0.6 / 2.355 mean = 3.3; mean = 3.3 random.seed(seed) random.seed(seed) energys = np.zeros(cr_event_size); energys = np.zeros(cr_event_size) for i in np.arange(cr_event_size): for i in np.arange(cr_event_size): energy_index = random.normalvariate(mean, sigma); energy_index = random.normalvariate(mean, sigma); energys[i] = pow(10, energy_index); energys[i] = pow(10, energy_index) return energys; return energys def convCR(CRmap=None, addPSF=None, sp_n=4): def convCR(CRmap=None, addPSF=None, sp_n=4): sh = CRmap.shape sh = CRmap.shape Loading Loading @@ -668,13 +668,13 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 # Return: an 2-D numpy array # Return: an 2-D numpy array # attachedSizes = np.loadtxt('./wfc-cr-attachpixel.dat'); # attachedSizes = np.loadtxt('./wfc-cr-attachpixel.dat'); np.random.seed(seed) np.random.seed(seed) CR_max_size = 20.0; CR_max_size = 20.0 cr_size = selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size); cr_size = selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size) cr_event_size = cr_size.shape[0]; cr_event_size = cr_size.shape[0] cr_energys = defineEnergyForCR(cr_event_size,seed = seed); cr_energys = defineEnergyForCR(cr_event_size, seed=seed) CRmap = np.zeros([yLen, xLen]); CRmap = np.zeros([yLen, xLen]) # produce conv kernel # produce conv kernel from astropy.modeling.models import Gaussian2D from astropy.modeling.models import Gaussian2D Loading @@ -690,27 +690,25 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 convKernel = addPSF/addPSF.sum() convKernel = addPSF/addPSF.sum() # --------------------------------- # --------------------------------- for i in np.arange(cr_event_size): for i in np.arange(cr_event_size): xPos = round((xLen - 1)* np.random.random()); xPos = round((xLen - 1)* np.random.random()); yPos = round((yLen - 1)* np.random.random()); yPos = round((yLen - 1)* np.random.random()); cr_lens = int(cr_size[i]); cr_lens = int(cr_size[i]); if cr_lens == 0: if cr_lens == 0: continue; continue pix_energy = cr_energys[i]/gain/cr_lens; pix_energy = cr_energys[i]/gain/cr_lens pos_angle = 1/2*math.pi*np.random.random(); pos_angle = 1/2*math.pi*np.random.random() crMatrix = np.zeros([cr_lens+1, cr_lens + 1]) crMatrix = np.zeros([cr_lens+1, cr_lens + 1]) for j in np.arange(cr_lens): for j in np.arange(cr_lens): x_n = int(np.cos(pos_angle)*j - np.sin(pos_angle)*0); x_n = int(np.cos(pos_angle)*j - np.sin(pos_angle)*0) if x_n < 0: if x_n < 0: x_n = x_n + cr_lens+1 x_n = x_n + cr_lens+1 y_n = int(np.sin(pos_angle)*j + np.cos(pos_angle)*0); y_n = int(np.sin(pos_angle)*j + np.cos(pos_angle)*0) if x_n < 0 or x_n > cr_lens or y_n < 0 or y_n > cr_lens: if x_n < 0 or x_n > cr_lens or y_n < 0 or y_n > cr_lens: continue; continue crMatrix[y_n, x_n] = pix_energy; crMatrix[y_n, x_n] = pix_energy crMatrix_n = convCR(crMatrix, convKernel, sp_n) crMatrix_n = convCR(crMatrix, convKernel, sp_n) # crMatrix_n = crMatrix # crMatrix_n = crMatrix Loading