From b8b5bed920ec9855906edb37d6ece2679c603739 Mon Sep 17 00:00:00 2001 From: Chengliang Date: Sat, 26 Oct 2024 09:36:35 +0800 Subject: [PATCH] update codestyle-PEP8 --- observation_sim/instruments/chip/effects.py | 326 ++++++++++---------- 1 file changed, 162 insertions(+), 164 deletions(-) diff --git a/observation_sim/instruments/chip/effects.py b/observation_sim/instruments/chip/effects.py index 21735ed..7620f85 100644 --- a/observation_sim/instruments/chip/effects.py +++ b/observation_sim/instruments/chip/effects.py @@ -22,7 +22,7 @@ def AddOverscan(GSImage, overscan=1000, gain=1, widthl=27, widthr=27, widtht=8, 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.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)) @@ -37,11 +37,11 @@ def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed= # Hot Pixel > 20e-/s # Dead Pixel < 70%*Mean rgf = Generator(PCG64(int(seed*1.1))) - if IfHotPix==True and IfDeadPix==True: + if IfHotPix is True and IfDeadPix is True: HotFraction = rgf.random() # fraction in total bad pixels - elif IfHotPix==False and IfDeadPix==False: + elif IfHotPix is False and IfDeadPix is False: return GSImage - elif IfHotPix==True: + elif IfHotPix is True: HotFraction = 1 else: HotFraction = 0 @@ -55,23 +55,23 @@ def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed= 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) + 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 + if IfHotPix is True: + GSImage.array[YPositHot, XPositHot] += rgh.gamma(2, 25*150, size=NPixHot) + if IfDeadPix is 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, logger=None): # Set bad column values - ysize,xsize = GSImage.array.shape + 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) @@ -82,8 +82,8 @@ def BadColumns(GSImage, seed=20240309, chipid=1, logger=None): 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)) + 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)) if logger is not None: logger.info(xposit+1) @@ -91,47 +91,47 @@ def BadColumns(GSImage, seed=20240309, chipid=1, logger=None): 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 + 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]) + 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]) + 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, logger=None): +def AddBiasNonUniform16(GSImage, bias_level=500, nsecy=2, nsecx=8, seed=202102, logger=None): # 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: + if int(bias_level) == 0: + BiasLevel = np.zeros((nsecy, nsecx)) + elif bias_level > 0: BiasLevel = Random16.reshape((nsecy,nsecx)) + bias_level if logger is not None: msg = str(" Biases of 16 channels: " + str(BiasLevel)) logger.info(msg) else: - print(" Biases of 16 channels:\n",BiasLevel) + 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] + 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, logger=None): # Start with 0 value bias GS-Image - ncombine=int(ncombine) + 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), - logger=logger) + BiasSngImg = AddBiasNonUniform16(BiasSngImg0, + bias_level=bias_level, + nsecy = 2, nsecx=8, + seed=int(seed), + logger=logger) BiasCombImg = BiasSngImg*ncombine rng = galsim.UniformDeviate() NoiseBias = galsim.GaussianNoise(rng=rng, sigma=read_noise*ncombine**0.5) @@ -139,7 +139,7 @@ def MakeBiasNcomb(npix_x, npix_y, bias_level=500, ncombine=1, read_noise=5, gain if ncombine == 1: BiasTag = 'Single' pass - elif ncombine >1: + elif ncombine > 1: BiasCombImg /= ncombine BiasTag = 'Combine' # BiasCombImg.replaceNegative(replace_value=0) @@ -147,32 +147,32 @@ def MakeBiasNcomb(npix_x, npix_y, bias_level=500, ncombine=1, read_noise=5, gain return BiasCombImg, BiasTag -def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=None): +def ApplyGainNonUniform16(GSImage, gain=1, nsecy=2, nsecx=8, seed=202102, logger=None): # 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 + Gain16 = Random16.reshape((nsecy, nsecx))/gain gain_array = np.ones(nsecy*nsecx)*gain if logger is not None: msg = str("Gain of 16 channels: " + str(Gain16)) logger.info(msg) else: - print("Gain of 16 channels: ",Gain16) + 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] - gain_array[rowi*nsecx+coli] = 1/Gain16[rowi,coli] + GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y, coli*secsize_x:(coli+1)*secsize_x] *= Gain16[rowi, coli] + gain_array[rowi*nsecx+coli] = 1/Gain16[rowi, coli] return GSImage, gain_array -def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=None): +def GainsNonUniform16(GSImage, gain=1, nsecy=2, nsecx=8, seed=202102, logger=None): # 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 + Gain16 = Random16.reshape((nsecy, nsecx))/gain if logger is not None: msg = str(seed-20210202, "Gains of 16 channels: " + str(Gain16)) logger.info(msg) @@ -190,22 +190,22 @@ def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=N def MakeFlatSmooth(GSBounds, seed): rg = Generator(PCG64(int(seed))) - r1,r2,r3,r4 = rg.random(4) + 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() + 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) + p1, p2, bg=rg.poisson(1000, 3) Fltz = 0.6*1e-7*(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, logger=None): - ncombine=int(ncombine) + ncombine = int(ncombine) FlatCombImg = flat_single_image*ncombine rng = galsim.UniformDeviate() NoiseFlatPoi = galsim.PoissonNoise(rng=rng, sky_level=0) @@ -215,15 +215,15 @@ def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan= # 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, + FlatCombImg, + bias_level=biaslevel, + nsecy=2, nsecx=8, seed=seed_bias, logger=logger) if ncombine == 1: FlatTag = 'Single' pass - elif ncombine >1: + elif ncombine > 1: FlatCombImg /= ncombine FlatTag = 'Combine' # FlatCombImg.replaceNegative(replace_value=0) @@ -232,7 +232,7 @@ def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan= 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, logger=None): - ncombine=int(ncombine) + ncombine = int(ncombine) darkpix = darkpsec*exptime DarkSngImg = galsim.Image(npix_x, npix_y, init_value=darkpix) rng = galsim.UniformDeviate() @@ -243,15 +243,15 @@ def MakeDarkNcomb(npix_x, npix_y, overscan=500, bias_level=500, seed_bias=202102 DarkCombImg.addNoise(NoiseReadN) for i in range(ncombine): DarkCombImg = AddBiasNonUniform16( - DarkCombImg, - bias_level=bias_level, - nsecy = 2, nsecx=8, + DarkCombImg, + bias_level=bias_level, + nsecy = 2, nsecx=8, seed=int(seed_bias), logger=logger) if ncombine == 1: DarkTag = 'Single' pass - elif ncombine >1: + elif ncombine > 1: DarkCombImg /= ncombine DarkTag = 'Combine' # DarkCombImg.replaceNegative(replace_value=0) @@ -261,7 +261,7 @@ def MakeDarkNcomb(npix_x, npix_y, overscan=500, bias_level=500, seed_bias=202102 def PRNU_Img(xsize, ysize, sigma=0.01, seed=202101): rg = Generator(PCG64(int(seed))) - prnuarr = rg.normal(1, sigma, (ysize,xsize)) + prnuarr = rg.normal(1, sigma, (ysize, xsize)) prnuimg = galsim.ImageF(prnuarr) return prnuimg @@ -272,11 +272,10 @@ def NonLinearity(GSImage, beta1=5E-7, beta2=0): return GSImage -######################################## Saturation & Bleeding Start ############################### - +#Saturation & Bleeding Start# def BleedingTrail(aa, yy): - if aa<0.2: - aa=0.2 + if aa < 0.2: + aa = 0.2 else: pass try: @@ -289,77 +288,78 @@ def BleedingTrail(aa, yy): 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 + 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: + 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 + 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 + 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 + 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: + 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]: + 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) + 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 + 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 + imgarr[yi, xi] += charge + if imgarr[yi, xi] > fullwell: + 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) + 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 @@ -371,12 +371,12 @@ def ChargeFlow(imgarr, fullwell=9E4): 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)]) + chargedict[yxidx] = imgarrorig[yi, xi]-fullwell - for yi,xi in zip(satupos_y,satupos_x): - yxidx = ''.join([str(yi),str(xi)]) + 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 @@ -384,11 +384,11 @@ def ChargeFlow(imgarr, fullwell=9E4): try: # Charge Clump moves up if yi>=0 and yi=4 and rowi<8 and coli==0: + 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: + 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 + 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 @@ -463,16 +463,16 @@ def readout16(GSImage, rowi=0, coli=0, overscan_value=0): 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_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) + 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) + 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 @@ -482,50 +482,50 @@ 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_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) + yidx, xidx = np.where(img >= threshold) if len(yidx) == 0: pass - elif len(yidx)>0: + elif len(yidx) > 0: # print(index) - for i, j in zip(yidx,xidx): - f = img[i,j] + 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) + 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 + 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 + 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: + if t_pow < 0: t_pow = 0 all_trail_pix += t_pow @@ -538,23 +538,23 @@ def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27 all_trail[0] = f - trail_f - for m in np.arange(0,xy_num,1): + 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: + if y_pos < 0 or y_pos >= sh1: break - n_img[y_pos,j] = n_img[y_pos,j] + all_trail[m] + 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: + if x_pos < 0 or x_pos >= sh2: break - n_img[i,x_pos] = n_img[i,x_pos] + all_trail[m] + n_img[i, x_pos] = n_img[i, x_pos] + all_trail[m] return n_img @@ -566,7 +566,7 @@ def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27 def getYValue(collection, x): index = 0; 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; if (index == collection.shape[0] or index == 0): return 0; @@ -578,11 +578,11 @@ def getYValue(collection, x): 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 a * (x - collection[index-1, 0]) + collection[index-1, 1]; return 0; -def selectCosmicRayCollection(attachedSizes, xLen, yLen,cr_pixelRatio,CR_max_size): +def selectCosmicRayCollection(attachedSizes, xLen, yLen, cr_pixelRatio, CR_max_size): normalRay = 0.90 nnormalRay = 1-normalRay @@ -591,7 +591,7 @@ def selectCosmicRayCollection(attachedSizes, xLen, yLen,cr_pixelRatio,CR_max_siz pixelNum_n = int(xLen * yLen * cr_pixelRatio * nnormalRay ) CRPixelNum = 0; - maxValue = max(attachedSizes[:,1]) + maxValue = max(attachedSizes[:, 1]) maxValue += 0.1; cr_event_num = 0; @@ -613,7 +613,7 @@ def selectCosmicRayCollection(attachedSizes, xLen, yLen,cr_pixelRatio,CR_max_siz return CRs[0:cr_event_num]; -def defineEnergyForCR(cr_event_size,seed = 12345): +def defineEnergyForCR(cr_event_size, seed = 12345): import random sigma = 0.6 / 2.355; mean = 3.3; @@ -637,7 +637,7 @@ def convCR(CRmap=None, addPSF=None, sp_n = 4): if CRmap[i,j] ==0: continue j_st = sp_n*j - pix_v1 = CRmap[i,j]*pix_v0 + 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 @@ -648,9 +648,9 @@ def convCR(CRmap=None, addPSF=None, sp_n = 4): 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 + 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 @@ -664,7 +664,7 @@ def convCR(CRmap=None, addPSF=None, sp_n = 4): for n in np.arange(sp_n): p_v += subCRmap_n[i_st+m, j_st + n] - CRmap_n[i,j] = p_v + CRmap_n[i, j] = p_v return CRmap_n @@ -681,7 +681,7 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 CRmap = np.zeros([yLen, xLen]); - ## produce conv kernel + # produce conv kernel from astropy.modeling.models import Gaussian2D o_size = 4 sp_n = 8 @@ -694,7 +694,7 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 addPSF = addPSF_(xp, yp) convKernel = addPSF/addPSF.sum() - ################################# + #--------------------------------- for i in np.arange(cr_event_size): @@ -713,9 +713,9 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 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: + 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[y_n, x_n] = pix_energy; crMatrix_n = convCR(crMatrix, convKernel, sp_n) # crMatrix_n = crMatrix @@ -729,9 +729,7 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 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] - - + CRmap[sly, slx] += crMatrix_n[oky, :][:, okx] return CRmap.astype(np.int32), cr_event_size @@ -770,9 +768,9 @@ def ShutterEffectArr(GSImage, t_exp=150, t_shutter=1.3, dist_bearing=735, dt=1E- 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]) & (idx s1idx[i]) & (idx < s2idx[i])] += dt - if t_exp>t_shutter*2: + if t_exp > t_shutter*2: brt = brt*2+(t_exp-t_shutter*2) else: brt = brt*2 @@ -785,16 +783,16 @@ def ShutterEffectArr(GSImage, t_exp=150, t_shutter=1.3, dist_bearing=735, dt=1E- xmax = GSImage.bounds.getXMax() ymin = GSImage.bounds.getYMin() ymax = GSImage.bounds.getYMax() - if xminnp.max(x): + if xmin < np.min(x) or xmax > np.max(x): raise LookupError("Out of focal-plane bounds in X-direction.") - if ymin<-25331 or ymax>25331: + 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 /= t_exp - exparrnormal = np.tile(expeffect, (sizey,1)) + exparrnormal = np.tile(expeffect, (sizey, 1)) # GSImage *= exparrnormal return exparrnormal -- GitLab