add CTI/1*16/prescan/overscan

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

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

ObservationSim/Instrument/Chip/libCTI/CTI_modeling.py

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

ObservationSim/Instrument/Chip/libCTI/Makefile

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

ObservationSim/Instrument/Chip/libCTI/readme

+CTI_lgl_v0.3/

ObservationSim/Instrument/Chip/libCTI/src/add_CTI1.c

+#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;
+  }*/