######################################################################
#                                                                    #
#     MAPLE V programs for Leaky Aquifer Function calculations       #
#                                                                    #
#   Frank E. Harris, Department of Chemistry, University of Utah     #
#    Salt Lake City, UT 84112 (e-mail harris@dirac.chem.utah.edu)    #
#                                                                    #
#   Reference:  paper submitted to Int. J. Quantum Chem. October     #
#    1996 by Frank E. Harris                                         #
#                                                                    #
#        These programs last modified October 23, 1996               #
#                                                                    #
######################################################################
#                                                                    #
#   Programs included:                                               #
#                                                                    #
#      FF(x,y,n): Computes integral (1 to infinity) exp(-xt-y/t)     #
#                  /t^(n+1) dt.  For n=0 this is the Leaky           #
#                  Aquifer function, denoted W(x,y);  non-zero n     #
#                  gives generalizations.  Computation is by Maple's #
#                  numerical integration procedures and is to the    #
#                  accuracy specified in the environmental variable  #
#                  "Digits".                                         #
#                                                                    #
#      Eint(n,x): Computes the generalized exponential integral      #
#                 E sub n (x) (as defined in Abramowitz and Stegun)  #
#                 by upward recursion from E sub 1, which is -Ei(-x).#
#                                                                    #
#      FF1(x,y,nmax): Computes W(x,y) as an expansion in powers of   #
#                  y; coefficients (which depend upon x) involve     #
#                  the generalized exponential integral functions.   #
#                  Expansion through y^nmax.                         #
#                                                                    #
#      FF2(x,y,nmax): Computes W(x,y) as an expansion whose leading  #
#                  term is a modified Bessel function and whose      #
#                  remaining terms contain successive powers of x,   #
#                  Expansion coefficients (depending on y) involve   #
#                  the generalized exponential integral functions.   #
#                  Expansion through x^nmax.                         #
#                                                                    #
#      FF3(x,y,nmax): Computes W(x,y) as a power series in v-1,      #
#                  where v = (x/y)^1/2.  Expansion through           #
#                  (v-1)^nmax.                                       #
#                                                                    #
#      FF4(x,y,nmax): Computes W(x,y) as a power series in v^-1 - 1, #
#                  where v^-1 = (y/x)^1/2.  Expansion through        #
#                  (v^-1 - 1)^nmax.                                  #
#                                                                    #
#      FFTST(x,y,nmax): Computes same function as FF3 (up to nmax    #
#                  = 14), but with hard-coded terms in place of      #
#                  dynamically generated ones (provides a check on   #
#                  Table I of the paper)                             #
#                                                                    #
######################################################################

Digits := 14;       # Change if desire different precision

FF := proc(x,y,n) local t;
evalf(int(E^(-x*t-y/t)/t^(n+1),t=1..infinity)) end;

Eint := proc(n,x);
if (n < 1) then RETURN('procname(args)') fi;
if (n = 1) then RETURN(evalf(-Ei(-x))) fi;
evalf((E^(-x) - x*Eint(n-1,x))/(n-1)) end;

FF1 := proc(x,y,nmax) local n;
evalf(sum('(-y)^n*Eint(n+1,x)/n!',n=0..nmax)) end;

FF2 := proc(x,y,nmax);
evalf(2*BesselK(0,2*sqrt(x*y)) - FF1(y,x,nmax)) end;

FF3 := proc(x,y,nmax) local F,G,t,u,v,m,n;
G := array(1..nmax);
u := sqrt(x*y);
v := sqrt(x/y);
F := -E^(-u*(t + 1/t))/t;
G[1] := subs(t=1,F);
m := 2;
while (m <= nmax) do
F := diff(F,t);
G[m] := subs(t=1,F);
m := m + 1; od;
evalf(BesselK(0,2*u) + sum('(v-1)^n*G[n]/n!',n=1..nmax)) end;

FF4 := proc(x,y,nmax);
evalf(2*BesselK(0,2*sqrt(x*y)) - FF3(y,x,nmax)) end;

FFTST := proc(x,y,nmax) local G;
if nmax > 14 then RETURN('procname(args)') fi;
G := array(1..14);
u := sqrt(x*y);
v := sqrt(x/y);
G[1] := -1;
G[2] := 1/2;
G[3] := (u-1)/3;
G[4] := (-2*u+1)/4;
G[5] := (-u^2+6*u-2)/10;
G[6] := (3*u^2-8*u+2)/12;
G[7] := (u^3-18*u^2+30*u-6)/42;
G[8] := (-2*u^3+15*u^2-18*u+3)/24;
G[9] := (-u^4+40*u^3-180*u^2+168*u-24)/216;
G[10] := (5*u^4-80*u^3+252*u^2-192*u+24)/240;
G[11] := (u^5-75*u^4+700*u^3-1680*u^2+1080*u-120)/1320;
G[12] := (-6*u^5+175*u^4-1120*u^3+2160*u^2-1200*u+120)/1440;
G[13] := (-u^6+126*u^5-2100*u^4+10080*u^3-16200*u^2+7920*u-720)/9360;
G[14] := (7*u^6-336*u^5+3780*u^4-14400*u^3+19800*u^2-8640*u+720)/10080;
evalf(BesselK(0,2*u) + E^(-2*u)*sum('(v-1)^n*G[n]',n=1..nmax)) end;

##########################That's all #################################