append unittest

            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

    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)

"""

from ObservationSim.PSF.PSFModel import PSFModel

LOG_DEBUG = False  #***#

NPSF      = 900    #***# 30*30

PixSizeInMicrons = 5.  #***# in microns

###define PSFInterp###

class PSFInterp(PSFModel):

    def __init__(self, chip, npsf=NPSF, PSF_data=None, PSF_data_file=None, PSF_data_prefix="", sigSpin=0, psfRa=0.15, HocBuild=False):

        if LOG_DEBUG:

    def __init__(self, chip, npsf=NPSF, PSF_data=None, PSF_data_file=None, PSF_data_prefix="", sigSpin=0, psfRa=0.15, HocBuild=False, LOG_DEBUG=False):

        self.LOG_DEBUG = LOG_DEBUG

        if self.LOG_DEBUG:

            print('===================================================')

            print('DEBUG: psf module for csstSim ' \

                  +time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True)

        self.npsf = npsf

        self.PSF_data = self._loadPSF(self.iccd, PSF_data_file, PSF_data_prefix)

        if LOG_DEBUG:

        if self.LOG_DEBUG:

            print('nwave-{:} on ccd-{:}::'.format(self.nwave, self.iccd), flush=True)

            print('self.PSF_data ... ok', flush=True)

            print('Preparing self.[psfMat,cen_col,cen_row] for psfMaker ... ', end='', flush=True)

                self.hoc.append(hoc)

                self.hoclist.append(hoclist)

        if LOG_DEBUG:

        if self.LOG_DEBUG:

            print('ok', flush=True)

            psfSet.append(psfWave)

        fq.close()

        if LOG_DEBUG:

        if self.LOG_DEBUG:

            print('psfSet has been loaded:', flush=True)

            print('psfSet[iwave][ipsf][keys]:', psfSet[0][0].keys(), flush=True)

        return psfSet

            assert(self.hoc != 0), 'hoclist should be built correctly!'

            imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True)

        '''

        ############TEST: START

        TestGaussian = False

        if TestGaussian:

            #pointing_pa = -23.433333

            imPSF= gsx.shear(g1=0.8, g2=0.).rotate(0.*galsim.degrees).drawImage(nx = 256, ny=256, scale=pixSize).array

        ############TEST: END

        '''

        if galsimGSObject:

            imPSFt = np.zeros([257,257])

            return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)

        return imPSF

    '''

    def PSFspin(self, x, y):

        """

        The PSF profile at a given image position relative to the axis center

        qr = np.sqrt((1.0+ell)/(1.0-ell))

        PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)

        return self.psf.shear(PSFshear), PSFshear

if __name__ == '__main__':

    pass

    '''

---

###############################################

#

#  Configuration file for CSST simulation

#           Overall settings

#  CSST-Sim Group, 2024/01/08

#

###############################################

# 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: "/public/home/chengliang/CSSOSDataProductsSims/csst-simulation/tests/UNIT_TEST_DATA/"

data_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/"

run_name: "testRun"

# Project cycle and run counter are used to name the outputs

project_cycle: 9

run_counter: 0

# Run options

run_option:

  use_mpi: YES

  # NOTE: "n_threads" paramters is currently not used in the backend

  # simulation codes. It should be implemented later in the web frontend

  # in order to config the number of threads to request from NAOC cluster

  n_threads: 80

  # Output catalog only?

  # If yes, no imaging simulation will run

  out_cat_only: NO

###############################################

# Catalog setting

###############################################

# Configure the input catalog: options should be implemented

# in the corresponding (user defined) 'Catalog' class

catalog_options:

  input_path:

    # cat_dir: "Catalog_C6_20221212"

    cat_dir: ""

    star_cat: "starcat/"

    galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/"

    # AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits"

  SED_templates_path:

    star_SED: "SpecLib.hdf5"

    galaxy_SED: "sedlibs/"

    AGN_SED: "qsocat/qsosed/"

    # AGN_SED_WAVE: "wave_ross13.npy"

  # Only simulate stars?

  star_only: NO

  # Only simulate galaxies?

  galaxy_only: NO

  # rotate galaxy ellipticity

  rotateEll: 0. # [degree]

  seed_Av: 121212    # Seed for generating random intrinsic extinction

###############################################

# Observation setting

###############################################

obs_setting:

  # (Optional) a file of point list 

  # if you just want to run default pointing:

  # - pointing_dir: null

  # - pointing_file: null

  pointing_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/pointing_gir25/"

  pointing_file: "pointing_50_1.dat"

  obs_config_file: "/public/home/chengliang/CSSOSDataProductsSims/csst-simulation/config/obs_config_SCI_WIDE_phot.yaml"

  # Run specific pointing(s):

  # - give a list of indexes of pointings: [ip_1, ip_2...]

  # - run all pointings: null

  # Note: only valid when a pointing list is specified

  run_pointings: [0]

  # Whether to enable astrometric modeling

  enable_astrometric_model: True

  # Cut by saturation magnitude in which band?

  cut_in_band: "z"

  # saturation magnitude margin

  mag_sat_margin: -2.5

  # mag_sat_margin: -15.

  # limiting magnitude margin

  mag_lim_margin: +1.0

###############################################

# PSF setting

###############################################

psf_setting:

  # Which PSF model to use:

  # "Gauss": simple gaussian profile

  # "Interp": Interpolated PSF from sampled ray-tracing data

  psf_model: "Interp"

  # PSF size [arcseconds]

  # radius of 80% energy encircled

  # NOTE: only valid for "Gauss" PSF

  # psf_rcont: 0.15

  # path to PSF data

  # NOTE: only valid for "Interp" PSF

  # PSF models for photometry survey simulation

  psf_pho_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/psfCube/set1_dynamic/"

  # PSF models for slitless spectrum survey simulation

  psf_sls_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/SLS_PSF_PCA_fp/"

###############################################

# Shear setting

###############################################

shear_setting:

  # Options to generate mock shear field:

  # "constant": all galaxies are assigned a constant reduced shear

  # "catalog": get shear values from catalog

  shear_type: "constant"

  # For constant shear field

  reduced_g1: 0.

  reduced_g2: 0.

###############################################

# Output options

###############################################

output_setting:

  output_format:  "channels" # Whether to export as 16 channels (subimages) with pre- and over-scan ("image"/"channels")

  shutter_output: OFF # Whether to export shutter effect 16-bit image

  prnu_output:    OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files

###############################################

# Random seeds

###############################################

random_seeds:

  seed_poisson:         20210601  # Seed for Poisson noise

  seed_CR:              20210317  # Seed for generating random cosmic ray maps

  seed_flat:            20210101  # Seed for generating random flat fields

  seed_prnu:            20210102  # Seed for photo-response non-uniformity

  seed_gainNonUniform:  20210202  # Seed for gain nonuniformity

  seed_biasNonUniform:  20210203  # Seed for bias nonuniformity

  seed_rnNonUniform:    20210204  # Seed for readout-noise nonuniformity

  seed_badcolumns:      20240309  # Seed for bad columns

  seed_defective:       20210304  # Seed for defective (bad) pixels

  seed_readout:         20210601  # Seed for read-out gaussian noise

...

import unittest

import sys,os,math

from itertools import islice

import numpy as np

import copy

import ctypes

import galsim

import yaml

from astropy.io import fits

from ObservationSim.Instrument import Chip, Filter, FilterParam, FocalPlane

from ObservationSim.Instrument.Chip import ChipUtils as chip_utils

#from ObservationSim.sim_steps import add_brighter_fatter_CTE

from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim

try:

    import importlib.resources as pkg_resources

except ImportError:

    # Try backported to PY<37 'importlib_resources'

    import importlib_resources as pkg_resources

### test FUNCTION --- START ###

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

### test FUNCTION --- END ###

def defineCCD(iccd, config_file):

    with open(config_file, "r") as stream:

        try:

            config = yaml.safe_load(stream)

            #for key, value in config.items():

            #    print (key + " : " + str(value))

        except yaml.YAMLError as exc:

            print(exc)

    chip = Chip(chipID=iccd, config=config)

    chip.img = galsim.ImageF(400, 200) #galsim.ImageF(chip.npix_x, chip.npix_y)

    focal_plane = FocalPlane(chip_list=[iccd])

    chip.img.wcs= focal_plane.getTanWCS(192.8595, 27.1283, -113.4333*galsim.degrees, chip.pix_scale)

    return chip

def defineFilt(chip):

    filter_param = FilterParam()

    filter_id, filter_type = chip.getChipFilter()

    filt = Filter(

        filter_id=filter_id,

        filter_type=filter_type,

        filter_param=filter_param,

        ccd_bandpass=chip.effCurve)

    bandpass_list = filt.bandpass_sub_list

    return filt

class detModule_coverage(unittest.TestCase):

    def __init__(self, methodName='runTest'):

        super(detModule_coverage, self).__init__(methodName)

        self.dataPath = "/public/home/chengliang/CSSOSDataProductsSims/csst-simulation/tests/UNIT_TEST_DATA" ##os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_fz_gc1')

        self.iccd = 1

    def test_add_brighter_fatter(self):

        config_file = os.path.join(self.dataPath, 'config_test.yaml')

        chip = defineCCD(self.iccd, config_file)

        filt = defineFilt(chip)

        print(chip.chipID)

        print(chip.cen_pix_x, chip.cen_pix_y)

        #objA-lowSFB

        obj = galsim.Gaussian(sigma=0.2, flux=1000)

        arr = obj.drawImage(nx=64, ny=64, scale=0.074).array

        chip.img.array[(100-32):(100+32),(200-32):(200+32)] = arr[:,:] 

        img_old = copy.deepcopy(chip.img)

        img_new = add_brighter_fatter(img=chip.img)

        arr1= img_old.array

        arr2= img_new.array

        deltaA_max = np.max(np.abs(arr2-arr1))

        print('deltaA-max:', np.max(np.abs(arr2-arr1))) 

        print('deltaA-min:', np.min(np.abs(arr2-arr1))) 

        #objB-highSFB

        obj = galsim.Gaussian(sigma=0.2, flux=10000)

        arr = obj.drawImage(nx=64, ny=64, scale=0.074).array

        chip.img.array[(100-32):(100+32),(200-32):(200+32)] = arr[:,:] 

        img_old = copy.deepcopy(chip.img)

        img_new = add_brighter_fatter(img=chip.img)

        arr3= img_old.array

        arr4= img_new.array

        deltaB_max = np.max(np.abs(arr4-arr3))

        print('deltaB-max:', np.max(np.abs(arr4-arr3))) 

        print('deltaB-min:', np.min(np.abs(arr4-arr3))) 

        self.assertTrue( deltaB_max > deltaA_max )

    def test_apply_CTE(self):

        config_file = os.path.join(self.dataPath, 'config_test.yaml')

        chip = defineCCD(self.iccd, config_file)

        filt = defineFilt(chip)

        print(chip.chipID)

        print(chip.cen_pix_x, chip.cen_pix_y)

        print("  Apply CTE Effect")

        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("UNIT_TEST_DATA/testCTE_image_before.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("UNIT_TEST_DATA/testCTE_image_after.fits",data=image_cti,overwrite=True)

if __name__ == '__main__':

    unittest.main()