Merge branch 'develop'

ObservationSim/Config/Pointing.py

-import numpy as np
+import os
+import shutil
+import yaml
 import galsim
+import numpy as np
 from astropy.time import Time
 
+from ObservationSim.Config._util import get_obs_id
+import ObservationSim.Instrument._util as _util
+
 class Pointing(object):
-    def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='MS'):
+    def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='SCI', pointing_type_code='101', pointing_id = '00000001', obs_config_file=None, t_shutter_open = 1.3, t_shutter_close = 1.3):
         self.id = id
         self.ra = ra
         self.dec = dec
@@ -14,15 +20,58 @@ class Pointing(object):
         self.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz
         self.exp_time = exp_time
         self.pointing_type = pointing_type
-        self.jdt = 0.
+        self.pointing_type_code = pointing_type_code
+        self.obs_id = pointing_id
+        self.survey_field_type = 'WIDE'
+        self.jdt = 0. 
+        self.obs_config_file = obs_config_file
+        self.t_shutter_open = t_shutter_open
+        self.t_shutter_close = t_shutter_close
+        self.output_dir = "."
+        if self.obs_config_file is not None:
+            with open(self.obs_config_file, "r") as stream:
+                try:
+                    self.obs_param = yaml.safe_load(stream)
+                except yaml.YAMLError as exc:
+                    print(exc)
+            if self.obs_param["obs_type"]:
+                self.pointing_type = self.obs_param["obs_type"]
+            if self.obs_param["obs_type_code"]:
+                self.pointing_type_code = self.obs_param["obs_type_code"]
+            if self.obs_param["obs_id"]:
+                self.obs_id = str(self.obs_param["obs_id"])
+
+    
+    def get_full_depth_exptime(self, filter_type):
+        if self.survey_field_type == 'WIDE':
+            if filter_type in _util.SPEC_FILTERS:
+                return 150. * 4
+            else:
+                if filter_type.lower() in ['nuv', 'y']:
+                    return 150. * 4
+                elif filter_type.lower() in ['u', 'g', 'r', 'i', 'z']:
+                    return 150. * 2
+                else:
+                    return max(150., self.exp_time) # [TODO] for FGS
+        elif self.survey_field_type == 'DEEP':
+            if filter_type in _util.SPEC_FILTERS:
+                return 250. * 4 * 4
+            else:
+                if filter_type.lower() in ['nuv', 'y']:
+                    return 250. * 4 * 4
+                elif filter_type.lower() in ['u', 'g', 'r', 'i', 'z']:
+                    return 250. * 2 * 4
+                else:
+                    return max(150., self.exp_time) # [TODO] for FGS
+
 
-    def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'):
+    def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='SCI'):
         self.id = id
         col_len = len(columns)
         self.ra = float(columns[0])
         self.dec = float(columns[1])
         self.img_pa = float(columns[4]) * galsim.degrees
-        self.pointing_type = pointing_type
+        # self.pointing_type = pointing_type
         if col_len > 5:
             jdt = np.double(columns[5])
             t_temp = Time(jdt, format='jd')
@@ -33,10 +82,55 @@ class Pointing(object):
             self.sat_z = float(columns[8])
             self.sun_x = float(columns[9])
             self.sun_y = float(columns[10])
-            self.sun_z = float(columns[1])
+            self.sun_z = float(columns[11])
             self.sat_vx = float(columns[15])
             self.sat_vy = float(columns[16])
             self.sat_vz = float(columns[17])
             self.exp_time = float(columns[18])
+            is_deep = float(columns[19])
+            # [TODO] Can also define other survey types
+            if is_deep != -1.0:
+                self.survey_field_type = "DEEP"
+            
+            if not self.obs_config_file:
+                self.obs_config_file = str(columns[20])
+                with open(self.obs_config_file, "r") as stream:
+                    try:
+                        self.obs_param = yaml.safe_load(stream)
+                    except yaml.YAMLError as exc:
+                        print(exc)
+            
+            self.pointing_type_code = columns[21][0:3]
+            self.obs_id = columns[21][3:]
+
+            self.pointing_type = self.obs_param["obs_type"]
         else:
             self.timestamp = t
+
+    def make_output_pointing_dir(self, overall_config, copy_obs_config=False):
+        run_dir = os.path.join(overall_config["work_dir"], overall_config["run_name"])
+        if not os.path.exists(run_dir):
+            try:
+                os.makedirs(run_dir, exist_ok=True)
+            except OSError:
+                pass
+        self.output_prefix = get_obs_id(
+            img_type=self.pointing_type, 
+            project_cycle=overall_config["project_cycle"], 
+            run_counter=overall_config["run_counter"], 
+            pointing_id=self.obs_id, 
+            pointing_type_code = self.pointing_type_code)
+        self.output_dir = os.path.join(run_dir, self.output_prefix)
+        if not os.path.exists(self.output_dir):
+            try:
+                os.makedirs(self.output_dir, exist_ok=True)
+            except OSError:
+                pass
+        if copy_obs_config and self.obs_config_file:
+            obs_config_output_path = os.path.join(self.output_dir, os.path.basename(self.obs_config_file))
+            if not os.path.exists(obs_config_output_path):
+                try:
+                    shutil.copy(self.obs_config_file, self.output_dir)
+                except OSError:
+                    pass
+

ObservationSim/Config/__init__.py

-from .Config import *
 from .ChipOutput import ChipOutput
 from .Pointing import Pointing
\ No newline at end of file

ObservationSim/Config/_util.py

+
+def get_obs_id(img_type='SCI', project_cycle=6, run_counter=0, pointing_id='00000001',pointing_type_code='101'):
+    # obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'}
+    # obs_type = {'SCIE': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'CAL': '01'}
+    obs_id = pointing_type_code + str(int(project_cycle)).rjust(2, '0') + str(int(run_counter)) + pointing_id
+    return obs_id
+
+# def get_obs_id(img_type='SCI', project_cycle=6, run_counter=0, pointing_num=0):
+#     # obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'}
+#     obs_type = {'SCIE': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'CAL': '01'}
+#     obs_id = '1'+ obs_type[img_type] + str(int(project_cycle)).rjust(2, '0') + str(int(run_counter)) + str(pointing_num).rjust(8,'0')
+#     return obs_id
+
+def get_file_type(img_type='SCI'):
+    file_type = {'SCI':'SCI', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'}
+    return file_type[img_type]
\ No newline at end of file

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 getChipSLSGratingID(chipID):
+    gratingID = ['','']
+    if chipID == 1: gratingID = ['GI2', 'GI1']
+    if chipID == 2: gratingID = ['GV4', 'GV3']
+    if chipID == 3: gratingID = ['GU2', 'GU1']
+    if chipID == 4: gratingID = ['GU4', 'GU3']
+    if chipID == 5: gratingID = ['GV2', 'GV1']
+    if chipID == 10: gratingID = ['GI4', 'GI3']
+    if chipID == 21: gratingID = ['GI6', 'GI5']
+    if chipID == 26: gratingID = ['GV8', 'GV7']
+    if chipID == 27: gratingID = ['GU6', 'GU5']
+    if chipID == 28: gratingID = ['GU8', 'GU7']
+    if chipID == 29: gratingID = ['GV6', 'GV5']
+    if chipID == 30: gratingID = ['GI8', 'GI7']
+    return gratingID
+
+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, pointing, img_type=None, img_type_code=None, project_cycle='9', run_counter='1'):
+    if (img_type is None) or (img_type_code is None):
+        img_type = pointing.pointing_type
+        img_type_code = pointing.pointing_type_code
+
+    h_prim = generatePrimaryHeader(
+        xlen=chip.npix_x, 
+        ylen=chip.npix_y, 
+        pointing_id = pointing.obs_id,
+        pointing_type_code = img_type_code, 
+        ra=pointing.ra, 
+        dec=pointing.dec, 
+        pixel_scale=chip.pix_scale,
+        time_pt = pointing.timestamp,
+        exptime=pointing.exp_time,
+        im_type=img_type,
+        sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
+        sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz],
+        project_cycle=project_cycle,
+        run_counter=run_counter,
+        chip_name=str(chip.chipID).rjust(2, '0'))
+    h_ext = generateExtensionHeader(
+        chip=chip,
+        xlen=chip.npix_x, 
+        ylen=chip.npix_y, 
+        ra=pointing.ra, 
+        dec=pointing.dec, 
+        pa=pointing.img_pa.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=img_type,
+        timestamp=pointing.timestamp,
+        exptime=pointing.exp_time,
+        readoutTime=chip.readout_time,
+        t_shutter_open=pointing.t_shutter_open, 
+        t_shutter_close=pointing.t_shutter_close)
+    return h_prim, h_ext
+
+def output_fits_image(chip, pointing, img, output_dir, img_type=None, img_type_code=None, project_cycle='9', run_counter='1'):
+    h_prim, h_ext = generateHeader(
+        chip=chip, 
+        pointing=pointing, 
+        img_type=img_type, 
+        img_type_code=img_type_code,
+        project_cycle=project_cycle,
+        run_counter=run_counter)
+    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, chip, read_noise, readout_time, dark_noise, exptime=150., InputDark=None):
+    if InputDark == 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)
+        base_img1 = base_level * np.ones_like(img.array)
+    else:
+        base_img1 = np.zeros_like(img.array)
+
+    ny = int(chip.npix_y/2)
+    nx = chip.npix_x
+    arr = np.arange(ny).reshape(ny, 1)
+    arr = np.broadcast_to(arr, (ny, nx))
+    base_img2 = np.zeros_like(img.array)
+    base_img2[:ny, :] = arr
+    base_img2[ny:, :] = arr[::-1,:]
+    base_img2[:,:]    = base_img2[:,:]*(readout_time/ny)*dark_noise
+    return base_img1+base_img2
+
+def add_poisson(img, chip, exptime=150., seed=0, sky_level=0., poisson_noise=None, dark_noise=None, InputDark=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, chip=chip, read_noise=read_noise, readout_time=chip.readout_time, dark_noise=dark_noise, exptime=exptime, InputDark=InputDark)
+    img += base_img
+    img.addNoise(poisson_noise)
+    # img -= read_noise**2
+
+    if InputDark != None:
+        hdu = fits.open(InputDark)  ##"Instrument/data/dark/dark_1000s_example_0.fits"
+        img += hdu[0].data/hdu[0].header['exptime']*exptime
+        hdu.close()
+    return img, base_img
+
+def add_brighter_fatter(img):
+    #Inital dynamic lib
+    try:
+        with pkg_resources.files('ObservationSim.Instrument.Chip.libBF').joinpath("libmoduleBF.so") as lib_path:
+            lib_bf = ctypes.CDLL(lib_path)
+    except AttributeError:
+        with pkg_resources.path('ObservationSim.Instrument.Chip.libBF', "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.reshape(arr_imc, [nx, ny])
+    del arr_ima, arr_imc
+    return img
+"""
+def add_inputdark(img, chip, exptime):
+    fname = "/share/home/weichengliang/CSST_git/test_new_sim/csst-simulation/ObservationSim/Instrument/data/dark/dark_1000s_example_0.fits"
+    hdu = fits.open(fname)
+    #ny, nx = img.array.shape
+    #inputdark = np.zeros([ny, nx])
+    img.array[:, :] += hdu[0].data/hdu[0].header['exptime']*exptime
+    hdu.close()
+    del inputdark
+    return img
+"""
+def AddPreScan(GSImage, pre1=27, pre2=4, over1=71, over2=80, nsecy = 2, nsecx=8):
+    img= GSImage.array
+    ny, nx = img.shape
+    dx = int(nx/nsecx)
+    dy = int(ny/nsecy)
+
+    imgt=np.zeros([int(nsecy*nsecx), int(ny/nsecy+pre2+over2), int(nx/nsecx+pre1+over1)])
+    for iy in range(nsecy):
+        for ix in range(nsecx):
+            if iy % 2 == 0:
+                tx = ix
+            else:
+                tx = (nsecx-1)-ix
+            ty = iy 
+            chunkidx = int(tx+ty*nsecx) #chunk1-[1,2,3,4], chunk2-[5,6,7,8], chunk3-[9,10,11,12], chunk4-[13,14,15,16]
+
+            imgtemp = np.zeros([int(ny/nsecy+pre2+over2), int(nx/nsecx+pre1+over1)])
+            if int(chunkidx/4) == 0:
+                imgtemp[pre2:pre2+dy, pre1:pre1+dx] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                imgt[chunkidx, :, :] = imgtemp
+            if int(chunkidx/4) == 1:
+                imgtemp[pre2:pre2+dy, over1:over1+dx] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                imgt[chunkidx, :, :] = imgtemp #[:, ::-1]
+            if int(chunkidx/4) == 2:
+                imgtemp[over2:over2+dy, over1:over1+dx] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                imgt[chunkidx, :, :] = imgtemp #[::-1, ::-1]
+            if int(chunkidx/4) == 3:
+                imgtemp[over2:over2+dy, pre1:pre1+dx] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                imgt[chunkidx, :, :] = imgtemp #[::-1, :]
+
+    imgtx1 = np.hstack(imgt[:nsecx:,       :, :])  #hstack chunk(1,2)-[1,2,3,4,5,6,7,8]
+    imgtx2 = np.hstack(imgt[:(nsecx-1):-1, :, :])  #hstack chunk(4,3)-[16,15,14,13,12,11,,10,9]
+
+    newimg = galsim.Image(int(nx+(pre1+over1)*nsecx), int(ny+(pre2+over2)*nsecy), init_value=0)
+    newimg.array[:, :] = np.concatenate([imgtx1, imgtx2]) #vstack chunk(1,2) & chunk(4,3)
+
+    newimg.wcs = GSImage.wcs
+    return newimg
+
+def AddPreScanFO(GSImage, pre1=27, pre2=4, over1=71, over2=80, nsecy = 1, nsecx=16):
+    img= GSImage.array
+    ny, nx = img.shape
+    dx = int(nx/nsecx)
+    dy = int(ny/nsecy)
+
+    newimg = galsim.Image(int(nx+(pre1+over1)*nsecx), int(ny+(pre2+over2)*nsecy), init_value=0)
+    for ix in range(nsecx):
+        newimg.array[pre2:pre2+dy, pre1+ix*(dx+pre1+over1):pre1+dx+ix*(dx+pre1+over1)] = img[0:dy, 0+ix*dx:dx+ix*dx]
+
+    newimg.wcs = GSImage.wcs
+    return newimg
+
+
+def formatOutput(GSImage, nsecy = 2, nsecx=8):
+    img = GSImage.array
+    ny, nx = img.shape
+    dx = int(nx/nsecx)
+    dy = int(ny/nsecy)
+    
+    imgt = np.zeros([int(nsecx*nsecy), dy, dx])
+    for iy in range(nsecy):
+        for ix in range(nsecx):
+            if iy % 2 == 0:
+                tx = ix
+            else:
+                tx = (nsecx-1)-ix
+            ty = iy
+            chunkidx = int(tx+ty*nsecx)
+            if int(chunkidx/4) == 0:
+                imgt[chunkidx, :, :] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+            if int(chunkidx/4) == 1:
+                imgt[chunkidx, :, :] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+            if int(chunkidx/4) == 2:
+                imgt[chunkidx, :, :] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+            if int(chunkidx/4) == 3:
+                imgt[chunkidx, :, :] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+    
+    imgttx0 = np.hstack(imgt[ 0:4:,    :,    :])
+    imgttx1 = np.hstack(imgt[ 4:8:,    :, ::-1])
+    imgttx2 = np.hstack(imgt[8:12:, ::-1, ::-1])
+    imgttx3 = np.hstack(imgt[12:16:,::-1,    :])
+    
+    newimg = galsim.Image(int(dx*nsecx*nsecy), dy, init_value=0)
+    newimg.array[:, :] = np.hstack([imgttx0, imgttx1, imgttx2, imgttx3])
+    return newimg
+
+def formatRevert(GSImage, nsecy = 1, nsecx=16):
+    img = GSImage.array
+    ny, nx = img.shape
+    dx = int(nx/nsecx)
+    dy = int(ny/nsecy)
+
+    newimg = galsim.Image(int(dx*8), int(dy*2), init_value=0)
+
+    for ix in range(0,4):
+        tx = ix
+        newimg.array[0:dy, 0+tx*dx:dx+tx*dx] = img[:, 0+ix*dx:dx+ix*dx]
+    for ix in range(4,8):
+        tx = ix
+        newimg.array[0:dy, 0+tx*dx:dx+tx*dx] = img[:, 0+ix*dx:dx+ix*dx][:, ::-1]
+    for ix in range(8,12):
+        tx = 7-(ix-8)
+        newimg.array[0+dy:dy+dy, 0+tx*dx:dx+tx*dx] = img[:, 0+ix*dx:dx+ix*dx][::-1, ::-1]
+    for ix in range(12,16):
+        tx = 7-(ix-8)
+        newimg.array[0+dy:dy+dy, 0+tx*dx:dx+tx*dx] = img[:, 0+ix*dx:dx+ix*dx][::-1, :]
+
+    return newimg

ObservationSim/Instrument/Chip/libBF/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/libBF/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/Instrument/Chip/libCTI/CTI_modeling.py

+from ctypes import CDLL, POINTER, c_int, c_double,c_float,c_long,c_char_p
+from numpy.ctypeslib import ndpointer
+import numpy.ctypeslib as clb
+import numpy as np
+from astropy.io import fits
+from scipy.stats import randint
+from glob import glob
+from datetime import datetime
+import os
+
+lib_path = os.path.dirname(os.path.realpath(__file__))
+
+#lib_path += "/add_CTI.so"
+lib_path += "/libmoduleCTI.so"
+lib = CDLL(lib_path)
+CTI_simul = lib.__getattr__('CTI_simul')
+CTI_simul.argtypes = [POINTER(POINTER(c_int)),c_int,c_int,c_int,c_int,POINTER(c_float),POINTER(c_float),\
+                    c_float,c_float,c_float,c_int,POINTER(c_int),c_int,POINTER(POINTER(c_int))]
+'''
+get_trap_h = lib.__getattr__('save_trap_map')
+get_trap_h.argtypes = [POINTER(c_int), c_int, c_int, c_int, c_int, POINTER(c_float), c_float, c_float, c_char_p]
+
+def get_trap_map(seeds,nx,ny,nmax,rho_trap,beta,c,out_dir):
+    hsp_result = np.zeros(ny*nx*nmax)
+    nsp = len(rho_trap)
+    seeds1 = seeds.astype(np.int32)
+    seeds_p = np.ctypeslib.as_ctypes(seeds1)
+    rho_trap1 = rho_trap.astype(np.float32)
+    rho_trap_p = np.ctypeslib.as_ctypes(rho_trap1)
+    filename = (out_dir+"/trap.bin").encode('utf-8')
+    get_trap_h(seeds_p,c_int(int(nsp)),c_int(int(nx)),c_int(int(ny)),\
+            c_int(int(nmax)),rho_trap_p,c_float(beta),\
+            c_float(c),filename)
+        
+def bin2fits(bin_file,fits_dir,nsp,nx,ny,nmax):
+    data = np.fromfile(bin_file,dtype=np.float32)
+    data = data.reshape(nx,nsp,ny,nmax).transpose(1,3,2,0)
+    for i in range(nsp):
+        print("transfering trap type "+str(i+1))
+        datai = data[i]
+        ntrap = datai[0,:,:]
+        for j in range(nmax-1):
+            h = datai[j+1,:,:]
+            h[np.where(ntrap<j+1)] = 0
+            datai[j+1,:,:] = h
+        fits.writeto(fits_dir+"/trap_"+str(i+1)+".fits",datai,overwrite=True)
+'''
+        
+def numpy_matrix_to_int_pointer(arr):
+    int_pointer_array = (POINTER(c_int)*arr.shape[0])()
+    for i in range(arr.shape[0]):
+        arr1 = np.array(arr[i].copy().tolist(),dtype=np.int32)
+        int_pointer_array[i] = np.ctypeslib.as_ctypes(arr1)
+    return int_pointer_array
+def pointer_to_numpy_matrix(arr_pointer,row,col):
+    arr = np.zeros((row,col))
+    for i in range(row):
+        for j in range(col):
+            arr[i,j] = arr_pointer[i][j]
+    return arr
+def CTI_sim(im,nx,ny,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed=0):
+    image = im.T
+    nx_c,ny_c,noverscan_c,nsp_c,nmax_c = c_int(nx),c_int(ny),c_int(noverscan),c_int(nsp),c_int(nmax)
+    ntotal = ny+noverscan
+    beta_c,w_c,c_c = c_float(beta),c_float(w),c_float(c)
+    t_p = np.ctypeslib.as_ctypes(t)
+    rho_trap_p = np.ctypeslib.as_ctypes(rho_trap)
+    image_p = numpy_matrix_to_int_pointer(image)
+    trap_seeds1 = trap_seeds.astype(np.int32)
+    trap_seeds_p = np.ctypeslib.as_ctypes(trap_seeds1)
+    release_seed_c = c_int(release_seed)
+    image_cti = np.zeros((nx,ntotal))
+    image_cti = image_cti.astype(np.int32)
+    image_cti_p = numpy_matrix_to_int_pointer(image_cti)
+    print(datetime.now())
+    CTI_simul(image_p,nx,ny,noverscan,nsp,rho_trap_p,t_p,beta,w,c,nmax,trap_seeds_p,release_seed_c,image_cti_p)
+    print(datetime.now())
+    image_cti_result = np.zeros((nx,ntotal))
+    for i in range(nx):
+        for j in range(ntotal):
+            image_cti_result[i,j] = image_cti_p[i][j]
+    return image_cti_result.T
+
+"""
+if __name__ =='__main__':
+    nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10
+    ntotal = 4700
+    beta,w,c = 0.478,84700,0
+    t = np.array([0.74,7.7,37],dtype=np.float32)
+    rho_trap = np.array([0.6,1.6,1.4],dtype=np.float32)
+    trap_seeds = np.array([0,100,1000],dtype=np.int32)
+    release_seed = 500
+    image = fits.getdata("inputdata/image.fits").astype(np.int32)
+    get_trap_map(trap_seeds,nx,ny,nmax,rho_trap,beta,c,".")
+    bin2fits("trap.bin",".",nsp,nx,ny,nmax)
+    image_cti = CTI_sim(image,nx,ny,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed)  
+    fits.writeto("output/image_CTI.fits",data=image_cti,overwrite=True)
+"""

ObservationSim/Instrument/Chip/libCTI/Makefile

+#include objects here:
+objects = src/add_CTI1.o src/nrutil.o src/ran1.o src/ran2.o src/poidev.o src/gammln.o src/gasdev.o src/sort.o src/creattraps.o
+
+add_CTI.so:  $(objects)
+	gcc -shared -fPIC -std=c99 -o $@ $(objects) -lm
+
+# general compilation rules 
+.SUFFIXES:  .c   .o
+.c.o:
+	cc -c $< -O3 -shared -fPIC -std=c99 -o $@
+.PHONY : clean
+clean:
+	rm -f src/*.o add_CTI.so

ObservationSim/Instrument/Chip/libCTI/readme

+CTI_lgl_v0.3/
+
+add_CTIfinal.c >> add_CTI.c