c   Although this program has been used by the USGS, no warranty,
c   expressed or implied, is made by the USGS as to the accuracy
c   and functioning of the program and related program material
c   nor shall the fact of distribution constitute any such warranty,
c   and no responsibility is assumed by the USGS in connection herewith.
c   
c   
c  recov.f  recovery   superposition of positive and negative theis
c    new variable is dimensionless time at which pump turned off.
c
c  Dan Goode                   djgoode@usgs.gov
c  U.S. Geological Survey
c  111 Great Valley Pkwy
c  Malver PA 19355  USA
c
c  see Goode, D.J., 1996, New type curves for estimation of aquifer 
c         properties from water-level recovery under Theis conditions: 
c         (abs.), p. S95 in 1996 Spring Meeting, supplement to Eos, 
c         23 April 1996, American Geophysical Union, Washington DC.
c
c  This program writes output to file recov.out
c  The program asks the user to enter the dimensionless pumping duration
c    from the keyboard.  The program then generates the recover 
c    solution in normalized time.
c
      parameter(small=1.e-20)
      iout=9
      open(iout,file='recov.out',status='unknown')
c
      write(*,*) ' recov    enter dimensionless pump duration'
      read(*,*) tpd
      write(*,*) ' dimensionless pumping duration =',tpd
c
      nstep=301
      cycles=6.d0
      stplog=cycles/float(nstep-1)
c
      write(iout,*) ' tpd=',tpd
      write(iout,*) ' i  tn  sd'
      write(*,*) ' i  tn  sd'
c
      stp1=-3.d0+alog10(tpd)
      write(*,*) ' stp1=',stp1
      stp=stp1-stplog
c
      do 20 i=1,nstep
        stp=stp+stplog
        td=10.0**stp
        tn=td/tpd
        u1=1.0/(4.0*(td+tpd))
        u2=1.0/(4.0*td)
        if(u1.lt.675.0) then
           sd1=E1(u1)
        else
           sd1=small
        end if
        if(u2.lt.675.0) then
           sd2=E1(u2)
        else
           sd2=small
        end if
        sd=sd1-sd2
        write(iout,1000) i,tn,sd
        write(*,1000) i,tn,sd
 20   continue
c
 1000 format(i10,1p2e12.4)
      stop
      end
c***********************************************
      function E1(x)
c***********************************************
c
c exint1.f   exponential integral E1(x)
c
c  Abramowitz, M., and I.A. Stegun (eds.), 1970, Handbook of mathematical
c    functions, National Bureau of Standards, Applied Mathematics Series
c    55, Washington, D.C.
c
c  page 231, eqn. 5.1.53 and 5.1.56
c
c  Dan Goode
c
      parameter (bignum=1.e30)
      double precision dterm,dx
c
      dx=dble(x)
c
      if(x.le.0.0) then
         E1=bignum
c
      else if(x.le.1.0) then
         dterm= -0.57721566d0 + dx*( 0.99999193d0 + 
     *     dx*( -0.24991055d0 + dx*( 0.05519968d0 +
     *     dx*( -0.00976004d0 + dx*  0.00107857d0 ))))
         E1=dterm-dlog(dx)
c
      else
         dterm=( 0.2677737343d0 + dx*(  8.6347608925d0 +
     *     dx*( 18.0590169730d0 + dx*(  8.5733287401d0 + dx))))/
     *      (    3.9584969228d0 + dx*( 21.0996530827d0 + 
     *     dx*( 25.6329561486d0 + dx*(  9.5733223454d0 + dx))))
         E1=dterm/(dx*dexp(dx))
c
      end if
c
      return
      end