module mie_theory

  implicit none

  private

  public :: CALLBH3

  contains

    !--------------------------------------------------------------------------

    subroutine CALLBH3(                &
      & Rmed, Rre, Rim, rad, wl,       &
      & QEXT, SSA, SSCoA, AF           &
      & )

      ! CALLBH calculates the size parameter X and relative 
      ! refreactive index refrel for a given sphere refractive 
      ! index, medium refractive index, radius, and free space 
      ! wavelenght.  
      ! It then calls BHMIE, the subroutine that computes 
      ! amplitude scattering matrix elements and efficiencies.  

      ! 1998.11.04
      ! Ver. 2

      ! -----CONSTANT-----
      ! number of angles between 0 and 90 degrees, matrix elements 
      ! calculated at 2*NANG-1 angles including 0, 90, and 180 
      ! degrees.  
      Integer, parameter :: NANG = 21

      ! -----INPUT-----
      ! (real) refractive index of surrounding medium
      real(8), intent(in) :: Rmed
      ! refractive index of sphere = Rre + i*Rim
      real(8), intent(in) :: Rre, Rim
      ! radius and wavelength  (same units)
      real(8), intent(in) :: rad, wl

      ! -----OUTPUT-----
      ! extinction efficiency
      real(8), intent(out) :: QEXT
      ! single scattering albedo and co-albedo
      real(8), intent(out) :: SSA
      real(8), intent(out) :: SSCoA
      ! asymmetry factor
      real(8), intent(out) :: AF

      ! -----LOCAL-----
      real(8) :: PI
      Integer :: NAN, i
      real(8) :: X
      real(8) :: REFRELr, S1r(200), S2r(200)
      real(8) :: REFRELi, S1i(200), S2i(200)
      real(8) :: S11(200), DANG, sum1, sum2, junk

      ! scattering efficiency
      real(8) :: QSCA
      ! (backscattered irradiance)/(total scattered irradiance)
      real(8) :: eBACK

      ! -----BEGIN-----
      ! set mathematical constant
      PI = acos(-1.0D0)

      ! calculation
      REFRELr = Rre/Rmed
      REFRELi = Rim/Rmed
      ! size parameter
      X=2.0D0*PI*rad*Rmed/wl
      ! calculate Mie scattering
      Call BHMIE2(                             &
        & X, REFRELr, REFRELi, NANG,           & ! (in)
        & S1r, S1i, S2r, S2i, QEXT, QSCA, junk & ! (out)
        & )
      ! calculate backscattering coefficient
      NAN=2*NANG-1
      DANG=PI/2.0D0/(NANG-1)
      ! phase function
      Do i = 1, NAN
        S11(i) =( (S1r(i)*S1r(i)+S1i(i)*S1i(i)) &
          &             +(S2r(i)*S2r(i)+S2i(i)*S2i(i)) ) / 2.0D0
      end Do

      ! integration 0-90 degree (forward)
      sum1 = 0.0D0
      Do i = 1, NANG-1
        sum1 = sum1 + (S11(i)+S11(i+1))*(sin(DANG*(i-1))+sin(DANG*i))
      end Do
      ! integration 90-180 degree (backward)
      sum2 = 0.0D0
      Do i = NANG, NAN-1
        sum2 = sum2 + (S11(i)+S11(i+1))*(sin(DANG*(i-1))+sin(DANG*i))
      end Do
      ! backscattering efficiency
      eBACK = (sum2/(sum1+sum2))

      ! Calculation of single scattering albedo and co-albedo
      !
      SSA   = QSCA/QEXT
      SSCoA = (QEXT-QSCA)/QEXT


      ! Normalization of phase function
      !   Integration over solid angle (4 \pi) divided by 4 pi has to be one.
      !
      !   Integration over solid angle (4 \pi) divided by 4 pi
      sum1 = 0.0D0
      i = 1
      sum1 = sum1                                                 &
        & - S11(i)                                                &
        &   * 2.0 * PI                                            &
        &   * (   ( cos(DANG*(i  )) + cos(DANG*(i+1)) ) / 2.0D0   &
        &       - cos(DANG*(i  ))                               )
      Do i = 2, NAN-1
        sum1 = sum1                                                 &
          & - S11(i)                                                &
          &   * 2.0 * PI                                            &
          &   * (   ( cos(DANG*(i  )) + cos(DANG*(i+1)) ) / 2.0D0   &
          &       - ( cos(DANG*(i-1)) + cos(DANG*(i  )) ) / 2.0D0 )
      end Do
      i = NAN
      sum1 = sum1                                                 &
        & - S11(i)                                                &
        &   * 2.0 * PI                                            &
        &   * (   cos(DANG*(i  ))                                 &
        &       - ( cos(DANG*(i-1)) + cos(DANG*(i  )) ) / 2.0D0 )
      sum1 = sum1 / ( 4.0D0 * PI )
      !   Normalization
      S11(1:NAN) = S11(1:NAN) / max( sum1, 1.0d-100 )


      !   Test output of phase function
!!$      Do i = 1, NAN
!!$        write( 60, * ) S11(i) * cos(DANG*(i-1)), S11(i) * sin(DANG*(i-1))
!!$      end Do


      ! Calculation of asymmetry factor
      !
      AF = 0.0D0
      i = 1
      AF = AF                                                            &
        & - S11(i)                                                       &
        &   * (   ( ( cos(DANG*(i  )) + cos(DANG*(i+1)) ) / 2.0D0 )**2   &
        &       - ( cos(DANG*(i  ))                               )**2 ) &
        &   / 2.0D0
      Do i = 2, NAN-1
        AF = AF                                                            &
          & - S11(i)                                                       &
          &   * (   ( ( cos(DANG*(i  )) + cos(DANG*(i+1)) ) / 2.0D0 )**2   &
          &       - ( ( cos(DANG*(i-1)) + cos(DANG*(i  )) ) / 2.0D0 )**2 ) &
          &   / 2.0D0
      end Do
      i = NAN
      AF = AF                                                            &
        & - S11(i)                                                       &
        &   * (   ( cos(DANG*(i  ))                               )**2   &
        &       - ( ( cos(DANG*(i-1)) + cos(DANG*(i  )) ) / 2.0D0 )**2 ) &
        &   / 2.0D0
      AF = AF / 2.0D0


      Return
    End subroutine CALLBH3

    ! -------------------------------------------------------------------------

    Subroutine BHMIE2(                         &
      & X, REFRELr, REFRELi, NANG,             & ! (in)
      & S1r, S1i, S2r, S2i, QEXT, QSCA, QBACK  & ! (out)
      & )

      ! Subroutine BHMIE calculates amplitude scattering matrix 
      ! elements and efficiencies for extinction, total scattering 
      ! and backscattering for a given size parameter and 
      ! relative refractive index.  

      ! 1998.10.25
      ! Ver. 1

      ! -----INPUT-----
      real(8), intent(in) :: X
      Integer, intent(in) :: NANG
      real(8), intent(in) :: REFRELr, REFRELi

      ! -----OUTPUT-----
      real(8), intent(out) :: S1r(200), S1i(200), S2r(200), S2i(200)
      real(8), intent(out) :: QEXT, QSCA, QBACK

      ! -----LOCAL-----
      real(8) :: XSTOP,YMOD,DANG,RN,CHI0,CHI1
      real(8) :: APSI0,APSI1,APSI,CHI,P,FN,T
      real(8) :: AMU(100),THETA(100),PI(100)
      real(8) :: TAU(100),PI0(100),PI1(100)
      integer, parameter :: NArrayMax = 10000
      real(8) :: Dr(NArrayMax),Yr,XIr,XI0r,XI1r,ANr,BNr
      real(8) :: Di(NArrayMax),Yi,XIi,XI0i,XI1i,ANi,BNi
      real(8) :: PSI0,PSI1,PSI,DN,DX
      Integer :: NSTOP,NMX,J,NN,N,JJ
      real(8) :: D1r,D1i,D2r,D2i
      real(8) :: AN1r,AN1i,AN2r,AN2i,AN3r,AN3i
      real(8) :: BN1r,BN1i,BN2r,BN2i,BN3r,BN3i

      ! -----BEGIN-----
      DX=X
      Yr=X*REFRELr
      Yi=X*REFRELi
      ! series terminated after NSTOP terms
      XSTOP=X+4.0D0*(X**(1.0D0/3.0D0))+2.0D0
      NSTOP=XSTOP
      YMOD=((Yr*Yr+Yi*Yi)**(0.5D0))
      If (XSTOP.ge.YMOD) then
        NMX=XSTOP+15
      Else
        NMX=YMOD+15
      End If
      DANG=3.1415927410126D0/2.0D0/FLOAT(NANG-1)
      DO J=1,NANG
        THETA(J)=(FLOAT(J)-1.0D0)*DANG
        AMU(J)=COS(THETA(J))
      end DO
      ! logarithmic derivative D(J) calculated by downward recurrence 
      ! beginning with initial value 0.0+i*0.0 at J=NMX
      if ( NMX > NArrayMax ) stop 'Array size is not enough.'
      Dr(NMX)=0.0D0
      Di(NMX)=0.0D0
      NN=NMX-1
      DO N=1,NN
        RN=NMX-N+1
        D1r=RN*Yr/(Yr*Yr+Yi*Yi)
        D1i=-RN*Yi/(Yr*Yr+Yi*Yi)
        D2r=Dr(NMX-N+1)+D1r
        D2i=Di(NMX-N+1)+D1i
        Dr(NMX-N)=D1r-D2r/(D2r*D2r+D2i*D2i)
        Di(NMX-N)=D1i+D2i/(D2r*D2r+D2i*D2i)
      end DO
      DO J=1,NANG
        PI0(J)=0.0D0
        PI1(J)=1.0D0
      end DO
      NN=2*NANG-1
      DO J=1,NN
        S1r(J)=0.0D0
        S1i(J)=0.0D0
        S2r(J)=0.0D0
        S2i(J)=0.0D0
      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
      XI0r=APSI0
      XI0i=-CHI0
      XI1r=APSI1
      XI1i=-CHI1
      QSCA=0.0D0
      N=1
200   CONTINUE
      DN=N
      RN=N
      FN=(2.0D0*RN+1.0D0)/(RN*(RN+1.0D0))
      PSI=(2.0D0*DN-1.0D0)*PSI1/DX-PSI0
      APSI=PSI
      CHI=(2.0D0*RN-1.0D0)*CHI1/X - CHI0
      XIr=APSI
      XIi=-CHI
      if ( N > NArrayMax ) stop 'Array size is not enough.'
      AN1r=(Dr(N)*REFRELr+Di(N)*REFRELi) &
        &          /(REFRELr*REFRELr+REFRELi*REFRELi) + RN/X
      AN1i=(Di(N)*REFRELr-Dr(N)*REFRELi) &
        &          /(REFRELr*REFRELr+REFRELi*REFRELi)
      AN2r=AN1r*APSI - APSI1
      AN2i=AN1i*APSI
      AN3r=AN1r*XIr-AN1i*XIi - XI1r
      AN3i=AN1r*XIi+AN1i*XIr - XI1i
      ANr=(AN2r*AN3r+AN2i*AN3i)/(AN3r*AN3r+AN3i*AN3i)
      ANi=(AN2i*AN3r-AN2r*AN3i)/(AN3r*AN3r+AN3i*AN3i)
      BN1r=REFRELr*Dr(N)-REFRELi*Di(N)+RN/X
      BN1i=REFRELr*Di(N)+REFRELi*Dr(N)
      BN2r=BN1r*APSI - APSI1
      BN2i=BN1i*APSI
      BN3r=BN1r*XIr-BN1i*XIi - XI1r
      BN3i=BN1i*XIr+BN1r*XIi - XI1i
      BNr=(BN2r*BN3r+BN2i*BN3i)/(BN3r*BN3r+BN3i*BN3i)
      BNi=(BN2i*BN3r-BN2r*BN3i)/(BN3r*BN3r+BN3i*BN3i)
      QSCA=QSCA+(2.0D0*RN+1.0D0) &
        &              *((ANr*ANr+ANi*ANi)+(BNr*BNr+BNi*BNi))
      DO J=1,NANG
        JJ=2*NANG-J
        PI(J)=PI1(J)
        TAU(J)=RN*AMU(J)*PI(J) - (RN+1.0D0)*PI0(J)
        P=(-1.0D0)**(N-1)
        S1r(J)=S1r(J)+FN*(ANr*PI(J)+BNr*TAU(J))
        S1i(J)=S1i(J)+FN*(ANi*PI(J)+BNi*TAU(J))
        T=(-1.0D0)**N
        S2r(J)=S2r(J)+FN*(ANr*TAU(J)+BNr*PI(J))
        S2i(J)=S2i(J)+FN*(ANi*TAU(J)+BNi*PI(J))
        IF(J.EQ.JJ) THEN
          cycle
        END IF
        S1r(JJ)=S1r(JJ)+FN*(ANr*PI(J)*P+BNr*TAU(J)*T)
        S1i(JJ)=S1i(JJ)+FN*(ANi*PI(J)*P+BNi*TAU(J)*T)
        S2r(JJ)=S2r(JJ)+FN*(ANr*TAU(J)*T+BNr*PI(J)*P)
        S2i(JJ)=S2i(JJ)+FN*(ANi*TAU(J)*T+BNi*PI(J)*P)
      end DO
      PSI0=PSI1
      PSI1=PSI
      APSI1=PSI1
      CHI0=CHI1
      CHI1=CHI
      XI1r=APSI1
      XI1i=-CHI1
      N=N+1
      RN=N
      DO J=1,NANG
        PI1(J)=((2.0D0*RN-1.0D0)/(RN-1.0D0))*AMU(J)*PI(J)
        PI1(J)=PI1(J)-RN*PI0(J)/(RN-1.0D0)
        PI0(J)=PI(J)
      end DO
      IF ((N-1-NSTOP).lt.0) then
        GO TO 200
      END IF
      QSCA=(2.0D0/(X*X))*QSCA
      QEXT=(4.0D0/(X*X))*S1r(1)
      QBACK=(4.0D0/(X*X)) &
        &    *(S1r(2*NANG-1)*S1r(2*NANG-1) &
        &      +S1i(2*NANG-1)*S1i(2*NANG-1))

      RETURN
    END Subroutine BHMIE2

    !--------------------------------------------------------------------------

  end module mie_theory
