mcmlgo.c

经典的蒙特卡罗程序

  /* compute array indices. */
  izd = Photon_Ptr->z/In_Ptr->dz;
  if(izd>In_Ptr->nz-1) iz=In_Ptr->nz-1;
  else iz = izd;
  
  ird = sqrt(x*x+y*y)/In_Ptr->dr;
  if(ird>In_Ptr->nr-1) ir=In_Ptr->nr-1;
  else ir = ird;
  
  /* update photon weight. */
  mua = In_Ptr->layerspecs[layer].mua;
  mus = In_Ptr->layerspecs[layer].mus;
  dwa = Photon_Ptr->w * mua/(mua+mus);
  Photon_Ptr->w -= dwa;
  
  /* assign dwa to the absorption array element. */
  Out_Ptr->A_rz[ir][iz]	+= dwa;
}

/***********************************************************
 *	The photon weight is small, and the photon packet tries 
 *	to survive a roulette.
 ****/
void Roulette(PhotonStruct * Photon_Ptr)
{
  if(Photon_Ptr->w == 0.0)	
    Photon_Ptr->dead = 1;
  else if(RandomNum() < CHANCE) /* survived the roulette.*/
    Photon_Ptr->w /= CHANCE;
  else 
    Photon_Ptr->dead = 1;
}

/***********************************************************
 *	Compute the Fresnel reflectance.
 *
 *	Make sure that the cosine of the incident angle a1
 *	is positive, and the case when the angle is greater 
 *	than the critical angle is ruled out.
 *
 * 	Avoid trigonometric function operations as much as
 *	possible, because they are computation-intensive.
 ****/
double RFresnel(double n1,	/* incident refractive index.*/
				double n2,	/* transmit refractive index.*/
				double ca1,	/* cosine of the incident */
							/* angle. 0<a1<90 degrees. */
				double * ca2_Ptr)  /* pointer to the */
							/* cosine of the transmission */
							/* angle. a2>0. */
{
  double r;
  
  if(n1==n2) {			  	/** matched boundary. **/
    *ca2_Ptr = ca1;
    r = 0.0;
  }
  else if(ca1>COSZERO) {	/** normal incident. **/
    *ca2_Ptr = ca1;
    r = (n2-n1)/(n2+n1);
    r *= r;
  }
  else if(ca1<COS90D)  {	/** very slant. **/
    *ca2_Ptr = 0.0;
    r = 1.0;
  }
  else  {			  		/** general. **/
    double sa1, sa2;	
	  /* sine of the incident and transmission angles. */
    double ca2;
    
    sa1 = sqrt(1-ca1*ca1);
    sa2 = n1*sa1/n2;
    if(sa2>=1.0) {	
	  /* double check for total internal reflection. */
      *ca2_Ptr = 0.0;
      r = 1.0;
    }
    else  {
      double cap, cam;	/* cosines of the sum ap or */
						/* difference am of the two */
						/* angles. ap = a1+a2 */
						/* am = a1 - a2. */
      double sap, sam;	/* sines. */
      
      *ca2_Ptr = ca2 = sqrt(1-sa2*sa2);
      
      cap = ca1*ca2 - sa1*sa2; /* c+ = cc - ss. */
      cam = ca1*ca2 + sa1*sa2; /* c- = cc + ss. */
      sap = sa1*ca2 + ca1*sa2; /* s+ = sc + cs. */
      sam = sa1*ca2 - ca1*sa2; /* s- = sc - cs. */
      r = 0.5*sam*sam*(cam*cam+cap*cap)/(sap*sap*cam*cam); 
		/* rearranged for speed. */
    }
  }
  return(r);
}

/***********************************************************
 *	Record the photon weight exiting the first layer(uz<0), 
 *	no matter whether the layer is glass or not, to the 
 *	reflection array.
 *
 *	Update the photon weight as well.
 ****/
void RecordR(double			Refl,	/* reflectance. */
			 InputStruct  *	In_Ptr,
			 PhotonStruct *	Photon_Ptr,
			 OutStruct *	Out_Ptr)
{
  double x = Photon_Ptr->x;
  double y = Photon_Ptr->y;
  short  ir, ia;	/* index to r & angle. */
  double ird, iad;	/* LW 5/20/98. To avoid out of short range.*/
  
  ird = sqrt(x*x+y*y)/In_Ptr->dr;
  if(ird>In_Ptr->nr-1) ir=In_Ptr->nr-1;
  else ir = ird;
  
  iad = acos(-Photon_Ptr->uz)/In_Ptr->da;
  if(iad>In_Ptr->na-1) ia=In_Ptr->na-1;
  else ia = iad;
  
  /* assign photon to the reflection array element. */
  Out_Ptr->Rd_ra[ir][ia] += Photon_Ptr->w*(1.0-Refl);
  
  Photon_Ptr->w *= Refl;
}

/***********************************************************
 *	Record the photon weight exiting the last layer(uz>0), 
 *	no matter whether the layer is glass or not, to the 
 *	transmittance array.
 *
 *	Update the photon weight as well.
 ****/
void RecordT(double 		Refl,
			 InputStruct  *	In_Ptr,
			 PhotonStruct *	Photon_Ptr,
			 OutStruct *	Out_Ptr)
{
  double x = Photon_Ptr->x;
  double y = Photon_Ptr->y;
  short  ir, ia;	/* index to r & angle. */
  double ird, iad;	/* LW 5/20/98. To avoid out of short range.*/
  
  ird = sqrt(x*x+y*y)/In_Ptr->dr;
  if(ird>In_Ptr->nr-1) ir=In_Ptr->nr-1;
  else ir = ird;
  
  iad = acos(Photon_Ptr->uz)/In_Ptr->da; /* LW 1/12/2000. Removed -. */
  if(iad>In_Ptr->na-1) ia=In_Ptr->na-1;
  else ia = iad;
  
  /* assign photon to the transmittance array element. */
  Out_Ptr->Tt_ra[ir][ia] += Photon_Ptr->w*(1.0-Refl);
  
  Photon_Ptr->w *= Refl;
}

/***********************************************************
 *	Decide whether the photon will be transmitted or 
 *	reflected on the upper boundary (uz<0) of the current 
 *	layer.
 *
 *	If "layer" is the first layer, the photon packet will 
 *	be partially transmitted and partially reflected if 
 *	PARTIALREFLECTION is set to 1,
 *	or the photon packet will be either transmitted or 
 *	reflected determined statistically if PARTIALREFLECTION 
 *	is set to 0.
 *
 *	Record the transmitted photon weight as reflection.  
 *
 *	If the "layer" is not the first layer and the photon 
 *	packet is transmitted, move the photon to "layer-1".
 *
 *	Update the photon parmameters.
 ****/
void CrossUpOrNot(InputStruct  *	In_Ptr, 
				  PhotonStruct *	Photon_Ptr,
				  OutStruct *		Out_Ptr)
{
  double uz = Photon_Ptr->uz; /* z directional cosine. */
  double uz1;	/* cosines of transmission alpha. always */
				/* positive. */
  double r=0.0;	/* reflectance */
  short  layer = Photon_Ptr->layer;
  double ni = In_Ptr->layerspecs[layer].n;
  double nt = In_Ptr->layerspecs[layer-1].n;
  
  /* Get r. */
  if( - uz <= In_Ptr->layerspecs[layer].cos_crit0) 
    r=1.0;		      /* total internal reflection. */
  else r = RFresnel(ni, nt, -uz, &uz1);
  
#if PARTIALREFLECTION
  if(layer == 1 && r<1.0) {	/* partially transmitted. */
    Photon_Ptr->uz = -uz1;	/* transmitted photon. */
    RecordR(r, In_Ptr, Photon_Ptr, Out_Ptr);
    Photon_Ptr->uz = -uz;	/* reflected photon. */
  }		
  else if(RandomNum() > r) {/* transmitted to layer-1. */
    Photon_Ptr->layer--;
    Photon_Ptr->ux *= ni/nt;
    Photon_Ptr->uy *= ni/nt;
    Photon_Ptr->uz = -uz1;
  }
  else			      		/* reflected. */
    Photon_Ptr->uz = -uz;
#else
  if(RandomNum() > r) {		/* transmitted to layer-1. */
    if(layer==1)  {
      Photon_Ptr->uz = -uz1;
      RecordR(0.0, In_Ptr, Photon_Ptr, Out_Ptr);
      Photon_Ptr->dead = 1;
    }
    else {
      Photon_Ptr->layer--;
      Photon_Ptr->ux *= ni/nt;
      Photon_Ptr->uy *= ni/nt;
      Photon_Ptr->uz = -uz1;
    }
  }
  else 						/* reflected. */
    Photon_Ptr->uz = -uz;
#endif
}

/***********************************************************
 *	Decide whether the photon will be transmitted  or be 
 *	reflected on the bottom boundary (uz>0) of the current 
 *	layer.
 *
 *	If the photon is transmitted, move the photon to 
 *	"layer+1". If "layer" is the last layer, record the 
 *	transmitted weight as transmittance. See comments for 
 *	CrossUpOrNot.
 *
 *	Update the photon parmameters.
 ****/
void CrossDnOrNot(InputStruct  *	In_Ptr, 
				  PhotonStruct *	Photon_Ptr,
				  OutStruct *		Out_Ptr)
{
  double uz = Photon_Ptr->uz; /* z directional cosine. */
  double uz1;	/* cosines of transmission alpha. */
  double r=0.0;	/* reflectance */
  short  layer = Photon_Ptr->layer;
  double ni = In_Ptr->layerspecs[layer].n;
  double nt = In_Ptr->layerspecs[layer+1].n;
  
  /* Get r. */
  if( uz <= In_Ptr->layerspecs[layer].cos_crit1) 
    r=1.0;		/* total internal reflection. */
  else r = RFresnel(ni, nt, uz, &uz1);
  
#if PARTIALREFLECTION	
  if(layer == In_Ptr->num_layers && r<1.0) {
    Photon_Ptr->uz = uz1;
    RecordT(r, In_Ptr, Photon_Ptr, Out_Ptr);
    Photon_Ptr->uz = -uz;
  }
  else if(RandomNum() > r) {/* transmitted to layer+1. */
    Photon_Ptr->layer++;
    Photon_Ptr->ux *= ni/nt;
    Photon_Ptr->uy *= ni/nt;
    Photon_Ptr->uz = uz1;
  }
  else 						/* reflected. */
    Photon_Ptr->uz = -uz;
#else
  if(RandomNum() > r) {		/* transmitted to layer+1. */
    if(layer == In_Ptr->num_layers) {
      Photon_Ptr->uz = uz1;
      RecordT(0.0, In_Ptr, Photon_Ptr, Out_Ptr);
      Photon_Ptr->dead = 1;
    }
    else {
      Photon_Ptr->layer++;
      Photon_Ptr->ux *= ni/nt;
      Photon_Ptr->uy *= ni/nt;
      Photon_Ptr->uz = uz1;
    }
  }
  else 						/* reflected. */
    Photon_Ptr->uz = -uz;
#endif
}

/***********************************************************
 ****/
void CrossOrNot(InputStruct  *	In_Ptr, 
				PhotonStruct *	Photon_Ptr,
				OutStruct    *	Out_Ptr)
{
  if(Photon_Ptr->uz < 0.0)
    CrossUpOrNot(In_Ptr, Photon_Ptr, Out_Ptr);
  else
    CrossDnOrNot(In_Ptr, Photon_Ptr, Out_Ptr);
}

/***********************************************************
 *	Move the photon packet in glass layer.
 *	Horizontal photons are killed because they will
 *	never interact with tissue again.
 ****/
void HopInGlass(InputStruct  * In_Ptr,
				PhotonStruct * Photon_Ptr,
				OutStruct    * Out_Ptr)
{
  double dl;     /* step size. 1/cm */
  
  if(Photon_Ptr->uz == 0.0) { 
	/* horizontal photon in glass is killed. */
    Photon_Ptr->dead = 1;
  }
  else {
    StepSizeInGlass(Photon_Ptr, In_Ptr);
    Hop(Photon_Ptr);
    CrossOrNot(In_Ptr, Photon_Ptr, Out_Ptr);
  }
}

/***********************************************************
 *	Set a step size, move the photon, drop some weight, 
 *	choose a new photon direction for propagation.  
 *
 *	When a step size is long enough for the photon to 
 *	hit an interface, this step is divided into two steps. 
 *	First, move the photon to the boundary free of 
 *	absorption or scattering, then decide whether the 
 *	photon is reflected or transmitted.
 *	Then move the photon in the current or transmission 
 *	medium with the unfinished stepsize to interaction 
 *	site.  If the unfinished stepsize is still too long, 
 *	repeat the above process.  
 ****/
void HopDropSpinInTissue(InputStruct  *  In_Ptr,
						 PhotonStruct *  Photon_Ptr,
						 OutStruct    *  Out_Ptr)
{
  StepSizeInTissue(Photon_Ptr, In_Ptr);
  
  if(HitBoundary(Photon_Ptr, In_Ptr)) {
    Hop(Photon_Ptr);	/* move to boundary plane. */
    CrossOrNot(In_Ptr, Photon_Ptr, Out_Ptr);
  }
  else {
    Hop(Photon_Ptr);
    Drop(In_Ptr, Photon_Ptr, Out_Ptr);
    Spin(In_Ptr->layerspecs[Photon_Ptr->layer].g, 
		Photon_Ptr);
  }
}

/***********************************************************
 ****/
void HopDropSpin(InputStruct  *  In_Ptr,
				 PhotonStruct *  Photon_Ptr,
				 OutStruct    *  Out_Ptr)
{
  short layer = Photon_Ptr->layer;

  if((In_Ptr->layerspecs[layer].mua == 0.0) 
  && (In_Ptr->layerspecs[layer].mus == 0.0)) 
	/* glass layer. */
    HopInGlass(In_Ptr, Photon_Ptr, Out_Ptr);
  else
    HopDropSpinInTissue(In_Ptr, Photon_Ptr, Out_Ptr);
  
  if( Photon_Ptr->w < In_Ptr->Wth && !Photon_Ptr->dead) 
    Roulette(Photon_Ptr);
}