diff --git a/observation_sim/instruments/chip/effects.py b/observation_sim/instruments/chip/effects.py index 7620f85ce1c2a31e86e83155928af9ea743fb584..f215a214a2579f9097d1bfb0115e33567c54a4e9 100644 --- a/observation_sim/instruments/chip/effects.py +++ b/observation_sim/instruments/chip/effects.py @@ -3,7 +3,8 @@ 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 +import math +import copy from numba import jit from astropy import stats @@ -50,11 +51,11 @@ def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed= NPixBad = int(NPix*fraction) NPixHot = int(NPix*fraction*HotFraction) NPixDead = NPixBad-NPixHot - - NPix_y,NPix_x = GSImage.array.shape + + NPix_y, NPix_x = GSImage.array.shape mean = np.mean(GSImage.array) rgp = Generator(PCG64(int(seed))) - yxposfrac = rgp.random((NPixBad,2)) + 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) @@ -108,7 +109,7 @@ def AddBiasNonUniform16(GSImage, bias_level=500, nsecy=2, nsecx=8, seed=202102, if int(bias_level) == 0: BiasLevel = np.zeros((nsecy, nsecx)) elif bias_level > 0: - BiasLevel = Random16.reshape((nsecy,nsecx)) + bias_level + BiasLevel = Random16.reshape((nsecy, nsecx)) + bias_level if logger is not None: msg = str(" Biases of 16 channels: " + str(BiasLevel)) logger.info(msg) @@ -129,7 +130,7 @@ def MakeBiasNcomb(npix_x, npix_y, bias_level=500, ncombine=1, read_noise=5, gain BiasSngImg0 = galsim.Image(npix_x, npix_y, init_value=0) BiasSngImg = AddBiasNonUniform16(BiasSngImg0, bias_level=bias_level, - nsecy = 2, nsecx=8, + nsecy=2, nsecx=8, seed=int(seed), logger=logger) BiasCombImg = BiasSngImg*ncombine @@ -198,7 +199,7 @@ def MakeFlatSmooth(GSBounds, seed): 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 @@ -245,7 +246,7 @@ def MakeDarkNcomb(npix_x, npix_y, overscan=500, bias_level=500, seed_bias=202102 DarkCombImg = AddBiasNonUniform16( DarkCombImg, bias_level=bias_level, - nsecy = 2, nsecx=8, + nsecy=2, nsecx=8, seed=int(seed_bias), logger=logger) if ncombine == 1: @@ -272,7 +273,7 @@ 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 @@ -280,7 +281,7 @@ def BleedingTrail(aa, yy): 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) + faa = 0.5*(math.e+1/math.e) trail_frac = 1-0.1*(fy-1)/(faa-1) except Exception as e: print(e) @@ -306,14 +307,14 @@ def MakeTrail(imgarr, satuyxtuple, charge, fullwell=9e4, direction='up', trailcu if direction == 'up': if imgarr[yi-1, xi] >= fullwell: imgarr[yi, xi] = fullwell - yi-=1 + yi -= 1 continue elif direction == 'down': if imgarr[yi+1, xi] >= fullwell: imgarr[yi, xi] = fullwell yi += 1 continue - if aa<=1: + if aa <= 1: while imgarr[yi, xi] >= fullwell: imgarr[yi, xi] = fullwell if direction == 'up': @@ -359,9 +360,9 @@ def MakeTrail(imgarr, satuyxtuple, charge, fullwell=9e4, direction='up', trailcu def ChargeFlow(imgarr, fullwell=9E4): size_y, size_x = imgarr.shape - satupos_y, satupos_x = np.where(imgarr>fullwell) + satupos_y, satupos_x = np.where(imgarr > fullwell) - if satupos_y.shape[0]==0: + if satupos_y.shape[0] == 0: # make no change for the image array return imgarr elif satupos_y.shape[0]/imgarr.size > 0.5: @@ -383,16 +384,16 @@ def ChargeFlow(imgarr, fullwell=9E4): try: # Charge Clump moves up - if yi>=0 and yi= 0 and yi < imgarr.shape[0]: imgarr = MakeTrail(imgarr, (yi, xi), chargeup, fullwell=9e4, direction='up', trailcutfrac=0.9) # Charge Clump moves down imgarr = MakeTrail(imgarr, (yi, xi), chargedn, fullwell=9e4, direction='down', trailcutfrac=0.9) except Exception as e: - print(e,'@pix ',(yi+1, xi+1)) + print(e, '@pix ', (yi+1, xi+1)) return imgarr - return imgarr + def SaturBloom(GSImage, nsect_x=1, nsect_y=1, fullwell=9e4): """ 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. @@ -424,7 +425,7 @@ def SaturBloom(GSImage, nsect_x=1, nsect_y=1, fullwell=9e4): 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. + # assuming image has 2 columns and 8 rows of output channels. # 00 01 # 10 11 # 20 21 @@ -453,7 +454,7 @@ def readout16(GSImage, rowi=0, coli=0, overscan_value=0): 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 @@ -469,19 +470,19 @@ def CTE_Effect(GSImage, threshold=27, direction='column'): 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[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 @jit() -def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27): +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 + # 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 @@ -491,7 +492,7 @@ def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27 idx = np.where(img < threshold) if len(idx[0]) == 0: pass - elif len(idx[0])>0: + elif len(idx[0]) > 0: n_img[idx] = img[idx] yidx, xidx = np.where(img >= threshold) @@ -505,13 +506,13 @@ def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27 # 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 + a1 = 12.97059491 b1 = 0.54286652 c1 = 0.69093105 a2 = 2.77298856 @@ -530,12 +531,8 @@ def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27 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): @@ -559,10 +556,8 @@ def CTEModelColRow(img, trail_direction = 'up', direction='column', threshold=27 return n_img - -#---------- For Cosmic-Ray Simulation ------------ -#---------- Zhang Xin ---------------------------- - +# ---------- For Cosmic-Ray Simulation ------------ +# ---------- Zhang Xin ---------------------------- def getYValue(collection, x): index = 0; if (collection.shape[1] == 2):