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 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=500): # 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=20210309, 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)] 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)) signs = 2*rgdn.integers(0,2,size=(nbadsecA+nbadsecD))-1 dn = rgdn.integers(low=meanimg*0.4, high=meanimg*0.8, size=(nbadsecA+nbadsecD))*signs for badcoli in range(nbadsecA): GSImage.array[(ysize-collen[badcoli]):ysize,xposit[badcoli]:(xposit[badcoli]+1)] += np.random.random((collen[badcoli],1))*stdimg*3+dn[badcoli] for badcoli in range(nbadsecD): GSImage.array[0:collen[badcoli+nbadsecA],xposit[badcoli+nbadsecA]:(xposit[badcoli+nbadsecA]+1)] += np.random.random((collen[badcoli+nbadsecA],1))*stdimg*3+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=1E-7, beta2=1E-10): 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); 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