import galsim from matplotlib.pyplot import flag import numpy as np from numpy.core.fromnumeric import mean, size from numpy.random import Generator, PCG64 import math from numba import jit from astropy import stats def AddOverscan(GSImage, overscan=1000, gain=1, widthl=27, widthr=27, widtht=8, widthb=8, read_noise=5): """ Add overscan/gain; gain=e-/ADU widthl: left pre-scan width widthr: right pre-scan width widtht: top over-scan width (the top of nd-array with small row-index) widthb: bottom over-scan width (the bottom of nd-array with large row-index) """ imgshape = GSImage.array.shape newimg = galsim.Image(imgshape[1]+widthl+widthr, imgshape[0]+widtht+widthb, init_value=0) rng = galsim.UniformDeviate() NoiseOS = galsim.GaussianNoise(rng, sigma=read_noise) newimg.addNoise(NoiseOS) newimg = (newimg+overscan)/gain newimg.array[widtht:(widtht+imgshape[0]),widthl:(widthl+imgshape[1])] = GSImage.array newimg.wcs = GSImage.wcs # if GSImage.wcs is not None: # newwcs = GSImage.wcs.withOrigin(galsim.PositionD(widthl,widtht)) # newimg.wcs = newwcs # else: # pass return newimg def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=20210304, biaslevel=0): # Also called bad pixels, including hot pixels and dead pixels # Hot Pixel > 20e-/s # Dead Pixel < 70%*Mean rgf = Generator(PCG64(int(seed*1.1))) if IfHotPix==True and IfDeadPix==True: HotFraction = rgf.random() # fraction in total bad pixels elif IfHotPix==False and IfDeadPix==False: return GSImage elif IfHotPix==True: HotFraction = 1 else: HotFraction = 0 NPix = GSImage.array.size NPixBad = int(NPix*fraction) NPixHot = int(NPix*fraction*HotFraction) NPixDead = NPixBad-NPixHot NPix_y,NPix_x = GSImage.array.shape mean = np.mean(GSImage.array) rgp = Generator(PCG64(int(seed))) yxposfrac = rgp.random((NPixBad,2)) YPositHot = np.array(NPix_y*yxposfrac[0:NPixHot,0]).astype(np.int32) XPositHot = np.array(NPix_x*yxposfrac[0:NPixHot,1]).astype(np.int32) YPositDead = np.array(NPix_y*yxposfrac[NPixHot:,0]).astype(np.int32) XPositDead = np.array(NPix_x*yxposfrac[NPixHot:,1]).astype(np.int32) rgh = Generator(PCG64(int(seed*1.2))) rgd = Generator(PCG64(int(seed*1.3))) if IfHotPix==True: GSImage.array[YPositHot,XPositHot] += rgh.gamma(2,25*150,size=NPixHot) if IfDeadPix==True: GSImage.array[YPositDead,XPositDead] = rgd.random(NPixDead)*(mean-biaslevel)*0.7+biaslevel+rgp.standard_normal()*5 return GSImage def BadColumns(GSImage, seed=20240309, chipid=1): # Set bad column values ysize,xsize = GSImage.array.shape subarr = GSImage.array[int(ysize*0.1):int(ysize*0.12), int(xsize*0.1):int(xsize*0.12)] subarr = stats.sigma_clip(subarr, sigma=4, cenfunc='median', maxiters=3, masked=False) meanimg = np.median(subarr) stdimg = np.std(subarr) seed += chipid rgn = Generator(PCG64(int(seed))) rgcollen = Generator(PCG64(int(seed*1.1))) rgxpos = Generator(PCG64(int(seed*1.2))) rgdn = Generator(PCG64(int(seed*1.3))) nbadsecA,nbadsecD = rgn.integers(low=1, high=5, size=2) collen = rgcollen.integers(low=int(ysize*0.1), high=int(ysize*0.7), size=(nbadsecA+nbadsecD)) xposit = rgxpos.integers(low=int(xsize*0.05), high=int(xsize*0.95), size=(nbadsecA+nbadsecD)) print(xposit+1) # signs = 2*rgdn.integers(0,2,size=(nbadsecA+nbadsecD))-1 if meanimg>0: dn = rgdn.integers(low=meanimg*1.3+50, high=meanimg*2+150, size=(nbadsecA+nbadsecD)) #*signs elif meanimg<0: dn = rgdn.integers(low=meanimg*2-150, high=meanimg*1.3-50, size=(nbadsecA+nbadsecD)) #*signs for badcoli in range(nbadsecA): GSImage.array[(ysize-collen[badcoli]):ysize,xposit[badcoli]:(xposit[badcoli]+1)] = (np.abs(np.random.normal(0, stdimg*2, (collen[badcoli],1)))+dn[badcoli]) for badcoli in range(nbadsecD): GSImage.array[0:collen[badcoli+nbadsecA],xposit[badcoli+nbadsecA]:(xposit[badcoli+nbadsecA]+1)] = (np.abs(np.random.normal(0, stdimg*2, (collen[badcoli+nbadsecA],1)))+dn[badcoli+nbadsecA]) return GSImage def AddBiasNonUniform16(GSImage, bias_level = 500, nsecy = 2, nsecx=8, seed=202102): # Generate Bias and its non-uniformity, and add the 16 bias values to the GS-Image rg = Generator(PCG64(int(seed))) Random16 = (rg.random(nsecy*nsecx)-0.5)*20 if int(bias_level)==0: BiasLevel = np.zeros((nsecy,nsecx)) elif bias_level>0: BiasLevel = Random16.reshape((nsecy,nsecx)) + bias_level arrshape = GSImage.array.shape secsize_x = int(arrshape[1]/nsecx) secsize_y = int(arrshape[0]/nsecy) for rowi in range(nsecy): for coli in range(nsecx): GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] += BiasLevel[rowi,coli] return GSImage def MakeBiasNcomb(npix_x, npix_y, bias_level=500, ncombine=1, read_noise=5, gain=1, seed=202102): # Start with 0 value bias GS-Image ncombine=int(ncombine) BiasSngImg0 = galsim.Image(npix_x, npix_y, init_value=0) BiasSngImg = AddBiasNonUniform16(BiasSngImg0, bias_level=bias_level, nsecy = 2, nsecx=8, seed=int(seed)) BiasCombImg = BiasSngImg*ncombine rng = galsim.UniformDeviate() NoiseBias = galsim.GaussianNoise(rng=rng, sigma=read_noise*ncombine**0.5) BiasCombImg.addNoise(NoiseBias) if ncombine == 1: BiasTag = 'Single' pass elif ncombine >1: BiasCombImg /= ncombine BiasTag = 'Combine' # BiasCombImg.replaceNegative(replace_value=0) # BiasCombImg.quantize() return BiasCombImg, BiasTag def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102): # Generate Gain non-uniformity, and multipy the different factors (mean~1 with sigma~1%) to the GS-Image rg = Generator(PCG64(int(seed))) Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1 # sigma~1% Gain16 = Random16.reshape((nsecy,nsecx))/gain print("Gain of 16 channels: ",Gain16) arrshape = GSImage.array.shape secsize_x = int(arrshape[1]/nsecx) secsize_y = int(arrshape[0]/nsecy) for rowi in range(nsecy): for coli in range(nsecx): GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] *= Gain16[rowi,coli] return GSImage def MakeFlatSmooth(GSBounds, seed): rg = Generator(PCG64(int(seed))) r1,r2,r3,r4 = rg.random(4) a1 = -0.5 + 0.2*r1 a2 = -0.5 + 0.2*r2 a3 = r3+5 a4 = r4+5 xmin,xmax,ymin,ymax = GSBounds.getXMin(), GSBounds.getXMax(), GSBounds.getYMin(), GSBounds.getYMax() Flty, Fltx = np.mgrid[ymin:(ymax+1), xmin:(xmax+1)] rg = Generator(PCG64(int(seed))) p1,p2,bg=rg.poisson(1000, 3) Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20 FlatImg = galsim.ImageF(Fltz) return FlatImg def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan=500, biaslevel=500, seed_bias=20210311): ncombine=int(ncombine) FlatCombImg = flat_single_image*ncombine rng = galsim.UniformDeviate() NoiseFlatPoi = galsim.PoissonNoise(rng=rng, sky_level=0) FlatCombImg.addNoise(NoiseFlatPoi) NoiseFlatReadN = galsim.GaussianNoise(rng=rng, sigma=read_noise*ncombine**0.5) FlatCombImg.addNoise(NoiseFlatReadN) # NoiseFlat = galsim.CCDNoise(rng, gain=gain, read_noise=read_noise*ncombine**0.5, sky_level=0) for i in range(ncombine): FlatCombImg = AddBiasNonUniform16( FlatCombImg, bias_level=biaslevel, nsecy=2, nsecx=8, seed=seed_bias) if ncombine == 1: FlatTag = 'Single' pass elif ncombine >1: FlatCombImg /= ncombine FlatTag = 'Combine' # FlatCombImg.replaceNegative(replace_value=0) # FlatCombImg.quantize() return FlatCombImg, FlatTag def MakeDarkNcomb(npix_x, npix_y, overscan=500, bias_level=500, seed_bias=202102, darkpsec=0.02, exptime=150, ncombine=10, read_noise=5, gain=1): ncombine=int(ncombine) darkpix = darkpsec*exptime DarkSngImg = galsim.Image(npix_x, npix_y, init_value=darkpix) rng = galsim.UniformDeviate() NoiseDarkPoi = galsim.PoissonNoise(rng=rng, sky_level=0) NoiseReadN = galsim.GaussianNoise(rng=rng, sigma=read_noise*ncombine**0.5) DarkCombImg = DarkSngImg*ncombine DarkCombImg.addNoise(NoiseDarkPoi) DarkCombImg.addNoise(NoiseReadN) for i in range(ncombine): DarkCombImg = AddBiasNonUniform16( DarkCombImg, bias_level=bias_level, nsecy = 2, nsecx=8, seed=int(seed_bias)) if ncombine == 1: DarkTag = 'Single' pass elif ncombine >1: DarkCombImg /= ncombine DarkTag = 'Combine' # DarkCombImg.replaceNegative(replace_value=0) # DarkCombImg.quantize() return DarkCombImg, DarkTag def PRNU_Img(xsize, ysize, sigma=0.01, seed=202101): rg = Generator(PCG64(int(seed))) prnuarr = rg.normal(1, sigma, (ysize,xsize)) prnuimg = galsim.ImageF(prnuarr) return prnuimg def NonLinearity(GSImage, beta1=5E-7, beta2=0): NonLinear_f = lambda x, beta_1, beta_2: x - beta_1*x*x + beta_2*x*x*x GSImage.applyNonlinearity(NonLinear_f, beta1, beta2) return GSImage def chargeflow(ndarr, fullwell=10E4): size_y,size_x = ndarr.shape satpos_y = np.where(ndarr>=fullwell)[0] satpos_x = np.where(ndarr>=fullwell)[1] Nsatpix = len(satpos_y) if Nsatpix==0: # make no change for the image array return ndarr else: for i in range(Nsatpix): ClumpFullwell=True satcharge = ndarr[satpos_y[i],satpos_x[i]]-fullwell ndarr[satpos_y[i],satpos_x[i]] = fullwell satpos_yi0 = satpos_y[i] # Define the x,y=0,0 element of image array is the lower-left corner. So y decreases being downward. # print('Charge Clump moves down') chargedn = ((np.random.random()-0.5)*0.05+0.5)*satcharge chargeup = satcharge - chargedn fwi = 1 aa = np.log(chargedn/fullwell)**3*0.9 # blooming length begin to has e- less than fullwell if aa < 0.05: ClumpFullwell=False # Test try: while ClumpFullwell==True: if satpos_y[i]<=0: break if ndarr[satpos_y[i]-1,satpos_x[i]]=fullwell: fx = 0.5*(math.exp(math.log(fwi)**3/aa)+np.exp(-1*math.log(fwi)**3/aa)) if fx>5: fx=5 faa= 0.5*(math.exp(aa/aa)+np.exp(-1*aa/aa)) rand_frac = 1-0.1*(fx-1)/(faa-1) chargedn = ndarr[satpos_y[i]-1,satpos_x[i]] - fullwell ndarr[satpos_y[i]-1,satpos_x[i]] = fullwell*rand_frac satpos_y[i] = satpos_y[i]-1 if satpos_y[i]<0: ClumpFullwell=False break else: ClumpFullwell=False fwi += 1 else: satpos_y[i] = satpos_y[i]-1 if satpos_y[i]<0: ClumpFullwell=False break # print('Charge Clump moves up') ClumpFullwell=True satpos_y[i] = satpos_yi0 fwi = 1 aa = np.log(chargeup/fullwell)**3*0.9 # blooming length at which it begins to have e- less than fullwell if aa < 0.05: ClumpFullwell=False while ClumpFullwell==True: if satpos_y[i]>=size_y-1: break if ndarr[satpos_y[i]+1,satpos_x[i]]=fullwell: fx = 0.5*(math.exp(math.log(fwi)**3/aa)+np.exp(-1*math.log(fwi)**3/aa)) if fx>5: fx=5 faa= 0.5*(math.exp(aa/aa)+np.exp(-1*aa/aa)) rand_frac = 1-0.1*(fx-1)/(faa-1) chargeup = ndarr[satpos_y[i]+1,satpos_x[i]] - fullwell ndarr[satpos_y[i]+1,satpos_x[i]] = fullwell*rand_frac satpos_y[i] = satpos_y[i]+1 if satpos_y[i]>=size_y-1: ClumpFullwell=False break else: ClumpFullwell=False fwi += 1 else: satpos_y[i] = satpos_y[i]+1 if satpos_y[i]>=size_y-1: ClumpFullwell=False break except Exception as e: print(e) print(fwi, aa) pass return ndarr def SaturBloom(GSImage, nsect_x=1, nsect_y=1, fullwell=10e4): """ To simulate digital detector's saturation and blooming effect. The blooming is along the read-out direction, perpendicular to the charge transfer direction. Charge clumpy overflows the pixel well will flow to two oposite directions with nearly same charges. Work together with chargeflow() function. Parameters: GSImage: a GalSim Image of a chip; nsect_x: number of sections divided along x direction; nsect_y: number of sections divided along y direction. Return: a saturated image array with the same size of input. """ imgarr = GSImage.array size_y, size_x = imgarr.shape subsize_y = int(size_y/nsect_y) subsize_x = int(size_x/nsect_x) for i in range(nsect_y): for j in range(nsect_x): subimg = imgarr[subsize_y*i:subsize_y*(i+1), subsize_x*j:subsize_x*(j+1)] subimg = chargeflow(subimg, fullwell=fullwell) imgarr[subsize_y*i:subsize_y*(i+1), subsize_x*j:subsize_x*(j+1)] = subimg return GSImage def readout16(GSImage, rowi=0, coli=0, overscan_value=0): # readout image as 16 outputs of sub-images plus prescan & overscan. # assuming image width and height sizes are both even. # assuming image has 2 columns and 8 rows of output channels. # 00 01 # 10 11 # 20 21 # ... # return: GS Image Object npix_y,npix_x = GSImage.array.shape subheight = int(8+npix_y/2+8) subwidth = int(16+npix_x/8+27) OutputSubimg = galsim.ImageUS(subwidth, subheight, init_value=overscan_value) if rowi<4 and coli==0: subbounds = galsim.BoundsI(1, int(npix_x/2), int(npix_y/8*rowi+1), int(npix_y/8*(rowi+1))) subbounds = subbounds.shift(galsim.PositionI(GSImage.bounds.getXMin()-1, GSImage.bounds.getYMin()-1)) subimg = GSImage[subbounds] OutputSubimg.array[27:int(npix_y/8)+27,8:int(npix_x/2)+8] = subimg.array elif rowi<4 and coli==1: subbounds = galsim.BoundsI(npix_x/2+1, npix_x, npix_y/8*rowi+1, npix_y/8*(rowi+1)) subbounds = subbounds.shift(galsim.PositionI(GSImage.bounds.getXMin()-1, GSImage.bounds.getYMin()-1)) subimg = GSImage[subbounds] OutputSubimg.array[27:int(npix_y/8)+27,8:int(npix_x/2)+8] = subimg.array elif rowi>=4 and rowi<8 and coli==0: subbounds = galsim.BoundsI(1, npix_x/2, npix_y/8*rowi+1, npix_y/8*(rowi+1)) subbounds = subbounds.shift(galsim.PositionI(GSImage.bounds.getXMin()-1, GSImage.bounds.getYMin()-1)) subimg = GSImage[subbounds] OutputSubimg.array[16:int(npix_y/8)+16,8:int(npix_x/2)+8] = subimg.array elif rowi>=4 and rowi<8 and coli==1: subbounds = galsim.BoundsI(npix_x/2+1, npix_x, npix_y/8*rowi+1, npix_y/8*(rowi+1)) subbounds = subbounds.shift(galsim.PositionI(GSImage.bounds.getXMin()-1, GSImage.bounds.getYMin()-1)) subimg = GSImage[subbounds] OutputSubimg.array[16 :int(npix_y/8)+16,8:int(npix_x/2)+8] = subimg.array else: print("\n\033[31mError: "+"Wrong rowi or coli assignment. Permitted: 0<=rowi<=7, 0<=coli<=1."+"\033[0m\n") return OutputSubimg return OutputSubimg def CTE_Effect(GSImage, threshold=27, direction='column'): # Devide the image into 4 sections and apply CTE effect with different trail directions. # GSImage: a GalSim Image object. size_y,size_x = GSImage.array.shape size_sect_y = int(size_y/2) size_sect_x = int(size_x/2) imgarr = GSImage.array if direction == 'column': imgarr[0:size_sect_y,:] = CTEModelColRow(imgarr[0:size_sect_y,:], trail_direction='down', direction='column', threshold=threshold) imgarr[size_sect_y:size_y,:] = CTEModelColRow(imgarr[size_sect_y:size_y,:], trail_direction='up', direction='column', threshold=threshold) elif direction == 'row': imgarr[:,0:size_sect_x] = CTEModelColRow(imgarr[:,0:size_sect_x], trail_direction='right', direction='row', threshold=threshold) imgarr[:,size_sect_x:size_x] = CTEModelColRow(imgarr[:,size_sect_x:size_x], trail_direction='left', direction='row', threshold=threshold) return GSImage @jit() def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27): #total trail flux vs (pixel flux)^1/2 is approximately linear #total trail flux = trail_a * (pixel flux)^1/2 + trail_b #trail pixel flux = pow(0.5,x)/0.5, normalize to 1 trail_a = 5.651803799619966 trail_b = -2.671933068990729 sh1 = img.shape[0] sh2 = img.shape[1] n_img = img*0 idx = np.where(img0: n_img[idx] = img[idx] yidx,xidx = np.where(img>=threshold) if len(yidx) == 0: pass elif len(yidx)>0: # print(index) for i, j in zip(yidx,xidx): f = img[i,j] trail_f = (np.sqrt(f)*trail_a + trail_b)*0.5 # trail_f=5E-5*f**1.5 xy_num = 10 all_trail = np.zeros(xy_num) xy_upstr = np.arange(1,xy_num,1) # all_trail_pix = np.sum(pow(0.5,xy_upstr)/0.5) all_trail_pix = 0 for m in xy_upstr: a1=12.97059491 b1=0.54286652 c1=0.69093105 a2=2.77298856 b2=0.11231055 c2=-0.01038675 # t_pow = 0 am=1 bm=1 t_pow = am*np.exp(-bm*m) # if m < 5: # t_pow = a1*np.exp(-b1*m)+c1 # else: # t_pow = a2*np.exp(-b2*m)+c2 if t_pow <0: t_pow = 0 all_trail_pix += t_pow all_trail[m] = t_pow trail_pix_eff = trail_f/all_trail_pix all_trail = trail_pix_eff*all_trail all_trail[0] = f - trail_f for m in np.arange(0,xy_num,1): if direction == 'column': if trail_direction == 'down': y_pos = i + m elif trail_direction == 'up': y_pos = i - m if y_pos < 0 or y_pos >=sh1: break n_img[y_pos,j] = n_img[y_pos,j] + all_trail[m] elif direction == 'row': if trail_direction == 'left': x_pos = j - m elif trail_direction == 'right': x_pos = j + m if x_pos < 0 or x_pos >=sh2: break n_img[i,x_pos] = n_img[i,x_pos] + all_trail[m] return n_img #---------- For Cosmic-Ray Simulation ------------ #---------- Zhang Xin ---------------------------- def getYValue(collection, x): index = 0; if (collection.shape[1] == 2): while(x>collection[index,0] and index < collection.shape[0]): index= index + 1; if (index == collection.shape[0] or index == 0): return 0; deltX = collection[index, 0] - collection[index-1, 0]; deltY = collection[index, 1] - collection[index-1, 1]; if deltX == 0: return (collection[index, 1] + collection[index-1, 1])/2.0 else: a = deltY/deltX; return a * (x - collection[index-1,0]) + collection[index-1, 1]; return 0; def selectCosmicRayCollection(attachedSizes, xLen, yLen,cr_pixelRatio,CR_max_size): normalRay = 0.90 nnormalRay = 1-normalRay max_nrayLen = 100 pixelNum = int(xLen * yLen * cr_pixelRatio * normalRay ); pixelNum_n = int(xLen * yLen * cr_pixelRatio * nnormalRay ) CRPixelNum = 0; maxValue = max(attachedSizes[:,1]) maxValue += 0.1; cr_event_num = 0; CRs = np.zeros(pixelNum); while (CRPixelNum < pixelNum): x = CR_max_size * np.random.random(); y = maxValue * np.random.random(); if (y <= getYValue(attachedSizes, x)): CRs[cr_event_num] = np.ceil(x); cr_event_num = cr_event_num + 1; CRPixelNum = CRPixelNum + round(x); while (CRPixelNum < pixelNum + pixelNum_n): nx = np.random.random()*(max_nrayLen-CR_max_size)+CR_max_size CRs[cr_event_num] = np.ceil(nx); cr_event_num = cr_event_num + 1; CRPixelNum = CRPixelNum + np.ceil(nx); return CRs[0:cr_event_num]; def defineEnergyForCR(cr_event_size): import random sigma = 0.6 / 2.355; mean = 3.3; energys = np.zeros(cr_event_size); for i in np.arange(cr_event_size): energy_index = random.normalvariate(mean,sigma); energys[i] = pow(10, energy_index); return energys; def convCR(CRmap=None, addPSF=None, sp_n = 4): sh = CRmap.shape # sp_n = 4 subCRmap = np.zeros(np.array(sh)*sp_n) pix_v0 = 1/(sp_n*sp_n) for i in np.arange(sh[0]): i_st = sp_n*i for j in np.arange(sh[1]): if CRmap[i,j] ==0: continue j_st = sp_n*j pix_v1 = CRmap[i,j]*pix_v0 for m in np.arange(sp_n): for n in np.arange(sp_n): subCRmap[i_st+m, j_st + n] = pix_v1 m_size = addPSF.shape[0] subCRmap_n = np.zeros(np.array(subCRmap.shape) + m_size -1) for i in np.arange(subCRmap.shape[0]): for j in np.arange(subCRmap.shape[1]): if subCRmap[i,j]>0: convPix = addPSF*subCRmap[i,j] subCRmap_n[i:i+m_size,j:j+m_size]+=convPix CRmap_n = np.zeros((np.array(subCRmap_n.shape)/sp_n).astype(np.int32)) sh_n = CRmap_n.shape for i in np.arange(sh_n[0]): i_st = sp_n*i for j in np.arange(sh_n[1]): p_v = 0 j_st=sp_n*j for m in np.arange(sp_n): for n in np.arange(sp_n): p_v += subCRmap_n[i_st+m, j_st + n] CRmap_n[i,j] = p_v return CRmap_n def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=20210317): # Return: an 2-D numpy array # attachedSizes = np.loadtxt('./wfc-cr-attachpixel.dat'); np.random.seed(seed) CR_max_size = 20.0; cr_size = selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size); cr_event_size = cr_size.shape[0]; cr_energys = defineEnergyForCR(cr_event_size); CRmap = np.zeros([yLen, xLen]); ## produce conv kernel from astropy.modeling.models import Gaussian2D o_size = 4 sp_n = 8 m_size = o_size*sp_n+1 m_cen = o_size*sp_n/2 sigma_psf = 0.2*sp_n addPSF_ = Gaussian2D(1, m_cen, m_cen, sigma_psf, sigma_psf) yp, xp = np.mgrid[0:m_size, 0:m_size] addPSF = addPSF_(xp, yp) convKernel = addPSF/addPSF.sum() ################################# for i in np.arange(cr_event_size): xPos = round((xLen - 1)* np.random.random()); yPos = round((yLen - 1)* np.random.random()); cr_lens = int(cr_size[i]); if cr_lens ==0: continue; pix_energy = cr_energys[i]/gain/cr_lens; pos_angle = 1/2*math.pi*np.random.random(); crMatrix = np.zeros([cr_lens+1, cr_lens + 1]) for j in np.arange(cr_lens): x_n = int(np.cos(pos_angle)*j - np.sin(pos_angle)*0); if x_n < 0: x_n = x_n + cr_lens+1 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: continue; crMatrix[y_n,x_n] = pix_energy; crMatrix_n = convCR(crMatrix, convKernel, sp_n) # crMatrix_n = crMatrix xpix = np.arange(crMatrix_n.shape[0]) + int(xPos) ypix = np.arange(crMatrix_n.shape[1]) + int(yPos) sh = CRmap.shape okx = (xpix >= 0) & (xpix < sh[1]) oky = (ypix >= 0) & (ypix < sh[0]) sly = slice(ypix[oky].min(), ypix[oky].max()+1) slx = slice(xpix[okx].min(), xpix[okx].max()+1) CRmap[sly, slx] += crMatrix_n[oky,:][:,okx] return CRmap.astype(np.int32), cr_event_size def ShutterEffectArr(GSImage, t_shutter=1.3, dist_bearing=735, dt=1E-3): # Generate Shutter-Effect normalized image # t_shutter: time of shutter movement # dist_bearing: distance between two bearings of shutter leaves # dt: delta_t of sampling from scipy import interpolate SampleNumb = int(t_shutter/dt+1) DistHalf = dist_bearing/2 t = np.arange(SampleNumb)*dt a_arr = 5.84*np.sin(2*math.pi/t_shutter*t) v = np.zeros(SampleNumb) theta = np.zeros(SampleNumb) x = np.arange(SampleNumb)/(SampleNumb-1)*dist_bearing s = np.zeros(SampleNumb) s1 = np.zeros(SampleNumb) s2 = np.zeros(SampleNumb) brt = np.zeros(SampleNumb) idx = np.arange(SampleNumb) sidx = np.zeros(SampleNumb) s1idx = np.zeros(SampleNumb) s2idx = np.zeros(SampleNumb) v[0] = 0 theta[0] = 0 for i in range(SampleNumb-1): v[i+1] = v[i]+a_arr[i]*dt theta[i+1] = theta[i]+v[i]*dt s1[i] = DistHalf*np.cos(theta[i]) s2[i] = dist_bearing-DistHalf*np.cos(theta[i]) s1idx[i] = int(s1[i]/dist_bearing*(SampleNumb)) s2idx[i] = int(s2[i]/dist_bearing*(SampleNumb)) brt[(idx>s1idx[i]) & (idxnp.max(x): raise LookupError("Out of focal-plane bounds in X-direction.") if ymin<-25331 or ymax>25331: raise LookupError("Out of focal-plane bounds in Y-direction.") sizex = xmax-xmin+1 sizey = ymax-ymin+1 xnewgrid = np.mgrid[xmin:(xmin+sizex)] expeffect = interpolate.splev(xnewgrid, intp, der=0) expeffect /= np.max(expeffect) exparrnormal = np.tile(expeffect, (sizey,1)) # GSImage *= exparrnormal return exparrnormal