Skip to content
GitLab
Commit b8b5bed9 authored by Wei Chengliang's avatar Wei Chengliang
Browse files

update codestyle-PEP8

parent ab1521fe
Pipeline #7101 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
observation_sim/instruments/chip/effects.py
View file @ b8b5bed9
......@@ -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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment