      subroutine parax2d(t,y,yp)
      implicit NONE
      double precision t,y(*),yp(*)
c--------------------------------------------------------------      
c Purpose:
c          Input subroutine, which defines system of differential
c          equations for dynamic and paraxial ray tracing in 2-D
c          both source and receiver position can be altered using
c          a propagator matrix method to compute the wavefield.
c--------------------------------------------------------------      
c Input    
c       t     - the independent variable (traveltime along ray). 
c       y(*)  - a vector defining the variables.
c       yp(*) - the first derivatives of the vector y.
c
c       y(1) - x-coordinate
c       y(2) - z-coordinate
c       y(3) - slowness x-coordinate
c       y(4) - slowness z-coordinate
c       y(5) - fundamental propagator matrix element (1,1) q1
c       y(6) - fundamental propagator matrix element (2,1) p1
c       y(7) - fundamental propagator matrix element (1,2) q2
c       y(8) - fundamental propagator matrix element (2,2) p2
c
c       See Cerveny, In "Shear waves" (1985)
c--------------------------------------------------------------      
      real*8 vel,dvdx,dvdz,d2vdq2,vel2,rvel
      integer ieq

      call findvel(y,vel,dvdx,dvdz,d2vdq2)

c      write (0,*) ' i parax2d etter findvel vel,dvdx,dvdz,d2vdq2'
c      write (0,*) vel,dvdx,dvdz,d2vdq2

      if (vel.gt.0.01) then

         vel2 = vel*vel
         rvel = 1.0/vel
         
c     definition of the equations.

         yp(1) = vel2*y(3)
         yp(2) = vel2*y(4)
         yp(3) = -dvdx*rvel
         yp(4) = -dvdz*rvel
         yp(5) = vel2*y(6)
         yp(6) = -d2vdq2*y(5)*rvel
         yp(7) = vel2*y(8)
         yp(8) = -d2vdq2*y(7)*rvel
      else
         do ieq=1,8
            yp(ieq) = 0.d0
         enddo
      endif

      return
      end

      
