new sim: modified ways to add detector effects: e.g. brighter-fatter etc.

ObservationSim/Instrument/Chip/ChipUtils.py

+import os
+import galsim
+import ctypes
+import numpy as np
+from astropy.io import fits
+from datetime import datetime
+
+from ObservationSim.Instrument.Chip import Effects as effects
+from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
+
+try:
+    import importlib.resources as pkg_resources
+except ImportError:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+
+def log_info(msg, logger=None):
+    if logger:
+        logger.info(msg)
+    else:
+        print(msg, flush=True)
+
+def getChipSLSConf(chipID):
+    confFile = ''
+    if chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf']
+    if chipID == 2: confFile = ['CSST_GV4.conf', 'CSST_GV3.conf']
+    if chipID == 3: confFile = ['CSST_GU2.conf', 'CSST_GU1.conf']
+    if chipID == 4: confFile = ['CSST_GU4.conf', 'CSST_GU3.conf']
+    if chipID == 5: confFile = ['CSST_GV2.conf', 'CSST_GV1.conf']
+    if chipID == 10: confFile = ['CSST_GI4.conf', 'CSST_GI3.conf']
+    if chipID == 21: confFile = ['CSST_GI6.conf', 'CSST_GI5.conf']
+    if chipID == 26: confFile = ['CSST_GV8.conf', 'CSST_GV7.conf']
+    if chipID == 27: confFile = ['CSST_GU6.conf', 'CSST_GU5.conf']
+    if chipID == 28: confFile = ['CSST_GU8.conf', 'CSST_GU7.conf']
+    if chipID == 29: confFile = ['CSST_GV6.conf', 'CSST_GV5.conf']
+    if chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
+    return confFile
+
+def generateHeader(chip, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., timestamp = 1621915200):
+    datetime_obs = datetime.utcfromtimestamp(timestamp)
+    date_obs = datetime_obs.strftime("%y%m%d")
+    time_obs = datetime_obs.strftime("%H%M%S")
+    h_prim = generatePrimaryHeader(
+        xlen=chip.npix_x, 
+        ylen=chip.npix_y, 
+        pointNum = str(pointing_ID),
+        ra=ra_cen, 
+        dec=dec_cen, 
+        pixel_scale=chip.pix_scale,
+        date=date_obs,
+        time_obs=time_obs,
+        im_type = im_type,
+        exptime=exptime,
+        chip_name=str(chip.chipID).rjust(2, '0')
+        )
+    h_ext = generateExtensionHeader(
+        chip=chip,
+        xlen=chip.npix_x, 
+        ylen=chip.npix_y, 
+        ra=ra_cen, 
+        dec=dec_cen,
+        pa=img_rot.deg, 
+        gain=chip.gain, 
+        readout=chip.read_noise, 
+        dark=chip.dark_noise, 
+        saturation=90000, 
+        pixel_scale=chip.pix_scale, 
+        pixel_size=chip.pix_size,
+        xcen=chip.x_cen,
+        ycen=chip.y_cen,
+        extName='SCI',
+        timestamp = timestamp,
+        exptime = exptime,
+        readoutTime = chip.readout_time)
+    return h_prim, h_ext
+
+def outputCal(chip, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., timestamp = 1621915200):
+    h_prim, h_ext = generateHeader(
+        chip=chip,
+        ra_cen=ra_cen,
+        dec_cen=dec_cen,
+        img_rot=img_rot,
+        im_type=im_type,
+        pointing_ID=pointing_ID,
+        exptime=exptime,
+        timestamp=timestamp)
+    hdu1 = fits.PrimaryHDU(header=h_prim)
+    hdu1.add_checksum()
+    hdu1.header.comments['CHECKSUM'] = 'HDU checksum'
+    hdu1.header.comments['DATASUM'] = 'data unit checksum'
+    hdu2 = fits.ImageHDU(img.array, header=h_ext)
+    hdu2.add_checksum()
+    hdu2.header.comments['XTENSION'] = 'extension type'
+    hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
+    hdu2.header.comments['DATASUM'] = 'data unit checksum'
+    hdu1 = fits.HDUList([hdu1, hdu2])
+    fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
+    hdu1.writeto(fname, output_verify='ignore', overwrite=True)
+
+def add_sky_background(img, filt, exptime, sky_map=None, tel=None):
+    # Add sky background
+    if sky_map is None:
+        sky_map = filt.getSkyNoise(exptime=exptime)
+        sky_map = sky_map * np.ones_like(img.array)
+        sky_map = galsim.Image(array=sky_map)
+        # Apply Poisson noise to the sky map
+        # # (NOTE): only for photometric chips if it utilizes the photon shooting to draw stamps
+        # if self.survey_type == "photometric":
+        #     sky_map.addNoise(poisson_noise)
+    elif img.array.shape != sky_map.shape:
+        raise ValueError("The shape img and sky_map must be equal.")
+    elif tel is not None: # If sky_map is given in flux
+        sky_map = sky_map * tel.pupil_area * exptime
+    img += sky_map
+    return img, sky_map
+
+def get_flat(img, seed):
+    flat_img = effects.MakeFlatSmooth(
+                GSBounds=img.bounds, 
+                seed=seed)
+    flat_normal = flat_img / np.mean(flat_img.array)
+    return flat_img, flat_normal
+
+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, 
+        exTime=exptime+0.5*chip.readout_time, 
+        cr_pixelRatio=0.003*(exptime+0.5*chip.readout_time)/600.,
+        gain=chip.gain, 
+        attachedSizes=chip.attachedSizes,
+        seed=seed)   # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
+    img += cr_map
+    cr_map[cr_map > 65535] = 65535
+    cr_map[cr_map < 0] = 0
+    crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
+    del cr_map
+    return img, crmap_gsimg, cr_event_num
+
+def add_PRNU(img, chip, seed=0):
+    prnu_img = effects.PRNU_Img(
+        xsize=chip.npix_x, 
+        ysize=chip.npix_y, 
+        sigma=0.01, 
+        seed=seed)
+    img *= prnu_img
+    return img, prnu_img
+
+def get_poisson(seed=0, sky_level=0.):
+    rng_poisson = galsim.BaseDeviate(seed)
+    poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=sky_level)
+    return rng_poisson, poisson_noise
+
+def get_base_img(img, read_noise, readout_time, dark_noise, exptime=150.):
+    base_level = read_noise**2 + dark_noise*(exptime+0.5*readout_time)
+    base_img = base_level * np.ones_like(img.array)
+    return base_img
+
+def add_poisson(img, chip, exptime=150., seed=0, sky_level=0., poisson_noise=None, dark_noise=None):
+    if poisson_noise is None:
+        _, poisson_noise = get_poisson(seed=seed, sky_level=sky_level)
+    read_noise = chip.read_noise
+    if dark_noise is None:
+        dark_noise = chip.dark_noise
+    base_img = get_base_img(img=img, read_noise=read_noise, readout_time=chip.readout_time, dark_noise=dark_noise, exptime=exptime)
+    img += base_img
+    img.addNoise(poisson_noise)
+    img -= read_noise**2
+    return img, base_img
+
+def add_brighter_fatter(img):
+    #Inital dynamic lib
+    try:
+        with pkg_resources.files('ObservationSim.Instrument.Chip.lib_bf').joinpath("libmoduleBF.so") as lib_path:
+            lib_bf = ctypes.CDLL(lib_path)
+    except AttributeError:
+        with pkg_resources.path('ObservationSim.Instrument.Chip.lib_bf', "libmoduleBF.so") as lib_path:
+            lib_bf = ctypes.CDLL(lib_path)
+    lib_bf.addEffects.argtypes = [ctypes.c_int, ctypes.c_int, ctypes.POINTER(ctypes.c_float), ctypes.POINTER(ctypes.c_float), ctypes.c_int]
+    
+    # Set bit flag
+    bit_flag = 1
+    bit_flag = bit_flag | (1 << 2)
+
+    nx, ny = img.array.shape
+    nn = nx * ny
+    arr_ima= (ctypes.c_float*nn)()
+    arr_imc= (ctypes.c_float*nn)()
+
+    arr_ima[:]= img.array.reshape(nn)
+    arr_imc[:]= np.zeros(nn)
+
+    lib_bf.addEffects(nx, ny, arr_ima, arr_imc, bit_flag)
+    img.array[:, :] = np.array(arr_imc[:]).reshape([nx, ny])
+    del arr_ima, arr_imc
+    return img
\ No newline at end of file

ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c

+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "nrutil.h"
+
+#define ISSETBITFLAG(x,b) ((x) & (1 << (b)))
+#define ADD_DIFFUSION    1
+#define ADD_BF_FILTER    2
+
+float linearInterp(float xp, float x0, float y0, float x1, float y1);
+
+void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bit_flag)
+{
+  int nx, ny, i,j,k,ks;
+  int it,jt,itt,jtt;
+  int diffuidx[26][2],diffuN,ilow,ih,im,dim[3];
+  float diffua[5][5],cdiffu[26],**bfa;
+  double mvar,mcov,tmp,ma,mb,mc;
+  char fname[100];
+ 
+  nx = ngx_ima;  //input-image size
+  ny = ngy_ima; 
+
+  //0. init. original image with an input array (arr_ima)
+  //1. Adding diffusion effect.
+  if(ISSETBITFLAG(bit_flag, ADD_DIFFUSION)) 
+  {
+    printf("adding diffusion.....\n");
+    printf("ERR: no diffusion filter ...");
+    exit(0);
+  }
+
+
+  //2. Adding BF effect
+  if(ISSETBITFLAG(bit_flag, ADD_BF_FILTER)) 
+  {
+    printf("Adding BF effect...\n");
+    //setup BF correlation fliter
+    float neX;
+    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;
+    float bfaP2[9]={0.9945288, 0.0003041936, 0.0007539311, 0.0002424675, 0.0001226098, 0.00009308617, 0.00008027447, 0.00006309676, 0.00006400052};
+
+    bfa=matrix(-2,2,-2,2);
+
+    // smooth with the BF filter
+    for(i=0;i<nx;i++)for(j=0;j<ny;j++) arr_imc[j+i*ny]=0;
+
+    for(i=0;i<nx;i++)
+    {
+      for(j=0;j<ny;j++)
+      {
+        //rescale BF filter with the local pix value
+        neX = arr_ima[j+i*ny];
+        if(neX >= 10000)
+        {
+          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;
+        }
+        else
+        {
+          neX=10000;
+          bfa[0][0]=0;
+          bfa[0][1]=bfa[0][-1]=bfaP2[1];
+          bfa[-1][0]=bfa[1][0]=bfaP2[2];
+          bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=bfaP2[3];
+          bfa[0][-2]=bfa[0][2]=bfaP2[4];
+          bfa[-2][0]=bfa[2][0]=bfaP2[5];
+          bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=bfaP2[6];
+          bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=bfaP2[7];
+          bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=bfaP2[8];
+        }
+        tmp = 0;
+        for(it=-2;it<=2;it++)
+          for(jt=-2;jt<=2;jt++)
+          {
+            bfa[it][jt] = bfa[it][jt]/neX*arr_ima[j+i*ny];
+            tmp += bfa[it][jt];
+          }
+        bfa[0][0]=1.-tmp;
+
+        // assign electrons according to the BF filter bfat
+        for(it=-2;it<=2;it++)
+        {
+          for(jt=-2;jt<=2;jt++)
+          {
+            itt=i+it;
+            jtt=j+jt;
+            if(itt>=0 && jtt>=0 && itt<nx && jtt<ny)
+              //c0[itt][jtt]+=bfa[it][jt]*b[i][j];
+              arr_imc[jtt+itt*ny] += bfa[it][jt]*arr_ima[j+i*ny];
+          }
+        }
+      }
+    }
+    free_matrix(bfa,-2,2,-2,2);
+  }
+  else
+  {
+    for(i=0;i<nx;i++) for(j=0;j<ny;j++) arr_imc[j+i*ny]=arr_ima[j+i*ny]; ////for ADD_BF False
+  }
+}
+
+
+float linearInterp(float xp, float x0, float y0, float x1, float y1)
+{
+  float yp;
+  yp = y0 + ((y1-y0)/(x1-x0)) * (xp - x0);
+  return yp;
+}
+

ObservationSim/Instrument/Chip/lib_bf/nrutil.h

+/* CAUTION: This is the ANSI C (only) version of the Numerical Recipes
+   utility file nrutil.h.  Do not confuse this file with the same-named
+   file nrutil.h that may be supplied in a 'misc' subdirectory.
+   *That* file is the one from the book, and contains both ANSI and
+   traditional K&R versions, along with #ifdef macros to select the
+   correct version.  *This* file contains only ANSI C.               */
+
+#ifndef _NR_UTILS_H_
+#define _NR_UTILS_H_
+
+static float sqrarg;
+#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
+
+static double dsqrarg;
+#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
+
+static double dmaxarg1,dmaxarg2;
+#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
+        (dmaxarg1) : (dmaxarg2))
+
+static double dminarg1,dminarg2;
+#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
+        (dminarg1) : (dminarg2))
+
+static float maxarg1,maxarg2;
+#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
+        (maxarg1) : (maxarg2))
+
+static float minarg1,minarg2;
+#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
+        (minarg1) : (minarg2))
+
+static long lmaxarg1,lmaxarg2;
+#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
+        (lmaxarg1) : (lmaxarg2))
+
+static long lminarg1,lminarg2;
+#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
+        (lminarg1) : (lminarg2))
+
+static int imaxarg1,imaxarg2;
+#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
+        (imaxarg1) : (imaxarg2))
+
+static int iminarg1,iminarg2;
+#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
+        (iminarg1) : (iminarg2))
+
+#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
+
+void nrerror(char error_text[]);
+float *vector(long nl, long nh);
+int *ivector(long nl, long nh);
+unsigned char *cvector(long nl, long nh);
+long *lvector(long nl, long nh);
+double *dvector(long nl, long nh);
+float **matrix(long nrl, long nrh, long ncl, long nch);
+double **dmatrix(long nrl, long nrh, long ncl, long nch);
+int **imatrix(long nrl, long nrh, long ncl, long nch);
+float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+	long newrl, long newcl);
+float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
+float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+void free_vector(float *v, long nl, long nh);
+void free_ivector(int *v, long nl, long nh);
+void free_cvector(unsigned char *v, long nl, long nh);
+void free_lvector(long *v, long nl, long nh);
+void free_dvector(double *v, long nl, long nh);
+void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
+void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
+void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
+void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
+void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
+void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+	long ndl, long ndh);
+
+
+int  ***i3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+void free_i3tensor(int ***t, long nrl, long nrh, long ncl, long nch,     
+	long ndl, long ndh);
+
+unsigned char  ***b3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+void free_b3tensor(unsigned char ***t, long nrl, long nrh, long ncl, long nch,     
+	long ndl, long ndh);
+
+double ***d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
+	long ndl, long ndh);
+
+
+char **cmatrix(long nrl, long nrh, long ncl, long nch);
+void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch);
+
+#endif /* _NR_UTILS_H_ */

ObservationSim/MockObject/Galaxy.py

@@ -14,18 +14,6 @@ from ObservationSim.MockObject.MockObject import MockObject
 class Galaxy(MockObject):
     def __init__(self, param, logger=None):
         super().__init__(param, logger=logger)
-        # self.thetaR = self.param["theta"]
-        # self.bfrac = self.param["bfrac"]
-        # self.hlr_disk = self.param["hlr_disk"]
-        # self.hlr_bulge = self.param["hlr_bulge"]
-
-        # Extract ellipticity components
-        # self.e_disk = galsim.Shear(g=self.param["ell_disk"], beta=self.thetaR*galsim.degrees)
-        # self.e_bulge = galsim.Shear(g=self.param["ell_bulge"], beta=self.thetaR*galsim.degrees)
-        # self.e_total = galsim.Shear(g=self.param["ell_tot"], beta=self.thetaR*galsim.degrees)
-        # self.e1_disk, self.e2_disk = self.e_disk.g1, self.e_disk.g2
-        # self.e1_bulge, self.e2_bulge = self.e_bulge.g1, self.e_bulge.g2
-        # self.e1_total, self.e2_total = self.e_total.g1, self.e_total.g2
 
         if not hasattr(self, "disk_sersic_idx"):
             self.disk_sersic_idx = 1.
@@ -101,10 +89,6 @@ class Galaxy(MockObject):
             if self.logger:
                 self.logger.error(e)
             return 2, None
-
-        nphotons_sum = 0
-        photons_list = []
-        xmax, ymax = 0, 0
         
         # # [C6 TEST]
         # print('hlr_disk = %.4f, hlr_bulge = %.4f'%(self.hlr_disk, self.hlr_bulge))
@@ -114,7 +98,7 @@ class Galaxy(MockObject):
         if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy
             big_galaxy = True
 
-        # (TEST) Galsim Parameters
+        # Set Galsim Parameters
         if self.getMagFilter(filt) <= 15 and (not big_galaxy):
             folding_threshold = 5.e-4
         else:
@@ -122,7 +106,7 @@ class Galaxy(MockObject):
         gsp = galsim.GSParams(folding_threshold=folding_threshold)
 
         self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
-                                        img_real_wcs=self.real_wcs)
+                                        img_real_wcs=self.chip_wcs)
 
         x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
         x_nominal = int(np.floor(x + 0.5))
@@ -130,9 +114,11 @@ class Galaxy(MockObject):
         dx = x - x_nominal
         dy = y - y_nominal
         offset = galsim.PositionD(dx, dy)
+        # Get real local wcs of object (deal with chip rotation w.r.t its center)
+        chip_wcs_local = self.chip_wcs.local(self.real_pos)
+        is_updated = 0
 
-        real_wcs_local = self.real_wcs.local(self.real_pos)
-
+        # Model the galaxy as disk + bulge
         disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
         disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
         disk = disk.shear(disk_shape)
@@ -140,33 +126,33 @@ class Galaxy(MockObject):
         bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
         bulge = bulge.shear(bulge_shape)
         
+        # Get shear and deal with shear induced by field distortion
         if fd_shear:
             g1 += fd_shear.g1
             g2 += fd_shear.g2
         gal_shear = galsim.Shear(g1=g1, g2=g2)
 
+        # Loop over all sub-bandpasses
         for i in range(len(bandpass_list)):
             bandpass = bandpass_list[i]
-
             try:
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
                 if self.logger:
                     self.logger.error(e)
-                # return False
                 continue
-            
             ratio = sub/full
             if not (ratio == -1 or (ratio != ratio)):
                 nphotons = ratio * nphotons_tot
             else:
                 continue
-            nphotons_sum += nphotons
-
+            
+            # nphotons_sum += nphotons
             # # [C6 TEST]
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
 
+            # Get PSF model
             psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
             
             gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk
@@ -194,39 +180,21 @@ class Galaxy(MockObject):
         # for stat in top_stats[:10]:
         #     print(stat)
 
-        stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
-        photons = stamp.photons
-        photons.x += x_nominal
-        photons.y += y_nominal
-        photons_list.append(photons)
-
-        stamp.wcs = real_wcs_local
+        # stamp = gal.drawImage(wcs=chip_wcs_local, method='phot', offset=offset, save_photons=True)
+        stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset)
+        if np.sum(np.isnan(stamp.array)) > 0:
+            # ERROR happens
+            return 2, pos_shear
         stamp.setCenter(x_nominal, y_nominal)
         bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
-
         if bounds.area() > 0:
             chip.img.setOrigin(0, 0)
-            stamp[bounds] = chip.img[bounds]
-
-            if not big_galaxy:
-                for i in range(len(photons_list)):
-                    if i == 0:
-                        chip.sensor.accumulate(photons_list[i], stamp)
-                    else:
-                        chip.sensor.accumulate(photons_list[i], stamp, resume=True)
-            else:
-                sensor = galsim.Sensor()
-                for i in range(len(photons_list)):
-                    if i == 0:
-                        sensor.accumulate(photons_list[i], stamp)
-                    else:
-                        sensor.accumulate(photons_list[i], stamp, resume=True)
-                del sensor
-
-            chip.img[bounds] = stamp[bounds]
-
+            chip.img[bounds] += stamp[bounds]
+            is_updated = 1
             chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
-        else:
+            del stamp
+        
+        if is_updated == 0:
             # Return code 0: object photons missed this detector
             print("obj %s missed"%(self.id))
             if self.logger:
@@ -235,8 +203,6 @@ class Galaxy(MockObject):
         
         # # [C6 TEST]
         # print("nphotons_sum = ", nphotons_sum)
-        del photons_list
-        del stamp
         return 1, pos_shear
 
     def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
@@ -255,7 +221,7 @@ class Galaxy(MockObject):
                           names=('WAVELENGTH', 'FLUX'))
 
         self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
-                                        img_real_wcs=self.real_wcs)
+                                        img_real_wcs=self.chip_wcs)
 
         x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
         x_nominal = int(np.floor(x + 0.5))
@@ -264,7 +230,7 @@ class Galaxy(MockObject):
         dy = y - y_nominal
         offset = galsim.PositionD(dx, dy)
 
-        real_wcs_local = self.real_wcs.local(self.real_pos)
+        chip_wcs_local = self.chip_wcs.local(self.real_pos)
 
 
         big_galaxy = False
@@ -313,7 +279,7 @@ class Galaxy(MockObject):
                 # if fd_shear is not None:
                 #     gal = gal.shear(fd_shear)
 
-            starImg = gal.drawImage(wcs=real_wcs_local, offset=offset)
+            starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset)
 
             origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
                            x_nominal - (starImg.center.x - starImg.xmin)]
@@ -340,7 +306,7 @@ class Galaxy(MockObject):
                                     isAlongY=0,
                                     flat_cube=flat_cube)
 
-                self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
+                self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
 
                 subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
                 star_p2 = galsim.Image(subImg_p2)
@@ -357,7 +323,7 @@ class Galaxy(MockObject):
                                        isAlongY=0,
                                        flat_cube=flat_cube)
 
-                self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
+                self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
 
                 del sdp_p1
                 del sdp_p2
@@ -369,7 +335,7 @@ class Galaxy(MockObject):
                                     conf=chip.sls_conf[1],
                                     isAlongY=0,
                                     flat_cube=flat_cube)
-                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
+                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
                 del sdp
             elif grating_split_pos_chip>=gal_end[1]:
                 sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
@@ -379,7 +345,7 @@ class Galaxy(MockObject):
                                     conf=chip.sls_conf[0],
                                     isAlongY=0,
                                     flat_cube=flat_cube)
-                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
+                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
                 del sdp
 
             # print(self.y_nominal, starImg.center.y, starImg.ymin)
@@ -404,46 +370,6 @@ class Galaxy(MockObject):
         final = galsim.Convolve(psf, gal)
         return final
 
-    def drawObject(self, img, final, noise_level=0.0, flux=None, filt=None, tel=None, exptime=150.):
-        """ Override the method in parent class 
-        Need to constrain the size of image stamp for extended objects
-        """
-        isUpdated = True
-        if flux == None:
-            flux = self.getElectronFluxFilt(filt, tel, exptime)
-        stamp = final.drawImage(wcs=self.localWCS, offset=self.offset)
-        stamp_arr = stamp.array
-        mask = (stamp_arr >= 0.001*noise_level) # why 0.001?
-        err = int(np.sqrt(mask.sum()))
-        if np.mod(err, 2) == 1:
-            err += 1
-        # if err == 1:
-        if err == 0:
-            subSize = 16 # why 16?
-        else:
-            subSize = max([err, 16])
-            fluxRatio = flux / stamp_arr[mask].sum()
-            final = final.withScaledFlux(fluxRatio)
-
-        imgSub = galsim.ImageF(subSize, subSize)
-
-        # Draw with FFT
-        # stamp = final.drawImage(image=imgSub, wcs=self.localWCS, offset=self.offset)
-
-        # Draw with Photon Shoot
-        stamp = final.drawImage(image=imgSub, wcs=self.localWCS, method='phot', offset=self.offset)
-        
-        stamp.setCenter(self.x_nominal, self.y_nominal)
-        if np.sum(np.isnan(stamp.array)) >= 1:
-            stamp.setZero()
-        bounds = stamp.bounds & img.bounds
-        if bounds.area() == 0:
-            isUpdated = False
-        else:
-            img[bounds] += stamp[bounds]
-        return img, stamp, isUpdated
-
-
     def getObservedEll(self, g1=0, g2=0):
         e1_obs, e2_obs, e_obs, theta = eObs(self.e1_total, self.e2_total, g1, g2)
         return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs

ObservationSim/MockObject/MockObject.py

@@ -13,7 +13,6 @@ from ObservationSim.MockObject.SpecDisperser import SpecDisperser
 class MockObject(object):
     def __init__(self, param, logger=None):
         self.param = param
-        
         for key in self.param:
             setattr(self, key, self.param[key])
         
@@ -27,14 +26,10 @@ class MockObject(object):
         elif self.param["star"] == 3:
             self.type = "stamp"
         ###mock_stamp_END
-
         self.sed = None
-
+        self.fd_shear = None
         # Place holder for outputs
         self.additional_output_str = ""
-        
-        self.fd_shear = None
-
         self.logger = logger
 
     def getMagFilter(self, filt):
@@ -65,6 +60,7 @@ class MockObject(object):
     def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, chip_wcs=None, img_header=None):
         self.posImg = img.wcs.toImage(self.getPosWorld())
         self.localWCS = img.wcs.local(self.posImg)
+        # Apply field distortion model
         if (fdmodel is not None) and (chip is not None):
             if verbose:
                 print("\n")
@@ -74,6 +70,7 @@ class MockObject(object):
             if verbose:
                 print("After field distortion:\n")
                 print("x = %.2f, y = %.2f\n" % (self.posImg.x, self.posImg.y), flush=True)
+        
         x, y = self.posImg.x + 0.5, self.posImg.y + 0.5
         self.x_nominal = int(np.floor(x + 0.5))
         self.y_nominal = int(np.floor(y + 0.5))
@@ -81,50 +78,22 @@ class MockObject(object):
         dy = y - self.y_nominal
         self.offset = galsim.PositionD(dx, dy)
         
+        # Deal with chip rotation
         if chip_wcs is not None:
-            self.real_wcs = chip_wcs
+            self.chip_wcs = chip_wcs
         elif img_header is not None:
-            self.real_wcs = galsim.FitsWCS(header=img_header)
+            self.chip_wcs = galsim.FitsWCS(header=img_header)
         else:
-            self.real_wcs = None
+            self.chip_wcs = None
 
-        return self.posImg, self.offset, self.localWCS, self.real_wcs, self.fd_shear
+        return self.posImg, self.offset, self.localWCS, self.chip_wcs, self.fd_shear
 
     def getRealPos(self, img, global_x=0., global_y=0., img_real_wcs=None):
         img_global_pos = galsim.PositionD(global_x, global_y)
         cel_pos = img.wcs.toWorld(img_global_pos)
         realPos = img_real_wcs.toImage(cel_pos)
-
         return realPos
 
-    def drawObject(self, img, final, flux=None, filt=None, tel=None, exptime=150.):
-        """ Draw (point like) object on img.
-        Should be overided for extended source, e.g. galaxy...
-        Paramter:
-            img: the "canvas"
-            final: final (after shear, PSF etc.) GSObject
-        Return:
-            img: the image with the GSObject added (or discarded)
-            isUpdated: is the "canvas" been updated? (a flag for updating statistcs)
-        """
-        isUpdated = True
-
-        # Draw with FFT
-        # stamp = final.drawImage(wcs=self.localWCS, offset=self.offset)
-
-        # Draw with Photon Shoot
-        stamp = final.drawImage(wcs=self.localWCS, method='phot', offset=self.offset)
-
-        stamp.setCenter(self.x_nominal, self.y_nominal)
-        if np.sum(np.isnan(stamp.array)) >= 1:
-            stamp.setZero()
-        bounds = stamp.bounds & img.bounds
-        if bounds.area() == 0:
-            isUpdated = False
-        else:
-            img[bounds] += stamp[bounds]
-        return img, stamp, isUpdated
-
     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:
@@ -139,29 +108,27 @@ class MockObject(object):
                 self.logger.error(e)
             return 2, None
 
-        nphotons_sum = 0
-        photons_list = []
-        xmax, ymax = 0, 0
-
-        # (TEST) Galsim Parameters
+        # Set Galsim Parameters
         if self.getMagFilter(filt) <= 15:
             folding_threshold = 5.e-4
         else:
             folding_threshold = 5.e-3
         gsp = galsim.GSParams(folding_threshold=folding_threshold)
-
+        
+        # Get real image position of object (deal with chip rotation w.r.t its center)
         self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
-                                        img_real_wcs=self.real_wcs)
-
+                                        img_real_wcs=self.chip_wcs)
         x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
         x_nominal = int(np.floor(x + 0.5))
         y_nominal = int(np.floor(y + 0.5))
         dx = x - x_nominal
         dy = y - y_nominal
         offset = galsim.PositionD(dx, dy)
+        # Get real local wcs of object (deal with chip rotation w.r.t its center)
+        chip_wcs_local = self.chip_wcs.local(self.real_pos)
+        is_updated = 0
 
-        real_wcs_local = self.real_wcs.local(self.real_pos)
-
+        # Loop over all sub-bandpasses
         for i in range(len(bandpass_list)):
             bandpass = bandpass_list[i]
             try:
@@ -170,61 +137,41 @@ class MockObject(object):
                 print(e)
                 if self.logger:
                     self.logger.error(e)
-                # return False
                 continue
-
             ratio = sub / full
-
             if not (ratio == -1 or (ratio != ratio)):
                 nphotons = ratio * nphotons_tot
             else:
-                # return False
                 continue
-            nphotons_sum += nphotons
+
+            # nphotons_sum += nphotons
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
+            
+            # Get PSF model
             psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
                                                folding_threshold=folding_threshold)
             star = galsim.DeltaFunction(gsparams=gsp)
             star = star.withFlux(nphotons)
             star = galsim.Convolve(psf, star)
 
-            stamp = star.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
-            xmax = max(xmax, stamp.xmax)
-            ymax = max(ymax, stamp.ymax)
-            photons = stamp.photons
-            photons.x += x_nominal
-            photons.y += y_nominal
-            photons_list.append(photons)
-
-        stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1))
-        stamp.wcs = real_wcs_local
-        stamp.setCenter(x_nominal, y_nominal)
-        bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
+            stamp = star.drawImage(wcs=chip_wcs_local, offset=offset)
+            if np.sum(np.isnan(stamp.array)) > 0:
+                continue
+            stamp.setCenter(x_nominal, y_nominal)
+            bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
+            if bounds.area() > 0:
+                chip.img.setOrigin(0, 0)
+                chip.img[bounds] += stamp[bounds]
+                is_updated = 1
+                chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
+                del stamp
         
-        # # (DEBUG)
-        # print("stamp bounds: ", stamp.bounds)
-        # print(bounds)
-        if bounds.area() > 0:
-            chip.img.setOrigin(0, 0)
-            stamp[bounds] = chip.img[bounds]
-            for i in range(len(photons_list)):
-                if i == 0:
-                    chip.sensor.accumulate(photons_list[i], stamp)
-                else:
-                    chip.sensor.accumulate(photons_list[i], stamp, resume=True)
-
-            chip.img[bounds] = stamp[bounds]
-
-            chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
-        else:
-            # Return code 0: object photons missed this detector
+        if is_updated == 0:
+            # Return code 0: object has missed this detector
             print("obj %s missed"%(self.id))
             if self.logger:
                 self.logger.info("obj %s missed"%(self.id))
             return 0, pos_shear
-
-        del photons_list
-        del stamp
         return 1, pos_shear # Return code 1: draw sucesss
 
     def addSLStoChipImage(self, sdp=None, chip=None, xOrderSigPlus=None, local_wcs=None):
@@ -299,7 +246,7 @@ class MockObject(object):
                           names=('WAVELENGTH', 'FLUX'))
 
         self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
-                                        img_real_wcs=self.real_wcs)
+                                        img_real_wcs=self.chip_wcs)
 
         x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
         x_nominal = int(np.floor(x + 0.5))
@@ -308,7 +255,7 @@ class MockObject(object):
         dy = y - y_nominal
         offset = galsim.PositionD(dx, dy)
 
-        real_wcs_local = self.real_wcs.local(self.real_pos)
+        chip_wcs_local = self.chip_wcs.local(self.real_pos)
 
         flat_cube = chip.flat_cube
 
@@ -323,7 +270,7 @@ class MockObject(object):
             star = star.withFlux(tel.pupil_area * exptime)
             star = galsim.Convolve(psf, star)
             
-            starImg = star.drawImage(nx=100, ny=100, wcs=real_wcs_local, offset=offset)
+            starImg = star.drawImage(nx=100, ny=100, wcs=chip_wcs_local, offset=offset)
 
             origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
                            x_nominal - (starImg.center.x - starImg.xmin)]
@@ -350,7 +297,7 @@ class MockObject(object):
                                        isAlongY=0,
                                        flat_cube=flat_cube)
 
-                self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local)
+                self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
 
                 subImg_p2 = starImg.array[:, subSlitPos + 1:starImg.array.shape[1]]
                 star_p2 = galsim.Image(subImg_p2)
@@ -367,7 +314,7 @@ class MockObject(object):
                                        isAlongY=0,
                                        flat_cube=flat_cube)
 
-                self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local)
+                self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
 
                 del sdp_p1
                 del sdp_p2
@@ -379,7 +326,7 @@ class MockObject(object):
                                     conf=chip.sls_conf[1],
                                     isAlongY=0,
                                     flat_cube=flat_cube)
-                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local)
+                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 del sdp
             elif grating_split_pos_chip >= gal_end[1]:
                 sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
@@ -389,7 +336,7 @@ class MockObject(object):
                                     conf=chip.sls_conf[0],
                                     isAlongY=0,
                                     flat_cube=flat_cube)
-                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=real_wcs_local)
+                self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus=xOrderSigPlus, local_wcs=chip_wcs_local)
                 del sdp
             del psf
         return 1, pos_shear

ObservationSim/ObservationSim.py

@@ -186,7 +186,7 @@ class Observation(object):
                 extName='SCI',
                 timestamp = pointing.timestamp,
                 exptime = pointing.exp_time,
-                readoutTime = 40.)
+                readoutTime = chip.readout_time)
             
             chip_wcs = galsim.FitsWCS(header=h_ext)
 
@@ -382,7 +382,7 @@ class Observation(object):
                     extName='SCI',
                     timestamp=pointing.timestamp,
                     exptime=pointing.exp_time,
-                    readoutTime=40.)
+                    readoutTime=chip.readout_time)
                 
                 chip.img = galsim.Image(chip.img.array, dtype=np.uint16)
                 hdu1 = fits.PrimaryHDU(header=h_prim)

config/config_50sqdeg.yaml

@@ -9,9 +9,9 @@
 # Base diretories and naming setup
 # Can add some of the command-line arguments here as well;
 # OK to pass either way or both, as long as they are consistent
-work_dir: "/share/home/fangyuedong/csst-simulation/workplace/"
+work_dir: "/share/home/fangyuedong/new_sim/workplace/"
 data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
-run_name: "rotate_0"
+run_name: "test_new_sim"
 
 # Whether to use MPI
 run_option:
@@ -45,10 +45,10 @@ catalog_options:
     # AGN_SED_WAVE: "wave_ross13.npy"
 
   # Only simulate stars?
-  star_only: NO
+  star_only: YES
 
   # Only simulate galaxies?
-  galaxy_only: YES
+  galaxy_only: NO
 
   # rotate galaxy ellipticity
   rotateEll: 0. # [degree]
@@ -171,7 +171,7 @@ ins_effects:
   non_linear:     ON  # Whether to add non-linearity
   cosmic_ray:     ON  # Whether to add cosmic-ray
   cray_differ:    ON  # Whether to generate different cosmic ray maps CAL and MS output
-  cte_trail:      ON  # Whether to simulate CTE trails
+  cte_trail:      OFF  # Whether to simulate CTE trails
   saturbloom:     ON  # Whether to simulate Saturation & Blooming
   add_badcolumns: ON  # Whether to add bad columns
   add_hotpixels:  ON  # Whether to add hot pixels

libmoduleDfBF/diffusion_X1.c

+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "nrutil.h"
+
+#define ISSETBITFLAG(x,b) ((x) & (1 << (b)))
+#define ADD_DIFFUSION    1
+#define ADD_BF_FILTER    2
+
+float linearInterp(float xp, float x0, float y0, float x1, float y1);
+
+void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bit_flag)
+{
+  int nx, ny, i,j,k,ks;
+  int it,jt,itt,jtt;
+  int diffuidx[26][2],diffuN,ilow,ih,im,dim[3];
+  float diffua[5][5],cdiffu[26],**bfa;
+  double mvar,mcov,tmp,ma,mb,mc;
+  char fname[100];
+ 
+  nx = ngx_ima;  //input-image size
+  ny = ngy_ima; 
+
+  //0. init. original image with an input array (arr_ima)
+  //1. Adding diffusion effect.
+  if(ISSETBITFLAG(bit_flag, ADD_DIFFUSION)) 
+  {
+    printf("adding diffusion.....\n");
+    printf("ERR: no diffusion filter ...");
+    exit(0);
+  }
+
+
+  //2. Adding BF effect
+  if(ISSETBITFLAG(bit_flag, ADD_BF_FILTER)) 
+  {
+    printf("Adding BF effect...\n");
+    //setup BF correlation fliter
+    float neX;
+    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;
+    float bfaP2[9]={0.9945288, 0.0003041936, 0.0007539311, 0.0002424675, 0.0001226098, 0.00009308617, 0.00008027447, 0.00006309676, 0.00006400052};
+
+    bfa=matrix(-2,2,-2,2);
+
+    // smooth with the BF filter
+    for(i=0;i<nx;i++)for(j=0;j<ny;j++) arr_imc[j+i*ny]=0;
+
+    for(i=0;i<nx;i++)
+    {
+      for(j=0;j<ny;j++)
+      {
+        //rescale BF filter with the local pix value
+        neX = arr_ima[j+i*ny];
+        if(neX >= 10000)
+        {
+          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;
+        }
+        else
+        {
+          neX=10000;
+          bfa[0][0]=0;
+          bfa[0][1]=bfa[0][-1]=bfaP2[1];
+          bfa[-1][0]=bfa[1][0]=bfaP2[2];
+          bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=bfaP2[3];
+          bfa[0][-2]=bfa[0][2]=bfaP2[4];
+          bfa[-2][0]=bfa[2][0]=bfaP2[5];
+          bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=bfaP2[6];
+          bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=bfaP2[7];
+          bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=bfaP2[8];
+        }
+        tmp = 0;
+        for(it=-2;it<=2;it++)
+          for(jt=-2;jt<=2;jt++)
+          {
+            bfa[it][jt] = bfa[it][jt]/neX*arr_ima[j+i*ny];
+            tmp += bfa[it][jt];
+          }
+        bfa[0][0]=1.-tmp;
+
+        // assign electrons according to the BF filter bfat
+        for(it=-2;it<=2;it++)
+        {
+          for(jt=-2;jt<=2;jt++)
+          {
+            itt=i+it;
+            jtt=j+jt;
+            if(itt>=0 && jtt>=0 && itt<nx && jtt<ny)
+              //c0[itt][jtt]+=bfa[it][jt]*b[i][j];
+              arr_imc[jtt+itt*ny] += bfa[it][jt]*arr_ima[j+i*ny];
+          }
+        }
+      }
+    }
+    free_matrix(bfa,-2,2,-2,2);
+  }
+  else
+  {
+    for(i=0;i<nx;i++) for(j=0;j<ny;j++) arr_imc[j+i*ny]=arr_ima[j+i*ny]; ////for ADD_BF False
+  }
+}
+
+
+float linearInterp(float xp, float x0, float y0, float x1, float y1)
+{
+  float yp;
+  yp = y0 + ((y1-y0)/(x1-x0)) * (xp - x0);
+  return yp;
+}
+

profile_C6.sh

@@ -2,13 +2,12 @@
 
 date
 
-python -m cProfile -o C6_profiler_test.pstats /share/home/fangyuedong/csst-simulation/run_sim.py \
-    --config_file config_fgs.yaml \
-    --catalog FGS_Catalog \
-    -c /share/home/fangyuedong/csst-simulation/config
-
-    # --config_file config_50sqdeg.yaml \
-    # --catalog C6_50sqdeg \
+python -m cProfile -o C6_profiler_test.pstats /share/home/fangyuedong/new_sim/csst-simulation/run_sim.py \
+    --config_file config_50sqdeg.yaml \
+    --catalog C6_50sqdeg \
+    -c /share/home/fangyuedong/new_sim/csst-simulation/config
+    # --config_file config_fgs.yaml \
+    # --catalog FGS_Catalog \
     # -c /share/home/fangyuedong/csst-simulation/config
 
     # --config_file config_C6.yaml \

run_C6.pbs

@@ -12,10 +12,10 @@ NP=96
 hostfile=wcl-1,wcl-2
 date
 
-mpirun --oversubscribe -H $hostfile -np $NP python /share/home/fangyuedong/csst-simulation/run_sim.py \
+mpirun --oversubscribe -H $hostfile -np $NP python /share/home/fangyuedong/new_sim/csst-simulation/run_sim.py \
     --config_file config_50sqdeg.yaml \
     --catalog C6_50sqdeg \
-    -c /share/home/fangyuedong/csst-simulation/config
+    -c /share/home/fangyuedong/new_sim/csst-simulation/config
     # --config_file config_C6.yaml \
     # --catalog C6_Catalog \
     # -c /share/home/fangyuedong/csst-simulation/config

setup.py

@@ -2,11 +2,30 @@ from setuptools import setup, find_packages
 
 from setuptools.extension import Extension
 from setuptools.config import read_configuration
+from distutils.command.build_ext import build_ext
 
 from Cython.Build import cythonize
 
 import numpy
 
+class CTypes(Extension): pass
+
+class build_ext(build_ext):
+
+    def build_extension(self, ext):
+        self._ctypes = isinstance(ext, CTypes)
+        return super().build_extension(ext)
+
+    def get_export_symbols(self, ext):
+        if self._ctypes:
+            return ext.export_symbols
+        return super().get_export_symbols(ext)
+
+    def get_ext_filename(self, ext_name):
+        if self._ctypes:
+            return ext_name + '.so'
+        return super().get_ext_filename(ext_name)
+
 extensions = [
     Extension("ObservationSim.MockObject.SpecDisperser.disperse_c.interp", ["ObservationSim/MockObject/SpecDisperser/disperse_c/interp.pyx"],
         include_dirs = [numpy.get_include()],
@@ -17,14 +36,17 @@ extensions = [
         libraries=["m"]),
 ]
 
+df_module = [CTypes('ObservationSim.Instrument.Chip.lib_bf.libmoduleBF',
+                    ['ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c', 'ObservationSim/Instrument/Chip/lib_bf/nrutil.c'],
+                    include_dirs=['ObservationSim/Instrument/Chip/lib_df_bf/', '/usr/include']
+)]
+
 
 # setup(
 #     name = "slssim_disperse",
 #     ext_modules = cythonize(extensions),
 # )
 
-
-
 setup(name='CSSTSim', 
     version='2.1.0', 
     packages=find_packages(),
@@ -43,6 +65,7 @@ setup(name='CSSTSim',
     ],
     package_data = {
         'ObservationSim.Astrometry.lib': ['libshao.so'], 
+        'ObservationSim.Instrument.Chip.lib_bf': ['libmoduleBF.so'],
         'ObservationSim.MockObject.data': ['*.dat'], 
         'ObservationSim.Instrument.data': ['*.txt', '*.dat', '*.json'],
         'ObservationSim.Instrument.data.field_distortion': ['*.pickle'],
@@ -57,5 +80,6 @@ setup(name='CSSTSim',
         'ObservationSim.Straylight.data.sky': ['*.dat'],
         'ObservationSim.Straylight.lib': ['*'],
     },
-    ext_modules = cythonize(extensions),
+    ext_modules = cythonize(extensions) + df_module,
+    cmdclass={'build_ext': build_ext}
 )