      program Raypath
c-----------------------------------------------
c
c     Raypath < in par='ray_par_file' > out  [outa='amp_front_file']
c
c-----------------------------------------------------------------
c     Purpose: Compute ray parameters along rays in a model
c              defined by Splintab and rays initially given in
c              ray_par_file. Only for 2D velocity model.
c              The program uses netlib's Runge Kutte program in
c              rksuite.f.
c--------------------------------------------------------------------------
c     Input file: 
c                 fastest dimension, z-direction
c                 Second dimension, x-direction
c                 Third direction:  n3=1 parameters as in input 
c                                   n3=2 1st derivatives in x-directionm
c                                   n3=3 1st derivatives in z-directionm
c                                   n3=4 cross derivatives. (dx dz)
c
c------------------------------------------------------------------------
c     Output file: Gives ray paths that can be plotted with Graph.
c
c------------------------------------------------------------------------
c%     
      implicit NONE

      include '../Src/raypar.inc'
      include '../Src/model.inc'

      logical in_bound,in_bound_intern
      integer n1,n2,esize,esize_cmp,it,uflag,ndum,idelt,inum,one,ibuf
      integer n1_ref,n2_ref,i1_ref,i2_ref
      integer auxin,soufld,psoufld,auxpar,zero_offset,angsoufld
      integer fetch,getch,n1a,sreed,j1
      real*4 d1,d2

      integer maxt,ndelt,nsou,method,isou,itres,iter
      
      integer OUTCH
      double precision MCHEPS,DWARF

      double precision tstart,ystart(maxeq),tend,tol,thres(maxeq),
     *                 hstart,work(maxeq*32),twant,tgot,
     *                 ygot(maxeq),ypgot(maxeq),ymax(maxeq)
      
      real*4 velocitn
      character task
      integer neq,lenwrk,idump,ampyes,check,array_sou,n1out
      logical errass,mesage
      real*4 delta,delt,angle,xso,zso,vel0,vel1,toler,KMAH,y7old
      real*4 dx_sou,dz_sou,rot_ang
      real*4 xsou(msou),zsou(msou),delmin,xre,zre,err
      real*4 src(2,msou,1),psrc(2,msou,1),ang_sou(msou,1)

      real*4 rtraj(2,mtim,mangle)
      real*4 rpar(2,mtim,mangle)
      real*4 rmod(2,mtim,mangle)
      real*4 bufout(4*mtim*mangle)
      real*4 toaux(4,mangle,mtim)


      external parax2d

			write(0,*) "alpha"
      call initpar()

      call doc('/homes/sep/bob/Src/velan/tomo/Raypath/Raypath.f')

c     open dump file.

      idump = 0

			write(0,*) "alpha"
c      open(idump,file='Raypath.dump',status='unknown')

c     fetch dimensions and esize.

      if (fetch('n1','d',n1).eq.0) call erexit('n1 must be supplied')
      if (fetch('n2','d',n2).eq.0) call erexit('n2 must be supplied')
      if (fetch('j1','d',j1).eq.0) j1=1
      if (fetch('d1','f',d1).eq.0) call erexit('d1 must be supplied')
      if (fetch('d2','f',d2).eq.0) call erexit('d2 must be supplied')
      if (fetch('o1','f',zminco).eq.0) zminco = 0.0
      if (fetch('o2','f',xminco).eq.0) xminco = 0.0
      if (fetch('esize','d',esize).eq.0) esize=4

c     define some values for common block model.inc

      delx  = d2
      delz  = d1
      rdelx = 1./delx
      rdelz = 1./delz
      numx = n2
      numz = n1

      zlength=zminco+(numz-1)*delz
      xlength=xminco+(numx-1)*delx

      if (maxnod.lt.(n1*n2))  then
				write(0,*) "maxnod ",maxnod, "n1*n2",n1*n2
        call erexit('model too large')
			endif


c     read input file
      
      if (sreed('in',velp,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading velocity file")
      end if
      if (sreed('in',dvdx,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading dvdx file")
      end if
      if (sreed('in',dvdz,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading dvdz file")
      end if
      if (sreed('in',d2vdxz,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading d2vdxz file")
      end if

c     read parameter filevalues.

      if (getch('maxt','d',maxt).le.0) 
     *                  call erexit('maxt must be supplied')
      if (maxt.gt.mtim) call erexit('to many time steps to be recorded')
      if (getch('delt','f',delt).le.0) 
     *                  call erexit('delt must be supplied')

      ndum=getch('delmin','f',delmin)
      ndum=getch('delta','f',delta)
      ndum=getch('ndelt','d',ndelt)
      ndum=getch('nsou','d',nsou)

      ndum=getch('tolerance','f',toler)
      ndum=getch('method','d',method)
      ndum=getch('amplitude','d',ampyes)
      
c      write (idump,*) 'toler =',toler
c      write (idump,*) 'method =',method

      if (getch('zero_offset','d',zero_offset).le.0) zero_offset=0
      call putch('zero_offset','d',zero_offset)
			write(0,*) "alph1a"

      if(abs(zero_offset) .ne. 1) then
c       {
c     read src positions.
			  if(-1 .ne. auxin("Source")) then 
					ndum=auxpar("n1","d",nsou,"Source")

c{
           if(nsou*8.ne.sreed("Source",src,nsou*8)) then
             call erexit("trouble reading sources ")
						end if
           do isou=1,nsou
             xsou(isou)=src(1,isou,1)
             zsou(isou)=src(2,isou,1)
					 enddo
				 angsoufld=auxin('Angsrc')
					if(angsoufld .ne. -1) then
				    if(4*nsou .ne. sreed('Angsrc',ang_sou(1,1),nsou*4) )then
             call erexit('Problems reading Angsrc')
            end if
          end if
c}
				else
c{
       	ndum=getch('xsou','f',xsou)
      	 ndum=getch('zsou','f',zsou)
      	 if(getch('array_sou','d',array_sou) .eq. 0) array_sou=0
      	 if(array_sou .ne. 1) then
c{
       		 do isou=1,nsou
C      		   xsou(isou) = 0.0
C      		   ndum=getch('xsou','f',xsou(isou))
C      		   ndum=getch('xsou','f',xsou(isou))
C      		   zsou(isou) = 0.0
C      			   ndum=getch('zsou','f',zsou(isou))
      			 write(0,*)'isou=',isou,' zsou,xsou=',zsou(isou),',',xsou(isou)
       		 enddo
c}
      	 else
c{
       	  if(getch('dx_sou','f',dx_sou).eq.0)call seperr('Need dx_sou')
       	  if(getch('dz_sou','f',dz_sou).eq.0)call seperr('Need dz_sou')
       	  call putch('dx_sou','f',dx_sou)
       	  call putch('dz_sou','f',dz_sou)
       	  call auxputch('dx_sou','f',dx_sou,'Raypar')
       	  call auxputch('dz_sou','f',dz_sou,'Raypar')
       	  isou=1
      		 write(0,*)'isou=',isou,'zsou,xsou=',zsou(isou),',',xsou(isou)
       	  do isou=2,nsou
       	    xsou(isou)=xsou(1)+dx_sou*(isou-1)
       	    zsou(isou)=zsou(1)+dz_sou*(isou-1)
      		   write(0,*)'isou=',isou,'zsou,xsou=',zsou(isou),',',xsou(isou)
 					end do
c }      	 
				end if
c}
        end if
c}
       else
c {
c if wants zero offset rays, read reflector locations

        soufld=auxin('Source')
        psoufld=auxin('Psrc')
        if(soufld .le. 0) call seperr("Cannot open tag Source") 
        if(auxpar('n1','d',n1_ref,'Source') .eq. 0) then
           call seperr("Need n1 in Source")
        end if
        if(auxpar('n2','d',n2_ref,'Source') .eq. 0) then
           call seperr("Need n1 in Source")
        end if
        if(n2_ref .ne. 1 ) then
           call seperr("n2_ref must be eq 1 (for the moment)")
        end if

        nsou=n1_ref
        ndelt=1
        if(zero_offset .gt. 0) then
	  rot_ang=deg2rad*180.
        else
	  rot_ang=0
        end if
       
        do i2_ref=1,n2_ref
c        {
          if(sreed('Source',src(1,1,i2_ref),n1_ref*8).ne.
     *                                         n1_ref*8) then
           call erexit('Problems reading Source')
          end if

          if(psoufld .gt. 0)  then
c          {
            write(0,*)'Read Psrc'
            if(sreed('Psrc',psrc(1,1,i2_ref),n1_ref*8).ne.
     *                                             n1_ref*8) then
             call erexit('Problems reading Psrc')
            end if
c          }
          else
c          {
            do i1_ref=2,n1_ref-1
              ang_sou(i1_ref,i2_ref) = 
     *          (src(2,i1_ref+1,i2_ref)-src(2,i1_ref-1,i2_ref))/
     *          (src(1,i1_ref+1,i2_ref)-src(1,i1_ref-1,i2_ref))
            end do
            ang_sou(1,i2_ref) = 
     *          (src(2,2,i2_ref)-src(2,1,i2_ref))/
     *          (src(1,2,i2_ref)-src(1,1,i2_ref))
            ang_sou(n1_ref,i2_ref) = 
     *          (src(2,n1_ref,i2_ref)-src(2,n1_ref-1,i2_ref))/
     *          (src(1,n1_ref,i2_ref)-src(1,n1_ref-1,i2_ref))
            do i1_ref=1,n1_ref
              ang_sou(i1_ref,i2_ref)=
     *           -atan(ang_sou(i1_ref,i2_ref)) + rot_ang
            end do
c          }
         end if
c        }
          end do

c       }
      end if


      if (getch('check','d',check).le.0) check=0
      call putch('check','d',check)

c check limits in shot points and angles
      if (ndelt.gt.mangle) call erexit('Too many shot directions')
      if (nsou.gt.msou) call erexit('Too many shot points')


c     write updated trace length to new header

      esize_cmp = 8
      one=1
			n1out=1.0* maxt/j1+.9999
			write(0,*) "bb",maxt, j1, maxt/j1, 1.0*maxt/j1+.9999
      call putch('n1','d',n1out)
      call putch('n2','d',ndelt*nsou)
C      call putch('n3','d',nsou)
      call putch('n3','d',one)
      call putch('nsou','d',nsou)
      call putch('esize','d',esize_cmp)
      call putch('d1','f',delt*j1)
      call putch('d2','f',delta)
      call putch('xsou','f',xsou)
      call putch('zsou','f',zsou)
      call putch('maxt','d',maxt)
      call putch('delt','f',delt)
      call putch('delmin','f',delmin)
      call putch('delta','f',delta)
      call putch('ndelt','d',ndelt)
      call putch('toler','f',toler)
      call putch('method','d',method)
      call putch('ampyes','d',ampyes)
      call putch('array_sou','d',array_sou)

c     close history file     

      call hclose()

      call auxputch('n1','d',n1out,'Raypar')
      call auxputch('n2','d',ndelt*nsou,'Raypar')
      call auxputch('n3','d',one,'Raypar')
      call auxputch('nsou','d',nsou,'Raypar')
      call auxputch('esize','d',esize_cmp,'Raypar')
      call auxputch('d1','f',delt*j1,'Raypar')
      call auxputch('d2','f',delta,'Raypar')
      call auxputch('xsou','f',xsou,'Raypar')
      call auxputch('zsou','f',zsou,'Raypar')
      call auxputch('maxt','d',maxt,'Raypar')
      call auxputch('delt','f',delt,'Raypar')
      call auxputch('delmin','f',delmin,'Raypar')
      call auxputch('delta','f',delta,'Raypar')
      call auxputch('ndelt','d',ndelt,'Raypar')
      call auxputch('toler','f',toler,'Raypar')
      call auxputch('method','d',method,'Raypar')
      call auxputch('ampyes','d',ampyes,'Raypar')
      call auxputch('array_sou','d',array_sou,'Raypar')

      call auxputch('n1','d',n1out,'Raymod')
      call auxputch('n2','d',ndelt*nsou,'Raymod')
      call auxputch('n3','d',one,'Raymod')
      call auxputch('nsou','d',nsou,'Raymod')
      call auxputch('esize','d',esize_cmp,'Raymod')
      call auxputch('d1','f',delt*j1,'Raymod')
      call auxputch('d2','f',delta,'Raymod')
      call auxputch('xsou','f',xsou,'Raymod')
      call auxputch('zsou','f',zsou,'Raymod')
      call auxputch('maxt','d',maxt,'Raymod')
      call auxputch('delt','f',delt,'Raymod')
      call auxputch('delmin','f',delmin,'Raymod')
      call auxputch('delta','f',delta,'Raymod')
      call auxputch('ndelt','d',ndelt,'Raymod')
      call auxputch('toler','f',toler,'Raymod')
      call auxputch('method','d',method,'Raymod')
      call auxputch('ampyes','d',ampyes,'Raymod')
      call auxputch('array_sou','d',array_sou,'Raymod')

C      call auxclose('Raypar')


      if (ampyes.eq.1) then

c     write to auxiliary output file

         n1a = 4         
         ndum=1
         call auxputch('n1','d',n1a,'outa')
         call auxputch('n2','d',ndelt,'outa')
         call auxputch('n3','d',maxt,'outa')
         call auxputch('n4','d',nsou,'outa')
         call auxputch('d2','f',delta,'outa')
         call auxputch('d3','f',delt,'outa')
         call auxputch('esize','d',esize,'outa')
         call auxputch('d1','f',delt,'outa')
         call auxputch('d2','f',delta,'outa')
         call auxputch('xsou','f',xsou,'outa')
         call auxputch('zsou','f',zsou,'outa')
         call auxputch('maxt','d',maxt,'outa')
         call auxputch('delt','f',delt,'outa')
         call auxputch('delmin','f',delmin,'outa')
         call auxputch('delta','f',delta,'outa')
         call auxputch('ndelt','d',ndelt,'outa')
         call auxputch('toler','f',toler,'outa')
         call auxputch('method','d',method,'outa')
         call auxputch('ampyes','d',ampyes,'outa')


c     close auxiliary output history file

         call auxclose('outa')

c     initialize toaux

            do it=1,maxt
               do idelt = 1,ndelt
                  toaux(1,idelt,it) = 0.0
                  toaux(2,idelt,it) = 0.0
                  toaux(3,idelt,it) = 1.0e+30
                  toaux(4,idelt,it) = 0.0
               enddo
            enddo

      endif
      
c     find machine dependent values.

c      call envirn(OUTCH,MCHEPS,DWARF)

      OUTCH = 0
      MCHEPS = 1.11D-15
      DWARF = 2.23D-100

c     Define src independent values.
c     these can later be set in a parameter file.

      neq = 8
      lenwrk = neq*32
      tol = toler
      task = 'U'
      errass = .false.
      hstart = 0.0d0
      mesage = .true.

c     Loop over src positions.

      do isou=1,nsou
c{

         if(abs(zero_offset) .ne. 1) then
            ystart(1) = xsou(isou)
            ystart(2) = zsou(isou)
         else
           ystart(1) = src(1,isou,1)
           ystart(2) = src(2,isou,1)
         end if
            
         xso = ystart(1) 
         zso = ystart(2)
         vel0 = velocitn(xso,zso)

C         write(idump,*) 'isou =',isou
C         write(idump,*) 'vel0 =',vel0

         do idelt=1,ndelt
c  {

c     set KMAH index.
            
            KMAH = 0.0

c     find next take-off angle.

            if(abs(zero_offset) .ne. 1) then
               if(angsoufld .gt. 0) then
 	               ystart(3) = sin(deg2rad*ang_sou(isou,1))/vel0
 	               ystart(4) = cos(deg2rad*ang_sou(isou,1))/vel0
               else 
	               angle = (delmin + real(idelt-1)*delta)*deg2rad
 	               ystart(3) = sin(angle)/vel0
 	               ystart(4) = cos(angle)/vel0
               end if
            else
              if(psoufld .gt. 0)  then
                ystart(3) = - psrc(1,isou,1)
C                ystart(4) = - psrc(2,isou,1)
                ystart(4)=sqrt((1./(vel0*vel0))-(ystart(3)*ystart(3)))
              else
                angle=ang_sou(isou,1)
                ystart(3) = sin(angle)/vel0
                ystart(4) = cos(angle)/vel0
              end if
C         write(idump,*) 'ystart(1) =',ystart(1)
C         write(idump,*) 'ystart(2) =',ystart(2)
C         write(idump,*) 'px =',psrc(1,isou,1)
C         write(idump,*) 'pv =',psrc(2,isou,1)
C         write(idump,*) 'angle =',acos(ystart(4)*vel0)/deg2rad

            end if

            
c     initialize values.

            tstart = 0.d0
            tend = tstart + real(maxt-1)*delt + 0.00001

c     define ystart and threshold values.

            ystart(5) = 1.d0
            ystart(6) = 0.d0
            ystart(7) = 0.d0
            ystart(8) = 1.d0

            thres(1) = 0.01
            thres(2) = 0.01
            thres(3) = 0.000001
            thres(4) = 0.000001
            thres(5) = 1.e-4
            thres(6) = 1.e-4
            thres(7) = 1.e-4
            thres(8) = 1.e-4

c     store  src point

            if(check .ne. 1) then
                  do inum=1,2
                     rtraj(inum,1,idelt) = ystart(inum)
                     rpar(inum,1,idelt) = ystart(inum+2)
                  enddo
                  rmod(1,1,idelt) = vel0
                  rmod(2,1,idelt) = 0.
            else
                  vel1 = vel0
                  err = 
     *       1. -
     *       (vel1**2)*((ystart(4)**2) + (ystart(3)**2))
                  rtraj(1,1,idelt) = err
                  rtraj(2,1,idelt) = err

             end if

c     initialize computation.
            
            call setup(neq,tstart,ystart,tend,tol,thres,method,task,
     *           errass,hstart,work,lenwrk,mesage)

C            write (idump,*) 'tend =',tend
            twant = tstart
            y7old = 1.0

            in_bound=.true.
            do it=2,maxt
c       {
               twant = twant + delt
               if(in_bound) then
c              {
c     proceed one step.
               
C               write (idump,*) ' foer ut twant,tgot =',twant,tgot
               call ut(parax2d,twant,tgot,ygot,ypgot,ymax,work,uflag)
               in_bound_intern=
     *      ((ymax(1).ge.xminco).and.(ymax(1).le.xlength).and.
     *       (ymax(2).ge.zminco).and.(ymax(2).le.zlength)) 

C               if(.not.in_bound_intern) then
C                 write(0,*)'it=',it
C                 write(0,*)'ymax(1)=',ymax(1)
C                 write(0,*)'ymax(2)=',ymax(2)
C               end if
c               write (idump,*)'ygot,ymax',ygot,ymax               
C               write (idump,*)'uflag',uflag
               if (uflag.ne.1) then
c              {
C                  write (0,*)
C     *                        'Problems in subroutine ut uflag=',uflag
C                  write (idump,*)'twant,tgot',twant,tgot
C                  write (idump,*)'ygot,ymax',ygot,ymax
                  if (uflag.ne.3) then
                     do itres=it-1,maxt
                        do inum=1,2
                         rtraj(inum,itres,idelt)=rtraj(inum,it-2,idelt) 
                         rpar(inum,itres,idelt)=rpar(inum,it-2,idelt) 
                        enddo               
                        rmod(1,itres,idelt)=rmod(1,it-2,idelt) 
                        rmod(2,itres,idelt)=0.
                     enddo
C                     write (idump,*)'uflag =',uflag
                     goto 888
                  endif
                  iter = 1
                  do while (uflag.ne.1)
                  call ut(parax2d,twant,tgot,ygot,ypgot,ymax,work,uflag)

                  if ((ygot(7)*y7old).lt.0.0) KMAH = KMAH + 1.0  
                  y7old = ygot(7)

                  if (iter.gt.10) then
                     do itres=it-1,maxt
                        do inum=1,2
                         rtraj(inum,itres,idelt)=rtraj(inum,it-2,idelt) 
                         rpar(inum,itres,idelt)=rpar(inum,it-2,idelt) 
                        enddo               
                        rmod(1,itres,idelt)=rmod(1,it-2,idelt) 
                        rmod(2,itres,idelt)=0.
                     enddo
                     goto 888
                  else if (uflag.ne.1.and.uflag.ne.3) then
                     do itres=it-1,maxt
                        do inum=1,2
                         rtraj(inum,itres,idelt)=rtraj(inum,it-2,idelt) 
                         rpar(inum,itres,idelt)=rpar(inum,it-2,idelt) 
                        enddo               
                        rmod(1,itres,idelt)=rmod(1,it-2,idelt) 
                        rmod(2,itres,idelt)=0.
                     enddo
                     goto 888
                  endif
                  iter = iter + 1
                  enddo
c              }
               endif

               xre = ygot(1)
               zre = ygot(2)

               in_bound=in_bound_intern .and. 
     *                  ((xre.ge.xminco).and.(xre.le.xlength).and.
     *                   (zre.ge.zminco).and.(zre.le.zlength)) 

               if(in_bound) then
c              {
                 vel1 = velocitn(xre,zre)

                 if(check .ne. 1) then
                       do inum=1,2
                          rtraj(inum,it,idelt) = ygot(inum)
                          rpar(inum,it,idelt) = ygot(inum+2)
                       enddo
                       rmod(1,it,idelt) = vel1
                       rmod(2,it,idelt) = 0.
                 else
                    err = 
     *       1. -
     *       (vel1**2)*((ygot(4)**2) + (ygot(3)**2))
                    rtraj(1,it,idelt) = err
                    rtraj(2,it,idelt) = err
                 end if
  
                 if ((ygot(7)*y7old).lt.0.0) KMAH = KMAH + 1.0 
                 y7old = ygot(7)
                 if (ampyes.eq.1) then
                    xre = ygot(1)
                    zre = ygot(2)
                    vel1 = velocitn(xre,zre)
                    toaux(1,idelt,it) = xre
                    toaux(2,idelt,it) = zre
                    vel1 = sqrt(vel0/vel1)
                    toaux(3,idelt,it) = ygot(7)*vel1
                    toaux(4,idelt,it) = KMAH
                 endif
c              }
            else
c              {
                do inum=1,2
                  rtraj(inum,it,idelt)=rtraj(inum,it-1,idelt) 
                  rpar(inum,it,idelt)=rpar(inum,it-1,idelt) 
                  rmod(inum,it,idelt)=rmod(inum,it-1,idelt) 
                end do
c              }
            endif
  
c          }
            else
c          {
              do inum=1,2
                rtraj(inum,it,idelt)=rtraj(inum,it-1,idelt) 
                rpar(inum,it,idelt)=rpar(inum,it-1,idelt) 
                rmod(inum,it,idelt)=rmod(inum,it-1,idelt) 
              end do
c          }
            end if

c       }
            enddo 

 888        continue
                        
c   }
         enddo

         ibuf=0
         do idelt=1,ndelt
            do it=1,maxt,j1
               do inum=1,2
                ibuf = ibuf+1
                  bufout(ibuf) = rtraj(inum,it,idelt)
               enddo
            enddo
         enddo

      call srite('out',bufout,n1out*ndelt*esize_cmp)

         ibuf=0
         do idelt=1,ndelt
            do it=1,maxt,j1
               do inum=1,2
                ibuf = ibuf+1
                  bufout(ibuf) = rpar(inum,it,idelt)
               enddo
            enddo
         enddo

      call srite('Raypar',bufout,n1out*ndelt*esize_cmp)

         ibuf=0
         do idelt=1,ndelt
            do it=1,maxt,j1
               do inum=1,2
                ibuf = ibuf+1
                  bufout(ibuf) = rmod(inum,it,idelt)
               enddo
            enddo
         enddo

      call srite('Raymod',bufout,n1out*ndelt*esize_cmp)

      if (ampyes.eq.1) then

c     write to auxiliary otput file
         ibuf=0
         do idelt=1,ndelt
            do it=1,maxt,j1
               do inum=1,4
                ibuf = ibuf+1
                  bufout(ibuf) = toaux(inum,it,idelt)
               enddo
            enddo
         enddo

         call srite('outa',bufout,4*ndelt*n1out*esize)

      endif

c}
      enddo
         
      end
