pro bhmie2,x,refrel,nang,s1,s2,qext,qsca,qback2

;----------------------------------------------------------------------------
;
; Purpose:
;
;   Subroutine bhmie2 calculates amplitude scattering matrix
;   elements and efficiencies for extinction, total scattering
;   and backscattering for a given size parameter and
;   refractive index,  according to Mie theory.
; 
; Parameters:
;
;  input:
;    x      = 2 pi r / wavelength  (r is particle radius,  
;             r and wavelength should be in the same units)
;    refrel = complex refractive index
;    nang   = # angles for integration (20 is OK)
;
;  output:
;    s1    =  the perpendicular "scattering functions", at each theta
;    s2    =  the parallel           "         "         "   "   " 
;    qext  = extinction efficiency
;    qsca  = scattering      "
;    qback2= backscattering efficiency (1/sr)
;
; Source: Mark Hervig
;
;-----------------------------------------------------------------------------

      amu = fltarr(nang+1) & theta=amu & pi=amu & tau=amu & pi0=amu & pi1=amu

      s1 = complexarr(2*nang) & s2=s1
 
      d = complexarr(3000)
      y = complex(0.,0.) & xi=y & xi0=y & xi1=y & an=y & bn=y
      ccan=y & ccbn=y & anmi1=y & bnmi1=y & sumbak=y     
      
;     double precision psi0,psi1,psi,dn,dx

      psi0=0. & psi1=0. & psi=0. & dn=0. & dx=0.

      psi0=0.d0 & psi1=0.d0 & psi=0.d0 & dn=0.d0 & dx=0.d0


;---- Mie x and y values.

      dx = x
      y  = x * refrel

;---- Series terminated after nstop terms

      xstop = x + 4. * x^0.3333 + 2.0
      nstop = xstop

;---- Will loop over nang angles.

      ymod = abs(y)
      nmx  = fix( max([xstop,ymod]) ) + 15
      dang = 1.570796327 / float(nang-1)

      for j = 1,nang do begin
        theta(j) = (float(j)-1.)*dang
        amu(j)   = cos(theta(j))
      endfor

;---- Logarithmic derivative d(j) calculated by downword
;     recurrence beginning with initial value 0.0 + i*0.0
;     at j = nmx

      d(nmx) = complex(0.0,0.0)
      nn     = nmx-1

      for n = 1,nn do begin
        rn      = nmx-n+1
        d(nmx-n)= (rn/y) - (1./ (d(nmx-n+1) + rn/y) )
      endfor

      pi0(*) = 0.0
      pi1(*) = 1.0
      s1(*)  = complex(0.0,0.0)
      s2(*)  = complex(0.0,0.0)

;---- Riccati-Bessel functions with real argument x
;     calculated by upward recurrence

      psi0 = cos(dx)
      psi1 = sin(dx)
      chi0 = -sin(x)
      chi1 = cos(x)
      apsi0= psi0
      apsi1= psi1
      xi0  = complex(apsi0,-chi0)
      xi1  = complex(apsi1,-chi1)
      qsca = 0.0
      g1   = 0.0
      g2   = 0.0
      n    = 1
      sumbak = complex(0.0,0.0)

;---- Loop over the terms n in the Mie series

      while ( (n-1-nstop) lt 0) do begin
 
        dn   = n
        rn   = n
        fn   = (2.*rn+1.)/(rn*(rn+1.))
        ffn  = (rn-1.)*(rn+1.)/rn
        psi  = (2.*dn-1.)*psi1/dx-psi0
        apsi = psi
        chi  = (2.*rn-1.)*chi1/x - chi0
        xi   = complex(apsi,-chi)
        an   = (d(n)/refrel+rn/x)*apsi - apsi1
        an   = an/((d(n)/refrel+rn/x)*xi-xi1)
        bn   = (refrel*d(n)+rn/x)*apsi - apsi1
        bn   = bn/((refrel*d(n)+rn/x)*xi - xi1)
        ccan = conj(an)
        ccbn = conj(bn)
        g2   = g2+fn*float(an*ccbn)

        if ((n-1) gt 0) then g1 = g1 + ffn*float(anmi1*ccan + bnmi1*ccbn)

        qsca   = qsca + (2.*rn+1.)*(abs(an)*abs(an)+abs(bn)*abs(bn))

;----   sum for back xsection

        sumbak = sumbak + (-1.0^n) *0.5* (2.0*rn + 1.0) * (an - bn) 

;----   loop over angles

        for j = 1,nang do begin
          jj     = 2*nang-j
          pi(j)  = pi1(j)
          tau(j) = rn *amu(j) *pi(j) - (rn +1.) *pi0(j)
          p      = (-1.)^(n-1)
          s1(j)  = s1(j) +fn *(an *pi(j) +bn *tau(j))
          t      = (-1.)^n
          s2(j)  = s2(j) + fn *(an *tau(j) +bn *pi(j))
          if (j eq jj) then goto, jump789
          s1(jj) = s1(jj) + fn * (an *pi(j) *p +bn *tau(j) *t)
          s2(jj) = s2(jj) + fn * (an *tau(j) *t +bn *pi(j) *p)
          jump789:
        endfor

        psi0 = psi1
        psi1 = psi
        apsi1= psi1
        chi0 = chi1
        chi1 = chi
        xi1  = complex(apsi1,-chi1)
        n    = n+1
        rn   = n

        pi1  = ((2.*rn-1.) / (rn-1.)) *amu *pi
        pi1  = pi1 - rn * pi0 / (rn-1.)
        pi0  = pi

        anmi1 = an
        bnmi1 = bn

      endwhile

      qsca  = (2./(x*x)) * qsca
      gfac  = (4./(x*x*qsca)) * (g1+g2)
      qext  = (4./(x*x)) * float(s1(1))
      qback = (4./(x*x)) * abs( s1(2*nang-1)) * abs(s1(2*nang-1) )

;---- qback2 is in 1/sr,  need to multiply by (wavl/2pi)^2 to get cm^2/sr

      qback2 = (abs(sumbak))^2                                       

      return
      end