diff --git a/observation_sim/instruments/chip/effects.py b/observation_sim/instruments/chip/effects.py
index 21735ede4ef8e5307fc645d75c4592b9e0c0e2d5..7620f85ce1c2a31e86e83155928af9ea743fb584 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_frac<trailcutfrac:
+ if trail_frac < trailcutfrac:
break
charge = fullwell*trail_frac
- imgarr[yi,xi] += charge
- if imgarr[yi,xi]>fullwell:
- 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<imgarr.shape[0]:
- imgarr = MakeTrail(imgarr, (yi,xi), chargeup, fullwell=9e4, direction='up', trailcutfrac=0.9)
+ 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)
+ 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
@@ -418,7 +418,7 @@ def SaturBloom(GSImage, nsect_x=1, nsect_y=1, fullwell=9e4):
return GSImage
-################################# Saturation & Bleeding End ####################################
+# Saturation & Bleeding End #
def readout16(GSImage, rowi=0, coli=0, overscan_value=0):
@@ -430,30 +430,30 @@ def readout16(GSImage, rowi=0, coli=0, overscan_value=0):
# 20 21
# ...
# return: GS Image Object
- npix_y,npix_x = GSImage.array.shape
+ 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:
+ 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:
+ 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:
+ 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(img<threshold)
+ idx = np.where(img < threshold)
if len(idx[0]) == 0:
pass
elif len(idx[0])>0:
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<s2idx[i])] += dt
+ brt[(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 xmin<np.min(x) or xmax>np.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