update add_CTI.c << add_CTIfinal.c

ObservationSim/Instrument/Chip/libCTI/CTI_modeling.py

@@ -16,7 +16,7 @@ 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]
 
@@ -44,6 +44,7 @@ def bin2fits(bin_file,fits_dir,nsp,nx,ny,nmax):
             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])()

ObservationSim/Instrument/Chip/libCTI/readme

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

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

@@ -5,11 +5,13 @@
 #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);
+float poidev(float x, long *idum);
+void creattraps(long *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,float *randv);
-float ran1(int *idum);
-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 read_fits_2D(const char *argv,float *galaxy,int imagesize);
+//int write_fits_2D(const char *argv,float **stars,int *dim);
+float ran1(long *idum);
+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, long *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]);
@@ -17,301 +19,219 @@ void CTI_simul(int **image, int nx, int ny, int noverscan, int nsp, float *rho_t
   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("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 tmp,*randv;
+  float *randv;
+  float tmp;
   int ntrap_tmp,nrandv=100000;
   int *a0,*acti,dim[2];
-  int i,j,k,l,iseed=-1;
+  int i,j,k,l;
+  long iseed=-1;
   ntotal=ny+noverscan;
   ntrap = f3tensor(0,nsp-1,0,ny,0,nmax+1); 
   a0=ivector(0,ny);
   acti=ivector(0,ntotal);
   randv=vector(0,nrandv);
   for(i=0;i<nrandv;i++)randv[i]=ran1(&iseed);
-  
+
+  printf("column k/%i, k=", nx);
   for(k=0;k<nx;k++){
-    printf("column k=%i/%i\r",k,nx);
+    if(k % 8 == 0) fprintf(stdout, "%i  ",k); fflush(stdout);
     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]);        
+        creattraps(&trap_seeds[l],ntrap_tmp,ny,nmax,c,beta,ntrap[l]);
+       
     }
     
-    
     for(i=0;i<ny;i++){a0[i]=image[k][i];}
-    
+
+    //for(i=0;i<ntotal;i++)acti[i]=0;
     addCTI(a0,ny,noverscan,nsp,beta,w,c,t,ntrap,acti,release_seed,randv);
     
     for(i=0;i<ntotal;i++){
-        imagecti[k][i]=acti[i];
+      imagecti[k][i]=acti[i];
     }
   }
   printf("\nCTI simulation finished!\n");
   free_ivector(a0,0,ny);
   free_ivector(acti,0,ntotal);
-  free_f3tensor(ntrap,0,nsp-1,0,ny,0,nmax+1);
+  free_f3tensor(ntrap,0,nsp-1,0,ny,0,nmax+1);  
   free_vector(randv,0,nrandv);
   return;
 }
 
 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 *randv)
 {
-  float ran1(int *idum);
-  float *f,*tmpntrap;
+  
+  void trap_release(int *flow, int *traped,int ntraped,float *randv,int *ncumran,float f);
+  //float ran1(long *idum);
+  float *tmpntrap,*f,tmph[50];
   float height,wre;
-  int *flow,**traped,flowt,ntraptot,ntraped,topoftrap,nrelease;
-  int ntotal,i,j,k,nmove,isum,ntmp,isp,*tmptraped,ncumran;
+  int flowt,ntraptot,ntraped,topoftrap,nrelease,ntraped0,tmptrapedt,ntrapsum;
+  int ntotal,i,j,k,nmove,isum,ntmp,isp,*tmptraped,**traped,*flow;
+  int ncumran,ncumran0,jc,*jcindx;
+  long hotseed=-1;
   wre=1/w;
   ntotal=ny+noverscan;
   f=vector(0,nsp-1);
   traped = imatrix(0,nsp-1,0,ntotal);
-  flow=ivector(0,ntotal);  
-  for(i=0;i<nsp;i++){
-    f[i] = 1-exp(-1/t[i]);
-  }
-  for(i=ny;i<ntotal;i++)flow[i]=0;
-  for(j=0;j<ntotal;j++){traped[0][j]=traped[1][j]=traped[2][j]=0;}
+  flow=ivector(0,ntotal); 
+  jcindx=ivector(0,ntotal);
   
+  for(i=0;i<nsp;i++){f[i]=1-exp(-1./t[i]);}
   
+  for(i=ny;i<ntotal;i++)flow[i]=0;
+  
+  for(i=0;i<nsp;i++){for(j=ny;j<ntotal;j++){traped[i][j] = 0;}}
   
-  //initialize
+ 
+  //initialize 
   for(i=0;i<ny;i++){
     flow[i]=a0[i];
-    //for(j=0;j<nsp;j++){
+    for(j=0;j<nsp;j++){
+      ntraptot = ntrap[j][i][0];
       height=flow[i]*wre;
-      //sp0
-      ntraptot = ntrap[0][i][0];
-      if(ntraptot>0){
-	ntraped=0;
-	tmpntrap=ntrap[0][i];
-	while(tmpntrap[ntraped+1]<=height){ntraped++;}
-	if(flow[i]<ntraped)ntraped=flow[i];
-	flow[i]-=ntraped; 
-	traped[0][i]=ntraped;
-      }
-      //sp1;
-      ntraptot = ntrap[1][i][0];
-      if(ntraptot>0){
-	ntraped=0;
-	tmpntrap=ntrap[1][i];
-	while(tmpntrap[ntraped+1]<=height){ntraped++;}
-	if(flow[i]<ntraped)ntraped=flow[i];
-	flow[i]-=ntraped; 
-	traped[1][i]=ntraped;
-      }
-      //sp2;
-      ntraptot = ntrap[2][i][0];
       if(ntraptot>0){
 	ntraped=0;
-	tmpntrap=ntrap[2][i];
+	tmpntrap=ntrap[j][i];
 	while(tmpntrap[ntraped+1]<=height){ntraped++;}
 	if(flow[i]<ntraped)ntraped=flow[i];
 	flow[i]-=ntraped; 
-	traped[2][i]=ntraped;
+	traped[j][i]=ntraped;
       }
+      else{traped[j][i]=0;}
+    }
   }
-  
+
+  // rearrange trapes
+  jc=-1;
+  for(i=0;i<ny;i++){
+    ntrapsum=ntrap[0][i][0]+ntrap[1][i][0]+ntrap[2][i][0];
+    //if((ntrapsum!=0 || hotpixflag){ if there is none-zero trap or it is a hot pix
+    if(ntrapsum!=0){
+      jc++;
+      for(j=0;j<nsp;j++){
+	traped[j][jc]=traped[j][i];
+	ntraptot=ntrap[j][i][0];
+	for(k=0;k<=ntraptot+1;k++)ntrap[j][jc][k]=ntrap[j][i][k];
+      }      
+    }
+    jcindx[i+1]=jc;
+  }
+  for(i=ny;i<ntotal;i++)jcindx[i]=jc;
+
   // add CTI effects
-  for(i=0;i<ntotal;i++){
-    acti[i] = flow[0];
-    nmove=ny;
-     ncumran=randv[i]*10000;
-    if(i>noverscan)nmove=ntotal-i;  
-    for(j=0;j<nmove;j++){
-      flow[j]=flow[j+1];
-      height=flow[j]*wre;
-      //sp0
-      ntraptot=ntrap[0][j][0];
-      if(ntraptot>0){
-	ntraped=0; 
-	//height=flow[j]*wre; 
-	tmpntrap=ntrap[0][j];
-	tmptraped=traped[0];
-	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++){
-	    ncumran++;
-	    if(randv[ncumran]<f[0]){
-	      flow[j]++;
-	      tmptraped[j]--;
-	    }
-	  }
-	}
-	else{     //trap
-	  ntmp = -topoftrap;
-	  if(flow[j]<ntmp)ntmp=flow[j];
-	  flow[j]-=ntmp; 
-	  tmptraped[j] += ntmp;
-	}
-      }
-      //sp1
-      ntraptot=ntrap[1][j][0];
-      if(ntraptot>0){
-	ntraped=0; 
-	//height=flow[j]*wre; 
-	tmpntrap=ntrap[1][j];
-	tmptraped=traped[1];
-	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++){
-	    ncumran++;
-	    if(randv[ncumran]<f[1]){
-	      flow[j]++;
-	      tmptraped[j]--;
-	    }
-	  }
-	}
-	else{     //trap
-	  ntmp = -topoftrap;
-	  if(flow[j]<ntmp)ntmp=flow[j];
-	  flow[j]-=ntmp; 
-	  tmptraped[j] += ntmp;
-	}
-      }
-      
-      //sp2
-      ntraptot=ntrap[2][j][0];
-      if(ntraptot>0){
-	ntraped=0; 
-	//height=flow[j]*wre; 
-	tmpntrap=ntrap[2][j];
-	tmptraped=traped[2];
-	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++){
-	    ncumran++;
-	    if(randv[ncumran]<f[2]){
-	      flow[j]++;
-	      tmptraped[j]--;
-	    }
-	  }
-	}
-	else{     //trap
-	  ntmp = -topoftrap;
-	  if(flow[j]<ntmp)ntmp=flow[j];
-	  flow[j]-=ntmp; 
-	  tmptraped[j] += ntmp;
-	}
-      }
+  acti[0]=flow[0];
+  for(i=1;i<ntotal;i++){
+    nmove=jcindx[i];
+    for(j=nmove;j>=0;j--){    
+      //flow[i]+=poidev(0.01, &hotseed);
+      ncumran=randv[j+i]*10000;
+      height=flow[i]*wre;
       
-      //if(nsp==3)continue;
+      isp=0;    
+      tmpntrap=ntrap[isp][j];
+      ntraped=0; 
+      tmptraped=traped[isp];
+      while(tmpntrap[ntraped+1]<=height){ntraped++;}
+      topoftrap=tmptraped[j]-ntraped;	
+      if(topoftrap!=0){trap_release(&flow[i], &tmptraped[j],topoftrap,randv,&ncumran,f[isp]);}    
       
-      /*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;
-	}
-	}*/
-    }
+      isp=1;
+      tmpntrap=ntrap[isp][j];
+      ntraped=0; 
+      //height=flow[j]*wre; ;
+      tmptraped=traped[isp];
+      while(tmpntrap[ntraped+1]<=height){ntraped++;}
+      topoftrap=tmptraped[j]-ntraped;
+      if(topoftrap!=0){trap_release(&flow[i], &tmptraped[j],topoftrap,randv,&ncumran,f[isp]);}
+	
+      isp=2;
+      tmpntrap=ntrap[isp][j];
+      ntraped=0; 
+      //height=flow[j]*wre; 
+      tmptraped=traped[isp];
+      while(tmpntrap[ntraped+1]<=height){ntraped++;}
+      topoftrap=tmptraped[j]-ntraped;
+      if(topoftrap!=0){trap_release(&flow[i], &tmptraped[j],topoftrap,randv,&ncumran,f[isp]);}
+
+      //isp=3...
+      }
+    acti[i]=flow[i];
   }
   free_ivector(flow,0,ntotal);
-  free_vector(f,0,nsp-1);
+  free_ivector(jcindx,0,ntotal);
   free_imatrix(traped,0,nsp-1,0,ntotal);
+  free_vector(f,0,nsp-1);
   return;
 }
 
-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);
+
+void trap_release(int *flow, int *traped,int topoftrap,float *randv,int *ncumran,float f){
+  int ntmp,k;
+  if (topoftrap>0){  // release
+    for(k=0;k<topoftrap;k++){
+      *ncumran+=1;
+      if(randv[*ncumran]<f){
+	*flow+=1;
+	*traped-=1;
+      }	    
+    }
+    return;
+  }
+  else{     //trap
+    ntmp = -topoftrap;
+    if(*flow<ntmp)ntmp=*flow;
+    *flow-=ntmp; 
+    *traped += ntmp;
+  }
+  return;
 }
 
 /*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;
+  int **image,nimg,kct,i,j;
   float *image1,**imagecti;
-  int nx=4608,ny=4616,noverscan=84,nmax=10,nsp=3;
-  int ntotal = ny+noverscan;
+  int nx=4608,ny=4616,noverscan=84,nmax=12,nsp=3;
+  int ntotal;
   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],i,j;
-  int trap_seeds[3] = {1,100,1000};
+  int dim[2];
+  long trap_seeds[3] = {1,100,1000};
   int release_seed = 10000;
   clock_t start, end;
-  
+  ntotal=ny+noverscan;
   nimg=nx*ny;
   image=imatrix(0,nx,0,ny);
-  
-  imagecti=matrix(0,nx,0,ntotal);
   image1=vector(0,nimg);
+  imagecti=matrix(0,nx,0,ntotal);
   strcpy(fname1,"inputdata/image.fits");
   read_fits_2D(fname1,image1,nimg);kct=0;    
-  for (j=0;j<ny;j++){for(i=0;i<nx;i++){
+  for (j=0; j<ny; j++){
+    for(i=0;i<nx;i++){
       image[i][j]=image1[kct];
       kct++;
     }}
   
   free_vector(image1,0,nimg);
-
-  
   start = clock();
   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\n",time);
+  printf("running time %f ms\n",time);
   
   // output final image
-  strcpy(fname1,"output/image_CTI.fits");
+  strcpy(fname1,"output/image_CTIfinal.fits");
   dim[0]=nx;dim[1]=ntotal;
   write_fits_2D(fname1,imagecti,dim);
-  free_matrix(imagecti,0,nx,0,ntotal); 
+  free_matrix(imagecti,0,nx,0,ntotal);
   free_imatrix(image,0,nx,0,ny);
   return 0;
   }*/