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,copy 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=np.abs(meanimg)*1.3+50, high=np.abs(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 print(" Biases of 16 channels:\n",BiasLevel) 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 GainsNonUniform16(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(seed-20210202, "Gains of 16 channels:\n", 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 return Gain16 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 ######################################## Saturation & Bleeding Start ############################### def BleedingTrail(aa, yy): if aa<0.2: aa=0.2 else: pass try: fy = 0.5*(math.exp(math.log(yy+1)**3/aa)+np.exp(-1*math.log(yy+1)**3/aa)) faa= 0.5*(math.e+1/math.e) trail_frac = 1-0.1*(fy-1)/(faa-1) except Exception as e: print(e) trail_frac = 1 return trail_frac def MakeTrail(imgarr, satuyxtuple, charge, fullwell=9e4, direction='up', trailcutfrac=0.9): ''' direction: "up" or "down". For "up", bleeds along Y-decreasing direction; for "down", bleeds along Y-increasing direction. ''' yi,xi = satuyxtuple aa = np.log(charge/fullwell)**3 # scale length of the bleeding trail yy = 1 while charge>0: if yi<0 or yi>imgarr.shape[0]-1: break if yi==0 or yi==imgarr.shape[0]-1: imgarr[yi,xi] = fullwell break if direction=='up': if imgarr[yi-1,xi]>=fullwell: imgarr[yi,xi] = fullwell yi-=1 continue elif direction=='down': if imgarr[yi+1,xi]>=fullwell: imgarr[yi,xi] = fullwell yi+=1 continue if aa<=1: while imgarr[yi,xi] >= fullwell: imgarr[yi,xi] = fullwell if direction=='up': imgarr[yi-1,xi] += charge charge = imgarr[yi-1,xi]-fullwell yi-=1 if yi<0: break elif direction=='down': imgarr[yi+1,xi] += charge charge = imgarr[yi+1,xi]-fullwell yi+=1 if yi>imgarr.shape[0]: break else: # calculate bleeding trail: trail_frac = BleedingTrail(aa,yy) # put charge upwards if trail_frac>=0.99: imgarr[yi,xi] = fullwell if direction=='up': yi-=1 elif direction=='down': yi+=1 yy += 1 else: if trail_fracfullwell: imgarr[yi,xi] = fullwell if direction=='up': yi-=1 elif direction=='down': yi+=1 yy += 1 return imgarr def ChargeFlow(imgarr, fullwell=9E4): size_y,size_x = imgarr.shape satupos_y,satupos_x = np.where(imgarr>fullwell) if satupos_y.shape[0]==0: # make no change for the image array return imgarr elif satupos_y.shape[0]/imgarr.size > 0.5: imgarr.fill(fullwell) return imgarr chargedict = {} imgarrorig = copy.deepcopy(imgarr) for yi,xi in zip(satupos_y,satupos_x): yxidx = ''.join([str(yi),str(xi)]) chargedict[yxidx] = imgarrorig[yi,xi]-fullwell for yi,xi in zip(satupos_y,satupos_x): yxidx = ''.join([str(yi),str(xi)]) satcharge = chargedict[yxidx] chargeup = ((np.random.random()-0.5)*0.05+0.5)*satcharge chargedn = satcharge - chargeup try: # Charge Clump moves up if yi>=0 and yi=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_exp=150, 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]) & (idxt_shutter*2: brt = brt*2+(t_exp-t_shutter*2) else: brt = brt*2 x = (x-dist_bearing/2)*100 intp = interpolate.splrep(x, brt, s=0) xmin = GSImage.bounds.getXMin() xmax = GSImage.bounds.getXMax() ymin = GSImage.bounds.getYMin() ymax = GSImage.bounds.getYMax() if xminnp.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