these Fortran90 source files comprise the streaming potential response
to a confined or unconfined (moving water table - aka Neuman 1972)
pumping test.

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),

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),

Code License
The code is available under the MIT open-source license, which allows
free use and redistribution, as long as you acknowledge the original
authors (by citing our original papers, please!)

Compilation Directions
The code requires either a recent version of the gfortran (>= 4.6)
compiler to properly compile (one that includes libquadmath), or a
somewhat recent version of the Intel Fortran compiler.

The code can be compiled, given you have the correct compiler set up 
already, by calling the compiler-specific makefile:

make -f Makefile.gfortran debug_driver
make -f Makefile.gfortran driver

these two makefile targets build the debug target (slower, checks for 
bounds overflows, etc) and the standard target (builds using OpenMP 
for a parallel build which is optimized and much faster. 

Input Instructions
The input.dat file has all the input parameters read by the program.
Parameters include physical properties of the formation, well geometry
parameters, numerical convergence parameters and the locations/times
to compute the solution at.

Anything to the right is a comment and is ignored by the program when
reading in the input.  Units are arbitrary but consistent.

----- example of input.dat -----

T  F  0  T  T  60            :: suppress output to screen?, dipole, BC type, hydrograph?, dist to bdry
4.10E-3                           :: pumping rates for all wells
5.0   20.0   10.0                  :: b1, b2, b3
4.0E-4  1.0  6.0E-5 3.756E-2        :: aquifer Kr, kappa, Ss, and Sy
1.0E-4  2.0E-3  8.0E-2            :: sigma for layers 1,2,3
1.0E-5  1.0E-5  1.0E-5            :: gamma*ell for layers 1,2,3
12                       :: number of times/radii to solve
6.633e+01  0.0  5.000e-02                  :: obs loc: x, y, depth
16                               :: Stehfest invlap N (must be even)
10  1  1  10  12  4             :: tanh-sinh 2^k-1 order, min/max split at, # zeros to integrate , order of G-Lquad, # extrapollation steps

----- end of example -----

Numerical Approach
If you are interested in the numerical approach used to invert the 
Laplace-Hankel transforms numerically, it is a slightly modified and 
updated version of the approach documented in Appendix B (repeated below 
citation) of the later of the two papers listed above: 

The numerical approach used to invert the double Laplace-Hankel
inversion is an improvement on the method in Malama et al. (2009). In
this study the Laplace transform is now inverted using the Stehfest
algorithm (Stehfest 1970) and a portion of the Hankel integral is
integrated using an accelerated form of tanh-sinh quadrature (Takahasi
& Mori 1974). The infinite integral in eq. (A2) for inverting the
Hankel transform is evalutated in two steps (Wieder 1999). One finite
(0 ≤ a ≤ j_{0,1}), the other approaching infinite length (j_{0,1} ≤ a ≤
j_{0,m}), where j_{0,n} is the nth zero of J_0(ar_D) and m is typically
10–12. The finite portion is evaluated using tanh-sinh quadrature
Takahasi & Mori (1974) where each level k of abcissa-weight pairs uses
a step size h = 4/2k and a number of abcissa N = 2k , similar to the
approach of Bailey et al. (2000). Evaluating the function at abcissa
related to level k, the function evalutations at lower levels k − j
are simply the 2^{k−j} th terms in the series (only requiring the
weights to be computed for each level, reusing the function
evaluations from the largest k). Richardson’s deferred approach to the
limit (Press et al. 2007) is used to extrapolate the approximate
solutions to the limit as h → 0 using a polynomial extrapolation. The
infinite portion of the Hankel integral is integrated between
successive j_{0,n} using Gauss–Lobatto quadrature and these alternating
sums are accelerated using Wynn’s ε-algorithm.

Contact Code Author
If you have any questions or require assistance using the code, please
contact me (Kris Kuhlman) at klkuhlm at sandia dot gov.  I will
try to help you.  I mostly work in Linux and on a Mac, but I have
gotten codes to work in Windows before.

July, 2009