mcmlgo.c

经典的蒙特卡罗程序

/***********************************************************
 *  Copyright Univ. of Texas M.D. Anderson Cancer Center
 *  1992.
 *
 *	Launch, move, and record photon weight.
 ****/

#include "mcml.h"

#define STANDARDTEST 0
  /* testing program using fixed rnd seed. */

#define PARTIALREFLECTION 0     
  /* 1=split photon, 0=statistical reflection. */

#define COSZERO (1.0-1.0E-12)	
  /* cosine of about 1e-6 rad. */

#define COS90D  1.0E-6		
  /* cosine of about 1.57 - 1e-6 rad. */


/***********************************************************
 *	A random number generator from Numerical Recipes in C.
 ****/
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC 1.0E-9

float ran3(int *idum)
{
  static int inext,inextp;
  static long ma[56];
  static int iff=0;
  long mj,mk;
  int i,ii,k;
  
  if (*idum < 0 || iff == 0) {
    iff=1;
    mj=MSEED-(*idum < 0 ? -*idum : *idum);
    mj %= MBIG;
    ma[55]=mj;
    mk=1;
    for (i=1;i<=54;i++) {
      ii=(21*i) % 55;
      ma[ii]=mk;
      mk=mj-mk;
      if (mk < MZ) mk += MBIG;
      mj=ma[ii];
    }
    for (k=1;k<=4;k++)
      for (i=1;i<=55;i++) {
	ma[i] -= ma[1+(i+30) % 55];
	if (ma[i] < MZ) ma[i] += MBIG;
      }
    inext=0;
    inextp=31;
    *idum=1;
  }
  if (++inext == 56) inext=1;
  if (++inextp == 56) inextp=1;
  mj=ma[inext]-ma[inextp];
  if (mj < MZ) mj += MBIG;
  ma[inext]=mj;
  return mj*FAC;
}

#undef MBIG
#undef MSEED
#undef MZ
#undef FAC


/***********************************************************
 *	Generate a random number between 0 and 1.  Take a 
 *	number as seed the first time entering the function.  
 *	The seed is limited to 1<<15.  
 *	We found that when idum is too large, ran3 may return 
 *	numbers beyond 0 and 1.
 ****/
double RandomNum(void)
{
  static Boolean first_time=1;
  static int idum;	/* seed for ran3. */
  
  if(first_time) {
#if STANDARDTEST /* Use fixed seed to test the program. */
    idum = - 1;
#else
    idum = -(int)time(NULL)%(1<<15);
	  /* use 16-bit integer as the seed. */
#endif
    ran3(&idum);
    first_time = 0;
    idum = 1;
  }
  
  return( (double)ran3(&idum) );
}

/***********************************************************
 *	Compute the specular reflection. 
 *
 *	If the first layer is a turbid medium, use the Fresnel
 *	reflection from the boundary of the first layer as the 
 *	specular reflectance.
 *
 *	If the first layer is glass, multiple reflections in
 *	the first layer is considered to get the specular
 *	reflectance.
 *
 *	The subroutine assumes the Layerspecs array is correctly 
 *	initialized.
 ****/
double Rspecular(LayerStruct * Layerspecs_Ptr)
{
  double r1, r2;
	/* direct reflections from the 1st and 2nd layers. */
  double temp;
  
  temp =(Layerspecs_Ptr[0].n - Layerspecs_Ptr[1].n)
	   /(Layerspecs_Ptr[0].n + Layerspecs_Ptr[1].n);
  r1 = temp*temp;
  
  if((Layerspecs_Ptr[1].mua == 0.0) 
  && (Layerspecs_Ptr[1].mus == 0.0))  { /* glass layer. */
    temp = (Layerspecs_Ptr[1].n - Layerspecs_Ptr[2].n)
		  /(Layerspecs_Ptr[1].n + Layerspecs_Ptr[2].n);
    r2 = temp*temp;
    r1 = r1 + (1-r1)*(1-r1)*r2/(1-r1*r2);
  }
  
  return (r1);	
}

/***********************************************************
 *	Initialize a photon packet.
 ****/
void LaunchPhoton(double Rspecular,
				  LayerStruct  * Layerspecs_Ptr,
				  PhotonStruct * Photon_Ptr)
{
  Photon_Ptr->w	 	= 1.0 - Rspecular;	
  Photon_Ptr->dead 	= 0;
  Photon_Ptr->layer = 1;
  Photon_Ptr->s	= 0;
  Photon_Ptr->sleft= 0;
  
  Photon_Ptr->x 	= 0.0;	
  Photon_Ptr->y	 	= 0.0;	
  Photon_Ptr->z	 	= 0.0;	
  Photon_Ptr->ux	= 0.0;	
  Photon_Ptr->uy	= 0.0;	
  Photon_Ptr->uz	= 1.0;	
  
  if((Layerspecs_Ptr[1].mua == 0.0) 
  && (Layerspecs_Ptr[1].mus == 0.0))  { /* glass layer. */
    Photon_Ptr->layer 	= 2;
    Photon_Ptr->z	= Layerspecs_Ptr[2].z0;	
  }
}

/***********************************************************
 *	Choose (sample) a new theta angle for photon propagation
 *	according to the anisotropy.
 *
 *	If anisotropy g is 0, then
 *		cos(theta) = 2*rand-1.
 *	otherwise
 *		sample according to the Henyey-Greenstein function.
 *
 *	Returns the cosine of the polar deflection angle theta.
 ****/
double SpinTheta(double g)
{
  double cost;
  
  if(g == 0.0) 
    cost = 2*RandomNum() -1;
  else {
    double temp = (1-g*g)/(1-g+2*g*RandomNum());
    cost = (1+g*g - temp*temp)/(2*g);
	if(cost < -1) cost = -1;
	else if(cost > 1) cost = 1;
  }
  return(cost);
}


/***********************************************************
 *	Choose a new direction for photon propagation by 
 *	sampling the polar deflection angle theta and the 
 *	azimuthal angle psi.
 *
 *	Note:
 *  	theta: 0 - pi so sin(theta) is always positive 
 *  	feel free to use sqrt() for cos(theta).
 * 
 *  	psi:   0 - 2pi 
 *  	for 0-pi  sin(psi) is + 
 *  	for pi-2pi sin(psi) is - 
 ****/
void Spin(double g,
		  PhotonStruct * Photon_Ptr)
{
  double cost, sint;	/* cosine and sine of the */
						/* polar deflection angle theta. */
  double cosp, sinp;	/* cosine and sine of the */
						/* azimuthal angle psi. */
  double ux = Photon_Ptr->ux;
  double uy = Photon_Ptr->uy;
  double uz = Photon_Ptr->uz;
  double psi;

  cost = SpinTheta(g);
  sint = sqrt(1.0 - cost*cost);	
	/* sqrt() is faster than sin(). */

  psi = 2.0*PI*RandomNum(); /* spin psi 0-2pi. */
  cosp = cos(psi);
  if(psi<PI)
    sinp = sqrt(1.0 - cosp*cosp);	
	  /* sqrt() is faster than sin(). */
  else
    sinp = - sqrt(1.0 - cosp*cosp);	
  
  if(fabs(uz) > COSZERO)  { 	/* normal incident. */
    Photon_Ptr->ux = sint*cosp;
    Photon_Ptr->uy = sint*sinp;
    Photon_Ptr->uz = cost*SIGN(uz);	
	  /* SIGN() is faster than division. */
  }
  else  {		/* regular incident. */
    double temp = sqrt(1.0 - uz*uz);
    Photon_Ptr->ux = sint*(ux*uz*cosp - uy*sinp)
					/temp + ux*cost;
    Photon_Ptr->uy = sint*(uy*uz*cosp + ux*sinp)
					/temp + uy*cost;
    Photon_Ptr->uz = -sint*cosp*temp + uz*cost;
  }
}

/***********************************************************
 *	Move the photon s away in the current layer of medium.  
 ****/
void Hop(PhotonStruct *	Photon_Ptr)
{
  double s = Photon_Ptr->s;

  Photon_Ptr->x += s*Photon_Ptr->ux;
  Photon_Ptr->y += s*Photon_Ptr->uy;
  Photon_Ptr->z += s*Photon_Ptr->uz;
}			

/***********************************************************
 *	If uz != 0, return the photon step size in glass, 
 *	Otherwise, return 0.
 *
 *	The step size is the distance between the current 
 *	position and the boundary in the photon direction.
 *
 *	Make sure uz !=0 before calling this function.
 ****/
void StepSizeInGlass(PhotonStruct *  Photon_Ptr,
					 InputStruct  *  In_Ptr)
{
  double dl_b;	/* step size to boundary. */
  short  layer = Photon_Ptr->layer;
  double uz = Photon_Ptr->uz;
  
  /* Stepsize to the boundary. */	
  if(uz>0.0)
    dl_b = (In_Ptr->layerspecs[layer].z1 - Photon_Ptr->z)
		   /uz;
  else if(uz<0.0)
    dl_b = (In_Ptr->layerspecs[layer].z0 - Photon_Ptr->z)
		   /uz;
  else
    dl_b = 0.0;
  
  Photon_Ptr->s = dl_b;
}

/***********************************************************
 *	Pick a step size for a photon packet when it is in 
 *	tissue.
 *	If the member sleft is zero, make a new step size 
 *	with: -log(rnd)/(mua+mus).
 *	Otherwise, pick up the leftover in sleft.
 *
 *	Layer is the index to layer.
 *	In_Ptr is the input parameters.
 ****/
void StepSizeInTissue(PhotonStruct * Photon_Ptr,
					  InputStruct  * In_Ptr)
{
  short  layer = Photon_Ptr->layer;
  double mua = In_Ptr->layerspecs[layer].mua;
  double mus = In_Ptr->layerspecs[layer].mus;
  
  if(Photon_Ptr->sleft == 0.0) {  /* make a new step. */
    double rnd;

    do rnd = RandomNum(); 
      while( rnd <= 0.0 );    /* avoid zero. */
	Photon_Ptr->s = -log(rnd)/(mua+mus);
  }
  else {	/* take the leftover. */
	Photon_Ptr->s = Photon_Ptr->sleft/(mua+mus);
	Photon_Ptr->sleft = 0.0;
  }
}

/***********************************************************
 *	Check if the step will hit the boundary.
 *	Return 1 if hit boundary.
 *	Return 0 otherwise.
 *
 * 	If the projected step hits the boundary, the members
 *	s and sleft of Photon_Ptr are updated.
 ****/
Boolean HitBoundary(PhotonStruct *  Photon_Ptr,
					InputStruct  *  In_Ptr)
{
  double dl_b;  /* length to boundary. */
  short  layer = Photon_Ptr->layer;
  double uz = Photon_Ptr->uz;
  Boolean hit;
  
  /* Distance to the boundary. */
  if(uz>0.0)
    dl_b = (In_Ptr->layerspecs[layer].z1 
			- Photon_Ptr->z)/uz;	/* dl_b>0. */
  else if(uz<0.0)
    dl_b = (In_Ptr->layerspecs[layer].z0 
			- Photon_Ptr->z)/uz;	/* dl_b>0. */
  
  if(uz != 0.0 && Photon_Ptr->s > dl_b) {
	  /* not horizontal & crossing. */
    double mut = In_Ptr->layerspecs[layer].mua 
				+ In_Ptr->layerspecs[layer].mus;

    Photon_Ptr->sleft = (Photon_Ptr->s - dl_b)*mut;
    Photon_Ptr->s    = dl_b;
    hit = 1;
  }
  else
    hit = 0;
  
  return(hit);
}

/***********************************************************
 *	Drop photon weight inside the tissue (not glass).
 *
 *  The photon is assumed not dead. 
 *
 *	The weight drop is dw = w*mua/(mua+mus).
 *
 *	The dropped weight is assigned to the absorption array 
 *	elements.
 ****/
void Drop(InputStruct  *	In_Ptr, 
		  PhotonStruct *	Photon_Ptr,
		  OutStruct *		Out_Ptr)
{
  double dwa;		/* absorbed weight.*/
  double x = Photon_Ptr->x;
  double y = Photon_Ptr->y;
  double izd, ird;	/* LW 5/20/98. To avoid out of short range.*/
  short  iz, ir;	/* index to z & r. */
  short  layer = Photon_Ptr->layer;
  double mua, mus;