diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py
index 33914c238e3eb16891675ddec3e1da63234c4942..4e143dd15453f097fca6660a68c9ac5f83706c57 100755
--- a/ObservationSim/Instrument/Chip/Chip.py
+++ b/ObservationSim/Instrument/Chip/Chip.py
@@ -15,6 +15,8 @@ from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensio
 from ObservationSim.Instrument._util import rotate_conterclockwise
 from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
 
+from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim
+
 try:
     import importlib.resources as pkg_resources
 except ImportError:
@@ -394,15 +396,58 @@ class Chip(FocalPlane):
             img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)
 
         # Apply CTE Effect
+        ###if config["ins_effects"]["cte_trail"] == True:
+        ###    chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
+        ###    img = effects.CTE_Effect(GSImage=img, threshold=27)
+
+        pre1 = self.prescan_x  #27
+        over1= self.overscan_x #71
+        pre2 = self.prescan_y  #0 #4
+        over2= self.overscan_y #84 #80
+
         if config["ins_effects"]["cte_trail"] == True:
             chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
-            img = effects.CTE_Effect(GSImage=img, threshold=27)
+            ###img = effects.CTE_Effect(GSImage=img, threshold=27)
+            ###CTI_modeling
+            ### 2*8 -> 1*16 img-layout
+            img = chip_utils.formatOutput(GSImage=img)
+            self.nsecy = 1
+            self.nsecx = 16
+
+            img_arr = img.array
+            ny, nx = img_arr.shape
+            dx = int(nx/self.nsecx)
+            dy = int(ny/self.nsecy)
+            newimg = galsim.Image(nx, int(ny+over2), init_value=0)
+            for ichannel in range(16):
+                #nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10
+                noverscan,nsp,nmax = over2,3,10
+                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,1000,10000],dtype=np.int32) + ichannel + self.chipID*16
+                release_seed = 50 + ichannel + pointing_ID*30  + self.chipID*16
+                newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(img_arr[:, 0+ichannel*dx:dx+ichannel*dx],dx,dy,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed)
+            newimg.wcs = img.wcs
+            del img
+            img = newimg
+
+            ### 1*16 -> 2*8 img-layout
+            img = chip_utils.formatRevert(GSImage=img)
+            self.nsecy = 2
+            self.nsecx = 8
 
         ### prescan & overscan
         if config["ins_effects"]["add_prescan"] == True:
-            img = chip_utils.AddPreScan(GSImage=img, pre1=27, pre2=4, over1=71, over2=80)
+            chip_utils.log_info(msg="  Apply pre/over-scan", logger=self.logger)
+            if config["ins_effects"]["cte_trail"] == False:
+                img = chip_utils.AddPreScan(GSImage=img, pre1=pre1, pre2=pre2, over1=over1, over2=over2)
+            if config["ins_effects"]["cte_trail"] == True:
+                img = chip_utils.AddPreScan(GSImage=img, pre1=pre1, pre2=pre2, over1=over1, over2=0)
+
         ### 1*16 output
         if config["ins_effects"]["format_output"] == True:
+            chip_utils.log_info(msg="  Apply 1*16 format", logger=self.logger)
             img = chip_utils.formatOutput(GSImage=img)
             self.nsecy = 1
             self.nsecx = 16
diff --git a/ObservationSim/Instrument/Chip/ChipUtils.py b/ObservationSim/Instrument/Chip/ChipUtils.py
index 381d310da490cf46704635b74f737b7d755074f4..8eb95f3cdbd61e5b2dddd50984b92cf6e3c4e03b 100644
--- a/ObservationSim/Instrument/Chip/ChipUtils.py
+++ b/ObservationSim/Instrument/Chip/ChipUtils.py
@@ -192,11 +192,10 @@ def add_poisson(img, chip, exptime=150., seed=0, sky_level=0., poisson_noise=Non
 def add_brighter_fatter(img):
     #Inital dynamic lib
     try:
-        with pkg_resources.files('ObservationSim.Instrument.Chip.lib_bf').joinpath("libmoduleBF.so") as lib_path:
-            print('--1', lib_path)
+        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.lib_bf', "libmoduleBF.so") as lib_path:
+        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]
     
@@ -245,27 +244,41 @@ def AddPreScan(GSImage, pre1=27, pre2=4, over1=71, over2=80, nsecy = 2, nsecx=8)
 
             imgtemp = np.zeros([int(ny/nsecy+pre2+over2), int(nx/nsecx+pre1+over1)])
             if int(chunkidx/4) == 0:
-                imgtemp[pre2:-over2, pre1:-over1] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                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:-over2, over1:-pre1] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                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:-pre2, over1:-pre1] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                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:-pre2, pre1:-over1] = img[iy*dy:(iy+1)*dy, ix*dx:(ix+1)*dx]
+                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:,       :, :])
-    imgtx2 = np.hstack(imgt[:(nsecx-1):-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])
+    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
@@ -299,3 +312,25 @@ def formatOutput(GSImage, nsecy = 2, nsecx=8):
     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
diff --git a/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..024f6e55eb27718a99d6153eca7714a75a1289f2
Binary files /dev/null and b/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc differ
diff --git a/ObservationSim/Instrument/Chip/__pycache__/ChipUtils.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/ChipUtils.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..54ec3424f546352d2327c7ecf9bb57982c696fe4
Binary files /dev/null and b/ObservationSim/Instrument/Chip/__pycache__/ChipUtils.cpython-38.pyc differ
diff --git a/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..661c5b51b54a66f76db46196c1aea5917dbe4314
Binary files /dev/null and b/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc differ
diff --git a/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..88fdd9681bd870a3648dd5635ee99fb7c01da9d6
Binary files /dev/null and b/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc differ
diff --git a/ObservationSim/Instrument/Chip/lib_bf/__init__.py b/ObservationSim/Instrument/Chip/libBF/__init__.py
similarity index 100%
rename from ObservationSim/Instrument/Chip/lib_bf/__init__.py
rename to ObservationSim/Instrument/Chip/libBF/__init__.py
diff --git a/ObservationSim/Instrument/Chip/libBF/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Instrument/Chip/libBF/__pycache__/__init__.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..d5b238d75bc180ae1222ae9e01f4f25e5c1dfd2e
Binary files /dev/null and b/ObservationSim/Instrument/Chip/libBF/__pycache__/__init__.cpython-38.pyc differ
diff --git a/ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c b/ObservationSim/Instrument/Chip/libBF/diffusion_X1.c
similarity index 100%
rename from ObservationSim/Instrument/Chip/lib_bf/diffusion_X1.c
rename to ObservationSim/Instrument/Chip/libBF/diffusion_X1.c
diff --git a/ObservationSim/Instrument/Chip/libBF/libmoduleBF.so b/ObservationSim/Instrument/Chip/libBF/libmoduleBF.so
new file mode 100755
index 0000000000000000000000000000000000000000..63400befa404b98806d231811541d12d6fa8f613
Binary files /dev/null and b/ObservationSim/Instrument/Chip/libBF/libmoduleBF.so differ
diff --git a/ObservationSim/Instrument/Chip/lib_bf/nrutil.c b/ObservationSim/Instrument/Chip/libBF/nrutil.c
similarity index 100%
rename from ObservationSim/Instrument/Chip/lib_bf/nrutil.c
rename to ObservationSim/Instrument/Chip/libBF/nrutil.c
diff --git a/ObservationSim/Instrument/Chip/lib_bf/nrutil.h b/ObservationSim/Instrument/Chip/libBF/nrutil.h
similarity index 100%
rename from ObservationSim/Instrument/Chip/lib_bf/nrutil.h
rename to ObservationSim/Instrument/Chip/libBF/nrutil.h
diff --git a/ObservationSim/Instrument/Chip/libCTI/CTI_modeling.py b/ObservationSim/Instrument/Chip/libCTI/CTI_modeling.py
new file mode 100644
index 0000000000000000000000000000000000000000..eced09ebdce098876b67c9de3c555241163fb6d5
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/CTI_modeling.py
@@ -0,0 +1,95 @@
+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)
diff --git a/ObservationSim/Instrument/Chip/libCTI/Makefile b/ObservationSim/Instrument/Chip/libCTI/Makefile
new file mode 100755
index 0000000000000000000000000000000000000000..f58d2d14e32afaf1bff8f834fb9281679ad0d75e
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/Makefile
@@ -0,0 +1,13 @@
+#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
diff --git a/ObservationSim/Instrument/Chip/libCTI/__init__.py b/ObservationSim/Instrument/Chip/libCTI/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/ObservationSim/Instrument/Chip/libCTI/__pycache__/CTI_modeling.cpython-38.pyc b/ObservationSim/Instrument/Chip/libCTI/__pycache__/CTI_modeling.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..6e94646d36b6bf60de4c8f04a75b5049a6d9c67d
Binary files /dev/null and b/ObservationSim/Instrument/Chip/libCTI/__pycache__/CTI_modeling.cpython-38.pyc differ
diff --git a/ObservationSim/Instrument/Chip/libCTI/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Instrument/Chip/libCTI/__pycache__/__init__.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..875c6baa7e8da652f4e6bb0ea811ec50c30a3c3f
Binary files /dev/null and b/ObservationSim/Instrument/Chip/libCTI/__pycache__/__init__.cpython-38.pyc differ
diff --git a/ObservationSim/Instrument/Chip/libCTI/libmoduleCTI.so b/ObservationSim/Instrument/Chip/libCTI/libmoduleCTI.so
new file mode 100755
index 0000000000000000000000000000000000000000..0421cf067b480f3b4c1967f6afe4cb8fe05a8e8d
Binary files /dev/null and b/ObservationSim/Instrument/Chip/libCTI/libmoduleCTI.so differ
diff --git a/ObservationSim/Instrument/Chip/libCTI/readme b/ObservationSim/Instrument/Chip/libCTI/readme
new file mode 100644
index 0000000000000000000000000000000000000000..8ab15085d2330d6dfe8a034894e896e1f91d802d
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/readme
@@ -0,0 +1 @@
+CTI_lgl_v0.3/
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/add_CTI1.c b/ObservationSim/Instrument/Chip/libCTI/src/add_CTI1.c
new file mode 100644
index 0000000000000000000000000000000000000000..0bfd1c54b555ecef9480131f15a5b1449391d2f9
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/add_CTI1.c
@@ -0,0 +1,313 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "nrutil.h"
+#include <time.h>
+
+float poidev(float x, int *idum);
+void creattraps(int *seed,int ntrap,int ny,int nmax,float c,float beta,float **sp);
+void addCTI(int *a0,int ny,int noverscan,int nsp,float beta,float w, float c,float *t, float ***ntrap, int *acti, int release_seed);
+//int read_fits_2D(const char *argv,float *galaxy,int imagesize);
+//int write_fits_2D(const char *argv,float **stars,int *dim);
+
+void CTI_simul(int **image, int nx, int ny, int noverscan, int nsp, float *rho_trap, float *t, float beta, float w, float c, int nmax, int *trap_seeds, int release_seed,int **imagecti){
+  int ntotal;
+  printf("image parameters: nx=%i ny=%i noverscan=%i\n",nx,ny,noverscan);
+  printf("input image: image[0][0]=%i image[50][60]=%i image[1000][20]=%i image[20][1000]=%i\n",image[0][0],image[50][60],image[1000][20],image[20][1000]);
+  printf("rho_trap rho1=%f rho2=%f rho3=%f\n",rho_trap[0],rho_trap[1],rho_trap[2]);
+  printf("t t1=%f t2=%f t3=%f\n",t[0],t[1],t[2]);
+  printf("volume parameter beta=%f,w=%f,c=%f\n",beta,w,c);
+  printf("nsp=%i,nmax=%i\n",nsp,nmax);
+  printf("trap_seeds = %i,%i,%i\n",trap_seeds[0],trap_seeds[1],trap_seeds[2]);
+  printf("release_seed = %i\n",release_seed);
+  printf("begin CTI simulation\n");
+  float ***ntrap;
+  float **sp;
+  float tmp;
+  int ntrap_tmp;
+  int *a0,*acti,dim[2];
+  int i,j,k,l;
+  ntotal=ny+noverscan;
+  ntrap = f3tensor(0,nsp,0,ny,0,nmax+1);  
+  sp = matrix(0,ny,0,nmax+1);
+  a0=ivector(0,ny);
+  acti=ivector(0,ntotal);
+
+  for(k=0;k<nx;k++){
+    printf("column k=%i/%i\r",k,nx);
+    for(l=0;l<nsp;l++){
+        tmp = rho_trap[l]*ny;
+        ntrap_tmp = poidev(tmp,&trap_seeds[l]);
+        creattraps(&trap_seeds[l],ntrap_tmp,ny,nmax,c,beta,ntrap[l]);
+        /*for(i=0;i<ny;i++){
+	  for(j=0;j<=nmax+1;j++){///!!   j<=nmax+1  LGL
+                ntrap[l][i][j] = sp[i][j];
+            }
+	    }*/
+    }
+    //printf("ntrap0 %f,%f,%f,%f\n",ntrap[0][0][0],ntrap[0][1000][0],ntrap[0][0][1],ntrap[0][1000][1]);
+    //printf("ntrap1 %f,%f,%f,%f\n",ntrap[1][0][0],ntrap[1][1000][0],ntrap[1][0][1],ntrap[1][1000][1]);
+    //printf("ntrap2 %f,%f,%f,%f\n",ntrap[2][0][0],ntrap[2][1000][0],ntrap[2][0][1],ntrap[2][1000][1]);
+    
+    for(i=0;i<ny;i++){
+        a0[i]=image[k][i];
+        //printf("image[k][i] %i\n",image[k][i]);
+    }
+
+    for(i=0;i<ntotal;i++)acti[i]=0;
+    //printf("trap_seed_before %i,%i,%i\n",trap_seeds[0],trap_seeds[1],trap_seeds[2]);
+    addCTI(a0,ny,noverscan,nsp,beta,w,c,t,ntrap,acti,release_seed);
+    //printf("trap_seed_after %i,%i,%i\n",trap_seeds[0],trap_seeds[1],trap_seeds[2]);
+    //printf("acti %i,%i,%i\n",acti[0],acti[1],acti[2]);
+    
+    for(i=0;i<ntotal;i++){
+        imagecti[k][i]=acti[i];
+        /*if (acti[i]>acti_max){
+            acti_max = acti[i];
+	    }*/
+    }
+  }
+  printf("\nCTI simulation finished!\n");
+  free_ivector(a0,0,ny);
+  free_ivector(acti,0,ntotal);
+  free_f3tensor(ntrap,0,nsp,0,ny,0,nmax+1);
+  free_matrix(sp,0,ny,0,nmax+1);
+}
+
+void addCTI(int *a0,int ny,int noverscan, int nsp, float beta,float w, float c,float *t, float ***ntrap, int *acti, int release_seed)
+{
+  float ran1(int *idum);
+  float *f,*tmpntrap;
+  float height,wre;
+  int *flow,**traped,flowt,ntraptot,ntraped,topoftrap,nrelease;
+  int ntotal,i,j,k,nmove,isum,ntmp,isp,*tmptraped;
+  wre=1/w;
+  ntotal=ny+noverscan;
+  f=vector(0,nsp);
+  traped = imatrix(0,nsp,0,ntotal);
+  flow=ivector(0,ntotal);  
+  for(i=0;i<nsp;i++){
+    f[i] = 1-exp(-1/t[i]);
+  }
+  for(i=0;i<ntotal;i++)flow[i]=0;
+  for(i=0;i<nsp;i++){
+    for(j=0;j<ntotal;j++){
+        traped[i][j] = 0;
+    }
+  }
+  
+  
+  for(i=0;i<ny;i++){
+    flow[i]=a0[i];
+    for(j=0;j<nsp;j++){
+      ntraptot = ntrap[j][i][0];
+      height=flow[i]*wre;
+      if(ntraptot>0){
+	ntraped=0;
+	tmpntrap=ntrap[j][i];///  LGL
+	while(tmpntrap[ntraped+1]<=height){ntraped++;}
+	if(flow[i]<ntraped)ntraped=flow[i];
+	flow[i]-=ntraped; // LGL
+	traped[j][i]=ntraped;
+      }
+    }
+  }
+  
+  for(i=0;i<ntotal;i++){
+    acti[i] = flow[0];
+    nmove=ny;
+    if(i>noverscan)nmove=ntotal-i;
+    for(j=0;j<nmove;j++){
+      flow[j]=flow[j+1];
+      height=flow[j]*wre;
+      //for(isp=0;isp<nsp;isp++){
+      isp=0;
+	  ntraptot=ntrap[isp][j][0];
+	  if(ntraptot>0){
+	    ntraped=0; 
+	    //height=flow[j]*wre; 
+	    tmpntrap=ntrap[isp][j];
+	    tmptraped=traped[isp];
+	    while(tmpntrap[ntraped+1]<=height){ntraped++;}
+	    topoftrap=tmptraped[j]-ntraped;
+	    if(topoftrap==0){} //do nothing
+	    else if (topoftrap>0){  // release
+	      for(k=0;k<topoftrap;k++){
+	        if(ran1(&release_seed)<f[isp]){
+		        flow[j]++;
+		        tmptraped[j]--;
+	        }
+	      }
+	    }
+	    else{     //trap
+	        ntmp = -topoftrap;
+	        if(flow[j]<ntmp)ntmp=flow[j];
+	        flow[j]-=ntmp; 
+	        tmptraped[j] += ntmp;
+	    }
+	  }
+	  if(nsp==1)continue;  
+	  
+	  isp=1;
+	  ntraptot=ntrap[isp][j][0];
+	  if(ntraptot>0){
+	    ntraped=0; 
+	    //height=flow[j]*wre; 
+	    tmpntrap=ntrap[isp][j];
+	    tmptraped=traped[isp];
+	    while(tmpntrap[ntraped+1]<=height){ntraped++;}
+	    topoftrap=tmptraped[j]-ntraped;
+	    if(topoftrap==0){} //do nothing
+	    else if (topoftrap>0){  // release
+	        for(k=0;k<topoftrap;k++){
+	            if(ran1(&release_seed)<f[isp]){
+		            flow[j]++;
+		            tmptraped[j]--;
+	            }
+	        }
+	    }
+	    else{     //trap
+	        ntmp = -topoftrap;
+	        if(flow[j]<ntmp)ntmp=flow[j];
+	        flow[j]-=ntmp; 
+	        tmptraped[j] += ntmp;
+	    }
+	  }
+	  if(nsp==2)continue; 
+	  
+	  isp=2;
+	  ntraptot=ntrap[isp][j][0];
+	  if(ntraptot>0){
+	    ntraped=0; 
+	    //height=flow[j]*wre; 
+	    tmpntrap=ntrap[isp][j];
+	    tmptraped=traped[isp];
+	    while(tmpntrap[ntraped+1]<=height){ntraped++;}
+	    topoftrap=tmptraped[j]-ntraped;
+	    if(topoftrap==0){} //do nothing
+	    else if (topoftrap>0){  // release
+	      for(k=0;k<topoftrap;k++){
+	        if(ran1(&release_seed)<f[isp]){
+		        flow[j]++;
+		        tmptraped[j]--;
+	        }
+	      }
+	    }
+	    else{     //trap
+	        ntmp = -topoftrap;
+	        if(flow[j]<ntmp)ntmp=flow[j];
+	        flow[j]-=ntmp; 
+	        tmptraped[j] += ntmp;
+	    }
+	  }
+	  if(nsp==3)continue;
+
+	  isp=3;
+	  ntraptot=ntrap[isp][j][0];
+	  if(ntraptot>0){
+	    ntraped=0; 
+	    //height=flow[j]*wre; 
+	    tmpntrap=ntrap[isp][j];
+	    tmptraped=traped[isp];
+	    while(tmpntrap[ntraped+1]<=height){ntraped++;}
+	    topoftrap=tmptraped[j]-ntraped;
+	    if(topoftrap==0){} //do nothing
+	    else if (topoftrap>0){  // release
+	      for(k=0;k<topoftrap;k++){
+	        if(ran1(&release_seed)<f[isp]){
+		        flow[j]++;
+		        tmptraped[j]--;
+	        }
+	      }
+	    }
+	    else{     //trap
+	        ntmp = -topoftrap;
+	        if(flow[j]<ntmp)ntmp=flow[j];
+	        flow[j]-=ntmp; 
+	        tmptraped[j] += ntmp;
+	    }
+	  }
+	  if(nsp==4)continue;
+	
+	//}
+    }
+  }
+  free_ivector(flow,0,ntotal);
+  free_vector(f,0,nsp-1);
+  free_imatrix(traped,0,nsp-1,0,ntotal);
+}
+
+void save_trap_map(int *trap_seeds, int nsp, int nx, int ny, int nmax, float *rho_trap, float beta, float c, char *filename){
+    printf("saving trap map rho_trap = %f,%f,%f, c = %f, beta = %f, nmax = %i, trap_seed = %i,%i,%i\n",rho_trap[0],rho_trap[1],rho_trap[2],c,beta,nmax,trap_seeds[0],trap_seeds[1],trap_seeds[2]);
+    float tmp;
+    int ntrap_tmp;
+    float **sp;
+    float hsp;
+    sp = matrix(0,ny,0,nmax);
+    FILE *file = fopen(filename,"wb");
+    for(int j=0;j<nx;j++){
+        for(int l=0;l<nsp;l++){
+            tmp = rho_trap[l]*ny;
+            ntrap_tmp = poidev(tmp,&trap_seeds[l]);
+            creattraps(&trap_seeds[l],ntrap_tmp,ny,nmax,c,beta,sp);
+            for(int i=0;i<ny;i++){
+                for(int k=0;k<nmax;k++){
+                    hsp = sp[i][k];
+                    unsigned char *hsp_char = (unsigned char *)&hsp;
+                    unsigned char bytes[4];
+                    bytes[0] = hsp_char[0];
+                    bytes[1] = hsp_char[1];
+                    bytes[2] = hsp_char[2];
+                    bytes[3] = hsp_char[3];
+                    fwrite(bytes,sizeof(unsigned char),4,file);
+                }
+            }
+        }
+    }    
+    fclose(file);
+    free_matrix(sp,0,ny,0,nmax);
+}
+
+/*int main(){
+  int read_fits_2D(const char *argv,float *galaxy,int imagesize);
+  int write_fits_2D(const char *argv,float **stars,int *dim);
+  int **image,nimg,kct;
+  float *image1,**imagecti;
+  int nx=4608,ny=4616,noverscan=84,nmax=10,nsp=3;
+  int ntotal = ny+noverscan;
+  float beta=0.478,w=84700,c=0;
+  float t[3] = {0.74,7.7,37};
+  float rho_trap[3] = {0.6,1.6,1.4};
+  char fname1[300];
+  int dim[2];
+  int trap_seeds[3] = {1,100,1000};
+  int release_seed = 10000;
+  clock_t start, end;
+  
+  nimg=nx*ny;
+  image=imatrix(0,nx,0,ny);
+  image1=vector(0,nimg);
+  strcpy(fname1,"inputdata/image.fits");
+  read_fits_2D(fname1,image1,nimg);kct=0;    
+  for (int j=0; j<ny; j++){for(int i=0;i<nx;i++){
+      image[i][j]=image1[kct];
+      kct++;
+    }}
+  
+  start = clock();
+  imagecti=matrix(0,nx,0,ntotal);
+  CTI_simul(image, nx, ny, noverscan, nsp, rho_trap, t, beta, w, c, nmax, trap_seeds, release_seed, imagecti);
+  end = clock();
+  float time;
+  time = (end-start)/CLOCKS_PER_SEC;
+  printf("running time %f",time);
+  
+  // output final image
+  strcpy(fname1,"output/image_CTI.fits");
+  dim[0]=nx;dim[1]=ntotal;
+  write_fits_2D(fname1,imagecti,dim);
+  free_vector(image1,0,nimg);
+  free_matrix(imagecti,0,nx,0,ntotal);
+  return 0;
+  }*/
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/creattraps.c b/ObservationSim/Instrument/Chip/libCTI/src/creattraps.c
new file mode 100644
index 0000000000000000000000000000000000000000..0922be11e00180a87a74a40a4262c0ca33e7e083
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/creattraps.c
@@ -0,0 +1,52 @@
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "nrutil.h"
+void creattraps(long *seed,int ntrap,int ny,int nmax,float c,float beta,float **sp)
+{
+  //creat ntrap traps along one column
+  //sp[i][0] the number of trap in the ith pixle;
+  //sp[i][j] (c,1+c) the height of the jth trap in the ith pixle;
+  
+  float ran1(long *idum);
+  void sort(unsigned long n, float arr[]);
+  float tmp,betare,height;
+  int i,nyt,ntmp,j;
+  float *tmpv;
+  tmpv = vector(1,100);
+
+  if(nmax>50){printf(" the upper limite of trap in each pixe is too large, nmax=%d\n",nmax); getchar();}
+
+  
+  betare=1./beta;
+  for(i=0;i<ny;i++)sp[i][0]=0;
+  for(i=0;i<ntrap;i++){
+    nyt=ran1(seed)*ny;
+    sp[nyt][0]++;
+  }
+  for(i=0;i<ny;i++){
+    ntmp=sp[i][0];
+    if(ntmp==1){
+      height=ran1(seed)-c;
+      if(height<=0){sp[i][1]=0;}
+      else{
+        sp[i][1]=pow(height,betare);
+      }
+    }
+    if(ntmp>1){
+      if(ntmp>nmax)sp[i][0]=ntmp=nmax; // upper limite of trap in each pixel is nmax
+      for(j=1;j<=ntmp;j++)tmpv[j]=ran1(seed)-c;
+      
+      sort(ntmp, tmpv);
+      for(j=1;j<=ntmp;j++){
+	    height=tmpv[j];
+	    if(height<=0){sp[i][j]=0;}
+	    else{sp[i][j]=pow(height,betare);}
+      }
+    }
+    sp[i][ntmp+1]=999999;
+  }
+  free_vector(tmpv,1,100);
+}
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/gammln.c b/ObservationSim/Instrument/Chip/libCTI/src/gammln.c
new file mode 100644
index 0000000000000000000000000000000000000000..0511563e06c29a34551e7c989fe93cf56cbf52d9
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/gammln.c
@@ -0,0 +1,18 @@
+#include <math.h>
+
+float gammln(float xx)
+{
+	double x,y,tmp,ser;
+	static double cof[6]={76.18009172947146,-86.50532032941677,
+		24.01409824083091,-1.231739572450155,
+		0.1208650973866179e-2,-0.5395239384953e-5};
+	int j;
+
+	y=x=xx;
+	tmp=x+5.5;
+	tmp -= (x+0.5)*log(tmp);
+	ser=1.000000000190015;
+	for (j=0;j<=5;j++) ser += cof[j]/++y;
+	return -tmp+log(2.5066282746310005*ser/x);
+}
+/* (C) Copr. 1986-92 Numerical Recipes Software )1!. */
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/gasdev.c b/ObservationSim/Instrument/Chip/libCTI/src/gasdev.c
new file mode 100644
index 0000000000000000000000000000000000000000..26926a38fd55bebaa3a756de9d1cca3ef9815cac
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/gasdev.c
@@ -0,0 +1,25 @@
+#include <math.h>
+
+float gasdev(long *idum)
+{
+	float ran1(long *idum);
+	static int iset=0;
+	static float gset;
+	float fac,rsq,v1,v2;
+
+	if  (iset == 0) {
+		do {
+			v1=2.0*ran1(idum)-1.0;
+			v2=2.0*ran1(idum)-1.0;
+			rsq=v1*v1+v2*v2;
+		} while (rsq >= 1.0 || rsq == 0.0);
+		fac=sqrt(-2.0*log(rsq)/rsq);
+		gset=v1*fac;
+		iset=1;
+		return v2*fac;
+	} else {
+		iset=0;
+		return gset;
+	}
+}
+/* (C) Copr. 1986-92 Numerical Recipes Software )1!. */
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/nrutil.c b/ObservationSim/Instrument/Chip/libCTI/src/nrutil.c
new file mode 100644
index 0000000000000000000000000000000000000000..8ecacf5d37563ab1848bf52a81d2f6ce955a708b
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/nrutil.c
@@ -0,0 +1,769 @@
+#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
+
+#include <stdio.h>
+#include <stddef.h>
+#include <stdlib.h>
+#define NR_END 1
+#define FREE_ARG char*
+
+void nrerror(char error_text[])
+/* Numerical Recipes standard error handler */
+{
+	fprintf(stderr,"Numerical Recipes run-time error...\n");
+	fprintf(stderr,"%s\n",error_text);
+	fprintf(stderr,"...now exiting to system...\n");
+	exit(1);
+}
+
+float *vector(long nl, long nh)
+/* allocate a float vector with subscript range v[nl..nh] */
+{
+	float *v;
+
+	v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
+	if (!v) nrerror("allocation failure in vector()");
+	return v-nl+NR_END;
+}
+
+int *ivector(long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+	int *v;
+
+	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
+	if (!v) nrerror("allocation failure in ivector()");
+	return v-nl+NR_END;
+}
+
+unsigned char *cvector(long nl, long nh)
+/* allocate an unsigned char vector with subscript range v[nl..nh] */
+{
+	unsigned char *v;
+
+	v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
+	if (!v) nrerror("allocation failure in cvector()");
+	return v-nl+NR_END;
+}
+
+long *lvector(long nl, long nh)
+/* allocate an long vector with subscript range v[nl..nh] */
+{
+	long *v;
+
+	v=(long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
+	if (!v) nrerror("allocation failure in lvector()");
+	return v-nl+NR_END;
+}
+
+double *dvector(long nl, long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+	double *v;
+
+	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
+	if (!v) nrerror("allocation failure in dvector()");
+	return v-nl+NR_END;
+}
+
+float **matrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+double **dmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+int **imatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	int **m;
+
+	/* allocate pointers to rows */
+	m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+	long newrl, long newcl)
+/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
+{
+	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
+	float **m;
+
+	/* allocate array of pointers to rows */
+	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in submatrix()");
+	m += NR_END;
+	m -= newrl;
+
+	/* set pointers to rows */
+	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
+declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
+and ncol=nch-ncl+1. The routine should be called with the address
+&a[0][0] as the first argument. */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in convert_matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* set pointers to rows */
+	m[nrl]=a-ncl;
+	for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+	float ***t;
+
+	/* allocate pointers to pointers to rows */
+	t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
+	if (!t) nrerror("allocation failure 1 in f3tensor()");
+	t += NR_END;
+	t -= nrl;
+
+	/* allocate pointers to rows and set pointers to them */
+	t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
+	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
+	t[nrl] += NR_END;
+	t[nrl] -= ncl;
+
+	/* allocate rows and set pointers to them */
+	t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
+	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
+	t[nrl][ncl] += NR_END;
+	t[nrl][ncl] -= ndl;
+
+	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+	for(i=nrl+1;i<=nrh;i++) {
+		t[i]=t[i-1]+ncol;
+		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+	}
+
+	/* return pointer to array of pointers to rows */
+	return t;
+}
+double ***d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+	double ***t;
+
+	/* allocate pointers to pointers to rows */
+	t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
+	if (!t) nrerror("allocation failure 1 in f3tensor()");
+	t += NR_END;
+	t -= nrl;
+
+	/* allocate pointers to rows and set pointers to them */
+	t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
+	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
+	t[nrl] += NR_END;
+	t[nrl] -= ncl;
+
+	/* allocate rows and set pointers to them */
+	t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
+	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
+	t[nrl][ncl] += NR_END;
+	t[nrl][ncl] -= ndl;
+
+	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+	for(i=nrl+1;i<=nrh;i++) {
+		t[i]=t[i-1]+ncol;
+		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+	}
+
+	/* return pointer to array of pointers to rows */
+	return t;
+}
+
+
+char **cmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+        char **m;
+
+	/* allocate pointers to rows */
+	m=(char **) malloc((size_t)((nrow+NR_END)*sizeof(char*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(char *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(char)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+void free_vector(float *v, long nl, long nh)
+/* free a float vector allocated with vector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_cvector(unsigned char *v, long nl, long nh)
+/* free an unsigned char vector allocated with cvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_lvector(long *v, long nl, long nh)
+/* free an long vector allocated with lvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_dvector(double *v, long nl, long nh)
+/* free a double vector allocated with dvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
+/* free a float matrix allocated by matrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by dmatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
+/* free an int matrix allocated by imatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a submatrix allocated by submatrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a matrix allocated by convert_matrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+	long ndl, long ndh)
+/* free a float f3tensor allocated by f3tensor() */
+{
+	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
+	free((FREE_ARG) (t[nrl]+ncl-NR_END));
+	free((FREE_ARG) (t+nrl-NR_END));
+}
+void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
+	long ndl, long ndh)
+/* free a double f3tensor allocated by f3tensor() */
+{
+	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
+	free((FREE_ARG) (t[nrl]+ncl-NR_END));
+	free((FREE_ARG) (t+nrl-NR_END));
+}
+void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch)
+/* free a character matrix allocated by matrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+
+#else /* ANSI */
+/* traditional - K&R */
+
+#include <stdio.h>
+#define NR_END 1
+#define FREE_ARG char*
+
+void nrerror(error_text)
+char error_text[];
+/* Numerical Recipes standard error handler */
+{
+	void exit();
+
+	fprintf(stderr,"Numerical Recipes run-time error...\n");
+	fprintf(stderr,"%s\n",error_text);
+	fprintf(stderr,"...now exiting to system...\n");
+	exit(1);
+}
+
+float *vector(nl,nh)
+long nh,nl;
+/* allocate a float vector with subscript range v[nl..nh] */
+{
+	float *v;
+
+	v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float)));
+	if (!v) nrerror("allocation failure in vector()");
+	return v-nl+NR_END;
+}
+
+int *ivector(nl,nh)
+long nh,nl;
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+	int *v;
+
+	v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int)));
+	if (!v) nrerror("allocation failure in ivector()");
+	return v-nl+NR_END;
+}
+
+unsigned char *cvector(nl,nh)
+long nh,nl;
+/* allocate an unsigned char vector with subscript range v[nl..nh] */
+{
+	unsigned char *v;
+
+	v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
+	if (!v) nrerror("allocation failure in cvector()");
+	return v-nl+NR_END;
+}
+
+long *lvector(nl,nh)
+long nh,nl;
+/* allocate an unsigned long vector with subscript range v[nl..nh] */
+{
+        long *v;
+
+	v=(long *)malloc((int) ((nh-nl+1+NR_END)*sizeof(long)));
+	if (!v) nrerror("allocation failure in lvector()");
+	return v-nl+NR_END;
+}
+
+double *dvector(nl,nh)
+long nh,nl;
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+	double *v;
+
+	v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
+	if (!v) nrerror("allocation failure in dvector()");
+	return v-nl+NR_END;
+}
+
+float **matrix(nrl,nrh,ncl,nch)
+long nch,ncl,nrh,nrl;
+/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((unsigned int)((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(float *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+double **dmatrix(nrl,nrh,ncl,nch)
+long nch,ncl,nrh,nrl;
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+int **imatrix(nrl,nrh,ncl,nch)
+long nch,ncl,nrh,nrl;
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	int **m;
+
+	/* allocate pointers to rows */
+	m=(int **) malloc((unsigned int)((nrow+NR_END)*sizeof(int*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(int *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(int)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
+float **a;
+long newcl,newrl,oldch,oldcl,oldrh,oldrl;
+/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
+{
+	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
+	float **m;
+
+	/* allocate array of pointers to rows */
+	m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in submatrix()");
+	m += NR_END;
+	m -= newrl;
+
+	/* set pointers to rows */
+	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float **convert_matrix(a,nrl,nrh,ncl,nch)
+float *a;
+long nch,ncl,nrh,nrl;
+/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
+declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
+and ncol=nch-ncl+1. The routine should be called with the address
+&a[0][0] as the first argument. */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
+	if (!m)	nrerror("allocation failure in convert_matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* set pointers to rows */
+	m[nrl]=a-ncl;
+	for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+double ***d3tensor(nrl,nrh,ncl,nch,ndl,ndh)
+long nch,ncl,ndh,ndl,nrh,nrl;
+/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+	double ***t;
+
+	/* allocate pointers to pointers to rows */
+	t=(double ***) malloc((unsigned int)((nrow+NR_END)*sizeof(double**)));
+	if (!t) nrerror("allocation failure 1 in f3tensor()");
+	t += NR_END;
+	t -= nrl;
+
+	/* allocate pointers to rows and set pointers to them */
+	t[nrl]=(double **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double*)));
+	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
+	t[nrl] += NR_END;
+	t[nrl] -= ncl;
+
+	/* allocate rows and set pointers to them */
+	t[nrl][ncl]=(double *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(double)));
+	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
+	t[nrl][ncl] += NR_END;
+	t[nrl][ncl] -= ndl;
+
+	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+	for(i=nrl+1;i<=nrh;i++) {
+		t[i]=t[i-1]+ncol;
+		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+	}
+
+	/* return pointer to array of pointers to rows */
+	return t;
+}
+
+float ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
+long nch,ncl,ndh,ndl,nrh,nrl;
+/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+	float ***t;
+
+	/* allocate pointers to pointers to rows */
+	t=(float ***) malloc((unsigned int)((nrow+NR_END)*sizeof(float**)));
+	if (!t) nrerror("allocation failure 1 in f3tensor()");
+	t += NR_END;
+	t -= nrl;
+
+	/* allocate pointers to rows and set pointers to them */
+	t[nrl]=(float **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float*)));
+	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
+	t[nrl] += NR_END;
+	t[nrl] -= ncl;
+
+	/* allocate rows and set pointers to them */
+	t[nrl][ncl]=(float *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(float)));
+	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
+	t[nrl][ncl] += NR_END;
+	t[nrl][ncl] -= ndl;
+
+	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+	for(i=nrl+1;i<=nrh;i++) {
+		t[i]=t[i-1]+ncol;
+		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+	}
+
+	/* return pointer to array of pointers to rows */
+	return t;
+}
+
+char **cmatrix(nrl,nrh,ncl,nch)
+long nch,ncl,nrh,nrl;
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	char **m;
+
+	/* allocate pointers to rows */
+	m=(char **) malloc((unsigned int)((nrow+NR_END)*sizeof(char*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(char *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(char)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+
+
+void free_vector(v,nl,nh)
+float *v;
+long nh,nl;
+/* free a float vector allocated with vector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_ivector(v,nl,nh)
+int *v;
+long nh,nl;
+/* free an int vector allocated with ivector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_cvector(v,nl,nh)
+long nh,nl;
+unsigned char *v;
+/* free an unsigned char vector allocated with cvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_lvector(v,nl,nh)
+long nh,nl;
+long *v;
+/* free an long vector allocated with lvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_dvector(v,nl,nh)
+double *v;
+long nh,nl;
+/* free a double vector allocated with dvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_matrix(m,nrl,nrh,ncl,nch)
+float **m;
+long nch,ncl,nrh,nrl;
+/* free a float matrix allocated by matrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_dmatrix(m,nrl,nrh,ncl,nch)
+double **m;
+long nch,ncl,nrh,nrl;
+/* free a double matrix allocated by dmatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_imatrix(m,nrl,nrh,ncl,nch)
+int **m;
+long nch,ncl,nrh,nrl;
+/* free an int matrix allocated by imatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_submatrix(b,nrl,nrh,ncl,nch)
+float **b;
+long nch,ncl,nrh,nrl;
+/* free a submatrix allocated by submatrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_convert_matrix(b,nrl,nrh,ncl,nch)
+float **b;
+long nch,ncl,nrh,nrl;
+/* free a matrix allocated by convert_matrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
+float ***t;
+long nch,ncl,ndh,ndl,nrh,nrl;
+/* free a float f3tensor allocated by f3tensor() */
+{
+	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
+	free((FREE_ARG) (t[nrl]+ncl-NR_END));
+	free((FREE_ARG) (t+nrl-NR_END));
+}
+void free_d3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
+double ***t;
+long nch,ncl,ndh,ndl,nrh,nrl;
+/* free a float f3tensor allocated by f3tensor() */
+{
+	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
+	free((FREE_ARG) (t[nrl]+ncl-NR_END));
+	free((FREE_ARG) (t+nrl-NR_END));
+}
+void free_cmatrix(m,nrl,nrh,ncl,nch)
+char **m;
+long nch,ncl,nrh,nrl;
+/* free a double matrix allocated by dmatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+#endif /* ANSI */
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/nrutil.h b/ObservationSim/Instrument/Chip/libCTI/src/nrutil.h
new file mode 100644
index 0000000000000000000000000000000000000000..48d23f345ebd1d8643ab9f82578a6c1d8c02415a
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/nrutil.h
@@ -0,0 +1,94 @@
+/* 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_ */
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/poidev.c b/ObservationSim/Instrument/Chip/libCTI/src/poidev.c
new file mode 100644
index 0000000000000000000000000000000000000000..cc3c50b888b50fbb65b5debd0a379870d8e7138f
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/poidev.c
@@ -0,0 +1,79 @@
+#include <math.h>
+#define PI 3.141592654
+
+float poidev(float xm, int *idum)
+{
+	float gammln(float xx);
+	float ran1(int *idum);
+	static float sq,alxm,g,oldm=(-1.0);
+	float em,t,y;
+
+	if (xm < 12.0) {
+		if (xm != oldm) {
+			oldm=xm;
+			g=exp(-xm);
+		}
+		em = -1;
+		t=1.0;
+		do {
+			++em;
+			t *= ran1(idum);
+		} while (t > g);
+	} else {
+		if (xm != oldm) {
+			oldm=xm;
+			sq=sqrt(2.0*xm);
+			alxm=log(xm);
+			g=xm*alxm-gammln(xm+1.0);
+		}
+		do {
+			do {
+				y=tan(PI*ran1(idum));
+				em=sq*y+xm;
+			} while (em < 0.0);
+			em=floor(em);
+			t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
+		} while (ran1(idum) > t);
+	}
+	return em;
+}
+
+float poidev_copy(float xm, int *idum)
+{
+	float gammln(float xx);
+	float ran1_copy(int *idum);
+	static float sq1,alxm1,g1,oldm1=(-1.0);
+	float em,t,y;
+
+	if (xm < 12.0) {
+		if (xm != oldm1) {
+			oldm1=xm;
+			g1=exp(-xm);
+		}
+		em = -1;
+		t=1.0;
+		do {
+			++em;
+			t *= ran1_copy(idum);
+		} while (t > g1);
+	} else {
+		if (xm != oldm1) {
+			oldm1=xm;
+			sq1=sqrt(2.0*xm);
+			alxm1=log(xm);
+			g1=xm*alxm1-gammln(xm+1.0);
+		}
+		do {
+			do {
+				y=tan(PI*ran1_copy(idum));
+				em=sq1*y+xm;
+			} while (em < 0.0);
+			em=floor(em);
+			t=0.9*(1.0+y*y)*exp(em*alxm1-gammln(em+1.0)-g1);
+		} while (ran1_copy(idum) > t);
+	}
+	return em;
+}
+
+#undef PI
+/* (C) Copr. 1986-92 Numerical Recipes Software )1!. */
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/ran1.c b/ObservationSim/Instrument/Chip/libCTI/src/ran1.c
new file mode 100644
index 0000000000000000000000000000000000000000..7f7bb3d302dde2fc87e06309b8fc54f6ad8d3248
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/ran1.c
@@ -0,0 +1,82 @@
+#define IA 16807
+#define IM 2147483647
+#define AM (1.0/IM)
+#define IQ 127773
+#define IR 2836
+#define NTAB 32
+#define NDIV (1+(IM-1)/NTAB)
+#define EPS 1.2e-7
+#define RNMX (1.0-EPS)
+
+float ran1(int *idum)
+{
+	int j;
+	int k;
+	//static
+	static int iy=0;
+	static int iv[NTAB];
+	
+	float temp;
+
+	if (*idum <= 0 || !iy) {
+		if (-(*idum) < 1) *idum=1;
+		else *idum = -(*idum);
+		for (j=NTAB+7;j>=0;j--) {
+			k=(*idum)/IQ;
+			*idum=IA*(*idum-k*IQ)-IR*k;
+			if (*idum < 0) *idum += IM;
+			if (j < NTAB) iv[j] = *idum;
+		}
+		iy=iv[0];
+	}
+	k=(*idum)/IQ;
+	*idum=IA*(*idum-k*IQ)-IR*k;
+	if (*idum < 0) *idum += IM;
+	j=iy/NDIV;
+	iy=iv[j];
+	iv[j] = *idum;
+	if ((temp=AM*iy) > RNMX) return RNMX;
+	else return temp;
+}
+
+float ran1_copy(int *idum)
+{
+	int j;
+	int k;
+	//static
+	static int iy1=0;
+	static int iv1[NTAB];
+	
+	float temp;
+
+	if (*idum <= 0 || !iy1) {
+		if (-(*idum) < 1) *idum=1;
+		else *idum = -(*idum);
+		for (j=NTAB+7;j>=0;j--) {
+			k=(*idum)/IQ;
+			*idum=IA*(*idum-k*IQ)-IR*k;
+			if (*idum < 0) *idum += IM;
+			if (j < NTAB) iv1[j] = *idum;
+		}
+		iy1=iv1[0];
+	}
+	k=(*idum)/IQ;
+	*idum=IA*(*idum-k*IQ)-IR*k;
+	if (*idum < 0) *idum += IM;
+	j=iy1/NDIV;
+	iy1=iv1[j];
+	iv1[j] = *idum;
+	if ((temp=AM*iy1) > RNMX) return RNMX;
+	else return temp;
+}
+
+#undef IA
+#undef IM
+#undef AM
+#undef IQ
+#undef IR
+#undef NTAB
+#undef NDIV
+#undef EPS
+#undef RNMX
+/* (C) Copr. 1986-92 Numerical Recipes Software )1!. */
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/ran2.c b/ObservationSim/Instrument/Chip/libCTI/src/ran2.c
new file mode 100644
index 0000000000000000000000000000000000000000..5dce1f0ae9f815eed4df7416c394bb6d97c1ca1e
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/ran2.c
@@ -0,0 +1,64 @@
+#define IM1 2147483563
+#define IM2 2147483399
+#define AM (1.0/IM1)
+#define IMM1 (IM1-1)
+#define IA1 40014
+#define IA2 40692
+#define IQ1 53668
+#define IQ2 52774
+#define IR1 12211
+#define IR2 3791
+#define NTAB 32
+#define NDIV (1+IMM1/NTAB)
+#define EPS 1.2e-7
+#define RNMX (1.0-EPS)
+
+float ran2(long *idum)
+{
+	int j;
+	long k;
+	static long idum2=123456789;
+	static long iy=0;
+	static long iv[NTAB];
+	float temp;
+
+	if (*idum <= 0) {
+		if (-(*idum) < 1) *idum=1;
+		else *idum = -(*idum);
+		idum2=(*idum);
+		for (j=NTAB+7;j>=0;j--) {
+			k=(*idum)/IQ1;
+			*idum=IA1*(*idum-k*IQ1)-k*IR1;
+			if (*idum < 0) *idum += IM1;
+			if (j < NTAB) iv[j] = *idum;
+		}
+		iy=iv[0];
+	}
+	k=(*idum)/IQ1;
+	*idum=IA1*(*idum-k*IQ1)-k*IR1;
+	if (*idum < 0) *idum += IM1;
+	k=idum2/IQ2;
+	idum2=IA2*(idum2-k*IQ2)-k*IR2;
+	if (idum2 < 0) idum2 += IM2;
+	j=iy/NDIV;
+	iy=iv[j]-idum2;
+	iv[j] = *idum;
+	if (iy < 1) iy += IMM1;
+	if ((temp=AM*iy) > RNMX) return RNMX;
+	else return temp;
+}
+#undef IM1
+#undef IM2
+#undef AM
+#undef IMM1
+#undef IA1
+#undef IA2
+#undef IQ1
+#undef IQ2
+#undef IR1
+#undef IR2
+#undef NTAB
+#undef NDIV
+#undef EPS
+#undef RNMX
+/* (C) Copr. 1986-92 Numerical Recipes Software )1!. */
diff --git a/ObservationSim/Instrument/Chip/libCTI/src/sort.c b/ObservationSim/Instrument/Chip/libCTI/src/sort.c
new file mode 100644
index 0000000000000000000000000000000000000000..4dfa8671f77ec2457e806367aa0449ede6583468
--- /dev/null
+++ b/ObservationSim/Instrument/Chip/libCTI/src/sort.c
@@ -0,0 +1,69 @@
+#define NRANSI
+#include "nrutil.h"
+#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
+#define M 7
+#define NSTACK 50
+
+void sort(unsigned long n, float arr[])
+{
+	unsigned long i,ir=n,j,k,l=1;
+	int jstack=0,*istack;
+	float a,temp;
+
+	istack=ivector(1,NSTACK);
+	for (;;) {
+		if (ir-l < M) {
+			for (j=l+1;j<=ir;j++) {
+				a=arr[j];
+				for (i=j-1;i>=1;i--) {
+					if (arr[i] <= a) break;
+					arr[i+1]=arr[i];
+				}
+				arr[i+1]=a;
+			}
+			if (jstack == 0) break;
+			ir=istack[jstack--];
+			l=istack[jstack--];
+		} else {
+			k=(l+ir) >> 1;
+			SWAP(arr[k],arr[l+1])
+			if (arr[l+1] > arr[ir]) {
+				SWAP(arr[l+1],arr[ir])
+			}
+			if (arr[l] > arr[ir]) {
+				SWAP(arr[l],arr[ir])
+			}
+			if (arr[l+1] > arr[l]) {
+				SWAP(arr[l+1],arr[l])
+			}
+			i=l+1;
+			j=ir;
+			a=arr[l];
+			for (;;) {
+				do i++; while (arr[i] < a);
+				do j--; while (arr[j] > a);
+				if (j < i) break;
+				SWAP(arr[i],arr[j]);
+			}
+			arr[l]=arr[j];
+			arr[j]=a;
+			jstack += 2;
+			if (jstack > NSTACK) nrerror("NSTACK too small in sort.");
+			if (ir-i+1 >= j-l) {
+				istack[jstack]=ir;
+				istack[jstack-1]=i;
+				ir=j-1;
+			} else {
+				istack[jstack]=j-1;
+				istack[jstack-1]=l;
+				l=i;
+			}
+		}
+	}
+	free_ivector(istack,1,NSTACK);
+}
+#undef M
+#undef NSTACK
+#undef SWAP
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software )1!. */
diff --git a/ObservationSim/Instrument/data/ccd/chip_definition.json b/ObservationSim/Instrument/data/ccd/chip_definition.json
index da04b866fc49054594648777aa2443f721ecee63..f50d90cf76f163c6a1eb851ff3983ec69757dc07 100644
--- a/ObservationSim/Instrument/data/ccd/chip_definition.json
+++ b/ObservationSim/Instrument/data/ccd/chip_definition.json
@@ -15,7 +15,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "32": {
         "chip_name": "FGS1A-D2",
@@ -33,7 +37,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "33": {
         "chip_name": "FGS1B-D1",
@@ -51,7 +59,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "34": {
         "chip_name": "FGS1B-D2",
@@ -69,7 +81,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "35": {
         "chip_name": "FGS2A-D1",
@@ -87,7 +103,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "36": {
         "chip_name": "FGS2A-D2",
@@ -105,7 +125,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "37": {
         "chip_name": "FGS2B-D1",
@@ -123,7 +147,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "38": {
         "chip_name": "FGS2B-D2",
@@ -141,7 +169,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "39": {
         "chip_name": "FGS3-D1",
@@ -159,7 +191,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "40": {
         "chip_name": "FGS3-D2",
@@ -177,7 +213,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "41": {
         "chip_name": "FGS4-D1",
@@ -195,7 +235,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "42": {
         "chip_name": "FGS4-D2",
@@ -213,7 +257,11 @@
         "df_strength": 2.3,
         "bias_level": 2000.0,
         "gain": 1.0,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "1": {
         "chip_name": "GI-1",
@@ -231,7 +279,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "2": {
         "chip_name": "GV-1",
@@ -249,7 +301,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "3": {
         "chip_name": "GU-1",
@@ -267,7 +323,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "4": {
         "chip_name": "GU-2",
@@ -285,7 +345,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "5": {
         "chip_name": "GV-2",
@@ -303,7 +367,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "6": {
         "chip_name": "y-1",
@@ -321,7 +389,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "7": {
         "chip_name": "i-1",
@@ -339,7 +411,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "8": {
         "chip_name": "g-1",
@@ -357,7 +433,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "9": {
         "chip_name": "r-1",
@@ -375,7 +455,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "10": {
         "chip_name": "GI-2",
@@ -393,7 +477,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "11": {
         "chip_name": "z-1",
@@ -411,7 +499,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "12": {
         "chip_name": "NUV-1",
@@ -429,7 +521,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "13": {
         "chip_name": "NUV-2",
@@ -447,7 +543,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "14": {
         "chip_name": "u-1",
@@ -465,7 +565,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "15": {
         "chip_name": "y-2",
@@ -483,7 +587,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "16": {
         "chip_name": "y-3",
@@ -501,7 +609,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "17": {
         "chip_name": "u-2",
@@ -519,7 +631,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "18": {
         "chip_name": "NUV-3",
@@ -537,7 +653,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "19": {
         "chip_name": "NUV-4",
@@ -555,7 +675,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "20": {
         "chip_name": "z-2",
@@ -573,7 +697,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "21": {
         "chip_name": "GI-3",
@@ -591,7 +719,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "22": {
         "chip_name": "r-2",
@@ -609,7 +741,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "23": {
         "chip_name": "g-2",
@@ -627,7 +763,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "24": {
         "chip_name": "i-2",
@@ -645,7 +785,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "25": {
         "chip_name": "y-4",
@@ -663,7 +807,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "26": {
         "chip_name": "GV-3",
@@ -681,7 +829,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "27": {
         "chip_name": "GU-3",
@@ -699,7 +851,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "28": {
         "chip_name": "GU-4",
@@ -717,7 +873,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "29": {
         "chip_name": "GV-4",
@@ -735,7 +895,11 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     },
     "30": {
         "chip_name": "GI-4",
@@ -753,6 +917,10 @@
         "df_strength": 2.3,
         "bias_level": 500,
         "gain": 1.1,
-        "full_well": 90000
+        "full_well": 90000,
+        "prescan_x": 27,
+        "overscan_x": 71,
+        "prescan_y": 0,
+        "overscan_y": 84
     }
-}
\ No newline at end of file
+}
diff --git a/config/config_C6_dev.yaml b/config/config_C6_dev.yaml
index 0619cb6cb80297200eecec248e94b792fe8649a7..c1b5ba1a74b2485a6776c07eac44805d80168b83 100644
--- a/config/config_C6_dev.yaml
+++ b/config/config_C6_dev.yaml
@@ -9,15 +9,14 @@
 # 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
-<<<<<<< HEAD
-work_dir: "/share/home/zhangxin/CSST_SIM/CSST_new_sim/csst-simulation/"
-=======
-work_dir: "/share/home/fangyuedong/new_sim/workplace/"
-# work_dir: "/share/C6_new_sim_2sq"
->>>>>>> new_sim
+##<<<<<<< HEAD
+##work_dir: "/share/home/zhangxin/CSST_SIM/CSST_new_sim/csst-simulation/"
+##=======
+work_dir: "/share/home/weichengliang/CSST_git/test_new_sim/outputs/"
+##>>>>>>> new_sim
 data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
-run_name: "C6_new_sim_2sq_run1"
-project_cycle: 6
+run_name: "testRun2"
+project_cycle: 8
 run_counter: 1
 
 # Whether to use MPI
@@ -71,7 +70,7 @@ obs_setting:
   # "Spectroscopic": simulate slitless spectroscopic chips only
   # "FGS": simulate FGS chips only (31-42)
   # "All": simulate full focal plane
-  survey_type: "Spectroscopic"
+  survey_type: "Photometric"
   
   # Exposure time [seconds]
   exp_time: 150.
@@ -107,7 +106,7 @@ obs_setting:
   # - give a list of indexes of chips: [ip_1, ip_2...]
   # - run all chips: null
   # Note: for all pointings
-  run_chips: [10]
+  run_chips: [8]
 
   # Whether to enable astrometric modeling
   enable_astrometric_model: True
@@ -178,12 +177,14 @@ ins_effects:
   non_linear:     YES  # Whether to add non-linearity
   cosmic_ray:     YES  # Whether to add cosmic-ray
   cray_differ:    YES  # Whether to generate different cosmic ray maps CAL and MS output
-  cte_trail:      YES  # Whether to simulate CTE trails
+  cte_trail:      YES  # Whether to simulate CTE trails, CTI_lgl_v0.3.tar.gz
   saturbloom:     YES  # Whether to simulate Saturation & Blooming
   add_badcolumns: YES  # Whether to add bad columns
   add_hotpixels:  YES  # Whether to add hot pixels
   add_deadpixels: YES  # Whether to add dead(dark) pixels
   bright_fatter:  YES  # Whether to simulate Brighter-Fatter (also diffusion) effect
+  add_prescan:    YES  # Whether to add pre/over-scan
+  format_output:  YES  ##1*16 output
 
   # Values: 
   # default values have been defined individually for each chip in: