! 
! Copyright (c) 2008 Kristopher L. Kuhlman (kuhlman at hwr dot arizona dot edu)
! 
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
! 
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
! 
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! 
! Malama, B., K.L. Kuhlman, and W. Barrash, 2007. Semi-analytical solution
!  for flow in leaky unconfined aquifer-aquitard systems, Journal of
!  Hydrology, 346(1-2), 59–68.
! http://dx.doi.org/10.1016/j.jhydrol.2007.08.018
!
! Malama, B., K.L. Kuhlman, and W. Barrash, 2008. Semi-analytical solution
!  for flow in a leaky unconfined aquifer toward a partially penetrating
!  pumping well, Journal of Hydrology, 356(1-2), 234–244.
! http://dx.doi.org/10.1016/j.jhydrol.2008.03.029
!

! $Id: invlap.f90,v 1.4 2007/01/06 16:55:03 kris Exp kris $
module invlap
implicit none

private
public :: dehoog_invlap, dehoog_invlap_ep

contains

  !! an implementation of the de Hoog method
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function deHoog_invLap(alpha,tol,t,tee,fp,M) result(ft)
    use constants, only : DP, PI, CZERO, CONE, EYE, HALF, RONE, RZERO
    real(DP), intent(in) :: alpha  ! abcissa of convergence
    real(DP), intent(in) :: tol    ! desired accuracy
    real(DP), intent(in) :: t,tee  ! time and scaling factor
    integer, intent(in) :: M
    complex(DP), intent(in), dimension(0:2*M) :: fp
    real(DP) :: ft ! output

    real(DP) :: gamma
    integer :: r, rq, n, max
    complex(DP), dimension(0:2*M,0:M) :: e
    complex(DP), dimension(0:2*M,1:M) :: q
    complex(DP), dimension(0:2*M) :: d
    complex(DP), dimension(-1:2*M) :: A,B
    complex(DP) :: brem,rem,z
    
    ! there will be problems is fp(:)=0
    if(maxval(abs(fp)) > epsilon(1.0_DP)) then

       ! Re(p) -- this is the de Hoog parameter c
       gamma = alpha - log(tol)/(2.0_DP*tee)

       ! initialize Q-D table 
       e(0:2*M,0) = CZERO
       q(0,1) = fp(1)/(fp(0)*HALF) ! half first term
       q(1:2*M-1,1) = fp(2:2*M)/fp(1:2*M-1)

       ! rhombus rule for filling in triangular Q-D table
       do r = 1,M
          ! start with e, column 1, 0:2*M-2
          max = 2*(M-r)
          e(0:max,r) = q(1:max+1,r) - q(0:max,r) + e(1:max+1,r-1)
          if (r /= M) then
             ! start with q, column 2, 0:2*M-3
             rq = r+1
             max = 2*(M-rq)+1
             q(0:max,rq) = q(1:max+1,rq-1) * e(1:max+1,rq-1) / e(0:max,rq-1)
          end if
       end do

       ! build up continued fraction coefficients
       d(0) = fp(0)*HALF ! half first term
       forall(r = 1:M)
          d(2*r-1) = -q(0,r) ! even terms
          d(2*r)   = -e(0,r) ! odd terms
       end forall

       ! seed A and B vectors for recurrence
       A(-1) = CZERO
       A(0) = d(0)
       B(-1:0) = CONE

       ! base of the power series
       z = exp(EYE*PI*t/tee)

       ! coefficients of Pade approximation
       ! using recurrence for all but last term
       do n = 1,2*M-1
          A(n) = A(n-1) + d(n)*A(n-2)*z
          B(n) = B(n-1) + d(n)*B(n-2)*z
       end do

       ! "improved remainder" to continued fraction
       brem = (CONE + (d(2*M-1) - d(2*M))*z)*HALF
       rem = -brem*(CONE - sqrt(CONE + d(2*M)*z/brem**2))

       ! last term of recurrence using new remainder
       A(2*M) = A(2*M-1) + rem*A(2*M-2)
       B(2*M) = B(2*M-1) + rem*B(2*M-2)

       ! diagonal Pade approximation
       ! F=A/B represents accelerated trapezoid rule
       ft =  exp(gamma*t)/tee * real(A(2*M)/B(2*M))

    else
       ft = RZERO
    end if

  end function deHoog_invLap

  !! an implementation of the de Hoog method (extended range)
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function deHoog_invLap_ep(alpha,tol,t,tee,fp,M) result(ft)
    use constants, only : DP, EP
    real(DP), intent(in) :: alpha  ! abcissa of convergence
    real(DP), intent(in) :: tol    ! desired accuracy
    real(DP), intent(in) :: t,tee  ! time and scaling factor
    integer, intent(in) :: M
    complex(EP), intent(in), dimension(0:2*M) :: fp
    real(EP) :: ft ! output

    real(EP) :: alphaEP,tolEP,tEP,teeEP

    real(EP) :: gamma
    integer :: r, rq, n, max
    complex(EP), dimension(0:2*M,0:M) :: e
    complex(EP), dimension(0:2*M,1:M) :: q
    complex(EP), dimension(0:2*M) :: d
    complex(EP), dimension(-1:2*M) :: A,B
    complex(EP) :: brem,rem,z
    
    real(EP), parameter :: PIEP = 3.141592653589793238462643383279503_EP

    alphaEP = real(alpha,EP); tolEP = real(tol,EP)
    tEP = real(t,EP); teeEP = real(tee,EP)

    ! there will be problems is fp(:)=0
    if(maxval(abs(fp)) > epsilon(1.0_EP)) then

       ! Re(p) -- this is the de Hoog parameter c
       gamma = alphaEP - log(tolEP)/(2.0_EP*teeEP)

       ! initialize Q-D table 
       e(0:2*M,0) = 0.0_EP
       q(0,1) = fp(1)/(fp(0)*0.5_EP) ! half first term
       q(1:2*M-1,1) = fp(2:2*M)/fp(1:2*M-1)

       ! rhombus rule for filling in triangular Q-D table
       do r = 1,M
          ! start with e, column 1, 0:2*M-2
          max = 2*(M-r)
          e(0:max,r) = q(1:max+1,r) - q(0:max,r) + e(1:max+1,r-1)
          if (r /= M) then
             ! start with q, column 2, 0:2*M-3
             rq = r+1
             max = 2*(M-rq)+1
             q(0:max,rq) = q(1:max+1,rq-1) * e(1:max+1,rq-1) / e(0:max,rq-1)
          end if
       end do

       ! build up continued fraction coefficients
       d(0) = fp(0)*0.5_EP ! half first term
       forall(r = 1:M)
          d(2*r-1) = -q(0,r) ! even terms
          d(2*r)   = -e(0,r) ! odd terms
       end forall

       ! seed A and B vectors for recurrence
       A(-1) = 0.0_EP
       A(0) = d(0)
       B(-1:0) = 1.0_EP

       ! base of the power series
       z = exp((0.0_EP,1.0_EP)*PIEP*tEP/teeEP)

       ! coefficients of Pade approximation
       ! using recurrence for all but last term
       do n = 1,2*M-1
          A(n) = A(n-1) + d(n)*A(n-2)*z
          B(n) = B(n-1) + d(n)*B(n-2)*z
       end do

       ! "improved remainder" to continued fraction
       brem = (1.0_EP + (d(2*M-1) - d(2*M))*z)*0.5_EP
       rem = -brem*(1.0_EP - sqrt(1.0_EP + d(2*M)*z/brem**2))

       ! last term of recurrence using new remainder
       A(2*M) = A(2*M-1) + rem*A(2*M-2)
       B(2*M) = B(2*M-1) + rem*B(2*M-2)

       ! diagonal Pade approximation
       ! F=A/B represents accelerated trapezoid rule
       ft =  exp(gamma*tEP)/teeEP * real(A(2*M)/B(2*M))

    else
       ft = 0.0_EP
    end if

  end function deHoog_invLap_ep
end module invlap