ECE532 Biomedical Optics
© 1998 Steven L. Jacques, Scott A. Prahl
Oregon Graduate Institute

Time-resolved Monte Carlo

Spin the photon trajectory due to scattering

The remaining photon weight is now scattered into a new trajectory. If the anisitropy (g) equals 0, then the scattering is isotropic. If g > 0, then the program uses the Henyey-Greenstein scattering function originally introduced by Henyey and Greenstein to mimic the scattering of light from distant stars by galactic dust. The same function approximately mimics the scattering function experimentally observed in biological tissues (see HG function).

The four sections of SPIN are:

/**** SPIN 
*   Scatter photon into new trajectory defined by theta and psi.
*   Theta is specified by cos(theta), which is determined 
*   based on the Henyey-Greenstein scattering function.
*   Convert theta and psi into cosines ux, uy, uz. 
*****/
  /* Sample for costheta */
  rnd = RandomNum;
     if (g == 0.0)
        costheta = 2.0*rnd - 1.0;
     else {
        double temp = (1.0 - g*g)/(1.0 - g + 2*g*rnd);
        costheta = (1.0 + g*g - temp*temp)/(2.0*g);
        }
  sintheta = sqrt(1.0 - costheta*costheta); /* sqrt() is faster than sin(). */

  /* Sample psi. */
  psi = 2.0*PI*RandomNum;
  cospsi = cos(psi);
  if (psi < PI)
    sinpsi = sqrt(1.0 - cospsi*cospsi);     /* sqrt() is faster than sin(). */
  else
    sinpsi = -sqrt(1.0 - cospsi*cospsi);

  /* New trajectory. */
  if (1 - fabs(uz) <= ONE_MINUS_COSZERO) {      /* close to perpendicular. */
    uxx = sintheta * cospsi;
    uyy = sintheta * sinpsi;
    uzz = costheta * SIGN(uz);   /* SIGN() is faster than division. */
    } 
  else {					/* usually use this option */
    temp = sqrt(1.0 - uz * uz);
    uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) / temp + ux * costheta;
    uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta;
    uzz = -sintheta * cospsi * temp + uz * costheta;
    }
    
  /* Update trajectory */
  ux = uxx;
  uy = uyy;
  uz = uzz;

NextMonte Carlo