

module mtrunc

  implicit none
  
  INTEGER, PUBLIC, PARAMETER        :: dp  = kind(1.0d0)
  
contains
  
  
  SUBROUTINE trunc_dm(nmat, ngauss, u, wth, F, nbetal, betal, coeftr_out)

    implicit none

    !--- Input Arguments
    integer, intent(in) :: nmat
    integer, intent(in) :: ngauss
    integer, intent(in)  :: nbetal
    real(kind=dp), dimension(ngauss), intent(in) :: u
    real(kind=dp), dimension(ngauss), intent(in) :: wth
    real(kind=dp), dimension(6,ngauss), intent(in) :: F

    !--- Output Arguments
    real(kind=dp), dimension(6,0:nbetal), intent(inout) :: betal
    real(kind=dp), dimension(1),          intent(inout) :: coeftr_out

    !--- Local Arguments
    real(kind=dp)  :: coeftr
    integer :: k, i, flag
    real(kind=dp), dimension(4,0:nbetal,ngauss)   :: apl
    real(kind=dp), dimension(ngauss)              :: F_trunc_recomp
    real(kind=dp), dimension(6,ngauss)            :: F_trunc
    real(kind=dp) :: renorm

    !**************************************************************************
    !                              DELTA-M
    !**************************************************************************
    !	  Compute the moments Bl (legendre) for F11 only
    ! *********************************************************************

    !print*, "(trunc_dm) coeftr =", coeftr

    flag = 0
    CALL comp_Bl(flag, ngauss, u, F, wth, nbetal, apl, betal)

    renorm     = betal(1,0)
    betal(1,:) = betal(1,:) / renorm 

    ! apply D-M
    coeftr     = betal(1,nbetal)
    !print*, "(trunc_dm) coeftr =", coeftr

    betal(1,:) = (betal(1,:) - coeftr) / (1.D0 - coeftr)

    !--- Now recompute F11 from the betals to recompute Fij accordingly
    F_trunc_recomp(:) = 0.0D0
    DO k = 0, nbetal
       DO i = 1, ngauss
          F_trunc_recomp(i) = F_trunc_recomp(i)  + DBLE(2*k+1)*betal(1,k)*apl(1,k,i) ! F11
       ENDDO
    ENDDO


    if (nmat.gt.1) then

       F_trunc(:,:) = 0.0D0
       !--- rescale Fij, Fii
       F_trunc(1,1:ngauss) = F_trunc_recomp(1:ngauss)
       F_trunc(2,1:ngauss) = F(2,:) * F_trunc_recomp(1:ngauss) / F(1,:)
       F_trunc(3,1:ngauss) = F(3,:) * F_trunc_recomp(1:ngauss) / F(1,:)
       F_trunc(4,1:ngauss) = F(4,:) * F_trunc_recomp(1:ngauss) / F(1,:) 
       !F_trunc(5,1:ngauss) = F(5,:) * F_trunc_recomp(1:ngauss) / F(1,:) 
       !F_trunc(6,1:ngauss) = F(6,:) * F_trunc_recomp(1:ngauss) / F(1,:) 
       F_trunc(5,1:ngauss) = F(5,:) / (1.0D0 - coeftr)
       F_trunc(6,1:ngauss) = F(6,:) / (1.0D0 - coeftr)

       flag = 1
       CALL comp_Bl(flag, ngauss, u, F_trunc, wth, nbetal, apl, betal)

    end if

    ! *********************************************************************
    !          Convert beta moment in legendre coefficient 
    !          (multiply by factor 2k+1)
    ! *********************************************************************
    DO k = 0,nbetal
       betal(:,k) = DBLE(2*k+1)*betal(:,k)
    ENDDO

    !print*, "(trunc_dm) coeftr =", coeftr
    coeftr_out(1) = coeftr

  END SUBROUTINE trunc_dm








  
  SUBROUTINE comp_Bl(flag, nangle, u, F, w, nbetal, apl, betal)

    !  Calculate the expansion coefficients of the scattering matrix in    
    !  generalized spherical functions by numerical integration over the   
    !  scattering angle.                                                   
    !  From the adding paper 
    !  de Haan, Bosma, Hovenier, 1987, A&A, 183, 371     

    implicit none

    INTEGER, INTENT(IN) :: flag
    INTEGER, INTENT(IN) :: nbetal
    INTEGER, INTENT(IN) :: nangle
    REAL(kind=dp), INTENT(IN) :: u(nangle)
    REAL(kind=dp), INTENT(IN) :: w(nangle)
    REAL(kind=dp), INTENT(IN) :: F(6,nangle)
    !                            F(1) = F11
    !                            F(2) = F22
    !                            F(3) = F33
    !                            F(4) = F44
    !                            F(5) = F12
    !                            F(6) = F34

    REAL(kind=dp), INTENT(INOUT) :: apl(4,0:nbetal,nangle) 
    REAL(kind=dp), INTENT(OUT)   :: betal(6,0:nbetal)
    !                             betal(1) = alpha1
    !                             betal(2) = alpha2
    !                             betal(3) = alpha3
    !                             betal(4) = alpha4
    !                             betal(5) = beta1
    !                             betal(6) = beta2

    ! local variables
    REAL(kind=dp) :: P00(nangle,2) 
    REAL(kind=dp) :: P22(nangle,2)
    REAL(kind=dp) :: P2m2(nangle,2)
    REAL(kind=dp) :: P02(nangle,2)
    REAL(kind=dp) :: qroot6
    REAL(kind=dp) :: fac1
    REAL(kind=dp) :: fac2
    REAL(kind=dp) :: fac3
    REAL(kind=dp) :: sql41
    REAL(kind=dp) :: sql4
    REAL(kind=dp) :: twol1
    REAL(kind=dp) :: tmp2
    REAL(kind=dp) :: tmp1
    REAL(kind=dp) :: denom
    REAL(kind=dp) :: alfap
    REAL(kind=dp) :: alfam

    INTEGER :: i, l
    INTEGER :: lnew
    INTEGER :: lold
    INTEGER :: itmp

    ! ----------------
    ! compute Legendre Polynomials

    if (flag .eq. 0) then

       apl(:,:,:) = -32768.0D0

       qroot6     = -0.25D0*dsqrt(6.D0)

       !  Start loop over the coefficient index l                             *
       !  first update generalized spherical functions, then calculate betal. *
       !  lold and lnew are pointer-like indices used in recurrence           *
       lnew = 1
       lold = 2

       do l= 0, nbetal

          if (l .eq. 0) then
             ! Adding paper Eq. (77) with m=0                           
             do i=1, nangle
                P00(i,lold) = 1.D0
                P00(i,lnew) = 0.D0
                P02(i,lold) = 0.D0
                P22(i,lold) = 0.D0
                P2m2(i,lold)= 0.D0
                P02(i,lnew) = 0.D0
                P22(i,lnew) = 0.D0
                P2m2(i,lnew)= 0.D0
             enddo
          else
             fac1 = (2.D0*l-1.D0)/dble(l)
             fac2 = dble(l-1.D0)/dble(l)
             ! Adding paper Eq. (81) with m=0                           *
             do i=1, nangle
                P00(i,lold) = fac1*u(i)*P00(i,lnew) - fac2*P00(i,lold)
             enddo
          endif

          if (l .eq. 2) then
             ! Adding paper Eqs. (78) and (80)  with m=2                *
             ! sql4 contains the factor dsqrt(l*l-4) needed in          *
             ! the recurrence Eqs. (81) and (82)                        *
             do i=1, nangle
                P02(i,lold) = qroot6*(1.D0-u(i)*u(i))
                P22(i,lold) = 0.25D0*(1.D0+u(i))*(1.D0+u(i))
                P2m2(i,lold)= 0.25D0*(1.D0-u(i))*(1.D0-u(i))
                P02(i,lnew) = 0.D0
                P22(i,lnew) = 0.D0
                P2m2(i,lnew)= 0.D0
             enddo
             sql41 = 0.D0
          endif

          if (l .gt. 2) then
             ! Adding paper Eq. (82) with m=0 and m=2                   *
             sql4  = sql41
             sql41 = dsqrt(dble(l*l)-4.D0)
             twol1 = 2.D0*dble(l)-1.D0
             tmp1  = twol1/sql41
             tmp2  = sql4/sql41
             denom = (dble(l)-1.D0)*(dble(l*l)-4.D0)
             fac1  = twol1*(dble(l)-1.D0)*dble(l)/denom
             fac2  = 4.D0*twol1/denom
             fac3  = dble(l)*((dble(l)-1.D0)*(dble(l)-1.D0)-4.D0)/denom
             do i=1, nangle
                P02(i,lold) = tmp1*u(i)*P02(i,lnew) - tmp2*P02(i,lold)
                P22(i,lold) = (fac1*u(i)-fac2)*P22(i,lnew) - fac3*P22(i,lold)
                P2m2(i,lold)= (fac1*u(i)+fac2)*P2m2(i,lnew) - fac3*P2m2(i,lold)
             enddo
          endif

          !  Switch indices so that lnew indicates the function with      *
          !  the present index value l, this mechanism prevents swapping  *
          !  of entire arrays.                                            *
          itmp = lnew
          lnew = lold
          lold = itmp

          apl(1,l,:) = P00(:,lnew)
          apl(2,l,:) = P02(:,lnew)
          apl(3,l,:) = P22(:,lnew)
          apl(4,l,:) = P2m2(:,lnew)

       ENDDO

    endif

    !--------------
    ! compute Betal

    if (flag .eq. 0) then

       betal(:,:)    = 0.0D0

       do l= 0, nbetal
          do i=1, nangle
             betal(1,l) = betal(1,l) + (apl(1,l,i) * w(i) * F(1,i) * 0.5D0) 
          enddo
       enddo ! End of loop over index l                                        

    else

       do l= 0, nbetal
          !  Now calculate the coefficients by integration over angle     *
          !  See de Haan et al. (1987) Eqs. (68)-(73).                    
          alfap = 0.0D0
          alfam = 0.0D0
          do i=1, nangle
             alfap = alfap           + apl(3,l,i) * w(i) * (F(2,i)+F(3,i)) * 0.5D0
             alfam = alfam           + apl(4,l,i) * w(i) * (F(2,i)-F(3,i)) * 0.5D0
             betal(4,l) = betal(4,l) + apl(1,l,i) * w(i) * F(4,i) *0.5D0
             betal(5,l) = betal(5,l) + apl(2,l,i) * w(i) * F(5,i) *0.5D0
             betal(6,l) = betal(6,l) + apl(2,l,i) * w(i) * F(6,i) *0.5D0
          enddo
          betal(2,l) = ( alfap + alfam ) / 2.0D0
          betal(3,l) = ( alfap - alfam ) / 2.0D0
       enddo ! End of loop over index l                                        

    endif

  end SUBROUTINE comp_Bl


  FUNCTION XINTEG2(imin, imax, n, xin, yin)
    ! computes integral of yin, variable step 
    ! make sure you have at least two points of integration
    ! Trapeze method

    IMPLICIT NONE

    REAL (KIND=dp) :: XINTEG2

    INTEGER, INTENT (IN)        :: imin
    INTEGER, INTENT (IN)        :: imax
    INTEGER, INTENT (IN)        :: n
    REAL (KIND=dp), INTENT (IN) :: xin(n)
    REAL (KIND=dp), INTENT (IN) :: yin(n)
    REAL (KIND=dp)              :: xa, ya, xb, yb
    REAL (KIND=dp)              :: primitaux
    INTEGER                     :: i

    xa = xin(imin)
    ya = yin(imin)
    primitaux = 0.0_dp
    DO i=imin+1,imax
       xb = xin(i)
       yb = yin(i)
       primitaux = primitaux + (xb - xa) * (ya + yb)
       xa = xb
       ya = yb
    ENDDO
    XINTEG2 = primitaux * 0.5_dp

  END FUNCTION XINTEG2

  
  subroutine get_f_recomp(nmat, nbetal, nangle, u, betal, F_recomp)

    !  Calculate the recomposed phase matrix from 
    !  Legendre expension coefficients
    !  From the adding paper 
    !  de Haan, Bosma, Hovenier, 1987, A&A, 183, 371     

    implicit none

    INTEGER, INTENT(IN) :: nmat
    INTEGER, INTENT(IN) :: nbetal
    INTEGER, INTENT(IN) :: nangle
    REAL(KIND=dp), INTENT(IN) :: u(nangle)
    REAL(KIND=dp), INTENT(IN) :: betal(6,0:nbetal)
    !                             betal(1) = alpha1
    !                             betal(2) = alpha2
    !                             betal(3) = alpha3
    !                             betal(4) = alpha4
    !                             betal(5) = beta1
    !                             betal(6) = beta2
    REAL(KIND=dp), INTENT(inout) :: F_recomp(6,nangle)

    ! local variables
    REAL(KIND=dp) :: F22p33(nangle)
    REAL(KIND=dp) :: F22m33(nangle)
    REAL(KIND=dp) :: P00(nangle,2) 
    REAL(KIND=dp) :: P22(nangle,2)
    REAL(KIND=dp) :: P2m2(nangle,2)
    REAL(KIND=dp) :: P02(nangle,2)
    REAL(KIND=dp) :: qroot6
    REAL(KIND=dp) :: fac1
    REAL(KIND=dp) :: fac2
    REAL(KIND=dp) :: fac3
    REAL(KIND=dp) :: sql41
    REAL(KIND=dp) :: sql4
    REAL(KIND=dp) :: twol1
    REAL(KIND=dp) :: tmp2
    REAL(KIND=dp) :: tmp1
    REAL(KIND=dp) :: denom

    INTEGER :: i, l
    INTEGER :: lnew
    INTEGER :: lold
    INTEGER :: itmp

    ! ====================
    !  Initialization    

    qroot6        = -0.25_DP*dsqrt(6._DP)
    F_recomp(:,:) = 0.0_DP
    F22p33(:)     = 0.0_DP
    F22m33(:)     = 0.0_DP

    !  Start loop over the coefficient index l                             *
    !  first update generalized spherical functions, then calculate betal. *
    !  lold and lnew are pointer-like indices used in recurrence           *
    lnew = 1
    lold = 2

    do l= 0, nbetal

       if (l .eq. 0) then
          ! Adding paper Eq. (77) with m=0                           
          do i=1, nangle
             P00(i,lold) = 1._DP
             P00(i,lnew) = 0._DP
             P02(i,lold) = 0._DP
             P22(i,lold) = 0._DP
             P2m2(i,lold)= 0._DP
             P02(i,lnew) = 0._DP
             P22(i,lnew) = 0._DP
             P2m2(i,lnew)= 0._DP
          enddo
       else
          fac1 = (2._DP*l-1._dp)/dble(l)
          fac2 = dble(l-1._dp)/dble(l)
          ! Adding paper Eq. (81) with m=0                           *
          do i=1, nangle
             P00(i,lold) = fac1*u(i)*P00(i,lnew) - fac2*P00(i,lold)
          enddo
       endif

       if (nmat.gt.1) then

          if (l .eq. 2) then
             ! Adding paper Eqs. (78) and (80)  with m=2                *
             ! sql4 contains the factor dsqrt(l*l-4) needed in          *
             ! the recurrence Eqs. (81) and (82)                        *
             do i=1, nangle
                P02(i,lold) = qroot6*(1._DP-u(i)*u(i))
                P22(i,lold) = 0.25_DP*(1._DP+u(i))*(1._DP+u(i))
                P2m2(i,lold)= 0.25_DP*(1._DP-u(i))*(1._DP-u(i))
                P02(i,lnew) = 0._DP
                P22(i,lnew) = 0._DP
                P2m2(i,lnew)= 0._DP
             enddo
             sql41 = 0._DP
          endif

          if (l .gt. 2) then
             ! Adding paper Eq. (82) with m=0 and m=2                   *
             sql4  = sql41
             sql41 = dsqrt(dble(l*l)-4._dp)
             twol1 = 2._DP*dble(l)-1._dp
             tmp1  = twol1/sql41
             tmp2  = sql4/sql41
             denom = (dble(l)-1._dp)*(dble(l*l)-4._dp)
             fac1  = twol1*(dble(l)-1._dp)*dble(l)/denom
             fac2  = 4._DP*twol1/denom
             fac3  = dble(l)*((dble(l)-1._dp)*(dble(l)-1._dp)-4._dp)/denom
             do i=1, nangle
                P02(i,lold) = tmp1*u(i)*P02(i,lnew) - tmp2*P02(i,lold)
                P22(i,lold) = (fac1*u(i)-fac2)*P22(i,lnew) - fac3*P22(i,lold)
                P2m2(i,lold)= (fac1*u(i)+fac2)*P2m2(i,lnew) - fac3*P2m2(i,lold)
             enddo
          endif

       end if
       
       !  Switch indices so that lnew indicates the function with      *
       !  the present index value l, this mechanism prevents swapping  *
       !  of entire arrays.                                            *
       itmp = lnew
       lnew = lold
       lold = itmp

       ! compute truncated phase matrix
       do i = 1, nangle
          F_recomp(1,i) = F_recomp(1,i) + betal(1,l)                * P00(i,lnew)
          if (nmat.gt.1) then
             F22p33(i)     = F22p33(i)     + (betal(2,l)+betal(3,l))   * P22(i,lnew)
             F22m33(i)     = F22m33(i)     + (betal(2,l)-betal(3,l))   * P2m2(i,lnew)
             F_recomp(4,i) = F_recomp(4,i) + betal(4,l)                * P00(i,lnew)
             F_recomp(5,i) = F_recomp(5,i) + betal(5,l)                * P02(i,lnew)
             F_recomp(6,i) = F_recomp(6,i) + betal(6,l)                * P02(i,lnew)
          end if
       enddo

    enddo ! End of loop over index l                                        
    
    if (nmat.gt.1) then
       F_recomp(2,:) = ( F22p33(:) + F22m33(:) ) / 2.0
       F_recomp(3,:) = ( F22p33(:) - F22m33(:) ) / 2.0
    end if
    
  end subroutine get_f_recomp








end module mtrunc
