/* ---------------------------------------------------------------

   The following function implements the Henyey-Greenstein phase
   function, so that a uniform random deviate r yields the value of
   the cosine of the scattering angle Phi.  The parameter g is the
   scattering asymmetry parameter, -1 < g < 1 */


double hg_phase_func(double g, double r) 
{
  double g2,a,b ;

  /* check for isotropic case */
  if (g == 0.0) return(2.0*r - 1.0) ;

  g2 = g*g ;
  a = (1.0-g2)/(1+g*(2.0*r-1.0)) ;
  b = a*a ;
  return( (1.0+g2-b)/(2.0*g) ) ;
}


/* ---------------------------------------------------------------
   The following fragment illustrates the process of converting the
   scattering direction in photon-relative coordinates to a new
   propagation direction in model coordinates. */


/* get cos and sin of scattering angle Theta  */

cosTheta = hg_phase_func(g,ran1()) ;     /* from Henyey-Greenstein phase function */
sinTheta = sqrt(1.0-cosTheta*cosTheta) ;

/* get cos(Phi) and sin(Phi) */
Phi = TWOPI*ran1() ;
cosPhi = cos(Phi) ;
sinPhi = sin(Phi) ;

/* define photon-relative x direction */
den = sqrt(k0*k0 + k1*k1) ;
if (den < 0.0001) {         /* deal with special case that k is vertical */
  xp0 = -sin(phi0) ;
  xp1 = cos(phi0) ;
  xp2 = 0.0 ;
} else {                    /* normal case (k not vertical) */
  xp0 = -k1/den ;
  xp1 = k0/den ;
  xp2 = 0.0 ;
}
kp0 = sinTheta*cosPhi ;     /* new direction in photon-relative coordinates */
kp1 = sinTheta*sinPhi ;
kp2 = cosTheta ;

/* new propagation vector in model coordinates */

k0n = xp0*kp0 - k2*xp1*kp1 + k0*kp2 ;
k1n = xp1*kp0 + k2*xp0*kp1 + k1*kp2 ;
k2n = kp1*k0*xp1 - kp1*k1*xp0 + k2*kp2 ;

/* ensure normalization to unit length */
den = sqrt(k0n*k0n + k1n*k1n + k2n*k2n) ;
k0 = k0n/den ; k1 = k1n/den ; k2 = k2n/den ;