! **********************************************************************
  subroutine bhmie(x,refrel,nang,s1,s2,qext,qsca,qback,gfac)
 
! Subroutine bhmie calculates amplitude scattering matrix
! elements and efficiencies for extinction, total scattering
! and backscattering for a given size parameter and
! relative refractive index.
 
! From the main program:
! refrel=cmplx(refre,refim)/refmed
 
  real :: amu(100),theta(100),pi(100),tau(100),pi0(100),pi1(100)
  complex :: d(3000),y,refrel,xi,xi0,xi1,an,bn,s1(200),s2(200)
  complex :: ccan,ccbn,anmi1,bnmi1
  double precision psi0,psi1,psi,dn,dx
   
! *****
! 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=cabs(y)
   nmx=amax1(xstop,ymod)+15
   dang=1.570796327/float(nang-1)
   do j=1,nang
    theta(j)=(float(j)-1.)*dang
    amu(j)=cos(theta(j))
   end do
 
! Logarithmic derivative d(j) calculated by downword
! recurrence beginning with initial value 0.0 + i*0.0
! at j = nmx
      d(nmx)=cmplx(0.0,0.0)
      nn=nmx-1
      do n=1,nn
       rn=nmx-n+1
       d(nmx-n)=(rn/y)-(1./(d(nmx-n+1)+rn/y))
      end do
      do j=1,nang
       pi0(j)=0.0
       pi1(j)=1.0
      end do
      nn=2*nang-1
      do j=1,nn
       s1(j)=cmplx(0.0,0.0)
       s2(j)=cmplx(0.0,0.0)
      end do
 
! Riccati-Bessel functions with real argument x
! calculated by upward recurrence
      psi0=dcos(dx)
      psi1=dsin(dx)
      chi0=-sin(x)
      chi1=cos(x)
      apsi0=psi0
      apsi1=psi1
      xi0=cmplx(apsi0,-chi0)
      xi1=cmplx(apsi1,-chi1)
      qsca=0.0
      g1=0.0
      g2=0.0
      n=1
 
! ******************
! Loop over the terms n in the Mie series
  200 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=cmplx(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=conjg(an)
      ccbn=conjg(bn)
      g2=g2+fn*real(an*ccbn)
      if (n-1) 55,55,56
  56  g1=g1+ffn*real(anmi1*ccan + bnmi1*ccbn)
  55  qsca=qsca+(2.*rn+1.)*(cabs(an)*cabs(an)+cabs(bn)*cabs(bn))
      do 789 j=1,nang
       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) go to 789
       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)
  789 continue
      psi0=psi1
      psi1=psi
      apsi1=psi1
      chi0=chi1
      chi1=chi
      xi1=cmplx(apsi1,-chi1)
      n=n+1
      rn=n
      do j=1,nang
       pi1(j)=((2.*rn-1.)/(rn-1.))*amu(j)*pi(j)
       pi1(j)=pi1(j)-rn*pi0(j)/(rn-1.)
       pi0(j)=pi(j)
      end do
      anmi1=an
      bnmi1=bn

      if (n-1-nstop) 200,300,300
 
! ******************
  300 qsca=(2./(x*x))*qsca
      gfac=(4./(x*x*qsca))*(g1+g2)
      qext=(4./(x*x))*real(s1(1))
      qback=(4./(x*x))*cabs(s1(2*nang-1))*cabs(s1(2*nang-1))
 
! *****
   return
   end
