*   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 */


      function hgphasefunc(g, r) 

      real g2,a,b, hgphasefunc

* check for isotropic case 
      if (g .eq. 0.0) then
         hgphasefunc = 2.0*r - 1.0
         return
      endif

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



* ---------------------------------------------------------------
*   The following code FRAGMENTS illustrates the process of converting the
*   scattering direction in photon-relative coordinates to a new
*   propagation direction in model coordinates.   You will need to paste these fragments
*   into your program at the appropriate locations.
*
*   Note that k0, k1, and k2 are the three (real) elements of the unit direction vector k
      

      real ran1, cosTheta,sinTheta,hgphasefunc,Phi,cosPhi,sinPhi
      real den,k0,k1,k2,xp0,xp1,xp2



* get cos and sin of scattering angle Theta  */
      cosTheta = hgphasefunc(g,ran1(idum)) 
      sinTheta = sqrt(1.0-cosTheta*cosTheta) 

* get cos(Phi) and sin(Phi) */

      Phi = TWOPI*ran1(idum)
      cosPhi = cos(Phi) 
      sinPhi = sin(Phi) 

* define photon-relative x direction */
      den = sqrt(k0*k0 + k1*k1) 


      if (den .lt. 0.0001) then
*    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 
      endif

* new direction in photon-relative coordinates */
      kp0 = sinTheta*cosPhi 
      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