!+ Subroutine to calculate Mie scattering cross-sections and asymmetry.
!
! Method:
!	This routine is closely based on that given in "Absorption 
!	and Scattering of Light by Small Particles" by  C. F. Bohren 
!	and Huffman. It uses downward recurrence.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE MIE_SCATTER(IERR, X, N_RELATIVE
     &   , Q_SCATTER, Q_EXT, ASYMMETRY
     &   , L_PHASE, N_ANGLE, MU_ANGLE, S1, S2)
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS.
      COMPLEX*16	!, INTENT(IN)
     &     N_RELATIVE
!             RELATIVE REFRACTIVE INDEX
      COMPLEX*16	!, INTENT(OUT)
     &     S1(NPD_SCATT_ANGLE)
!             ELEMENT OF SCATTERING MATRIX
     &   , S2(NPD_SCATT_ANGLE)
!             ELEMENT OF SCATTERING MATRIX
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     N_ANGLE
!             NUMBER OF SCATTERING ANGLES
      LOGICAL 	!, INTENT(IN)
     &     L_PHASE
!             CALCULATE PHASE FUNCTION
      REAL*8	!, INTENT(IN)
     &     X
!             SIZE PARAMETER
     &   , MU_ANGLE(NPD_SCATT_ANGLE)
!             COSINES OF SCATTERING ANGLES
      REAL*8	!, INTENT(OUT)
     &     Q_SCATTER
!             SCATTERING EFFICIENCY
     &   , Q_EXT
!             EXTINCTION EFFICIENCY
     &   , ASYMMETRY
!             ASYMMETRY FACTOR
!
!     LOCAL ARGUMENTS.
      COMPLEX*16
     &     D(NPD_LOG_DERIV)
!             LOGARITHMIC DERIVATIVES
     &   , C
!             COMPLEX NUMBER OF TERM
     &   , Y
!             SIZE PARAMETER*REFRACTIVE INDEX
     &   , KSI_N
!             RICCATI-BESSEL FUNCTION AT N
     &   , KSI_NM1
!             RICCATI-BESSEL FUNCTION AT N-1
     &   , KSI_NM2
!             RICCATI-BESSEL FUNCTION AT N-2
     &   , A
!             COEFFICIENT IN SERIES
     &   , B
!             COEFFICIENT IN SERIES
     &   , A_OLD
!             PREVIOUS COEFFICIENT
     &   , B_OLD
!             PREVIOUS COEFFICIENT
      INTEGER
     &     N_TERM
!             NUMBER OF TERMS COMPUTED
     &   , N_DOWN
!             START OF DOWNWARD RECURRENCE	
     &   , N_Y
!             INTEGRAL PART OF Y
     &   , I_TERM
!             TERM INDEX
     &   , I_ANGLE
!             ANGLE INDEX
      REAL*8
     &     R
!             NUMERICAL FACTOR
     &   , RX
!             "
     &   , FACTOR
!             "
     &   , X2
!             "
     &   , PSI_N
!             RICCATI-BESSEL FUNCTION AT N
     &   , PSI_NM1
!             RICCATI-BESSEL FUNCTION AT N-1
     &   , PSI_NM2
!             RICCATI-BESSEL FUNCTION AT N-2
     &   , CHI_N
!             RICCATI-BESSEL FUNCTION AT N
     &   , CHI_NM1
!             RICCATI-BESSEL FUNCTION AT N-1
     &   , CHI_NM2
!             RICCATI-BESSEL FUNCTION AT N-2
     &   , PI_N(NPD_SCATT_ANGLE)
!             RICCATI-BESSEL FUNCTION AT N
     &   , PI_NM1(NPD_SCATT_ANGLE)
!             RICCATI-BESSEL FUNCTION AT N-1
     &   , PI_NM2(NPD_SCATT_ANGLE)
!             RICCATI-BESSEL FUNCTION AT N-2
     &   , TAU_N(NPD_SCATT_ANGLE)
!             RICCATI-BESSEL FUNCTION AT N
!
!
      Y=X*N_RELATIVE
!
!     SET THE NUMBER OF TERMS REQUIRED IN THE SCATTERING SERIES.
      N_TERM=INT(2.0D+00+4.0D+00*X**3.333D-01+X)
!     SET THE TERM AT WHICH DOWNWARD RECURRENCE BEGINS TO BE
!     LARGE ENOUGH.
      N_Y=1.1D+00*ABS(Y)
      N_DOWN=MAX(N_TERM, N_Y)+15
!
!     CHECK ARRAY BOUNDS.
      IF (N_DOWN.GT.NPD_LOG_DERIV) THEN
         WRITE(IU_ERR, '(/A)') '*** ERROR: TOO MANY LOGARITHMIC '
     &      //'DERIVATIVES ARE REQUIRED.'
         WRITE(IU_ERR, '(A, 1PE12.5)') 'SIZE PARAMETER = ', X
         WRITE(IU_ERR, '(A, I8)') 'RECURRENCE STARTS AT ', N_DOWN
         WRITE(IU_ERR, '(A)')
     &      'INCREASE NPD_LOG_DERIV AND RE-COMPILE.'
         IERR=I_ERR_FATAL
         RETURN
      ENDIF

!
!     CALCULATE THE LOGARITHMIC DERIVATIVES, D, BY DOWNWARD RECURRENCE.
!     D(N_DOWN) IS INITIALIZED TO 0.
      D(N_DOWN)=DCMPLX(0.0D+00, 0.0D+00)
      DO I_TERM=N_DOWN-1, 1, -1
         C=(DBLE(I_TERM)+1.0D+00)/Y
         D(I_TERM)=C-1.0D+00/(D(I_TERM+1)+C)
      ENDDO
!
!     MAIN LOOP OF SERIES FOLLOWS: RICCATI-BESSEL FUNCTIONS 
!     ARE CALCULATED BY UPWARD RECURRENCE.
!
!     INITIALIZTION: THIS STARTS WITH AT -1, I.E. N-2 FOR N=1 BELOW.
      PSI_NM2=COS(X)
      PSI_NM1=SIN(X)
      CHI_NM2=-SIN(X)
      CHI_NM1=COS(X)
      KSI_NM2=DCMPLX(PSI_NM2, -CHI_NM2)
      KSI_NM1=DCMPLX(PSI_NM1, -CHI_NM1)
      Q_SCATTER=0.0D+00
      Q_EXT=0.0D+00
      ASYMMETRY=0.0D+00
!     INITIALIZATION OF ANGULARLY-DEPENDENT TERMS.
      IF (L_PHASE) THEN
         DO I_ANGLE=1, N_ANGLE
            PI_NM2(I_ANGLE)=0.0D+00
            PI_NM1(I_ANGLE)=0.0D+00
            S1(I_ANGLE)=0.0D+00
            S2(I_ANGLE)=0.0D+00
         ENDDO
      ENDIF
!
      I_TERM=0
!
100      I_TERM=I_TERM+1
         R=DBLE(I_TERM)
         RX=R/X
         FACTOR=(2.0D+00*R+1.0D+00)/(R*(R+1.0D+00))
!        COMPUTE NEW PSI AND KSI BY RECURRENCE.
         PSI_N=(2.0D+00*R-1.0D+00)*PSI_NM1/X-PSI_NM2
         CHI_N=(2.0D+00*R-1.0D+00)*CHI_NM1/X-CHI_NM2
         KSI_N=DCMPLX(PSI_N, -CHI_N)
         A=((D(I_TERM)/N_RELATIVE+RX)*PSI_N-PSI_NM1)
     &      /((D(I_TERM)/N_RELATIVE+RX)*KSI_N-KSI_NM1)
         B=((N_RELATIVE*D(I_TERM)+RX)*PSI_N-PSI_NM1)
     &      /((N_RELATIVE*D(I_TERM)+RX)*KSI_N-KSI_NM1)
!
         Q_SCATTER=Q_SCATTER+(2.0D+00*R+1.0D+00)
     &      *((ABS(A))**2+(ABS(B))**2)
         Q_EXT=Q_EXT+(2.0D+00*R+1.0D+00)*DBLE(A+B)
!
         IF (I_TERM.GT.1) THEN
            ASYMMETRY=ASYMMETRY+(R*R-1.0D+00)
     &         *DBLE(A_OLD*CONJG(A)+B_OLD*CONJG(B))
     &         /R+(2.0D+00*R-1.0D+00)*DBLE(A_OLD*CONJG(B_OLD))
     &         /(R*(R-1.0D+00))
         ENDIF
!
         IF (L_PHASE) THEN
!           FUNCTIONS PI AND TAU ARE REQUIRED FOR ELEMENTS OF SCATTERING
!           MATRIX.
            DO I_ANGLE=1, N_ANGLE
               IF (I_TERM.EQ.1) THEN
                  PI_N(I_ANGLE)=1.0D+00
               ELSE
                  PI_N(I_ANGLE)=((2.0D+00*R-1.0D+00)/(R-1.0D+00))
     &               *MU_ANGLE(I_ANGLE)
     &               *PI_NM1(I_ANGLE)-(R/(R-1.0D+00))*PI_NM2(I_ANGLE)
               ENDIF
               TAU_N(I_ANGLE)=R*MU_ANGLE(I_ANGLE)*PI_N(I_ANGLE)
     &            -(R+1.0D+00)*PI_NM1(I_ANGLE)
               S1(I_ANGLE)=S1(I_ANGLE)+FACTOR*(A*PI_N(I_ANGLE)
     &            +B*TAU_N(I_ANGLE))
               S2(I_ANGLE)=S2(I_ANGLE)+FACTOR*(A*TAU_N(I_ANGLE)
     &            +B*PI_N(I_ANGLE))
            ENDDO
         ENDIF
!
!        TEST FOR NEXT TERM OF SERIES AND LOOP AROUND IF REQUIRED.
         IF (I_TERM.LT.N_TERM) THEN
!           ADVANCE STORED TERMS FOR RECURRENCE.
            A_OLD=A
            B_OLD=B
            PSI_NM2=PSI_NM1
            PSI_NM1=PSI_N
            CHI_NM2=CHI_NM1
            CHI_NM1=CHI_N
            KSI_NM1=KSI_N
            IF (L_PHASE) THEN
               DO I_ANGLE=1, N_ANGLE
                  PI_NM2(I_ANGLE)=PI_NM1(I_ANGLE)
                  PI_NM1(I_ANGLE)=PI_N(I_ANGLE)
               ENDDO
            ENDIF
            GOTO 100
         ENDIF
!
          
      X2=1.0D+00/X**2
      Q_SCATTER=2.0D+00*Q_SCATTER*X2
      Q_EXT=2.0D+00*Q_EXT*X2
      ASYMMETRY=4.0D+00*ASYMMETRY*X2/Q_SCATTER
    
!
!
      RETURN
      END
