Merge branch 'master' into 'release_v3.0'

Release version v3.2.0

See merge request !33

observation_sim/instruments/chip/Chip.py

@@ -191,7 +191,7 @@ class Chip(FocalPlane):
         """Return the filter index and type for a given chip #(chipID)
         """
         filter_type_list = _util.ALL_FILTERS
-        if chipID == None:
+        if chipID is None:
             chipID = self.chipID
 
         # updated configurations

observation_sim/instruments/chip/chip_utils.py

@@ -173,6 +173,13 @@ def get_flat(img, seed):
     return flat_img, flat_normal
 
 
+def get_innerflat(chip=None, filt=None):
+    from observation_sim.mock_objects import FlatLED
+    led_obj = FlatLED(chip, filt)
+    flat_img = led_obj.getInnerFlat()
+    return flat_img
+
+
 def add_cosmic_rays(img, chip, exptime=150, seed=0):
     cr_map, cr_event_num = effects.produceCR_Map(
         xLen=chip.npix_x, yLen=chip.npix_y,
@@ -206,7 +213,7 @@ def get_poisson(seed=0, sky_level=0.):
 
 
 def get_base_img(img, chip, read_noise, readout_time, dark_noise, exptime=150., InputDark=None):
-    if InputDark == None:
+    if InputDark is None:
         # base_level = read_noise**2 + dark_noise*(exptime+0.5*readout_time)
         # base_level = dark_noise*(exptime+0.5*readout_time)
         base_level = dark_noise*(exptime)
@@ -237,7 +244,7 @@ def add_poisson(img, chip, exptime=150., seed=0, sky_level=0., poisson_noise=Non
     img.addNoise(poisson_noise)
     # img -= read_noise**2
 
-    if InputDark != None:
+    if InputDark is not None:
         # "Instrument/data/dark/dark_1000s_example_0.fits"
         hdu = fits.open(InputDark)
         img += hdu[0].data/hdu[0].header['exptime']*exptime

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
 
@@ -22,7 +23,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 +38,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
@@ -51,27 +52,27 @@ def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=
     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))
-    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)
+    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
+    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 +83,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.5), 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,45 +92,45 @@ 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:
-        BiasLevel = Random16.reshape((nsecy,nsecx)) + bias_level
+    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, 
+                                     nsecy=2, nsecx=8,
                                      seed=int(seed),
                                      logger=logger)
     BiasCombImg = BiasSngImg*ncombine
@@ -139,7 +140,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 +148,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 +191,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)
@@ -223,7 +224,7 @@ def MakeFlatNcomb(flat_single_image, ncombine=1, read_noise=5, gain=1, overscan=
     if ncombine == 1:
         FlatTag = 'Single'
         pass
-    elif ncombine >1:
+    elif ncombine > 1:
         FlatCombImg /= ncombine
         FlatTag = 'Combine'
     # FlatCombImg.replaceNegative(replace_value=0)
@@ -232,7 +233,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()
@@ -245,13 +246,13 @@ 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:
         DarkTag = 'Single'
         pass
-    elif ncombine >1:
+    elif ncombine > 1:
         DarkCombImg /= ncombine
         DarkTag = 'Combine'
     # DarkCombImg.replaceNegative(replace_value=0)
@@ -261,27 +262,30 @@ 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
 
 
+def NonLinear_f(x, beta_1, beta_2):
+    return x - beta_1 * x * x + beta_2 * x * x * x
+
+
 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
+    # 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    ###############################
-
+# Saturation & Bleeding Start#
 def BleedingTrail(aa, yy):
-    if aa<0.2:
-        aa=0.2
+    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)                
+        faa = 0.5*(math.e+1/math.e)
         trail_frac = 1-0.1*(fy-1)/(faa-1)
     except Exception as e:
         print(e)
@@ -289,79 +293,96 @@ 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
+                yi -= 1
+                # [TEST] charge in the middle
+                if yi == (imgarr.shape[0] // 2 - 1):
                     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
+        elif direction == 'down':
+            if imgarr[yi+1, xi] >= fullwell:
+                imgarr[yi, xi] = fullwell
+                yi += 1
+                if yi == (imgarr.shape[0] // 2):
+                    break
                 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:
+        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:
+                    if yi < 0 or yi == (imgarr.shape[0]//2 - 1):
                         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]:
+                    if yi > imgarr.shape[0] or yi == (imgarr.shape[0]//2):
                         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
+                    if yi == (imgarr.shape[0]//2 - 1):
+                        break
+                elif direction == 'down':
+                    yi += 1
+                    if yi == (imgarr.shape[0]//2):
+                        break
                 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
+                    if yi == (imgarr.shape[0]//2 - 1):
+                        break
+                elif direction == 'down':
+                    yi += 1
+                    if yi == (imgarr.shape[0]//2):
+                        break
                 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:
+    if satupos_y.shape[0] == 0:
         # make no change for the image array
         return imgarr
     elif satupos_y.shape[0]/imgarr.size > 0.5:
@@ -371,28 +392,28 @@ 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
 
         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)
+            if 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)
+                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.
@@ -418,7 +439,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 +451,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,169 +484,163 @@ 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
 
 
 @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
 
     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:
+    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
                 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):
+            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
 
 
-
-#---------- For Cosmic-Ray Simulation ------------
-#---------- Zhang Xin ----------------------------
-
+# ---------- For Cosmic-Ray Simulation ------------
+# ---------- Zhang Xin ----------------------------
 def getYValue(collection, x):
-    index = 0;
+    index = 0
     if (collection.shape[1] == 2):
-        while(x>collection[index,0] and index < collection.shape[0]):
-            index= index + 1;
+        while (x > collection[index, 0] and index < collection.shape[0]):
+            index = index + 1
         if (index == collection.shape[0] or index == 0):
-            return 0;
+            return 0
 
-        deltX = collection[index, 0] - collection[index-1, 0];
-        deltY = collection[index, 1] - collection[index-1, 1];
+        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;
+            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):
+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;
+    pixelNum = int(xLen * yLen * cr_pixelRatio * normalRay)
+    pixelNum_n = int(xLen * yLen * cr_pixelRatio * nnormalRay)
+    CRPixelNum = 0
 
-    maxValue = max(attachedSizes[:,1])
-    maxValue += 0.1;
+    maxValue = max(attachedSizes[:, 1])
+    maxValue += 0.1
 
-    cr_event_num = 0;
-    CRs = np.zeros(pixelNum);
+    cr_event_num = 0
+    CRs = np.zeros(pixelNum)
     while (CRPixelNum < pixelNum):
-        x = CR_max_size * np.random.random();
-        y = maxValue * np.random.random();
+        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);
+            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);
+        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];
+    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;
+    sigma = 0.6 / 2.355
+    mean = 3.3
     random.seed(seed)
-    energys = np.zeros(cr_event_size);
+    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);
+        energy_index = random.normalvariate(mean, sigma)
+        energys[i] = pow(10, energy_index)
+    return energys
 
-    return energys;
 
-def convCR(CRmap=None, addPSF=None, sp_n = 4):
+def convCR(CRmap=None, addPSF=None, sp_n=4):
     sh = CRmap.shape
 
     # sp_n = 4
@@ -634,23 +649,23 @@ def convCR(CRmap=None, addPSF=None, sp_n = 4):
     for i in np.arange(sh[0]):
         i_st = sp_n*i
         for j in np.arange(sh[1]):
-            if CRmap[i,j] ==0:
+            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
 
     m_size = addPSF.shape[0]
 
-    subCRmap_n = np.zeros(np.array(subCRmap.shape) + m_size -1)
+    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
+            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
@@ -659,12 +674,12 @@ def convCR(CRmap=None, addPSF=None, sp_n = 4):
         i_st = sp_n*i
         for j in np.arange(sh_n[1]):
             p_v = 0
-            j_st=sp_n*j
+            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
+            CRmap_n[i, j] = p_v
 
     return CRmap_n
 
@@ -673,15 +688,15 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2
     # 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_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,seed = seed);
+    cr_event_size = cr_size.shape[0]
+    cr_energys = defineEnergyForCR(cr_event_size, seed=seed)
 
-    CRmap = np.zeros([yLen, xLen]);
+    CRmap = np.zeros([yLen, xLen])
 
-    ## produce conv kernel
+    # produce conv kernel
     from astropy.modeling.models import Gaussian2D
     o_size = 4
     sp_n = 8
@@ -694,28 +709,26 @@ 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):
-        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();
+        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);
+            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;
+            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
@@ -729,9 +742,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 +781,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 +796,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

observation_sim/instruments/chip/libBF/diffusion_X1.c

@@ -37,7 +37,7 @@ void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bi
   {
     printf("Adding BF effect...\n");
     //setup BF correlation fliter
-    float neX;
+    float neX, neXtemp;
     float neP1 = 50000;
     float bfaP1[9]={0.9707182, 0.002143905, 0.004131103, 0.001149542, 0.0005501739, 0.0005469659, 0.0003726081, 0.0003795207, 0.0001633302};
     float neP2 = 10000;
@@ -56,15 +56,18 @@ void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bi
         neX = arr_ima[j+i*ny];
         if(neX >= 10000)
         {
+          neXtemp = neX;
+          if(neXtemp > 100000)
+              neXtemp = 100000;
           bfa[0][0]=0; //linearInterp(neX, neP1, bfaP1[0], neP2, bfaP2[0]); //0;
-          bfa[0][1]=bfa[0][-1]=linearInterp(neX, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575;
-          bfa[-1][0]=bfa[1][0]=linearInterp(neX, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652;
-          bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neX, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335;
-          bfa[0][-2]=bfa[0][2]=linearInterp(neX, neP1, bfaP1[4], neP2, bfaP2[4]);
-          bfa[-2][0]=bfa[2][0]=linearInterp(neX, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118;
-          bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neX, neP1, bfaP1[6], neP2, bfaP2[6]);
-          bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neX, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083;
-          bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neX, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043;
+          bfa[0][1]=bfa[0][-1]=linearInterp(neXtemp, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575;
+          bfa[-1][0]=bfa[1][0]=linearInterp(neXtemp, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652;
+          bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neXtemp, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335;
+          bfa[0][-2]=bfa[0][2]=linearInterp(neXtemp, neP1, bfaP1[4], neP2, bfaP2[4]);
+          bfa[-2][0]=bfa[2][0]=linearInterp(neXtemp, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118;
+          bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neXtemp, neP1, bfaP1[6], neP2, bfaP2[6]);
+          bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neXtemp, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083;
+          bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neXtemp, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043;
         }
         else
         {

observation_sim/mock_objects/CatalogBase.py

@@ -12,6 +12,9 @@ class CatalogBase(metaclass=ABCMeta):
     def __init__(self):
         pass
 
+    def free_mem(self, **kward):
+        pass
+
     @abstractmethod
     def load_sed(self, obj, **kward):
         pass
@@ -29,8 +32,8 @@ class CatalogBase(metaclass=ABCMeta):
         param = {
             "star": -1,
             "id": -1,
-            "ra": 0,
-            "dec": 0.,
+            "ra": -999.,
+            "dec": -999.,
             "ra_orig": 0.,
             "dec_orig": 0.,
             "z": 0.,
@@ -44,15 +47,15 @@ class CatalogBase(metaclass=ABCMeta):
             "bfrac": 0.,
             "av": 0.,
             "redden": 0.,
-            "hlr_bulge": 0.,
-            "hlr_disk": 0.,
+            "hlr_bulge": -999.,
+            "hlr_disk": -999.,
             "ell_bulge": 0.,
             "ell_disk": 0.,
             "ell_tot": 0.,
-            "e1_disk": 0.,
-            "e2_disk": 0.,
-            "e1_bulge": 0.,
-            "e2_bulge": 0.,
+            "e1_disk": -999.,
+            "e2_disk": -999.,
+            "e1_bulge": -999.,
+            "e2_bulge": -999.,
             "teff": 0.,
             "logg": 0.,
             "feh": 0.,
@@ -91,6 +94,8 @@ class CatalogBase(metaclass=ABCMeta):
         bandpass = target_filt.bandpass_full
 
         if norm_filt is not None:
+            if hasattr(norm_filt, 'bandpass_full'):
+                norm_filt = Table(np.array(np.array([norm_filt.bandpass_full.wave_list*10.0, norm_filt.bandpass_full.func(norm_filt.bandpass_full.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY']))
             norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
         else:
             norm_filt = Table(

observation_sim/mock_objects/FlatLED.py

@@ -50,6 +50,8 @@ fluxLED = {'LED1': 15, 'LED2': 15, 'LED3': 12.5, 'LED4': 9, 'LED5': 9,
 #            'LED14': 10}
 
 mirro_eff = {'GU': 0.61, 'GV': 0.8, 'GI': 0.8}
+
+bandtoLed = {'NUV': ['LED1', 'LED2'], 'u': ['LED13', 'LED14'], 'g': ['LED3', 'LED4', 'LED5'], 'r': ['LED6', 'LED7'], 'i': ['LED8'], 'z': ['LED9', 'LED10'], 'y': ['LED10'], 'GU': ['LED1', 'LED2', 'LED13', 'LED14'], 'GV': ['LED3', 'LED4', 'LED5', 'LED6'], 'GI': ['LED7', 'LED8', 'LED9', 'LED10']}
 # mirro_eff = {'GU':1, 'GV':1, 'GI':1}
 
 
@@ -69,10 +71,18 @@ class FlatLED(MockObject):
                 with pkg_resources.path('observation_sim.mock_objects.data.led', "") as ledDir:
                     self.flatDir = ledDir.as_posix()
 
+    def getInnerFlat(self):
+        ledflats = bandtoLed[self.chip.filter_type]
+        iFlat = np.zeros([self.chip.npix_y, self.chip.npix_x])
+        for nled in ledflats:
+            iFlat = iFlat + self.getLEDImage1(led_type=nled, LED_Img_flag=False)
+        iFlat = iFlat/len(ledflats)
+        return iFlat
+
     ###
     # return LED flat, e/s
     ###
-    def getLEDImage(self, led_type='LED1'):
+    def getLEDImage(self, led_type='LED1', LED_Img_flag=True):
         # cwave = cwaves[led_type]
         flat = fits.open(os.path.join(self.flatDir, 'model_' +
                          cwaves_name[led_type] + 'nm.fits'))
@@ -102,7 +112,10 @@ class FlatLED(MockObject):
             N[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)]),
             method='linear')
         U = U/np.mean(U)
-        flatImage = U*fluxLED[led_type]*1000
+
+        flatImage = U
+        if LED_Img_flag:
+            flatImage = flatImage*fluxLED[led_type]*1000
         gc.collect()
         return flatImage
 
@@ -110,38 +123,41 @@ class FlatLED(MockObject):
     # return LED flat, e/s
     ###
 
-    def getLEDImage1(self, led_type='LED1'):
+    def getLEDImage1(self, led_type='LED1', LED_Img_flag=True):
         # cwave = cwaves[led_type]
         flat = fits.open(os.path.join(self.flatDir, 'model_' +
                          cwaves_name[led_type] + 'nm.fits'))
         xlen = flat[0].header['NAXIS1']
         ylen = 601
-
         i = self.chip.rowID - 1
         j = self.chip.colID - 1
-
         x = np.linspace(0, self.chip.npix_x, int(xlen/6.))
         y = np.linspace(0, self.chip.npix_y, int(ylen/5.))
         xx, yy = np.meshgrid(x, y)
-
         a1 = flat[0].data[int(ylen*i/5.):int(ylen*i/5.)+int(ylen/5.),
                           int(xlen*j/6.):int(xlen*j/6.)+int(xlen/6.)]
         # z = np.sin((xx+yy+xx**2+yy**2))
         # fInterp = interp2d(xx, yy, z, kind='linear')
-
         X_ = np.hstack((xx.flatten()[:, None], yy.flatten()[:, None]))
         Z_ = a1.flatten()
-
         n_x = np.arange(0, self.chip.npix_x, 1)
         n_y = np.arange(0, self.chip.npix_y, 1)
-
         M, N = np.meshgrid(n_x, n_y)
-
-        U = griddata(X_, Z_, (
-            M[0:self.chip.npix_y, 0:self.chip.npix_x],
-            N[0:self.chip.npix_y, 0:self.chip.npix_x]),
-            method='linear')
+        x_seg_len = 4
+        y_seg_len = 8
+        x_seg = int(self.chip.npix_x/x_seg_len)
+        y_seg = int(self.chip.npix_y/y_seg_len)
+        U = np.zeros([self.chip.npix_y, self.chip.npix_x], dtype=np.float32)
+        for y_seg_i in np.arange(y_seg_len):
+            for x_seg_i in np.arange(x_seg_len):
+                U[y_seg_i*y_seg:(y_seg_i+1)*y_seg, x_seg_i*x_seg:(x_seg_i+1)*x_seg] = griddata(X_, Z_, (M[y_seg_i*y_seg:(y_seg_i+1)*y_seg, x_seg_i*x_seg:(x_seg_i+1)*x_seg], N[y_seg_i*y_seg:(y_seg_i+1)*y_seg, x_seg_i*x_seg:(x_seg_i+1)*x_seg]), method='linear')
+        # U = griddata(X_, Z_, (
+        #     M[0:self.chip.npix_y, 0:self.chip.npix_x],
+        #     N[0:self.chip.npix_y, 0:self.chip.npix_x]),
+        #     method='nearest').astype(np.float32)
         U = U/np.mean(U)
+        flatImage = U
+        if LED_Img_flag:
             flatImage = U*fluxLED[led_type]*1000
         gc.collect()
         return flatImage

observation_sim/mock_objects/Galaxy.py

@@ -33,7 +33,7 @@ class Galaxy(MockObject):
             raise ValueError(
                 "!!!The number of PSF profiles and the number of bandpasses must be equal.")
         objs = []
-        if nphotons_tot == None:
+        if nphotons_tot is None:
             nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
         # print("nphotons_tot = ", nphotons_tot)
 
@@ -92,7 +92,7 @@ class Galaxy(MockObject):
         return final
 
     def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
-        if nphotons_tot == None:
+        if nphotons_tot is None:
             nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
         # print("nphotons_tot = ", nphotons_tot)
 
@@ -115,7 +115,7 @@ class Galaxy(MockObject):
 
         # Set Galsim Parameters
         if self.getMagFilter(filt) <= 15 and (not big_galaxy):
-            folding_threshold = 5.e-4
+            folding_threshold = 5.e-8
         else:
             folding_threshold = 5.e-3
         gsp = galsim.GSParams(folding_threshold=folding_threshold)
@@ -170,8 +170,11 @@ class Galaxy(MockObject):
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
 
             # Get PSF model
+            EXTRA = False
+            if self.getMagFilter(filt) <= filt.mag_saturation-1.:
+                EXTRA = True
             psf, pos_shear = psf_model.get_PSF(
-                chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
+                chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold, extrapolate=EXTRA)
 
             if self.bfrac == 0:
                 gal_temp = disk
@@ -325,23 +328,28 @@ class Galaxy(MockObject):
 
             galImg_List = []
             try:
-                pos_img_local = [0,0]
+                pos_img_local = [0, 0]
                 x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
                 y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
                 pos_img_local[0] = pos_img.x - x_start
                 pos_img_local[1] = pos_img.y - y_start
                 nnx = 0
                 nny = 0
-                for order in ["A","B"]:
+                for order in ["A", "B"]:
+                    EXTRA = False
+                    if self.getMagFilter(filt) <= filt.mag_saturation-2.:
+                        EXTRA = True
                     psf, pos_shear = psf_model.get_PSF(
-                        chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos)
+                        chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos, extrapolate=EXTRA, ngg=3072)
                     star_p = galsim.Convolve(psf, gal)
                     if nnx == 0:
-                        galImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
+                        galImg = star_p.drawImage(
+                            wcs=chip_wcs_local, offset=offset, method='no_pixel')
                         nnx = galImg.xmax - galImg.xmin + 1
                         nny = galImg.ymax - galImg.ymin + 1
                     else:
-                        galImg = star_p.drawImage(nx = nnx, ny = nny, wcs=chip_wcs_local, offset=offset)
+                        galImg = star_p.drawImage(
+                            nx=nnx, ny=nny, wcs=chip_wcs_local, offset=offset, method='no_pixel')
                     galImg.setOrigin(0, 0)
                     # n1 = np.sum(np.isinf(galImg.array))
                     # n2 = np.sum(np.isnan(galImg.array))
@@ -351,19 +359,23 @@ class Galaxy(MockObject):
                         # ERROR happens
                         return 2, pos_shear
                     galImg_List.append(galImg)
-                for order in ["C","D","E"]:
+                for order in ["C", "D", "E"]:
                     galImg_List.append(galImg)
             except:
-                psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
+                try:
+                    psf, pos_shear = psf_model.get_PSF(
+                        chip=chip, pos_img=pos_img)
                     star_p = galsim.Convolve(psf, gal)
-                galImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
+                    galImg = star_p.drawImage(
+                        wcs=chip_wcs_local, offset=offset)
                     galImg.setOrigin(0, 0)
                     if np.sum(np.isnan(galImg.array)) > 0:
                         # ERROR happens
                         return 2, pos_shear
-                for order in ["A","B","C","D","E"]:
+                    for order in ["A", "B", "C", "D", "E"]:
                         galImg_List.append(galImg)
-
+                except Exception as e:
+                    continue
             # starImg = gal.drawImage(
             #     wcs=chip_wcs_local, offset=offset, method='real_space')
 
@@ -378,7 +390,7 @@ class Galaxy(MockObject):
                 subSlitPos = int(grating_split_pos_chip - gal_origin[1])
                 # part img disperse
 
-                star_p1s=[]
+                star_p1s = []
                 for galImg in galImg_List:
 
                     subImg_p1 = galImg.array[:, 0:subSlitPos]
@@ -397,13 +409,14 @@ class Galaxy(MockObject):
                                        isAlongY=0,
                                        flat_cube=flat_cube)
 
-                self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+                self.addSLStoChipImage(
+                    sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
                 #                                           psf_model=psf_model, bandNo=i + 1,
                 #                                           grating_split_pos=grating_split_pos,
                 #                                           local_wcs=chip_wcs_local, pos_img=pos_img)
 
-                star_p2s=[]
+                star_p2s = []
                 for galImg in galImg_List:
 
                     subImg_p2 = galImg.array[:,
@@ -423,7 +436,8 @@ class Galaxy(MockObject):
                                        isAlongY=0,
                                        flat_cube=flat_cube)
 
-                self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+                self.addSLStoChipImage(
+                    sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
                 #                                           psf_model=psf_model, bandNo=i + 1,
                 #                                           grating_split_pos=grating_split_pos,
@@ -439,7 +453,8 @@ class Galaxy(MockObject):
                                     conf=chip.sls_conf[1],
                                     isAlongY=0,
                                     flat_cube=flat_cube)
-                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+                self.addSLStoChipImage(
+                    sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
                 #                                           psf_model=psf_model, bandNo=i + 1,
                 #                                           grating_split_pos=grating_split_pos,
@@ -453,7 +468,8 @@ class Galaxy(MockObject):
                                     conf=chip.sls_conf[0],
                                     isAlongY=0,
                                     flat_cube=flat_cube)
-                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
+                self.addSLStoChipImage(
+                    sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
                 #                                           psf_model=psf_model, bandNo=i + 1,
                 #                                           grating_split_pos=grating_split_pos,
@@ -465,7 +481,7 @@ class Galaxy(MockObject):
         return 1, pos_shear
 
     def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
-        if flux == None:
+        if flux is None:
             flux = self.getElectronFluxFilt(filt, tel, exptime)
         disk = galsim.Sersic(n=self.disk_sersic_idx,
                              half_light_radius=self.hlr_disk, flux=1.0)

observation_sim/mock_objects/MockObject.py

@@ -78,6 +78,8 @@ class MockObject(object):
                       (self.posImg.x, self.posImg.y), flush=True)
             self.posImg, self.fd_shear = fdmodel.get_distorted(
                 chip=chip, pos_img=self.posImg)
+            if self.posImg is None:
+                return None, None, None, None, None
             if verbose:
                 print("After field distortion:\n")
                 print("x = %.2f, y = %.2f\n" %
@@ -108,7 +110,7 @@ class MockObject(object):
 
     def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
                           exptime=150., fd_shear=None):
-        if nphotons_tot == None:
+        if nphotons_tot is None:
             nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
         # print("nphotons_tot = ", nphotons_tot)
 
@@ -122,7 +124,7 @@ class MockObject(object):
             return 2, None
         # Set Galsim Parameters
         if self.getMagFilter(filt) <= 15:
-            folding_threshold = 5.e-4
+            folding_threshold = 5.e-8
         else:
             folding_threshold = 5.e-3
         gsp = galsim.GSParams(folding_threshold=folding_threshold)
@@ -160,14 +162,17 @@ class MockObject(object):
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
 
             # Get PSF model
+            EXTRA = False
+            if self.getMagFilter(filt) <= filt.mag_saturation-1.:
+                EXTRA = True
             psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
-                                               folding_threshold=folding_threshold)
+                                               folding_threshold=folding_threshold, extrapolate=EXTRA)
             # star = galsim.DeltaFunction(gsparams=gsp)
             # star = star.withFlux(nphotons)
             # star = galsim.Convolve(psf, star)
             star = psf.withFlux(nphotons)
 
-            stamp = star.drawImage(wcs=chip_wcs_local, offset=offset)
+            stamp = star.drawImage(method='no_pixel', wcs=chip_wcs_local, offset=offset)
             if np.sum(np.isnan(stamp.array)) > 0:
                 continue
             stamp.setCenter(x_nominal, y_nominal)
@@ -242,7 +247,7 @@ class MockObject(object):
                     print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
                 #########################################################
                 # img_s, orig_off = convolveImg(img_s, psf_img_m)
-                orig_off = [0,0]
+                orig_off = [0, 0]
                 origin_order_x = v[1] - orig_off[0]
                 origin_order_y = v[2] - orig_off[1]
 
@@ -257,15 +262,18 @@ class MockObject(object):
                     continue
 
                 # orders = {'A': 'order1', 'B': 'order0', 'C': 'order2', 'D': 'order-1', 'E': 'order-2'}
-                orders = {'A': 'order1', 'B': 'order0', 'C': 'order0', 'D': 'order0', 'E': 'order0'}
+                orders = {'A': 'order1', 'B': 'order0',
+                          'C': 'order0', 'D': 'order0', 'E': 'order0'}
                 gratingN = chip_utils.getChipSLSGratingID(chip.chipID)[1]
                 if pos_img_local[0] < grating_split_pos:
                     gratingN = chip_utils.getChipSLSGratingID(chip.chipID)[0]
 
-
-                chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)].setOrigin(0, 0)
-                chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] = chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] + specImg[bounds]
-                chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)].setOrigin(chip.bound.xmin, chip.bound.ymin)
+                chip.img_stack[gratingN][orders[k]
+                                         ]['w' + str(bandNo)].setOrigin(0, 0)
+                chip.img_stack[gratingN][orders[k]]['w' + str(
+                    bandNo)][bounds] = chip.img_stack[gratingN][orders[k]]['w' + str(bandNo)][bounds] + specImg[bounds]
+                chip.img_stack[gratingN][orders[k]]['w' +
+                                                    str(bandNo)].setOrigin(chip.bound.xmin, chip.bound.ymin)
 
         else:
             for k, v in spec_orders.items():
@@ -330,9 +338,11 @@ class MockObject(object):
                 specImg.wcs = local_wcs
                 specImg.setOrigin(origin_order_x, origin_order_y)
                 try:
-                    specImg = psf_model.get_PSF_AND_convolve_withsubImg(chip, cutImg=specImg, pos_img_local=pos_img_local, bandNo=bandNo, g_order=k, grating_split_pos=grating_split_pos)
+                    specImg = psf_model.get_PSF_AND_convolve_withsubImg(
+                        chip, cutImg=specImg, pos_img_local=pos_img_local, bandNo=bandNo, g_order=k, grating_split_pos=grating_split_pos)
                 except:
-                    psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
+                    psf, pos_shear = psf_model.get_PSF(
+                        chip=chip, pos_img=pos_img)
 
                     psf_img = psf.drawImage(nx=100, ny=100, wcs=local_wcs)
 
@@ -376,7 +386,7 @@ class MockObject(object):
             sedNormFactor = 1.
 
         if self.getMagFilter(filt) <= 15:
-            folding_threshold = 5.e-4
+            folding_threshold = 5.e-8
         else:
             folding_threshold = 5.e-3
         gsp = galsim.GSParams(folding_threshold=folding_threshold)
@@ -425,45 +435,51 @@ class MockObject(object):
             # star = galsim.DeltaFunction(gsparams=gsp)
             # star = star.withFlux(tel.pupil_area * exptime)
 
-            #psf list :["A","B","C","D","E"]
+            # psf list :["A","B","C","D","E"]
             starImg_List = []
             try:
-                pos_img_local = [0,0]
+                pos_img_local = [0, 0]
                 x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
                 y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
                 pos_img_local[0] = pos_img.x - x_start
                 pos_img_local[1] = pos_img.y - y_start
                 nnx = 0
                 nny = 0
-                for order in ["A","B"]:
+                for order in ["A", "B"]:
+                    EXTRA = False
+                    if self.getMagFilter(filt) <= filt.mag_saturation-2.:
+                        EXTRA = True
+                        nnx = 2048
+                        nny = 2048
                     psf, pos_shear = psf_model.get_PSF(
-                        chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos)
+                        chip, pos_img_local=pos_img_local, bandNo=i+1, galsimGSObject=True, g_order=order, grating_split_pos=grating_split_pos, extrapolate=EXTRA, ngg=3072)
                     # star_p = galsim.Convolve(psf, star)
                     star_p = psf.withFlux(tel.pupil_area * exptime)
                     if nnx == 0:
-                        starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
+                        starImg = star_p.drawImage(
+                            wcs=chip_wcs_local, offset=offset, method='no_pixel')
                         nnx = starImg.xmax - starImg.xmin + 1
                         nny = starImg.ymax - starImg.ymin + 1
                     else:
-                        starImg = star_p.drawImage(nx = nnx, ny = nny, wcs=chip_wcs_local, offset=offset)
+                        starImg = star_p.drawImage(
+                            nx=nnx, ny=nny, wcs=chip_wcs_local, offset=offset, method='no_pixel')
                     # n1 = np.sum(np.isinf(starImg.array))
                     # n2 = np.sum(np.isnan(starImg.array))
                     # if n1>0 or n2 > 0:
                     #     print("DEBUG: MockObject, inf:%d, nan:%d"%(n1, n2))
                     starImg.setOrigin(0, 0)
                     starImg_List.append(starImg)
-                for order in ["C","D","E"]:
+                for order in ["C", "D", "E"]:
                     starImg_List.append(starImg)
             except:
                 psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
                 # star_p = galsim.Convolve(psf, star)
                 star_p = psf.withFlux(tel.pupil_area * exptime)
-                starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset)
+                starImg = star_p.drawImage(wcs=chip_wcs_local, offset=offset, method='no_pixel')
                 starImg.setOrigin(0, 0)
-                for order in ["A","B","C","D","E"]:
+                for order in ["A", "B", "C", "D", "E"]:
                     starImg_List.append(starImg)
 
-            
             # psf_tmp = galsim.Gaussian(sigma=0.002)
             # star = galsim.Convolve(psf_tmp, star)
 
@@ -479,7 +495,7 @@ class MockObject(object):
             if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
                 subSlitPos = int(grating_split_pos_chip - gal_origin[1])
                 # part img disperse
-                star_p1s=[]
+                star_p1s = []
                 for starImg in starImg_List:
 
                     subImg_p1 = starImg.array[:, 0:subSlitPos]
@@ -498,13 +514,13 @@ class MockObject(object):
                                        isAlongY=0,
                                        flat_cube=flat_cube)
 
-                self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+                self.addSLStoChipImage(
+                    sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
                 #                                           psf_model=psf_model, bandNo=i+1, grating_split_pos=grating_split_pos,
                 #                                           local_wcs=chip_wcs_local, pos_img=pos_img)
 
-
-                star_p2s=[]
+                star_p2s = []
                 for starImg in starImg_List:
 
                     subImg_p2 = starImg.array[:,
@@ -524,7 +540,8 @@ class MockObject(object):
                                        isAlongY=0,
                                        flat_cube=flat_cube)
 
-                self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+                self.addSLStoChipImage(
+                    sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
                 #                                           psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
                 #                                           local_wcs=chip_wcs_local, pos_img=pos_img)
@@ -539,7 +556,8 @@ class MockObject(object):
                                     conf=chip.sls_conf[1],
                                     isAlongY=0,
                                     flat_cube=flat_cube)
-                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+                self.addSLStoChipImage(
+                    sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
                 #                                           psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
                 #                                           local_wcs=chip_wcs_local, pos_img=pos_img)
@@ -552,7 +570,8 @@ class MockObject(object):
                                     conf=chip.sls_conf[0],
                                     isAlongY=0,
                                     flat_cube=flat_cube)
-                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
+                self.addSLStoChipImage(
+                    sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 # pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
                 #                                           psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos,
                 #                                           local_wcs=chip_wcs_local, pos_img=pos_img)
@@ -572,7 +591,7 @@ class MockObject(object):
 
     def drawObj_PSF(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
                     exptime=150., fd_shear=None, chip_output=None):
-        if nphotons_tot == None:
+        if nphotons_tot is None:
             nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
         # print("nphotons_tot = ", nphotons_tot)
 
@@ -641,7 +660,7 @@ class MockObject(object):
         os.makedirs(fn, exist_ok=True)
         fn = fn + "/ccd_{:}".format(chip.chipID) + \
             "_psf_"+str(self.param['id'])+".fits"
-        if fn != None:
+        if fn is not None:
             if os.path.exists(fn):
                 os.remove(fn)
         hdu = fitsio.PrimaryHDU()

observation_sim/mock_objects/Quasar.py

@@ -98,7 +98,7 @@ class Quasar(MockObject):
         del self.sed
 
     def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
-        if flux == None:
+        if flux is None:
             flux = self.getElectronFluxFilt(filt, tel, exptime)
         qso = galsim.Gaussian(sigma=1.e-8, flux=1.)
         qso = qso.withFlux(flux)

observation_sim/mock_objects/SpecDisperser/SpecDisperser.py

@@ -50,8 +50,21 @@ def rotate90(array_orig=None, xc=0, yc=0, isClockwise=0):
 
 
 class SpecDisperser(object):
-    def __init__(self, orig_img=None, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=None, band_start=6200,
-                 band_end=10000, isAlongY=0, conf='../param/CONF/csst.conf', gid=0, flat_cube=None, ignoreBeam=[]):
+    def __init__(
+        self,
+        orig_img=None,
+        xcenter=0,
+        ycenter=0,
+        origin=[100, 100],
+        tar_spec=None,
+        band_start=6200,
+        band_end=10000,
+        isAlongY=0,
+        conf="../param/CONF/csst.conf",
+        gid=0,
+        flat_cube=None,
+        ignoreBeam=[],
+    ):
         """
         orig_img: normal image,galsim image
         xcenter, ycenter: the center of galaxy in orig_img
@@ -66,23 +79,21 @@ class SpecDisperser(object):
         # self.img_y = orig_img.shape[0]
 
         # 5 orders, A, B ,
-        orderName=["A","B","C","D","E"]
+        orderName = ["A", "B", "C", "D", "E"]
         self.orig_img_orders = OrderedDict()
         if isinstance(orig_img, list):
             orig_img_list = orig_img
             list_len = len(orig_img_list)
             if list_len < 5:
-                for i in np.arange(5-list_len):
-                    orig_img_list.append(orig_img_list[list_len-1])
+                for i in np.arange(5 - list_len):
+                    orig_img_list.append(orig_img_list[list_len - 1])
             for i, k in enumerate(orig_img_list):
                 self.orig_img_orders[orderName[i]] = k
 
-        
         if isinstance(orig_img, galsim.Image):
             for i in np.arange(5):
                 self.orig_img_orders[orderName[i]] = orig_img
 
-
         orig_img_one = self.orig_img_orders["A"]
 
         self.thumb_img = np.abs(orig_img_one.array)
@@ -99,8 +110,12 @@ class SpecDisperser(object):
         self.flat_cube = flat_cube
         if self.isAlongY == 1:
             for order in orderName:
-                self.orig_img_orders[order], self.thumb_x, self.thumb_y = rotate90(array_orig=self.orig_img_orders[order], xc=orig_img_one.center.x,
-                                                                  yc=orig_img_one.center.y, isClockwise=1)
+                self.orig_img_orders[order], self.thumb_x, self.thumb_y = rotate90(
+                    array_orig=self.orig_img_orders[order],
+                    xc=orig_img_one.center.x,
+                    yc=orig_img_one.center.y,
+                    isClockwise=1,
+                )
             # self.thumb_img, self.thumb_x, self.thumb_y = rotate90(array_orig=self.thumb_img, xc=orig_img_one.center.x,
             #                                                       yc=orig_img_one.center.y, isClockwise=1)
 
@@ -139,6 +154,7 @@ class SpecDisperser(object):
 
         from .disperse_c import interp
         from .disperse_c import disperse
+
         # from MockObject.disperse_c import disperse
         self.thumb_img = np.abs(self.orig_img_orders[beam].array)
         self.thumb_x = self.orig_img_orders[beam].center.x
@@ -146,11 +162,12 @@ class SpecDisperser(object):
         self.img_sh = self.orig_img_orders[beam].array.shape
         dx = self.grating_conf.dxlam[beam]
         xoff = 0
-        ytrace_beam, lam_beam = self.grating_conf.get_beam_trace(x=self.xcenter, y=self.ycenter, dx=(dx + xoff),
-                                                                 beam=beam)
+        ytrace_beam, lam_beam = self.grating_conf.get_beam_trace(
+            x=self.xcenter, y=self.ycenter, dx=(dx + xoff), beam=beam
+        )
 
         # Account for pixel centering of the trace
-        yfrac_beam = ytrace_beam - floor(ytrace_beam+0.5)
+        yfrac_beam = ytrace_beam - floor(ytrace_beam + 0.5)
 
         ysens = lam_beam * 0
         lam_index = argsort(lam_beam)
@@ -160,8 +177,8 @@ class SpecDisperser(object):
             (self.band_end - self.band_start) / 0.1))
 
         thri = interpolate.interp1d(
-            conf_sens['WAVELENGTH'], conf_sens['SENSITIVITY'])
-        spci = interpolate.interp1d(self.spec['WAVELENGTH'], self.spec['FLUX'])
+            conf_sens["WAVELENGTH"], conf_sens["SENSITIVITY"])
+        spci = interpolate.interp1d(self.spec["WAVELENGTH"], self.spec["FLUX"])
 
         beam_thr = thri(lam_intep)
         spec_sample = spci(lam_intep)
@@ -179,8 +196,9 @@ class SpecDisperser(object):
         #
         # self.writerSensitivityFile(conffile = self.grating_conf_file,  beam = beam, w = lam_beam[lam_index], sens = ysensitivity[lam_index])
 
-        ysens[lam_index] = interp.interp_conserve_c(lam_beam[lam_index], lam_intep, bean_thr_spec, integrate=1, left=0,
-                                                    right=0)
+        ysens[lam_index] = interp.interp_conserve_c(
+            lam_beam[lam_index], lam_intep, bean_thr_spec, integrate=1, left=0, right=0
+        )
 
         sensitivity_beam = ysens
 
@@ -196,7 +214,7 @@ class SpecDisperser(object):
 
         dxpix = dx - dx[0] + x0[1]
 
-        dyc = cast[int](np.floor(ytrace_beam+0.5))
+        dyc = cast[int](np.floor(ytrace_beam + 0.5))
 
         # dypix = cast[int](np.floor(ytrace_beam - dyc[0] + x0[0] + 0.5))
         dypix = dyc - dyc[0] + x0[0]
@@ -204,7 +222,7 @@ class SpecDisperser(object):
         frac_ids = yfrac_beam < 0
 
         dypix[frac_ids] = dypix[frac_ids] - 1
-        yfrac_beam[frac_ids] = 1+yfrac_beam[frac_ids]
+        yfrac_beam[frac_ids] = 1 + yfrac_beam[frac_ids]
 
         flat_index = idx[dypix, dxpix]
 
@@ -233,7 +251,7 @@ class SpecDisperser(object):
 
             sub_flat_cube = zeros(
                 [len(self.flat_cube), beam_sh[0], beam_sh[1]])
-            sub_flat_cube[0] = sub_flat_cube[0] + 1.
+            sub_flat_cube[0] = sub_flat_cube[0] + 1.0
 
             overlap_flag = 1
 
@@ -260,33 +278,38 @@ class SpecDisperser(object):
                     overlap_flag = 0
 
             if overlap_flag == 1:
-                sub_flat_cube[:, beam_y_s-originOut_y:beam_y_e-originOut_y+1, beam_x_s-originOut_x:beam_x_e -
-                              originOut_x+1] = self.flat_cube[:, beam_y_s:beam_y_e+1, beam_x_s:beam_x_e+1]
+                sub_flat_cube[
+                    :,
+                    beam_y_s - originOut_y: beam_y_e - originOut_y + 1,
+                    beam_x_s - originOut_x: beam_x_e - originOut_x + 1,
+                ] = self.flat_cube[:, beam_y_s: beam_y_e + 1, beam_x_s: beam_x_e + 1]
 
             for i in arange(0, len(self.flat_cube), 1):
                 beam_flat[:, i] = sub_flat_cube[i].flatten()
-#            beam_flat = zeros([len(modelf), len(self.flat_cube)])
-#            flat_sh = self.flat_cube[0].shape
-#            for i in arange(0, beam_sh[0], 1):
-#                for j in arange(0, beam_sh[1], 1):
-#                    k = i * beam_sh[1] + j
-#                    if originOut_y + i >= flat_sh[0] or originOut_y + i < 0 or originOut_x + j>= flat_sh[1] or originOut_x+j < 0:
-#                        temp_bf = np.zeros_like(self.flat_cube[:, 0, 0])
-#                        temp_bf[0] = 1.0
-#                        beam_flat[k] = temp_bf
-#                    else:
-#                        beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j]
-
-        status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32),
+        #            beam_flat = zeros([len(modelf), len(self.flat_cube)])
+        #            flat_sh = self.flat_cube[0].shape
+        #            for i in arange(0, beam_sh[0], 1):
+        #                for j in arange(0, beam_sh[1], 1):
+        #                    k = i * beam_sh[1] + j
+        #                    if originOut_y + i >= flat_sh[0] or originOut_y + i < 0 or originOut_x + j>= flat_sh[1] or originOut_x+j < 0:
+        #                        temp_bf = np.zeros_like(self.flat_cube[:, 0, 0])
+        #                        temp_bf[0] = 1.0
+        #                        beam_flat[k] = temp_bf
+        #                    else:
+        #                        beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j]
+
+        status = disperse.disperse_grism_object(
+            self.thumb_img.astype(np.float32),
             flat_index[nonz],
             yfrac_beam[nonz],
             sensitivity_beam[nonz],
-                                                modelf, x0,
-                                                array(self.img_sh,
-                                                      dtype=int64),
+            modelf,
+            x0,
+            array(self.img_sh, dtype=int64),
             array(beam_sh, dtype=int64),
             beam_flat,
-                                                lam_beam[lam_index][nonz])
+            lam_beam[lam_index][nonz],
+        )
 
         model = modelf.reshape(beam_sh)
         # n1 = np.sum(np.isinf(model))
@@ -309,14 +332,14 @@ class SpecDisperser(object):
             model, _, _ = rotate90(array_orig=model, isClockwise=0)
         return model, originOut_x, originOut_y, dxpix, dypix, lam_beam, ysens
 
-    def writerSensitivityFile(self, conffile='', beam='', w=None, sens=None):
-        orders = {'A': '1st', 'B': '0st', 'C': '2st', 'D': '-1st', 'E': '-2st'}
+    def writerSensitivityFile(self, conffile="", beam="", w=None, sens=None):
+        orders = {"A": "1st", "B": "0st", "C": "2st", "D": "-1st", "E": "-2st"}
         sens_file_name = conffile[0:-5] + \
-            '_sensitivity_' + orders[beam] + '.fits'
-        if not os.path.exists(sens_file_name) == True:
+            "_sensitivity_" + orders[beam] + ".fits"
+        if os.path.exists(sens_file_name) is False:
             senstivity_out = Table(
-                array([w, sens]).T, names=('WAVELENGTH', 'SENSITIVITY'))
-            senstivity_out.write(sens_file_name, format='fits')
+                array([w, sens]).T, names=("WAVELENGTH", "SENSITIVITY"))
+            senstivity_out.write(sens_file_name, format="fits")
 
 
 """
@@ -324,8 +347,8 @@ Demonstrate aXe trace polynomials.
 """
 
 
-class aXeConf():
-    def __init__(self, conf_file='WFC3.IR.G141.V2.5.conf'):
+class aXeConf:
+    def __init__(self, conf_file="WFC3.IR.G141.V2.5.conf"):
         """Read an aXe-compatible configuration file
 
         Parameters
@@ -340,17 +363,17 @@ class aXeConf():
             self.count_beam_orders()
 
             # Global XOFF/YOFF offsets
-            if 'XOFF' in self.conf.keys():
-                self.xoff = np.float(self.conf['XOFF'])
+            if "XOFF" in self.conf.keys():
+                self.xoff = np.float(self.conf["XOFF"])
             else:
-                self.xoff = 0.
+                self.xoff = 0.0
 
-            if 'YOFF' in self.conf.keys():
-                self.yoff = np.float(self.conf['YOFF'])
+            if "YOFF" in self.conf.keys():
+                self.yoff = np.float(self.conf["YOFF"])
             else:
-                self.yoff = 0.
+                self.yoff = 0.0
 
-    def read_conf_file(self, conf_file='WFC3.IR.G141.V2.5.conf'):
+    def read_conf_file(self, conf_file="WFC3.IR.G141.V2.5.conf"):
         """Read an aXe config file, convert floats and arrays
 
         Parameters
@@ -366,11 +389,11 @@ class aXeConf():
         lines = open(conf_file).readlines()
         for line in lines:
             # empty / commented lines
-            if (line.startswith('#')) | (line.strip() == '') | ('"' in line):
+            if (line.startswith("#")) | (line.strip() == "") | ('"' in line):
                 continue
 
             # split the line, taking out ; and # comments
-            spl = line.split(';')[0].split('#')[0].split()
+            spl = line.split(";")[0].split("#")[0].split()
             param = spl[0]
             if len(spl) > 2:
                 value = np.cast[float](spl[1:])
@@ -385,22 +408,20 @@ class aXeConf():
         return conf
 
     def count_beam_orders(self):
-        """Get the maximum polynomial order in DYDX or DLDP for each beam
-        """
+        """Get the maximum polynomial order in DYDX or DLDP for each beam"""
         self.orders = {}
-        for beam in ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J']:
+        for beam in ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J"]:
             order = 0
-            while 'DYDX_{0:s}_{1:d}'.format(beam, order) in self.conf.keys():
+            while "DYDX_{0:s}_{1:d}".format(beam, order) in self.conf.keys():
                 order += 1
 
-            while 'DLDP_{0:s}_{1:d}'.format(beam, order) in self.conf.keys():
+            while "DLDP_{0:s}_{1:d}".format(beam, order) in self.conf.keys():
                 order += 1
 
             self.orders[beam] = order - 1
 
     def get_beams(self):
-        """Get beam parameters and read sensitivity curves
-        """
+        """Get beam parameters and read sensitivity curves"""
         import os
         from collections import OrderedDict
         from astropy.table import Table, Column
@@ -413,12 +434,15 @@ class aXeConf():
         for beam in self.orders:
             if self.orders[beam] > 0:
                 self.beams.append(beam)
-                self.dxlam[beam] = np.arange(self.conf['BEAM{0}'.format(beam)].min(),
-                                             self.conf['BEAM{0}'.format(beam)].max(), dtype=int)
+                self.dxlam[beam] = np.arange(
+                    self.conf["BEAM{0}".format(beam)].min(), self.conf["BEAM{0}".format(beam)].max(), dtype=int
+                )
                 self.nx[beam] = int(self.dxlam[beam].max() -
                                     self.dxlam[beam].min()) + 1
                 self.sens[beam] = Table.read(
-                    '{0}/{1}'.format(os.path.dirname(self.conf_file), self.conf['SENSITIVITY_{0}'.format(beam)]))
+                    "{0}/{1}".format(os.path.dirname(self.conf_file),
+                                     self.conf["SENSITIVITY_{0}".format(beam)])
+                )
                 # self.sens[beam].wave = np.cast[np.double](self.sens[beam]['WAVELENGTH'])
                 # self.sens[beam].sens = np.cast[np.double](self.sens[beam]['SENSITIVITY'])
 
@@ -471,41 +495,41 @@ class aXeConf():
         return a
 
     def evaluate_dp(self, dx, dydx):
-        """Evalate arc length along the trace given trace polynomial coefficients
+        # """Evalate arc length along the trace given trace polynomial coefficients
 
-        Parameters
-        ----------
-        dx : array-like
-            x pixel to evaluate
+        # Parameters
+        # ----------
+        # dx : array-like
+        #     x pixel to evaluate
 
-        dydx : array-like
-            Coefficients of the trace polynomial
+        # dydx : array-like
+        #     Coefficients of the trace polynomial
 
-        Returns
-        -------
-        dp : array-like
-            Arc length along the trace at position `dx`.
+        # Returns
+        # -------
+        # dp : array-like
+        #     Arc length along the trace at position `dx`.
 
-        For `dydx` polynomial orders 0, 1 or 2, integrate analytically.
-        Higher orders must be integrated numerically.
+        # For `dydx` polynomial orders 0, 1 or 2, integrate analytically.
+        # Higher orders must be integrated numerically.
 
-        **Constant:**
-            .. math:: dp = dx
+        # **Constant:**
+        #     .. math:: dp = dx
 
-        **Linear:**
-            .. math:: dp = \sqrt{1+\mathrm{DYDX}[1]}\cdot dx
+        # **Linear:**
+        #     .. math:: dp = \sqrt{1+\mathrm{DYDX}[1]}\cdot dx
 
-        **Quadratic:**
-            .. math:: u = \mathrm{DYDX}[1] + 2\ \mathrm{DYDX}[2]\cdot dx
+        # **Quadratic:**
+        #     .. math:: u = \mathrm{DYDX}[1] + 2\ \mathrm{DYDX}[2]\cdot dx
 
-            .. math:: dp = (u \sqrt{1+u^2} + \mathrm{arcsinh}\ u) / (4\cdot \mathrm{DYDX}[2])
+        #     .. math:: dp = (u \sqrt{1+u^2} + \mathrm{arcsinh}\ u) / (4\cdot \mathrm{DYDX}[2])
 
-        """
+        # """
         # dp is the arc length along the trace
         # $\lambda = dldp_0 + dldp_1 dp + dldp_2 dp^2$ ...
 
         poly_order = len(dydx) - 1
-        if (poly_order == 2):
+        if poly_order == 2:
             if np.abs(np.unique(dydx[2])).max() == 0:
                 poly_order = 1
 
@@ -515,10 +539,9 @@ class aXeConf():
             dp = np.sqrt(1 + dydx[1] ** 2) * (dx)
         elif poly_order == 2:  # quadratic trace
             u0 = dydx[1] + 2 * dydx[2] * (0)
-            dp0 = (u0 * np.sqrt(1 + u0 ** 2) + np.arcsinh(u0)) / (4 * dydx[2])
+            dp0 = (u0 * np.sqrt(1 + u0**2) + np.arcsinh(u0)) / (4 * dydx[2])
             u = dydx[1] + 2 * dydx[2] * (dx)
-            dp = (u * np.sqrt(1 + u ** 2) + np.arcsinh(u)) / \
-                (4 * dydx[2]) - dp0
+            dp = (u * np.sqrt(1 + u**2) + np.arcsinh(u)) / (4 * dydx[2]) - dp0
         else:
             # high order shape, numerical integration along trace
             # (this can be slow)
@@ -530,7 +553,7 @@ class aXeConf():
                 dyfull += i * dydx[i] * (xfull - 0.5) ** (i - 1)
 
             # Integrate from 0 to dx / -dx
-            dpfull = xfull * 0.
+            dpfull = xfull * 0.0
             lt0 = xfull < 0
             if lt0.sum() > 1:
                 dpfull[lt0] = np.cumsum(
@@ -548,7 +571,7 @@ class aXeConf():
 
         return dp
 
-    def get_beam_trace(self, x=507, y=507, dx=0., beam='A'):
+    def get_beam_trace(self, x=507, y=507, dx=0.0, beam="A"):
         """Get an aXe beam trace for an input reference pixel and list of output x pixels `dx`
 
         Parameters
@@ -579,17 +602,17 @@ class aXeConf():
 
         xi, yi = x - self.xoff, y - self.yoff
         xoff_beam = self.field_dependent(
-            xi, yi, self.conf['XOFF_{0}'.format(beam)])
+            xi, yi, self.conf["XOFF_{0}".format(beam)])
         yoff_beam = self.field_dependent(
-            xi, yi, self.conf['YOFF_{0}'.format(beam)])
+            xi, yi, self.conf["YOFF_{0}".format(beam)])
 
         # y offset of trace (DYDX)
         dydx = np.zeros(NORDER)  # 0 #+1.e-80
         dydx = [0] * NORDER
 
         for i in range(NORDER):
-            if 'DYDX_{0:s}_{1:d}'.format(beam, i) in self.conf.keys():
-                coeffs = self.conf['DYDX_{0:s}_{1:d}'.format(beam, i)]
+            if "DYDX_{0:s}_{1:d}".format(beam, i) in self.conf.keys():
+                coeffs = self.conf["DYDX_{0:s}_{1:d}".format(beam, i)]
                 dydx[i] = self.field_dependent(xi, yi, coeffs)
 
         # $dy = dydx_0+dydx_1 dx+dydx_2 dx^2+$ ...
@@ -603,14 +626,20 @@ class aXeConf():
         dldp = [0] * NORDER
 
         for i in range(NORDER):
-            if 'DLDP_{0:s}_{1:d}'.format(beam, i) in self.conf.keys():
-                coeffs = self.conf['DLDP_{0:s}_{1:d}'.format(beam, i)]
+            if "DLDP_{0:s}_{1:d}".format(beam, i) in self.conf.keys():
+                coeffs = self.conf["DLDP_{0:s}_{1:d}".format(beam, i)]
                 dldp[i] = self.field_dependent(xi, yi, coeffs)
 
-        self.eval_input = {'x': x, 'y': y, 'beam': beam, 'dx': dx}
-        self.eval_output = {'xi': xi, 'yi': yi, 'dldp': dldp, 'dydx': dydx,
-                            'xoff_beam': xoff_beam, 'yoff_beam': yoff_beam,
-                            'dy': dy}
+        self.eval_input = {"x": x, "y": y, "beam": beam, "dx": dx}
+        self.eval_output = {
+            "xi": xi,
+            "yi": yi,
+            "dldp": dldp,
+            "dydx": dydx,
+            "xoff_beam": xoff_beam,
+            "yoff_beam": yoff_beam,
+            "dy": dy,
+        }
 
         dp = self.evaluate_dp(dx - xoff_beam, dydx)
         # ## dp is the arc length along the trace
@@ -648,13 +677,13 @@ class aXeConf():
         #     dp = np.interp(dx-xoff_beam, xfull, dpfull)
 
         # Evaluate dldp
-        lam = dp * 0.
+        lam = dp * 0.0
         for i in range(NORDER):
-            lam += dldp[i] * dp ** i
+            lam += dldp[i] * dp**i
 
         return dy, lam
 
-    def show_beams(self, beams=['E', 'D', 'C', 'B', 'A']):
+    def show_beams(self, beams=["E", "D", "C", "B", "A"]):
         """
         Make a demo plot of the beams of a given configuration file
         """
@@ -663,45 +692,47 @@ class aXeConf():
         x0, x1 = 507, 507
         dx = np.arange(-800, 1200)
 
-        if 'WFC3.UV' in self.conf_file:
+        if "WFC3.UV" in self.conf_file:
             x0, x1 = 2073, 250
             dx = np.arange(-1200, 1200)
-        if 'G800L' in self.conf_file:
+        if "G800L" in self.conf_file:
             x0, x1 = 2124, 1024
             dx = np.arange(-1200, 1200)
 
         s = 200  # marker size
         fig = plt.figure(figsize=[10, 3])
-        plt.scatter(0, 0, marker='s', s=s, color='black', edgecolor='0.8',
-                    label='Direct')
+        plt.scatter(0, 0, marker="s", s=s, color="black",
+                    edgecolor="0.8", label="Direct")
 
         for beam in beams:
-            if 'XOFF_{0}'.format(beam) not in self.conf.keys():
+            if "XOFF_{0}".format(beam) not in self.conf.keys():
                 continue
 
             xoff = self.field_dependent(
-                x0, x1, self.conf['XOFF_{0}'.format(beam)])
+                x0, x1, self.conf["XOFF_{0}".format(beam)])
             dy, lam = self.get_beam_trace(x0, x1, dx=dx, beam=beam)
-            xlim = self.conf['BEAM{0}'.format(beam)]
+            xlim = self.conf["BEAM{0}".format(beam)]
             ok = (dx >= xlim[0]) & (dx <= xlim[1])
-            plt.scatter(dx[ok] + xoff, dy[ok], c=lam[ok] / 1.e4, marker='s', s=s,
-                        alpha=0.5, edgecolor='None')
-            plt.text(np.median(dx[ok]), np.median(dy[ok]) + 1, beam,
-                     ha='center', va='center', fontsize=14)
-            print('Beam {0}, lambda=({1:.1f} - {2:.1f})'.format(beam,
+            plt.scatter(dx[ok] + xoff, dy[ok], c=lam[ok] / 1.0e4,
+                        marker="s", s=s, alpha=0.5, edgecolor="None")
+            plt.text(np.median(dx[ok]), np.median(
+                dy[ok]) + 1, beam, ha="center", va="center", fontsize=14)
+            print("Beam {0}, lambda=({1:.1f} - {2:.1f})".format(beam,
                   lam[ok].min(), lam[ok].max()))
 
         plt.grid()
-        plt.xlabel(r'$\Delta x$')
-        plt.ylabel(r'$\Delta y$')
+        plt.xlabel(r"$\Delta x$")
+        plt.ylabel(r"$\Delta y$")
 
         cb = plt.colorbar(pad=0.01, fraction=0.05)
-        cb.set_label(r'$\lambda\,(\mu\mathrm{m})$')
+        cb.set_label(r"$\lambda\,(\mu\mathrm{m})$")
         plt.title(self.conf_file)
         plt.tight_layout()
-        plt.savefig('{0}.pdf'.format(self.conf_file))
+        plt.savefig("{0}.pdf".format(self.conf_file))
 
     # def load_grism_config(conf_file):
+
+
 #     """Load parameters from an aXe configuration file
 
 #     Parameters

observation_sim/mock_objects/SpecDisperser/setup_c.py

@@ -8,16 +8,16 @@ import numpy
 
 extensions = [
     Extension("disperse_c.interp", ["disperse_c/interp.pyx"],
-        include_dirs = [numpy.get_include()],
+              include_dirs=[numpy.get_include()],
               libraries=["m"]),
 
     Extension("disperse_c.disperse", ["disperse_c/disperse.pyx"],
-        include_dirs = [numpy.get_include()],
+              include_dirs=[numpy.get_include()],
               libraries=["m"]),
 ]
 
 
 setup(
-    name = "slssim_disperse",
-    ext_modules = cythonize(extensions),
+    name="slssim_disperse",
+    ext_modules=cythonize(extensions),
 )

observation_sim/mock_objects/Stamp.py

@@ -21,7 +21,7 @@ class Stamp(MockObject):
         del self.sed
 
     def drawObj_multiband(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
-        if nphotons_tot == None:
+        if nphotons_tot is None:
             nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
 
         try:
@@ -92,7 +92,7 @@ class Stamp(MockObject):
             else:
                 gal = gal + gal_temp
 
-        stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset)
+        stamp = gal.drawImage(method='no_pixel', wcs=chip_wcs_local, offset=offset)
         if np.sum(np.isnan(stamp.array)) > 0:
             # ERROR happens
             return 2, pos_shear

observation_sim/mock_objects/Star.py

@@ -22,7 +22,7 @@ class Star(MockObject):
         del self.sed
 
     def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
-        if flux == None:
+        if flux is None:
             flux = self.getElectronFluxFilt(filt, tel, exptime)
         # star = galsim.Gaussian(sigma=1.e-8, flux=1.)
         star = galsim.DeltaFunction()
@@ -35,7 +35,7 @@ class Star(MockObject):
             raise ValueError(
                 "!!!The number of PSF profiles and the number of bandpasses must be equal.")
         objs = []
-        if nphotons_tot == None:
+        if nphotons_tot is None:
             nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
 
         try:

observation_sim/psf/FieldDistortion.py

@@ -62,7 +62,8 @@ class FieldDistortion(object):
                                 the distored position.
         """
         if not self.isContainObj_FD(chip=chip, pos_img=pos_img):
-            return galsim.PositionD(-1, -1), None
+            # return galsim.PositionD(-1, -1), None
+            return None, None
         if not img_rot:
             img_rot = np.radians(self.img_rot)
         else:

observation_sim/psf/PSFGauss.py

@@ -50,7 +50,7 @@ class PSFGauss(PSFModel):
 
         Return: the flux ratio
         """
-        if pscale == None:
+        if pscale is None:
             pscale = self.pix_size
         gaussx = galsim.Gaussian(flux=1.0, sigma=sig)
         gaussImg = gaussx.drawImage(scale=pscale, method='no_pixel')
@@ -68,7 +68,7 @@ class PSFGauss(PSFModel):
 
         return the fwhm in arcsec
         """
-        if pscale == None:
+        if pscale is None:
             pscale = self.pix_size
         err = 1.0e-3
         nxx = 100