      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, tau-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 tau-directionm
c                                   n3=4 cross derivatives. (dx dtau)
c                                   n3=5 sigma
c                                   n3=6 d_sigma/dx
c                                   n3=5 d_sigma/dtau
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
      integer fetch,getch,n1a,sreed
      real*4 d1,d2

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

      double precision tstart,tend,tol,thres(maxeq),
     *                 hstart,work(maxeq*32),twant,tgot,
     *                 ymax(maxeq)
      double precision ygot(maxeq),ypgot(maxeq),ystart(maxeq)
      
      real*4 velocitn,sigmatn
      character task
      integer neq,lenwrk,idump,ampyes,zmap,check,array_sou
      integer ixray,itauray,auxin,zero_offset,auxpar
      real*4 xray,tauray,wx,wtau,ztaumap,zmax,err,px,ptau
      logical errass,mesage
      real*4 delta,delt,angle,xso,tauso,vel0,vel1,sig0,sig1
      real*4 dx_sou,dtau_sou
      real*4 toler,KMAH,y7old
      real*4 xsou(msou),tausou(msou),delmin,xre,taure,c_source(2,msou)

      real*4 rtraj(2,mtim,mangle)
      real*4 rpar(2,mtim,mangle)
      real*4 rmod(2,mtim,mangle)
      real*4 toaux(4,mangle,mtim),bufout(4*mangle*mtim)
      real*4 ztau(maxnod)
      external parax2d

      call initpar()

      call doc(
     * '/home/kana/sep/biondo/Code/Focus/Raytrace/Raytau/Raytau.f')

c     open dump file.

      idump = 0

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('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',tauminco).eq.0) tauminco = 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
      deltau  = d1
      rdelx = 1./delx
      rdeltau = 1./deltau
      numx = n2
      numtau = n1

      taulength=tauminco+(numtau-1)*deltau
      xlength=xminco+(numx-1)*delx

      if (maxnod.lt.(n1*n2)) 
     *     call erexit('input model too large')


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',dvdtau,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading dvdtau file")
      end if
      if (sreed('in',d2vdxtau,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading d2vdxtau file")
      end if
      if (sreed('in',sigmap,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading sigmap file")
      end if
      if (sreed('in',dsigdx,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading dsigdx file")
      end if
      if (sreed('in',dsigdtau,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading dsigdtau file")
      end if
      if (sreed('in',d2sigdxtau,n1*n2*esize) .ne. n1*n2*esize) then
	 call erexit("Problems reading d2sigdxtau 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('too 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)
      if (ndelt.gt.mangle) call erexit('Too many shot directions')
      ndum=getch('nsou','d',nsou)
      if (nsou.gt.msou) call erexit('Too many shot points')

      ndum=getch('tolerance','f',toler)
      ndum=getch('method','d',method)
      ndum=getch('amplitude','d',ampyes)
      ndum=getch('zero_offset','d',zero_offset)

      
c      write (idump,*) 'toler =',toler
c      write (idump,*) 'method =',method

c     read source positions.
	
		   if(zero_offset .eq. 1) then
				write(0,*) "setting stuff",nsou
			  if(0 .eq. auxpar("n1","d",nsou,"Source"))  then
					call erexit("trouble reading in number of source locations")
				end if
				if(8*nsou .ne. sreed("Source",c_source,8*nsou)) then
					call erexit("trouble reading source locations")
				end if

				do isou=1,nsou
					xsou(isou)=c_source(1,isou)
					tausou(isou)=c_source(2,isou)
				end do


		   else

       ndum=getch('xsou','f',xsou)
       ndum=getch('tausou','f',tausou)
       if(getch('array_sou','d',array_sou) .eq. 0) array_sou=0
       if(array_sou .ne. 1) then
        do isou=1,nsou
C         xsou(isou) = 0.0
C         ndum=getch('xsou','f',xsou(isou))
C         ndum=getch('xsou','f',xsou(isou))
C         tausou(isou) = 0.0
C         ndum=getch('tausou','f',tausou(isou))
	write(0,*)'isou=',isou,'tausou,xsou=',
     *             tausou(isou),',',xsou(isou)
        enddo
       else
         if(getch('dx_sou','f',dx_sou).eq.0)
     *             call seperr('Need dx_sou')
         if(getch('dtau_sou','f',dtau_sou).eq.0)
     *             call seperr('Need dtau_sou')
         call putch('dx_sou','f',dx_sou)
         call putch('dtau_sou','f',dtau_sou)
         call auxputch('dx_sou','f',dx_sou,'Raypar')
         call auxputch('dtau_sou','f',dtau_sou,'Raypar')
         call auxputch('dx_sou','f',dx_sou,'Raymod')
         call auxputch('dtau_sou','f',dtau_sou,'Raymod')
         isou=1
	 write(0,*)'isou=',isou,'tausou,xsou=',
     *             tausou(isou),',',xsou(isou)
         do isou=2,nsou
           xsou(isou)=xsou(1)+dx_sou*(isou-1)
           tausou(isou)=tausou(1)+dtau_sou*(isou-1)
	   write(0,*)'isou=',isou,'tausou,xsou=',
     *             tausou(isou),',',xsou(isou)
         end do

       end if
	     end if

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

      if (getch('zmap','d',zmap).le.0) zmap=0
      call putch('zmap','d',zmap)
      if(zmap .eq. 1) then

C         if(auxin('Zmap' .eq. 0)) call seperr('Zmap file not found')

         if (sreed('Zmap',ztau,n1*n2*esize) .ne. n1*n2*esize) then
	     call erexit("Problems reading Zmap file")
         end if
C         zmax=-100000.
C         do inum=1,n1*n2
C           if(ztau(inum) .gt. zmax) zmax=ztau(inum)
C         end do
      end if
 

c     write updated trace length to new header

      esize_cmp = 8
      one=1
      call putch('n1','d',maxt)
      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)
      call putch('d2','f',delta)
      call putch('o1','f',0.)
      call putch('o2','f',delmin)
      call putch('xsou','f',xsou)
      call putch('tausou','f',tausou)
      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)

c     close history file     

      call hclose()

      call auxputch('n1','d',maxt,'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,'Raypar')
      call auxputch('d2','f',delta,'Raypar')
      call auxputch('xsou','f',xsou,'Raypar')
      call auxputch('tausou','f',tausou,'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',maxt,'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,'Raymod')
      call auxputch('d2','f',delta,'Raymod')
      call auxputch('xsou','f',xsou,'Raymod')
      call auxputch('tausou','f',tausou,'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')


      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('tausou','f',tausou,'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.

			write(0,*) "here i am one"

c      call envirn(OUTCH,MCHEPS,DWARF)

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

c     Define source 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 source positions.
			write(0,*) "here i am two",nsou

      do isou=1,nsou
c{
				write(0,*) "here i am three",isou

         xso = xsou(isou)
         tauso = tausou(isou)
         vel0 = velocitn(xso,tauso)

         xso = xsou(isou)
         tauso = tausou(isou)
         sig0 = sigmatn(xso,tauso)

         write(idump,*) 'xso =',xso
	write(0,*)'tauso=',tauso
         write(idump,*) 'vel0 =',vel0
	write(0,*)'sig0=',sig0
         do idelt=1,ndelt
c  {

c     set KMAH index.
            
            KMAH = 0.0

c     find next take-off angle.

            angle = (delmin + real(idelt-1)*delta)*deg2rad
            
c     initialize values.

            tstart = 0.d0
            tend = tstart + real(maxt-1)*delt + 0.00001
C            write (idump,*) 'tstart tend,delt'
C            write (idump,*) tstart,tend,delt

c     define ystart and threshold values.
          
            ystart(1) = xsou(isou)
            ystart(2) = tausou(isou)

C            ystart(3) = sin(angle)/vel0
C            ystart(4) = cos(angle)/vel0

            ystart(4) = .5*cos(angle)
            ystart(3) = (sin(angle)/vel0 - sig0*ystart(4))

C            ystart(3) = (sin(angle)/vel0 - sig0*.5*cos(angle))

            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  source 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) = sig0
            else
                  vel1 = vel0
                  sig1 = sig0
                  px=ystart(3)
                  ptau=ystart(4)

                  err =  1.- 
     *       (4.*(ptau**2) + (vel1**2)*(px+sig1*ptau)**2)

C            ystart(4) = .5*cos(angle)
C            ystart(3) = (sin(angle)/vel0 - sig0*ystart(4))

                  rtraj(1,1,idelt) = err
                  rtraj(2,1,idelt) = 
     *              1.-(cos(angle)**2+sin(angle)**2)

             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.tauminco).and.(ymax(2).le.taulength)) 

c               write (idump,*)'ygot,ymax',ygot,ymax               
C               write (idump,*)'uflag',uflag
               if (uflag.ne.1) then
c {
C                  write (idump,*)
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
c	{
                     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) 
                           rmod(inum,itres,idelt) =
     *                          rmod(inum,it-2,idelt) 
                        enddo               
                     enddo
C                     write (idump,*)'uflag =',uflag
                     goto 888
c	}
                  endif
                  iter = 1
                  do while (uflag.ne.1)
c	{
                  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) 
                           rmod(inum,itres,idelt) =
     *                          rmod(inum,it-2,idelt) 
                        enddo               
                     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) 
                           rmod(inum,itres,idelt) =
     *                          rmod(inum,it-2,idelt) 
                        enddo               
                     enddo
                     goto 888
                  endif
                  iter = iter + 1
c	}
                  enddo
c }
               endif

               xre = ygot(1)
               taure = ygot(2)

               in_bound=in_bound_intern .and. 
     *                  ((xre.ge.xminco).and.(xre.le.xlength).and.
     *                   (taure.ge.tauminco).and.(taure.le.taulength)) 

               if(in_bound) then
c              {

               vel1 = velocitn(xre,taure)
               sig1 = sigmatn(xre,taure)
               
               if(check .ne. 1) then
                     do inum=1,2
C                  write (idump,*) 'ygot(inum)',ygot(inum)
                        rtraj(inum,it,idelt) = ygot(inum)
                        rpar(inum,it,idelt) = ygot(inum+2)
                     enddo
                     rmod(1,it,idelt) = vel1
                     rmod(2,it,idelt) = sig1
               else

                  px=ygot(3)
                  ptau=ygot(4)

                  err =  1.- 
     *       (4.*(ptau**2) + (vel1**2)*(px+sig1*ptau)**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)
                  taure = ygot(2)
                  vel1 = velocitn(xre,taure)
                  sig1 = sigmatn(xre,taure)
                  toaux(1,idelt,it) = xre
                  toaux(2,idelt,it) = taure
                  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

c     write to file

      if(zmap .ne. 1) then
            ibuf = 0
            do idelt=1,ndelt
               do it=1,maxt
                  do inum=1,2
                     ibuf = ibuf+1
                     bufout(ibuf) = rtraj(inum,it,idelt)
                  enddo
               enddo
             enddo
       else
            ibuf = 0
            do idelt=1,ndelt
               do it=1,maxt
                  tauray=rtraj(2,it,idelt)
                  xray=rtraj(1,it,idelt)

                  tauray=(tauray-tauminco)/d1+1.
                  xray=(xray-xminco)/d2+1.

                  itauray=int(tauray)
                  ixray=int(xray)
                  if((ixray   .ge. 1).and.(ixray   .lt. n2) 
     *            .and. (itauray .ge. 1).and.(itauray .lt. n1)) then
                        wtau=tauray-itauray
                        wx=xray-ixray
                        ztaumap=
     *         ztau(itauray   + (ixray-1)*n1)*(1.-wtau)*(1.-wx)+
     *         ztau(itauray+1 + (ixray-1)*n1)*(wtau   )*(1.-wx)+
     *         ztau(itauray   + (ixray  )*n1)*(1.-wtau)*(wx   )+
     *         ztau(itauray+1 + (ixray  )*n1)*(wtau   )*(wx   )
                        ibuf = ibuf+1
                        bufout(ibuf) = rtraj(1,it,idelt)
                        ibuf = ibuf+1
                        bufout(ibuf) = ztaumap
                   else
                         ibuf = ibuf+1
                         bufout(ibuf) = bufout(ibuf-2)
                         ibuf = ibuf+1
                         bufout(ibuf) = bufout(ibuf-2)
                   end if
                enddo
             enddo
       end if

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

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

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

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

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

      if (ampyes.eq.1) then

         ibuf = 0
            do it=1,maxt
               do idelt = 1,ndelt
                  do inum=1,4
                     ibuf = ibuf+1
                     bufout(ibuf) = toaux(inum,idelt,it)
                  enddo
               enddo
            enddo

c     write to auxiliary otput file

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

      endif

c}
      enddo
         
      end


