creattraps.c

#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);
}