update add_CTI.c << add_CTI3.c

ObservationSim/Instrument/Chip/libCTI/readme

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

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

@@ -7,10 +7,8 @@
 
 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 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 ntotal;
   printf("image parameters: nx=%i ny=%i noverscan=%i\n",nx,ny,noverscan);
@@ -23,16 +21,16 @@ void CTI_simul(int **image, int nx, int ny, int noverscan, int nsp, float *rho_t
   printf("release_seed = %i\n",release_seed);
   printf("begin CTI simulation\n");
   float ***ntrap;
-  float **sp;
-  float tmp;
-  int ntrap_tmp;
+  float tmp,*randv;
+  int ntrap_tmp,nrandv=100000;
   int *a0,*acti,dim[2];
-  int i,j,k,l;
+  int i,j,k,l,iseed=-1;
   ntotal=ny+noverscan;
-  ntrap = f3tensor(0,nsp,0,ny,0,nmax+1);  
-  sp = matrix(0,ny,0,nmax+1);
+  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);
   
   for(k=0;k<nx;k++){
     printf("column k=%i/%i\r",k,nx);
@@ -40,101 +38,105 @@ void CTI_simul(int **image, int nx, int ny, int noverscan, int nsp, float *rho_t
         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<ny;i++){a0[i]=image[k][i];}
+    
+    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];
-        /*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);
+  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)
+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;
   float height,wre;
   int *flow,**traped,flowt,ntraptot,ntraped,topoftrap,nrelease;
-  int ntotal,i,j,k,nmove,isum,ntmp,isp,*tmptraped;
+  int ntotal,i,j,k,nmove,isum,ntmp,isp,*tmptraped,ncumran;
   wre=1/w;
   ntotal=ny+noverscan;
-  f=vector(0,nsp);
-  traped = imatrix(0,nsp,0,ntotal);
+  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=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=ny;i<ntotal;i++)flow[i]=0;
+  for(j=0;j<ntotal;j++){traped[0][j]=traped[1][j]=traped[2][j]=0;}
+  
   
   
+  //initialize
   for(i=0;i<ny;i++){
     flow[i]=a0[i];
-    for(j=0;j<nsp;j++){
-      ntraptot = ntrap[j][i][0];
+    //for(j=0;j<nsp;j++){
       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[j][i];///  LGL
+	tmpntrap=ntrap[1][i];
 	while(tmpntrap[ntraped+1]<=height){ntraped++;}
 	if(flow[i]<ntraped)ntraped=flow[i];
-	flow[i]-=ntraped; // LGL
-	traped[j][i]=ntraped;
+	flow[i]-=ntraped; 
+	traped[1][i]=ntraped;
       }
+      //sp2;
+      ntraptot = ntrap[2][i][0];
+      if(ntraptot>0){
+	ntraped=0;
+	tmpntrap=ntrap[2][i];
+	while(tmpntrap[ntraped+1]<=height){ntraped++;}
+	if(flow[i]<ntraped)ntraped=flow[i];
+	flow[i]-=ntraped; 
+	traped[2][i]=ntraped;
       }
   }
   
+  // 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;
-      //for(isp=0;isp<nsp;isp++){
-      isp=0;
-	  ntraptot=ntrap[isp][j][0];
+      //sp0
+      ntraptot=ntrap[0][j][0];
       if(ntraptot>0){
 	ntraped=0; 
 	//height=flow[j]*wre; 
-	    tmpntrap=ntrap[isp][j];
-	    tmptraped=traped[isp];
+	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++){
-	        if(ran1(&release_seed)<f[isp]){
+	    ncumran++;
+	    if(randv[ncumran]<f[0]){
 	      flow[j]++;
 	      tmptraped[j]--;
 	    }
@@ -147,21 +149,20 @@ void addCTI(int *a0,int ny,int noverscan, int nsp, float beta,float w, float c,f
 	  tmptraped[j] += ntmp;
 	}
       }
-	  if(nsp==1)continue;  
-	  
-	  isp=1;
-	  ntraptot=ntrap[isp][j][0];
+      //sp1
+      ntraptot=ntrap[1][j][0];
       if(ntraptot>0){
 	ntraped=0; 
 	//height=flow[j]*wre; 
-	    tmpntrap=ntrap[isp][j];
-	    tmptraped=traped[isp];
+	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++){
-	            if(ran1(&release_seed)<f[isp]){
+	    ncumran++;
+	    if(randv[ncumran]<f[1]){
 	      flow[j]++;
 	      tmptraped[j]--;
 	    }
@@ -174,21 +175,21 @@ void addCTI(int *a0,int ny,int noverscan, int nsp, float beta,float w, float c,f
 	  tmptraped[j] += ntmp;
 	}
       }
-	  if(nsp==2)continue; 
       
-	  isp=2;
-	  ntraptot=ntrap[isp][j][0];
+      //sp2
+      ntraptot=ntrap[2][j][0];
       if(ntraptot>0){
 	ntraped=0; 
 	//height=flow[j]*wre; 
-	    tmpntrap=ntrap[isp][j];
-	    tmptraped=traped[isp];
+	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++){
-	        if(ran1(&release_seed)<f[isp]){
+	    ncumran++;
+	    if(randv[ncumran]<f[2]){
 	      flow[j]++;
 	      tmptraped[j]--;
 	    }
@@ -201,9 +202,10 @@ void addCTI(int *a0,int ny,int noverscan, int nsp, float beta,float w, float c,f
 	  tmptraped[j] += ntmp;
 	}
       }
-	  if(nsp==3)continue;
       
-	  isp=3;
+      //if(nsp==3)continue;
+      
+      /*isp=3;
 	ntraptot=ntrap[isp][j][0];
 	if(ntraptot>0){
 	ntraped=0; 
@@ -227,15 +229,13 @@ void addCTI(int *a0,int ny,int noverscan, int nsp, float beta,float w, float c,f
 	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);
+  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){
@@ -280,34 +280,38 @@ void save_trap_map(int *trap_seeds, int nsp, int nx, int ny, int nmax, float *rh
   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 dim[2],i,j;
   int trap_seeds[3] = {1,100,1000};
   int release_seed = 10000;
   clock_t start, end;
   
   nimg=nx*ny;
   image=imatrix(0,nx,0,ny);
+  
+  imagecti=matrix(0,nx,0,ntotal);
   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++){
+  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();
-  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);
+  printf("running time %f\n",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); 
+  free_imatrix(image,0,nx,0,ny);
   return 0;
   }*/

setup.py

@@ -41,7 +41,7 @@ df_module = [CTypes('ObservationSim.Instrument.Chip.libBF.libmoduleBF',
                     include_dirs=['ObservationSim/Instrument/Chip/libBF/', '/usr/include']
 )]
 cti_module = [CTypes('ObservationSim.Instrument.Chip.libCTI.libmoduleCTI',
-                    ['ObservationSim/Instrument/Chip/libCTI/src/add_CTI1.c', 'ObservationSim/Instrument/Chip/libCTI/src/nrutil.c', 'ObservationSim/Instrument/Chip/libCTI/src/ran1.c', 'ObservationSim/Instrument/Chip/libCTI/src/ran2.c', 'ObservationSim/Instrument/Chip/libCTI/src/poidev.c', 'ObservationSim/Instrument/Chip/libCTI/src/gammln.c', 'ObservationSim/Instrument/Chip/libCTI/src/gasdev.c', 'ObservationSim/Instrument/Chip/libCTI/src/sort.c', 'ObservationSim/Instrument/Chip/libCTI/src/creattraps.c'],
+                    ['ObservationSim/Instrument/Chip/libCTI/src/add_CTI.c', 'ObservationSim/Instrument/Chip/libCTI/src/nrutil.c', 'ObservationSim/Instrument/Chip/libCTI/src/ran1.c', 'ObservationSim/Instrument/Chip/libCTI/src/ran2.c', 'ObservationSim/Instrument/Chip/libCTI/src/poidev.c', 'ObservationSim/Instrument/Chip/libCTI/src/gammln.c', 'ObservationSim/Instrument/Chip/libCTI/src/gasdev.c', 'ObservationSim/Instrument/Chip/libCTI/src/sort.c', 'ObservationSim/Instrument/Chip/libCTI/src/creattraps.c'],
                     include_dirs=['ObservationSim/Instrument/Chip/libCTI/src/', '/usr/include']
 )]