pro bhmie, x,refrel,nang,s1,s2,qext,qsca,qback,gsca
;***********************************************************************
; Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine
;    to calculate scattering and absorption by a homogenous isotropic
;    sphere.
; Given:
;    X = 2*pi*a/lambda
;    REFREL = (complex refr. index of sphere)/(real index of medium)
;           ** This should be a COMPLEX() value in IDL ! **
;    NANG = number of angles between 0 and 90 degrees
;           (will calculate 2*NANG-1 directions from 0 to 180 deg.)
;           if called with NANG<2, will set NANG=2 and will compute
;           scattering for theta=0,90,180.
; Returns:
;    S1(0 : 2*(NANG-1)) = -i*f_22 (incid. E perp. to scatt. plane,
;                                scatt. E perp. to scatt. plane)
;    S2(0 : 2*(NANG-1)) = -i*f_11 (incid. E parr. to scatt. plane,
;                                scatt. E parr. to scatt. plane)
;    QEXT = C_ext/pi*a**2 = efficiency factor for extinction
;    QSCA = C_sca/pi*a**2 = efficiency factor for scattering
;    QBACK = (dC_sca/domega)/pi*a**2
;          = backscattering efficiency [NB: this is (1/4*pi) smaller
;            than the "radar backscattering efficiency"; see Bohren &
;            Huffman 1983 pp. 120-123]
;    GSCA = <cos(theta)> for scattering
;
; Original program taken from Bohren and Huffman (1983), Appendix A
; Modified by B.T.Draine, Princeton Univ. Obs., 90/10/26
; in order to compute <cos(theta)>
; Converted to IDL by P. J. Flatau, Scripps Inst. Oceanography. UCSD 96/11/14
; Optimized by Henry Throop, U. Colorado, 99/04/10
;
; 91/05/07 (BTD): Modified to allow NANG=1
; 91/08/15 (BTD): Corrected error (failure to initialize P)
; 91/08/15 (BTD): Modified to enhance vectorizability.
; 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1
; 91/08/15 (BTD): Changed definition of QBACK.
; 92/01/08 (BTD): Converted to full double precision and double complex
;                 eliminated 2 unneed lines of code
;                 eliminated redundant variables (e.g. APSI,APSI0)
;                 renamed RN -> EN = double precision N
;                 Note that DOUBLE COMPLEX and DCMPLX are not part
;                 of f77 standard, so this version may not be fully
;                 portable.  In event that portable version is
;                 needed, use src/bhmie_f77.f
; 93/06/01 (BTD): Changed AMAX1 to generic function MAX
; 96/11/14 (PJF): Converted to double precision (again)
;                 added  declaration standarization (strong typing)
;                 some code "polishing" using Toolpack
;                 Converted to IDL
; 99/04/10 (HBT): Vectorized several (not all) for-next loops.
;                 For X=1000, speed-up is ~ factor of 10.
;                 NB: This IDL code successfully reproduces Mie scattering
;                 in the geometric optics regime; e.g., fig. 5, X=600, of
;                 Hansen & Travis 1974.
; 00/03/15 (MTC): Dimension of output arrays S1 and S2 changed to 2*NANG-1,
;                 and indexing of these arrays now starts at zero
;                 Dimension of local array D changed to NMX+1
;                 Dimension of other local arrays changed to NANG, and
;                 indexing of these arrays now starts at zero
;                 NMX redefined as long integer instead of short integer
;                 Several Fortran style assignment statements corrected
;                 to have intended meaning in IDL
;                 Simplified parts of code, removing redundant statements
;                 Corrected error in calculation of NMX: replaced
;                 MAX(xstop,ymod) with (xstop > ymod)
; 00/03/21 (MTC): Made the following changes to speed up the code by
;                 a factor of 1.8 (NANG=2) and 3.0 (NANG=91):
;                 Added /NOZERO keyword to definition of logarithmic
;                 derivative array D, and reduced loop for calculation of D
;                 into one statement
;                 Reduced calculation of AN and BN to one statement each
;                 Simplified calculation of GSCA by using IMAGINARY function
;                 Improved vectorization of scattering intensity pattern
;                 calculations
;                 Removed unnecessary array definitions for PI and TAU
; 00/03/21 (MTC): Completely vectorized code
;***********************************************************************
;IDL example test code for "bhmie"
;pro test
; read,x,x1,x2
; refrel=dcomplex(x1,x2)
; bhmie, x,refrel,nang,s1,s2,qext,qsca,qback,gsca
; print,x, refrel, qext,qsca,qback,gsca
;stop
;end
;
;     .. Parameters ..
;      INTEGER mxnang, nmxx
;      PARAMETER (mxnang=1000,nmxx=150000)
;  The constants mxnang and nmxx are obsolete and are no longer used - MTC
;
;     .. Scalar Arguments ..
;      DOUBLE COMPLEX refrel
;      DOUBLE PRECISION gsca, qback, qext, qsca, x
;      INTEGER nang
;
;     .. Local Scalars ..
;      DOUBLE COMPLEX an, an1, bn, bn1, drefrl, xi, xi1, y
;      DOUBLE PRECISION chi, chi0, chi1, dang, dx, en, fn, p, pii, psi,$
;                      psi0, psi1, theta, xstop, ymod
;      INTEGER j, jj, n, nmx, nn, nstop
;
;     .. Intrinsic Functions ..
;      INTRINSIC abs, atan, cos, dble, dcmplx, max, sin
;
      nang = fix(nang) > 2	;redefined - MTC
;
;     .. Array Arguments ..
;      DOUBLE COMPLEX s1(2*mxnang-1), s2(2*mxnang-1)
;  s1 and s2 are defined here not in the subroutine calling  "bhmie"
      s1=dcomplexarr(2*nang-1,/nozero)
      s2=dcomplexarr(2*nang-1,/nozero)
;
;*** Obtain pi:
      pii = 4.D0*atan(1.D0)
      dx = double(x)
      drefrl = dcomplex(refrel)
      y = dx*drefrl
      ymod = abs(y)
;
;*** Series expansion terminated after NSTOP terms
;    Logarithmic derivatives calculated from NMX on down
      xstop = dx + 4.D0*dx^0.3333D0 + 2.D0
      nstop = long(xstop)
      nmx = long(xstop > ymod) + 15
;      print, ' nmx ', nmx		; commented out by HT
;
;     .. Local Arrays ..
;     DOUBLE COMPLEX d(nmxx)
      d=dcomplexarr(nmx+1,/nozero)
;
;     DOUBLE PRECISION amu(mxnang), pi(mxnang), pi0(mxnang),
;    &                 pi1(mxnang), tau(mxnang)
;
;*** Compute coefficient arrays PI and TAU needed for
;    scattering intensity pattern calculations
      dang = .5D0*pii/double(nang-1)
      amu = cos(dindgen(nang)*dang)				; ht
      pi = dblarr(nang,nstop+1,/nozero)
      pi(0,0) = dblarr(nang)
      pi(0,1) = dblarr(nang) + 1.D0
      for n=2,nstop do pi(0,n) = (double(2*n-1)*amu*pi(*,n-1) - $
                                  double(n)*pi(*,n-2))/double(n-1)
      tau = dblarr(nang,nstop+1,/nozero)
      for n=1,nstop do tau(0,n) = double(n)*amu*pi(*,n) - $
                                  double(n+1)*pi(*,n-1)
      pi = transpose(pi)
      tau = transpose(tau)
;
; BTD experiment 91/1/15: add one more term to series and compare results
;      NMX=AMAX1(XSTOP,YMOD)+16
; test: compute 7001 wavelengths between .0001 and 1000 micron
; for a=1.0micron SiC grain.  When NMX increased by 1, only a single
; computed number changed (out of 4*7001) and it only changed by 1/8387
; conclusion: we are indeed retaining enough terms in series!
;
;*** Logarithmic derivative D(J) calculated by downward recurrence
;    beginning with initial value (0.,0.) at J=NMX
;
      d(nmx) = dcomplex(0.D0,0.D0)
      for n=nmx-1,1,-1 do d(n) = double(n+1)/y - $
                                 1.D0/(d(n+1)+double(n+1)/y)
;
;*** Riccati-Bessel functions with real argument X
;    calculated by upward recurrence
;
      psi = dblarr(nstop+2,/nozero)
      psi(0) = cos(dx)
      psi(1) = sin(dx)
      for n=2,nstop+1 do psi(n) = double(2*n-3)*psi(n-1)/dx - psi(n-2)
      chi = dblarr(nstop+2,/nozero)
      chi(0) = -sin(dx)
      chi(1) = cos(dx)
      for n=2,nstop+1 do chi(n) = double(2*n-3)*chi(n-1)/dx - chi(n-2)
      xi = dcomplex(psi,-chi)
      p = double(lonarr(nstop)+1 - 2*(lindgen(nstop) mod 2)) ; 1.D0,-1.D0,...
;
;*** Compute coefficient arrays A and B needed for calculation
;    of scattering efficiency QSCA and asymmetry parameter <cos(theta)>,
;    or GSCA
      e = dindgen(nstop) + 1.D0
      da = e/dx
      db = d(1:nstop)*drefrl + da
      da = d(1:nstop)/drefrl + da
      a = (da*psi(2:nstop+1) - psi(1:nstop))/ $
          (da*xi(2:nstop+1) - xi(1:nstop))
      b = (db*psi(2:nstop+1) - psi(1:nstop))/ $
          (db*xi(2:nstop+1) - xi(1:nstop))
;
;*** Sum terms for QSCA and GSCA
      g = 2.D0*e + 1.D0
      f = g/(e*(e + 1.D0))
      qsca = total(g*(abs(a)^2 + abs(b)^2))
      gsca = total(f*(double(a)*double(b) + imaginary(a)*imaginary(b))) + $
             total((e(1:nstop-1) - 1.D0/e(1:nstop-1))* $
                   (double(a(0:nstop-2))*double(a(1:nstop-1)) + $
                    imaginary(a(0:nstop-2))*imaginary(a(1:nstop-1)) + $
                    double(b(0:nstop-2))*double(b(1:nstop-1)) + $
                    imaginary(b(0:nstop-2))*imaginary(b(1:nstop-1))))
;
;*** Now calculate scattering intensity pattern
;    First do angles from 0 to 90
      i = dblarr(nang)+1.D0
      fa = f*a
      fb = f*b
      pfa = (p*fa)#i
      pfb = (p*fb)#i
      fa = fa#i
      fb = fb#i
      s1(0) = total(fa*pi(1:nstop,*) + fb*tau(1:nstop,*),1)
      s2(0) = total(fa*tau(1:nstop,*) + fb*pi(1:nstop,*),1)
;
;*** Now do angles greater than 90 using PI and TAU from
;    angles less than 90.
;    P=1 for N=1,3,...; P=-1 for N=2,4,...
      s1(nang) = rotate([total(pfa*pi(1:nstop,0:nang-2) - $
                                    pfb*tau(1:nstop,0:nang-2),1)],2)
      s2(nang) = rotate([total(pfb*pi(1:nstop,0:nang-2) - $
                                    pfa*tau(1:nstop,0:nang-2),1)],2)
;
;*** Have summed sufficient terms.
;    Now compute QSCA,QEXT,QBACK,and GSCA
      gsca = 2.D0*gsca/qsca
      qsca = (2.D0/ (dx*dx))*qsca
      qext = (4.D0/ (dx*dx))*double(s1(0))
      qback = (abs(s1(2*nang-2))/dx)^2/pii
      RETURN
      END