! 
! 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., A. Revil, and K.L. Kuhlman, 2009. A semi-analytical solution 
!  for transient streaming potentials associated with confined aquifer 
!  pumping tests, Geophysical Journal International, 176(3), 1007–1016.
! http://dx.doi.org/10.1111/j.1365-246X.2008.04014.x
!
! Malama, B., K.L. Kuhlman, and A. Revil, 2009. Theory of transient 
!  streaming potentials associated with axial-symmetric flow in unconfined 
!  aquifers, Geophysical Journal International, 179(2), 990–1003.
! http://dx.doi.org/10.1111/j.1365-246X.2009.04336.x
!

! $Id: utility.f90,v 1.3 2009/03/15 03:38:31 kris Exp kris $

!! this file has _three_ modules in it 

!! ######################################################################
!! ######################################################################

module constants

  ! integers with range of 9 and 18 orders of magnitude
  integer, parameter :: I4B = selected_int_kind(r=9)   ! std 32-bit integer

  ! real with range 300 orders of mag, 15 sig figs (8 on both g95 & ifort)
  integer, parameter :: DP = selected_real_kind (p=15,r=300)

  !! extended range internal variables (10 on g95, 16 on ifort)
  integer, parameter :: EP = selected_real_kind(r=3000)

  !! full quad precision
  integer, parameter :: QP = selected_real_kind(p=33,r=3000)

  ! useful? constants related to pi and ln
  ! calculated to precision=33 using ifort or mathematica
  real(DP), parameter :: PI =        3.141592653589793238462643383279503_DP
  real(DP), parameter :: TWOPI =     6.28318530717958647692528676655901_DP
  real(DP), parameter :: PIOV2 =     1.57079632679489661923132169163975_DP
  real(DP), parameter :: PIOV3 =     1.04719755119659774615421446109317_DP
  real(DP), parameter :: PIOV4 =     0.785398163397448309615660845819876_DP
  real(DP), parameter :: ONEOVPI =   0.318309886183790671537767526745029_DP
  real(DP), parameter :: LN10  =     2.30258509299404568401799145468436_DP
  real(DP), parameter :: LN2 =       0.693147180559945309417232121458177_DP
  real(DP), parameter :: LNPIOV2 =   0.451582705289454864726195229894882_DP 
  real(DP), parameter :: EULER =     0.5772156649015328606065120900824025_DP
  real(DP), parameter :: PISQ =      9.86960440108935861883449099987615_DP
  real(DP), parameter :: SQRTPI =    1.77245385090551602729816748334115_DP
  real(DP), parameter :: SQRTPIOV2 = 1.25331413731550025120788264240552_DP
  real(DP), parameter :: INVLOG10 =  0.434294481903251827651128918916605_DP

  real(EP), parameter :: LN2EP = 0.693147180559945309417232121458177_EP
  real(QP), parameter :: LN2QP = 0.693147180559945309417232121458177_QP
  real(EP), parameter :: PIEP = 3.141592653589793238462643383279503_EP

  ! a very small number (if < SMALL, it is effectively zero, compared to one)
  real(DP), parameter :: SMALL = epsilon(1.0_DP) ! approx 2.22E-16

  ! the largest number which can be represented accurately
  real(DP), parameter :: LARGE = huge(1.0_DP)  ! approx 1.79E+308

  ! commonly(?) used constants
  real(EP), parameter :: EONE = 1.0_EP
  real(EP), parameter :: EZERO = 0.0_EP
  real(EP), parameter :: EHALF = 0.5_EP
  real(EP), parameter :: ETWO = 2.0_EP
  complex(DP), parameter :: CZERO = (0.0_DP, 0.0_DP)
  complex(DP), parameter :: CONE  = (1.0_DP, 0.0_DP)
  complex(DP), parameter :: CTWO =  (2.0_DP, 0.0_DP)
  complex(DP), parameter :: EYE = (0.0_DP, 1.0_DP)
  real(DP),  parameter :: RZERO = 0.0_DP
  real(DP),  parameter :: RONE  = 1.0_DP
  real(DP),  parameter :: RTWO =  2.0_DP
  real(DP),  parameter :: HALF =  0.5_DP

end module constants

!! ######################################################################
!! ######################################################################

module shared_data
  use constants, only : DP,EP,QP
  implicit none

  integer, parameter :: MAXNLAYERS = 3

  ! variables to pass using modules

  ! parameters to change between calculations
  real(DP), save :: tsval   ! time - passed through hankel integral to Lap xform

  ! aquifer parameters (primary - read from input)
  real(DP), save, dimension(MAXNLAYERS) :: b,sigma,gamell
  real(DP), save :: Kr,Kz,Ss,Sy,zd,rd

  ! aquifer parameters (secondary - computed from above)
  real(DP), save :: alphaR, alphaY, alphaD, kappa

  real(DP), save, dimension(1:3) :: sigmaD

  real(DP), save :: bd1, bd3 
  logical, save, dimension(MAXNLAYERS) :: layer

  integer, save :: stehfestN
  real(QP), save, allocatable :: stehfestV(:)


end module shared_data

!! ######################################################################
!! ######################################################################

module invlap_fcns
implicit none
private
public :: stehfest_setup, stehfest, stehfest_coeff

contains
  
  function stehfest_coeff(N) result(V)
    use constants, only : QP
    integer, intent(in) :: N
    real(QP), dimension(N) :: V
    integer :: num
    real :: x
    
    if (mod(N,2) == 1) then
       print *, 'order of stehfest coefficient must be even,',N
       stop 
    end if

    select case(N)
       case(14)
          V = (/ 0.002777777777777777777777777777777777777778_QP,&
               & -6.402777777777777777777777777777777777778_QP, &
               & 924.05_QP, &
               & -34597.92777777777777777777777777777777778_QP, &
               & 540321.1111111111111111111111111111111111_QP, &
               & -4398346.366666666666666666666666666666667_QP, &
               & 21087591.77777777777777777777777777777778_QP, &
               & -63944913.04444444444444444444444444444444_QP, &
               & 127597579.55_QP, &
               & -170137188.0833333333333333333333333333333_QP,&
               & 150327467.0333333333333333333333333333333_QP, &
               & -84592161.5_QP, &
               & 27478884.76666666666666666666666666666667_QP, &
               & -3925554.966666666666666666666666666666667_QP /)
       case(16)
          V = (/ -0.0003968253968253968253968253968253968253968_QP, &
               & 2.133730158730158730158730158730158730159_QP, &
               & -551.0166666666666666666666666666666666667_QP, &
               & 33500.16111111111111111111111111111111111_QP, &
               & -812665.1111111111111111111111111111111111_QP, &
               & 10076183.76666666666666666666666666666667_QP, &
               & -73241382.97777777777777777777777777777778_QP, &
               & 339059632.0730158730158730158730158730159_QP, &
               & -1052539536.278571428571428571428571428571_QP, &
               & 2259013328.583333333333333333333333333333_QP, &
               & -3399701984.433333333333333333333333333333_QP, &
               & 3582450461.7_QP,  &
               & -2591494081.366666666666666666666666666667_QP, &
               & 1227049828.766666666666666666666666666667_QP, &
               & -342734555.4285714285714285714285714285714_QP, &
               & 42841819.42857142857142857142857142857143_QP /)
       case(18)
          V = (/ 0.0000496031746031746031746031746031746031746_QP, &
               & -0.6095734126984126984126984126984126984127_QP, &
               & 274.5940476190476190476190476190476190476_QP, &
               & -26306.95674603174603174603174603174603175_QP, &
               & 957257.2013888888888888888888888888888889_QP, &
               & -17358694.84583333333333333333333333333333_QP, &
               & 182421222.6472222222222222222222222222222_QP, &
               & -1218533288.309126984126984126984126984127_QP, &
               & 5491680025.283035714285714285714285714286_QP, &
               & -17362131115.2068452380952380952380952381_QP, &
               & 39455096903.52738095238095238095238095238_QP, &
               & -65266516985.175_QP,  &
               & 78730068328.22083333333333333333333333333_QP, &
               & -68556444196.12083333333333333333333333333_QP, &
               & 41984343475.05357142857142857142857142857_QP, &
               & -17160934711.83928571428571428571428571429_QP, &
               & 4204550039.102678571428571428571428571429_QP, &
               & -467172226.5669642857142857142857142857143_QP /)
       case default

          open(unit=27,file='stehfest_coefficients.dat',action='read',status='old')
          readV: do
             read(27,*) num
             if (num == N) then
                read(27,*) V
                exit readV
             elseif (num > N) then
                stop 'stehfestN used not supported'
             else
                read(27,*)  x ! throw away this data, onto 
             end if
          end do readV
          close(27)

    end select
    
  end function stehfest_coeff
  
  subroutine stehfest_setup(t,p)
    use constants, only : DP,EP,LN2EP
    use shared_data, only : stehfestN
    real(EP), dimension(stehfestN), intent(out) :: p
    real(DP), intent(in) :: t
    integer :: i

    p = LN2EP*real((/(i, i=1,stehfestN)/),EP)/t

  end subroutine stehfest_setup
  
  function stehfest(fp,t) result(f)
    use constants, only : DP,EP,LN2EP
    use shared_data, only : stehfestV
    real(EP), intent(in), dimension(:) :: fp
    real(DP), intent(in) :: t
    real(DP) :: f

    f = real(LN2EP/t*sum(stehfestV*fp(:)),DP)

  end function stehfest
  
end module invlap_fcns

