program laptest
c	
c  This program tests the numerical inverse laplace transform of a
c  function with an analytically-known inverse.  The function to
c  be transformed is exp(a*t).  

c  Use IMSL Numerical Libraries (comment out on most UNIX machines,
c  but needed if using a PC version of compiler.
c      use numerical_libraries

c  Allow only variables that are explcitly declared   
      implicit none

c  Set up all variables, vectors, and common blocks needed
	integer i, n, kmax
	parameter (n=30)
      double precision t(n),ft(n),ftrue(n),a,relerr,alpha
      complex*16 func
      external func
      common /prm/ a

c  Open a file to which results will be printed
      open (unit=16,file= 'laptest.out')

c  Define necessary parameters.  
c  a = parameter in exp(a*t)
c  alpha = An estimate for the maximum of the real parts of the 
c          singularities of F. If unknown, set ALPHA = 0.
c  kmax = The number of function evaluations allowed for each t(i).
c  relerr = The relative accuracy desired.
      a = 1.d0
      alpha  = 0.0d0
      kmax   = 1000
      relerr = 1.0d-6

c  Set up the vector of times to invert for, and calculate the true
c  value (ftrue) of the function exp(a*t)
      do 50,i=1,n
        t(i) = 10.d0**(dble(i)*0.5d0-13.d0)
        ftrue(i) = exp(a*t(i))
   50 continue

c  Call the inverse laplace transform.
      call dinlap (func,n,t,alpha,relerr,kmax,ft)

c  Write out the values to the file. The last number is the true relative
c  error.
      do 60,i=1,n
        write(16,'(4e15.8)') t(i),ft(i),ftrue(i),(ft(i)/ftrue(i)-1.d0)
   60 continue

c  Close up the file and end the program.
      close (16)
      stop
      end

c***********************************************************************
      complex*16 function func(p)
c  This function contains the Laplace-transform of the function exp(a*t).
c  This transform is 1/(p-a).

c  Set up the function parameters, etc.
      implicit none
      complex*16 p
      double precision a
      common /prm/ a

c  Define the function to be inverted.
      func = 1.d0/(p-a)

c  Return to the routine that called func(p).
      return
      end