      Program MieScat
      Implicit NONE


c -----CONSTANT-----
c mathematical constant
      Real*8  PI

c -----INPUT-----
c (real) refractive index of surrounding medium
      Real*8  Rmed
c refractive index of sphere = Rre + i*Rim
      Real*8  Rre, Rim
c radius and wavelength  [m]
      Real*8  rad, wl

c -----OUTPUT-----
c extinction efficiency
      Real*8  QEXT
c scattering efficiency
      Real*8  QSCA
c (backscattered irradiance)/(total scattered irradiance)
      Real*8  eBACK

c size parameter
      Real*8  X
c extinction cross section  [m2]
      Real*8  SEXT

c single scattering albedo and co-albedo
      Real*8  ssa, cssa

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

c set input parameters
c   (real) refractive index of surrounding medium
      Rmed = 1.00045D0
c   refractive index of sphere = Rre + i*Rim
      Rre = 1.427D0
      Rim = 1.19D-6
c   radius and wavelength  (same units)
      rad = 1.0D-6
      wl  = 1.0D-6

c output parameters
      Write(*,*) '-----input-----'
      Write(*,*) 'refractive index of surrounding medium (real)',
     &           Rmed
      Write(*,*) 'refractive index of particle (real)          ',
     &           Rre
      Write(*,*) 'refractive index of particle (imaginary)     ',
     &           Rim
      Write(*,*)
      Write(*,*) 'particle radius [m]                          ',
     &           rad
      Write(*,*) 'wavelength [m]                               ',
     &           wl
      Write(*,*)

c calculation
      X = 2.0D0*PI*rad*Rmed/wl
      Call CALLBH3(Rmed,Rre,Rim,rad,wl,QEXT,QSCA,eBACK)
      SEXT = PI*(rad*rad)*QEXT
      ssa  = QSCA/QEXT
      cssa = (QEXT-QSCA)/QEXT

c output results
      Write(*,*) '-----output-----'
      Write(*,*) 'size parameter               ', X
      Write(*,*)
      Write(*,*) 'extinction efficiency        ', QEXT
      Write(*,*) 'scattering efficiency        ', QSCA
      Write(*,*) 'extinction cross section [m2]', SEXT
      Write(*,*)
      Write(*,*) 'single scattering albedo     ', ssa
      Write(*,*) 'co-single scattering albedo  ', cssa


      Stop
      End





      Subroutine CALLBH3(Rmed,Rre,Rim,rad,wl,QEXT,QSCA,eBACK)
      Implicit NONE

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

c       1998.11.04
c                     Ver. 2

c -----CONSTANT-----
c   number of angles between 0 and 90 degrees, matrix elements 
c   calculated at 2*NANG-1 angles including 0, 90, and 180 
c   degrees.  
      Integer NANG
      Parameter (NANG=21)

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

c -----OUTPUT-----
c extinction efficiency
      Real*8  QEXT
c scattering efficiency
      Real*8  QSCA
c (backscattered irradiance)/(total scattered irradiance)
      Real*8  eBACK

c -----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

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

c calculation
        REFRELr = Rre/Rmed
        REFRELi = Rim/Rmed
c   size parameter
        X=2.0D0*PI*rad*Rmed/wl
c   calculate Mie scattering
        Call BHMIE2(X,REFRELr,REFRELi,NANG,S1r,S1i,S2r,S2i,
     &              QEXT,QSCA,junk)
c   calculate backscattering coefficient
        NAN=2*NANG-1
        DANG=PI/2.0D0/(NANG-1)
c     phase function
        Do 110 i = 1, NAN
          S11(i) =( (S1r(i)*S1r(i)+S1i(i)*S1i(i))
     &             +(S2r(i)*S2r(i)+S2i(i)*S2i(i)) ) / 2.0D0
 110    Continue
c     integration 0-90 degree (forward)
        sum1 = 0.0D0
        Do 120 i = 1, NANG-1
          sum1 = sum1
     &         + (S11(i)+S11(i+1))*(sin(DANG*(i-1))+sin(DANG*i))
 120    Continue
c     integration 90-180 degree (backward)
        sum2 = 0.0D0
        Do 130 i = NANG, NAN-1
          sum2 = sum2
     &         + (S11(i)+S11(i+1))*(sin(DANG*(i-1))+sin(DANG*i))
 130    Continue
c     backscattering efficiency
        eBACK = (sum2/(sum1+sum2))


      Return
      End





      Subroutine BHMIE2(X,REFRELr,REFRELi,NANG,S1r,S1i,S2r,S2i,
     &                  QEXT,QSCA,QBACK)
      Implicit NONE

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

c       1998.10.25
c                     Ver. 1

c -----INPUT-----
      Real*8  X
      Integer NANG
      Real*8  REFRELr, REFRELi

c -----OUTPUT-----
      Real*8  S1r(200), S1i(200), S2r(200), S2i(200)
      Real*8  QEXT, QSCA, QBACK

c -----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)
      Real*8  Dr(3000),Yr,XIr,XI0r,XI1r,ANr,BNr
      Real*8  Di(3000),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

c -----BEGIN-----
        DX=X
        Yr=X*REFRELr
        Yi=X*REFRELi
c 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 555 J=1,NANG
          THETA(J)=(FLOAT(J)-1.0D0)*DANG
          AMU(J)=COS(THETA(J))
 555    CONTINUE
c logarithmic derivative D(J) calculated by downward recurrence 
c beginning with initial value 0.0+i*0.0 at J=NMX
        Dr(NMX)=0.0D0
        Di(NMX)=0.0D0
        NN=NMX-1
        DO 120 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)
 120    CONTINUE
        DO 666 J=1,NANG
          PI0(J)=0.0D0
          PI1(J)=1.0D0
 666    CONTINUE
        NN=2*NANG-1
        DO 777 J=1,NN
          S1r(J)=0.0D0
          S1i(J)=0.0D0
          S2r(J)=0.0D0
          S2i(J)=0.0D0
 777    CONTINUE
c Riccati-Bessel functions with real argument X calculated by 
c 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
        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 789 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
              GO TO 789
          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)
 789    CONTINUE
        PSI0=PSI1
        PSI1=PSI
        APSI1=PSI1
        CHI0=CHI1
        CHI1=CHI
        XI1r=APSI1
        XI1i=-CHI1
        N=N+1
        RN=N
        DO 999 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)
 999    CONTINUE
        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
