!! Original program was written by George L. Hashimoto
!! program name was 'clascat6.f'
!! This program was modified from 'clascat6.f' by M. Onishi, 2017/01/18
!!

module mod_clscat6

  implicit none

  private

  public :: miescat
  public :: clscat6

contains

  subroutine miescat(                              & ! 2017/01/23 onishi
    & WaveNum, Rmean, Rsigma, Rmed, Rre, Rim,      & ! (in )
    & ExtCrossSec, SingleScatA, AsymFact           ) ! (out)

! -----constant-----
    Include 'clscat6.con'

! -----variables-----
    real(8), intent(in ) :: WaveNum
    real(8), intent(in ) :: Rmean
    real(8), intent(in ) :: Rsigma
    real(8), intent(in ) :: Rmed
    real(8), intent(in ) :: Rre
    real(8), intent(in ) :: Rim
    real(8), intent(out) :: ExtCrossSec
    real(8), intent(out) :: SingleScatA
    real(8), intent(out) :: AsymFact

! number of angles between 0 and 90 degrees, matix elements
! calculated at 2*NANG-1 angles including 0, 90, and 180 degrees.
    integer :: NANG
! number of calculation for change direction ratio
    integer :: NDCR
! accuracy
    real(8) :: EPS
! wavelength
    real(8) :: wl
! single scattering co-albedo
    Real*8  cSA
! direction change ratio
! (direction(cosine of zenith angle), ratio)
    Real*8  mu(NDCRMAX), bs(NDCRMAX)
! phase function
    Real*8  ANG(NXMAX), SP(NXMAX)

    integer :: k


    NANG = 91!NXMAX
    NDCR = NDCRMAX
    EPS = 1.0D-6

    if( WaveNum .gt. 0.0_8 ) then
      wl = 1.0d+0/WaveNum
    else
      print *, 'WaveNum is invalid value:,', WaveNum
      stop 
    end if

    call clscat6(                                           &
      & wl, Rmean, Rsigma, NANG, NDCR, EPS, Rmed, Rre, Rim, & ! (in )
      & ExtCrossSec, SingleScatA, cSA, mu, bs, ANG, SP      ) ! (out)

    !print *, 'OptDep', OptDep
    !print *, 'SingleScatA', SingleScatA
    !print *, '==SP'
    AsymFact = 0.0_8
    do k = 1, 2*NANG-1-1
      !print *, cos(ANG(k))!,SP(k)
      AsymFact = AsymFact + &
        & (SP(k+1)*cos(ANG(k+1))+ SP(k)*cos(ANG(k)))* &
        & (cos(ANG(k+1))-cos(ANG(k)))
    end do
    AsymFact = -AsymFact * 0.5_8 * 0.5_8
    !print *, 'AsymFact', AsymFact

  end subroutine miescat

  subroutine clscat6(                              &
    & wl, RM, RS, NANG, NDCR, EPS, Rmed, Rre, Rim, & ! (in )
    & SE, SA, cSA, mu, bs, ANG, SP                 ) ! (out)
!    & SE, SA, cSA, mu(1:NDCR), bs(1:NDCR),         & ! (out)
!    & ANG(1:2*NANG-1),SP(1:2*NANG-1)               ) ! (out)

!      Program ClScat
      Implicit NONE

! George L. Hashimoto

! 040729    Ver. 6
!   accuracy is set by input file
!   input and output is modified
!     input:  data.cl
!     output: out.cl

! 040712    Ver. 5
!   integration region is automatically determined
!   [from -5*sigma to +5*sigma]

! 040622    Ver. 4
!   integration region is automatically determined

! 040622    Ver. 3
!   data output is modified

! 020206    Ver. 2

! 020128
!   data length output

! 020124

! Calculate scattering propety of clouds.  
!   input (fort.32)
!       wavelength
!       refractive index of surrounding medium
!       refractive index of cloud particle
!       size distribution parameter (log-normal)
!   output (fort.52)
!       extinction cross section
!       single scattering albedo
!       single scattering co-albedo
!       direction change ratio
!       phase function

! -----constant-----
      Include 'clscat6.con'

! -----variables-----
! number of angles between 0 and 90 degrees, matix elements
! calculated at 2*NANG-1 angles including 0, 90, and 180 degrees.  
      Integer NANG
! number of calculation for change direction ratio
      Integer NDCR
! accuracy
      Real*8  EPS
! (real) refractive index of surrounding medium
      Real*8  Rmed
! refractive index of sphere = Rre + i*Rim
      Real*8  Rre, Rim
! wavelength
      Real*8  wl

! size distribution parameter
! [ log normal: v=exp((ln(R/RM)/ln(RS))**2 / 2) ]
      Real*8  RM, RS

! lower and upper bound of size distribution
      Real*8  A, B

! extinction cross section
      Real*8  SE
! single scattering albedo
      Real*8  SA
! single scattering co-albedo
      Real*8  cSA
! direction change ratio
! (direction(cosine of zenith angle), ratio)
      Real*8  mu(NDCRMAX), bs(NDCRMAX)
! phase function
      Real*8  ANG(NXMAX), SP(NXMAX)

! local
      Integer IO
      Real*8  PI, DANG
      Integer NAN
      Real*8  gz0, gz, pfb1, pfb2

! loop control
      Integer i, j

! -----begin-----
! set pi
      PI = acos(-1.0D0)

! read data
!      IO = 33
!      Open(IO,FILE='data.cl',STATUS='OLD',ERR=910)
!      Read(IO,*)
!      Read(IO,*) wl
!      Read(IO,*)
!      Read(IO,*) 
!      Read(IO,*)
!      Read(IO,*) RM, RS
!      Read(IO,*)
!      Read(IO,*)
!      Read(IO,*) NANG
!      Read(IO,*)
!      Read(IO,*)
!      Read(IO,*) NDCR
!      Read(IO,*)
!      Read(IO,*)
!      Read(IO,*) EPS
!      Read(IO,*)
!      Read(IO,*)
!      Read(IO,*) Rmed
!      Read(IO,*)
!      Read(IO,*)
!      Read(IO,*) Rre, Rim
!      Close(IO)
!   check data
      If (wl.le.(0.0D0)) then
          Write(*,*) 'invalid wavelength', wl
          Stop
      End If
      If (RM.le.(0.0D0)) then
          Write(*,*) 'invalid size parameter (RM)', RM
          Stop
      End If
      If (RS.le.(1.0D0)) then
          Write(*,*) 'invalid size parameter (RS)', RS
          Stop
      End If
      If (NANG.le.2) then
          Write(*,*) 'invalid (NANG)', NANG
          Stop
      End If
      NAN = 2*NANG-1
      If (NAN.gt.NXMAX) then
          Write(*,*) 'invalid (NANG)', NANG, NXMAX
          Stop
      End If
      If (NDCR.le.2) then
          Write(*,*) 'invalid (NDCR)', NDCR
          Stop
      End If
      If (NDCR.gt.NDCRMAX) then
          Write(*,*) 'invalid (NDCR)', NDCR, NDCRMAX
          Stop
      End If
      If (EPS.le.(0.0D0)) then
          Write(*,*) 'invalid (EPS)', EPS
          Stop
      End If
      If (Rmed.lt.(1.0D0)) then
          Write(*,*) 'invalid (Rmed)', Rmed
          Stop
      End If
      If (Rre.lt.(1.0D0)) then
          Write(*,*) 'invalid (Rre)', Rre
          Stop
      End If
      If (Rim.le.(0.0D0)) then
          Write(*,*) 'invalid (Rim)', Rim
          Stop
      End If


! set integration region
      A = RM * (RS**(-ta))
      B = RM * (RS**(+ta))
      Write(*,*) A, B

! calculation
      Call QROMB(A,B,NANG,Rmed,Rre,Rim,wl,RM,RS,SE,SA,ANG,SP,EPS)
      cSA = 1.0D0 - SA
!   direction change ratio
      NAN = 2*NANG-1
      DANG=PI/2.0D0/(NANG-1)
      Do 120 i = 1, NDCR
        mu(i) = 1.0D0*(i-1)/(NDCR-1)
        gz0   = mu(i)
        bs(i) = 0.0D0
        gz = cos(ANG(NANG))
        Call Pfbar(gz0,gz,pfb1,NAN,ANG,SP)
        Do 110 j = NANG+1, NAN
          gz = cos(ANG(j))
          Call Pfbar(gz0,gz,pfb2,NAN,ANG,SP)
          bs(i) = bs(i) &
     & + ( ((pfb1*sin(ANG(j-1)))+(pfb2*sin(ANG(j))))/2.0D0 &
     &     *(2.0D0*PI)*DANG )
          pfb1 = pfb2
 110    Continue
 120  Continue

! output
  if(.false.) then
      IO = 52
      Open(IO,FILE='out.cl',STATUS='UNKNOWN')
      Write(IO,*) 'output of clscat6.f'
      Write(IO,*)
      Write(IO,*) '--input parameters--'
      Write(IO,*) 'wavelength [m]'
      Write(IO,*) wl
      Write(IO,*) 'size distribution parameter (RM RS)'
      Write(IO,*) '[log normal: v=exp((-ln(R/RM)/ln(RS))**2 /2)]'
      Write(IO,*) RM, RS
      Write(IO,*) 'lower and upper bound of size distribution [m]'
      Write(IO,*) A, B
      Write(IO,*) 'refractive index of surrounding medium (real)'
      Write(IO,*) Rmed
      Write(IO,*) 'refractive index of particle (Rre + i*Rim)'
      Write(IO,*) Rre, Rim
      Write(IO,*)
      Write(IO,*) 'accuracy'
      Write(IO,*) EPS
      Write(IO,*)
      Write(IO,*) '--result--'
      Write(IO,*) 'extinction cross section'
      Write(IO,*) SE
      Write(IO,*) 'single scattering albedo'
      Write(IO,*) SA
      Write(IO,*) 'single scattering co-albedo'
      Write(IO,*) cSA
      Write(IO,*)
      Write(IO,*) 'data length of direction change ratio'
      Write(IO,*) NDCR
      Write(IO,*) 'direction change ratio'
      Write(IO,*) '(cosine of zenith angle, ratio)'
      Do 130 i = 1, NDCR
        Write(IO,*) mu(i), bs(i)
 130  Continue
      Write(IO,*)
      Write(IO,*) 'data length of phase function'
      Write(IO,*) NAN
      Write(IO,*) 'phase function (angle[deg], phase function)'
      Do 140 i = 1, NAN
        ANG(i) = ANG(i)*(180.0D0/PI)
        Write(IO,*) ANG(i), SP(i)
 140  Continue
      Close(IO)


      Stop
  end if

 if(.false.) then
 910  Continue
      Write(*,*) 'cannot open file (data)'
      Stop
  end if

  End subroutine



      SUBROUTINE QROMB(A,B,NANG,Rmed,Rre,Rim,wl,RM,RS,SE,SA,ANG,SP,EPS)
      Implicit NONE

!        Return as S the integral of the function FUNC from A to B.
!        Integration is performed by Romberg's method of order 2K,
!        where, e.g., K=2 is Simpson's rule.
!        Here EPS is the fractional accuracy desired, as determined
!        by the extrapolation error estimate; JMAX limits the total
!        number of steps; K is the number of steps; K is the number
!        of points used in the extrapolation.

!       2001.01.31      Ver. 1

! -----CONSTANTS-----
      Include 'clscat6.con'

! number of points used in the extrapolation
      Integer K0
      Parameter (K0=5)

! -----INPUT-----
! start and end point of integration
      Real*8  A, B

      Integer NANG
      Real*8  Rmed, Rre, Rim, wl, RM, RS
      Real*8  EPS

! -----OUTPUT-----
      Real*8  SE, SA, ANG(NXMAX), SP(NXMAX)

! -----LOCAL-----
      Integer JMAX0
      Parameter (JMAX0=20)
      Integer JMAXP
      Parameter (JMAXP=JMAX0+1)
      Integer KM
      Parameter (KM=K0-1)

      Real*8  S(JMAXP), H(JMAXP), DSS, SS, SD
      Integer NP, J, k

! -----BEGIN-----
        If (A.eq.B) then
            Write(*,*) 'A=B !', A, B
            Stop
        End If

        H(1)=1.0D0
        DO 11 J=1,JMAX0
          Call TRAPZD(A,B,J,SD,SE,SA,ANG,SP,      &
     &                NANG,Rmed,Rre,Rim,wl,RM,RS)
! totyu-keika (1)
          Write(*,*) J, SE
! 
          S(J)=SE
          IF (J.GE.K0) THEN
              NP=J-KM
              CALL POLINT(H,S,NP,K0,0.0D0,SS,DSS)
! totyu-keika (2)
              Write(*,*) SS, DSS
              IF (ABS(DSS).LT.EPS*ABS(SS)) then
                  SA=SA/SE
                  Do 210 k=1,2*NANG-1
                    SP(k)=SP(k)/SE
 210              Continue
                  SE=SE/SD
                  RETURN
              End If
          END IF
          H(J+1)=0.25D0*H(J)
 11     CONTINUE

!x        Write(*,*) 'QROMB  Too many steps.'
!x        Stop

        SA=SA/SE
        Do 310 k=1,2*NANG-1
          SP(k)=SP(k)/SE
 310    Continue
        SE=SE/SD
        RETURN

      END





      SUBROUTINE TRAPZD(A,B,N,SD,SE,SA,ANG,SP,      &
     &                  NANG,Rmed,Rre,Rim,wl,RM,RS)
      Implicit NONE

!        This routine computes the N'th stage of refinement of an
!        extended trapezoidal rule.
!        The function to be integrated between limits A and B, also
!        input. When called with N=1, the routine returns as S the
!        crudest estimate of \int_a^b f(x) dx. Subsequent calls
!        with N=2,3,... (in that sequential order) will improve the
!        accuracy of S by adding 2^(N-2) additional interior points.
!        S should not be modified between sequential calls.

!       2001.01.31      Ver. 1

! -----constant-----
      Include 'clscat6.con'

! -----INPUT-----
      Real*8  A, B
      Integer N

      Integer NANG
      Real*8  Rmed, Rre, Rim, wl, RM, RS

! -----OUTPUT-----
      Real*8  SD, SE, SA, ANG(NXMAX), SP(NXMAX)

! -----LOCAL-----
      Real*8  DEL, X
      Integer TNM, J

      Real*8  PI
      Real*8  WA, QextA, ssaA, pfA(NXMAX), SWA
      Real*8  WB, QextB, ssaB, pfB(NXMAX), SWB
      Real*8  WX, QextX, ssaX, pfX(NXMAX), SWX
      Real*8  SUMd, SUMe, SUMa, SUMp(NXMAX)
      Integer k

! -----BEGIN-----
        PI=acos(-1.0D0)

        IF (N.EQ.1) THEN
            Call FUNC(A,QextA,ssaA,ANG,pfA,NANG,Rmed,Rre,Rim,wl)
            Call Sdist(A,WA,RM,RS)
            SWA=QextA*(PI*A*A)*WA
            Call FUNC(B,QextB,ssaB,ANG,pfB,NANG,Rmed,Rre,Rim,wl)
            Call Sdist(B,WB,RM,RS)
            SWB=QextB*(PI*B*B)*WB
            SD=0.5D0*(WA+WB)
            SE=0.5D0*(SWA+SWB)
            SA=0.5D0*((ssaA*SWA)+(ssaB*SWB))
            Do 210 k = 1, 2*NANG-1
              SP(k)=0.5D0*((pfA(k)*SWA)+(pfB(k)*SWB))
 210        Continue
          ELSE
            TNM=2**(N-2)
            DEL=(B-A)/TNM
            X=A+(0.5D0*DEL)
            SUMd=0.0D0
            SUMe=0.0D0
            SUMa=0.0D0
            Do 220 k = 1, 2*NANG-1
              SUMp(k)=0.0D0
 220        Continue
            DO 11 J=1,TNM
              Call FUNC(X,QextX,ssaX,ANG,pfX,NANG,Rmed,Rre,Rim,wl)
              Call Sdist(X,WX,RM,RS)
              SWX=QextX*(PI*X*X)*WX
              SUMd=SUMd+WX
              SUMe=SUMe+SWX
              SUMa=SUMa+(ssaX*SWX)
              Do 230 k = 1, 2*NANG-1
                SUMp(k)=SUMp(k)+(pfX(k)*SWX)
 230          Continue
              X=X+DEL
 11         CONTINUE
            SD=0.5D0*(SD+SUMd/TNM)
            SE=0.5D0*(SE+SUMe/TNM)
            SA=0.5D0*(SA+SUMa/TNM)
            Do 240 k = 1, 2*NANG-1
              SP(k)=0.5D0*(SP(k)+SUMp(k)/TNM)
 240        Continue
        ENDIF

      RETURN
      END





      SUBROUTINE POLINT(XA,YA,NP,N,X,Y,DY)
      Implicit NONE

!        Given arrays XA and YA, each of length N, and given a value
!        X, this routine returns a value Y, and an error estimate DY.
!        If P(x) is the polynomial of degree N-1 such that
!        P(XAj)=YAj, i=1,...,N, then the returned value Y=P(X).

!       1997.11.04      Ver. 1

! -----CONSTANTS-----
      Integer JMAX0
      Parameter (JMAX0=20)
      Integer JMAXP
      Parameter (JMAXP=JMAX0+1)

! -----INPUT-----
      Real*8  XA(JMAXP), YA(JMAXP)
      Integer NP, N
      Real*8  X

! -----OUTPUT-----
      Real*8  Y, DY

! -----LOCAL-----
      Real*8  C(JMAXP), D(JMAXP)
      Integer NS
      Real*8  DIF, DIFT
      Real*8  HO, HP
      Real*8  W
      Real*8  DEN

      Integer I, M

! -----BEGIN-----
        NS=1
        DIF=ABS(X-XA(NP))
        DO 11 I=1,N
          DIFT=ABS(X-XA(NP-1+I))
          IF (DIFT.LT.DIF) THEN
              NS=I
              DIF=DIFT
          END IF
          C(I)=YA(NP-1+I)
          D(I)=YA(NP-1+I)
 11     CONTINUE
        Y=YA(NP-1+NS)
        NS=NS-1
        DO 13 M=1,N-1
          DO 12 I=1,N-M
            HO=XA(NP-1+I)-X
            HP=XA(NP-1+I+M)-X
            W=C(I+1)-D(I)
            DEN=HO-HP
            IF (DEN.EQ.(0.0D0)) then
                Write(*,*) 'POLINT  DEN=0.0'
                Stop
            End If
            DEN=W/DEN
            D(I)=HP*DEN
            C(I)=HO*DEN
 12       CONTINUE
          IF (2*NS.LT.N-M) THEN
              DY=C(NS+1)
            ELSE
              DY=D(NS)
              NS=NS-1
          END IF
          Y=Y+DY
 13     CONTINUE


      RETURN
      END





      Subroutine Sdist(X,FX,RM,RS)
      Implicit NONE

! -----input-----
      Real*8  X
      Real*8  RM, RS

! -----output-----
      Real*8  FX

! -----begin-----
        FX = exp(-((log(X/RM)/log(RS))**2)/2.0D0)

      Return
      End





      Subroutine FUNC(X,Qext,ssa,ANG,pf,NANG,Rmed,Rre,Rim,wl)
      Implicit NONE

! -----constant-----
      Include 'clscat6.con'

! -----INPUT-----
! number of angles between 0 and 90 degrees, matix elements
! calculated at 2*NANG-1 angles including 0, 90, and 180 degrees.  
      Integer NANG
! radius
      Real*8  X
! (real) refractive index of surrounding medium
      Real*8  Rmed
! refractive index of sphere = Rre + i*Rim
      Real*8  Rre, Rim
! wavelength
      Real*8  wl

! -----OUTPUT-----
! extinction efficiency
      Real*8  Qext
! single scattering albedo
      Real*8  ssa
! phase function
      Real*8  ANG(NXMAX), pf(NXMAX)

! -----LOCAL-----
      Real*8  PI
      Integer NAN, i
      Real*8  sp
      Real*8  REFRELr, S1r(200), S2r(200)
      Real*8  REFRELi, S1i(200), S2i(200)
      Real*8  Qsca, S11(200), DANG, xxx
      Real*8  sum, gz0, gz, pfb1, pfb2

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

! set ratio of refractive index
        REFRELr = Rre/Rmed
        REFRELi = Rim/Rmed
! set size parameter
        sp=2.0D0*PI*X*Rmed/wl
! calculate Mie scattering
        Call BHMIE2(sp,REFRELr,REFRELi,NANG,S1r,S1i,S2r,S2i, &
     &              Qext,Qsca,xxx)
! calculate single scattering albedo
        ssa = QSCA / Qext
! calculate phase function
        NAN=2*NANG-1
        DANG=PI/2.0D0/(NANG-1)
!   set angle
        Do 110 i = 1, NAN
          ANG(i) = DANG*(i-1)
 110    Continue
!   phase function
        Do 120 i = 1, NAN
          S11(i) = ( (S1r(i)*S1r(i)+S1i(i)*S1i(i))           &
     &              +(S2r(i)*S2r(i)+S2i(i)*S2i(i)) ) / 2.0D0
 120  Continue
!  normalization
! (rotate 90 degree to include 0 and 180 degree scattering)
        gz0 = 0.0D0
        sum = 0.0D0
        gz  = cos(ANG(1))
        Call Pfbar(gz0,gz,pfb1,NAN,ANG,S11)
        Do 130 i = 2, NAN
          gz = cos(ANG(i))
          Call Pfbar(gz0,gz,pfb2,NAN,ANG,S11)
          sum = sum                                         &
     & + ( ((pfb1*sin(ANG(i-1)))+(pfb2*sin(ANG(i))))/2.0D0  &
     &     *(2.0D0*PI)*DANG )
          pfb1 = pfb2
 130    Continue
        Do 140 i = 1, NAN
          pf(i) = S11(i)/sum
 140    Continue

      Return
      End





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

! 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  X
      Integer NANG
      Real*8  REFRELr, REFRELi

! -----OUTPUT-----
      Real*8  S1r(200), S1i(200), S2r(200), S2i(200)
      Real*8  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)
      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

! -----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 555 J=1,NANG
          THETA(J)=(FLOAT(J)-1.0D0)*DANG
          AMU(J)=COS(THETA(J))
 555    CONTINUE
! logarithmic derivative D(J) calculated by downward recurrence 
! 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
! 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
        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





      Subroutine Pfbar(gz0,gz,pfb,NAN,ANG,SP)
      Implicit NONE

! -----constant-----
      Include 'clscat6.con'

! -----input-----
      Real*8  gz0, gz
      Integer NAN
      Real*8  ANG(NXMAX), SP(NXMAX)

! -----output-----
      Real*8  pfb

! -----local-----
      Real*8  PI
      Real*8  DANG, PHI, theta, pfp
      Integer i

! -----begin-----
! set pi
        PI = acos(-1.0D0)

! set angle
        DANG = PI/(NAN-1)

! calculation
        pfb = 0.0D0
        Do 110 i = 1, (NAN-1)*2
          PHI = DANG*i
          theta = acos( &
     & gz0*gz+(((1.0D0-gz0**2)*(1.0D0-gz**2))**(0.5D0))*cos(PHI) )
          Call FindPF(theta,pfp,NAN,ANG,SP)
          pfb = pfb + pfp
 110    Continue
        pfb = pfb/((NAN-1)*2)

      Return
      End





      Subroutine FindPF(theta,pfp,NAN,ANG,SP)
      Implicit NONE

! -----constant-----
      Include 'clscat6.con'

! -----input-----
      Real*8  theta
      Integer NAN
      Real*8  ANG(NXMAX), SP(NXMAX)

! -----output-----
      Real*8  pfp

! -----local-----
      Integer i

! -----begin-----
        i = 1
        If (theta.lt.(ANG(i))) then
            Go to 120
        End If
        Do 110 i = 2, NAN
          If (theta.le.(ANG(i))) then
              pfp = ( (SP(i-1)*(ANG(i)-theta))    &
     &               +(SP(i)*(theta-ANG(i-1))) )  &
     &                  / (ANG(i)-ANG(i-1))
              Go to 210
          End If
 110    Continue
 120    Continue
        Write(*,*) 'Cannot find phase function'
        Write(*,*) theta, ANG(1), ANG(NAN)
        Stop
 210    Continue

      Return
      End

end module
