c LAYERTRACE calculates the travel time and range offset
c for ray tracing through a single layer.
c
c Input:    p     =  horizontal slowness
c           h     =  layer thickness
c           utop  =  slowness at top of layer
c           ubot  =  slowness at bottom of layer
c           imth  =  interpolation method
c                    imth = 1,  v(z) = 1/sqrt(a - 2*b*z)
c                         = 2,  v(z) = a - b*z
c                         = 3,  v(z) = a*exp(-b*z)
c
c Returns:  dx    =  range offset
c           dt    =  travel time
c           irtr  =  return code
c                 = -1, zero thickness layer
c                 =  0,  ray turned above layer
c                 =  1,  ray passed through layer
c                 =  2,  ray turned within layer, 1 segment counted
c
c Note:  This version does calculation in double precision,
c        but all i/o is still single precision
c
      subroutine LAYERTRACE(p1,h1,utop1,ubot1,imth,dx1,dt1,irtr)
      implicit real*8 (a-h,o-z)
      real*4 p1,h1,utop1,ubot1,dx1,dt1
      p=dble(p1)
      h=dble(h1)
      utop=dble(utop1)
      ubot=dble(ubot1)
c
      if (h.eq.0.) then      !check for zero thickness layer
         dx1=0.
         dt1=0.
         irtr=-1
         return         
      end if
c
      u=utop
      y=u-p
      if (y.le.0.) then   !complex vertical slowness
         dx1=0.
         dt1=0.
         irtr=0
         return
      end if
c
      q=y*(u+p)
      qs=dsqrt(q)
c
c special function needed for integral at top of layer
      if (imth.eq.2) then
         y=u+qs
         if (p.ne.0.) y=y/p
         qr=dlog(y)
      else if (imth.eq.3) then
         qr=atan2(qs,p)
      end if      
c
      if (imth.eq.1) then
          b=-(utop**2-ubot**2)/(2.*h)
      else if (imth.eq.2) then
          vtop=1./utop
          vbot=1./ubot
          b=-(vtop-vbot)/h
      else
          b=-dlog(ubot/utop)/h
      end if  
c
      if (b.eq.0.) then     !constant velocity layer
         b=1./h
         etau=qs
         ex=p/qs
         irtr=1
         go to 160
      end if
c
c integral at upper limit, 1/b factor omitted until end
      if (imth.eq.1) then
         etau=-q*qs/3.
         ex=-qs*p
      else if (imth.eq.2) then
         ex=qs/u                    !*** - in some versions (wrongly)
         etau=qr-ex
         if (p.ne.0.) ex=ex/p
      else
         etau=qs-p*qr
         ex=qr
      end if
c
c check lower limit to see if we have turning point
      u=ubot
      if (u.le.p) then   !if turning point,
         irtr=2          !then no contribution
         go to 160       !from bottom point
      end if 
      irtr=1
      q=(u-p)*(u+p)
      qs=dsqrt(q)
c
      if (imth.eq.1) then
         etau=etau+q*qs/3.
         ex=ex+qs*p
      else if (imth.eq.2) then
         y=u+qs
         z=qs/u
         etau=etau+z
         if (p.ne.0.) then
            y=y/p
            z=z/p
         end if
         qr=dlog(y)
         etau=etau-qr
         ex=ex-z
      else
         qr=atan2(qs,p)
         etau=etau-qs+p*qr
         ex=ex-qr
      end if      
c
160   dx=ex/b
      dtau=etau/b
      dt=dtau+p*dx     !convert tau to t
c
      dx1=sngl(dx)
      dt1=sngl(dt)
      return
      end