      program MatchRay
c-----------------------------------------------
c
c     MatchRay 
c
c-----------------------------------------------------------------
c     Purpose: Compute ray parameters along rays in a model
c--------------------------------------------------------------------------
c     Input file: 
c
c------------------------------------------------------------------------
c     Output file: Gives ray paths that can be plotted with Graph.
c
c------------------------------------------------------------------------
c%     
      implicit none

      include '../Src/raypar.inc'

      logical ray_in_tau,found,found_multi,end_ray,ray_far_ref
      integer esize,esize4,it,idelt,iref,one,count_match
      integer fetch,getch,auxpar,sreed,sseek,srite,auxin,tempi

      integer maxt,ndelt,nsou,nx_ref,n_ref,isou,ix_ref,n_mid,n_h
      integer nby_ray,nby_seek,nby_sou,nby_mod
      integer idelt_ref,idelt_inc,idelt_ref_min,idelt_inc_min
      integer isou_ref,isou_inc,imid,ih
      integer nt_inc,nt_ref,zmap,n_zmap
      integer want_multi,cross_time,cross_only,cross_dir,bad
      integer it_min,ix_ref_min,ix_ref_max,dipfd,middle
      
      real xray1,xray2,vray1,vray2,xref1,xref2,vref1,vref2,vref0
      real wxref,xref,vref,pxref,pvref,velref,sigref
      real xa1,xa2,xb1,xb2,pa1,pa2,pb1,pb2,va1,vb1
      real dx_ref,idx_ref,tmatch_ref,tmatch_inc,tmatch_tot
      real mid_0,h_0,d_mid,d_h,sou_inc,sou_ref
      real wt,w_inc,w_ref
      real sin_cross,cos_cross,tandip_cross
      real tmin_ref,vmin_ref,vmax_ref

      real delta,delt,delmin
      real x0_ref,dip_ang
      real xsou,zsou,tausou,dx_sou,dz_sou,dtau_sou
      double precision depsi,sepsi,det,a,b,a_int,b_int
      double precision Dxa,Dva,Dpa,Dxb,Dvb,Dpb
      double precision Dxab,Dvab,Dpab

      logical crossed(mangle,msou)

      real rtraj(2,mtim,mangle),rpar(2,mtim,mangle)
      real rmod(2,mtim,mangle)
      real rtraj_inc(2,mtim,2),rpar_inc(2,mtim,2)
      real rmod_inc(2,mtim,2)
      real rtraj_ref(2,mtim,2),rpar_ref(2,mtim,2)
      real rmod_ref(2,mtim,2)
      real rtraj_zmap_inc(2,mtim,2),rtraj_zmap_ref(2,mtim,2)
      real rtraj_match(2,2*mtim),rpar_match(2,2)
      real rmod_match(2,2*mtim)

      real refl(2,mnx_ref)
      real cosrefl(mnx_ref),sinrefl(mnx_ref),tandip(mnx_ref)

      real tcross(mangle,msou)
      real xvcross(2,mangle,msou)
      real pcross_inc(2,mangle,msou)
      real dipcross(mangle,msou)
      real pcross_ref(2,mangle,msou)
      real velcross_inc(mangle,msou)
      real velcross_ref(mangle,msou)
      real sigcross_inc(mangle,msou)
      real sigcross_ref(mangle,msou)
      real location(3)


      call initpar()

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

c set constants
      depsi= 1.11D-15
      sepsi= 1.0d-4

c     fetch dimensions and esize.

      if (fetch('n1','d',maxt).eq.0) 
     *    call erexit('n1 must be supplied')
      if (fetch('n2','d',ndelt).eq.0) 
     *    call erexit('n2 must be supplied')
      if (fetch('d1','f',delt).eq.0) 
     *    call erexit('d1 must be supplied')
      if (fetch('d2','f',delta).eq.0) 
     *    call erexit('d2 must be supplied')
      if (fetch('o2','f',delmin).eq.0) 
     *    call erexit('o2 must be supplied')

      if (getch('middle','d',middle).eq.0) then
         middle=0
			end if

      if (fetch('nsou','d',nsou).eq.0)  then
			if(fetch('n3','d',nsou).eq.0) 
     *   call erexit('nsou must be supplied')
			else
      	ndelt=ndelt/nsou
			end if

      if (fetch('esize','d',esize).eq.0) 
     *   call erexit('esize must be supplied')

c get parametes from parameter file
      if (getch('tmin_ref','f',tmin_ref).le.0) tmin_ref=0.
      call putch('tmin_ref','f',tmin_ref)
      it_min=int(tmin_ref/delt)
      write(0,*)'it_min=',it_min

      if (getch('zmap','d',zmap).le.0) zmap=0
      call putch('zmap','d',zmap)
      if(zmap .eq. 1) then
       if (auxpar('n1','d',n_zmap,'Rayzmap').eq.0) 
     *   call erexit('n1 must be supplied in Rayzmap')
      end if
      if (getch('want_multi','d',want_multi).le.0) want_multi=0
      call putch('want_multi','d',want_multi)

      if (getch('cross_time','d',cross_time).le.0) cross_time=0
      call putch('cross_time','d',cross_time)
      if (getch('cross_only','d',cross_only).le.0) cross_only=0
      call putch('cross_only','d',cross_only)
      if (getch('cross_dir','d',cross_dir).le.0) cross_dir=1
      call putch('cross_dir','d',cross_dir)

      if(cross_only .ne. 1) then 
c     {
      if (fetch('xsou','f',xsou).eq.0) 
     *    call erexit('xsou must be supplied')
      if (fetch('dx_sou','f',dx_sou).eq.0) 
     *    call erexit('dx_sou must be supplied')
      if (fetch('dz_sou','f',dz_sou).eq.0) then
        if (fetch('dtau_sou','f',dtau_sou).eq.0) then 
					if(fetch('output_tau','d',tempi).eq.0) then
        	 call erexit('Either dz_sou or dtau_sou or output_tau must be supplied')
					else
						if(tempi.eq.1) ray_in_tau=.true.
					end if
				else
        	ray_in_tau=.true.
				end if
			end if
c        if (fetch('tausou','f',tausou).eq.0) 
c     *      call erexit('tausou must be supplied')
c      else
c        if (fetch('zsou','f',zsou).eq.0) 
c     *      call erexit('zsou must be supplied')
c        ray_in_tau=.false.
c      end if
      write(0,*)'ray_in_tau=',ray_in_tau

      if (getch('n_mid','d',n_mid).eq.0)  n_mid=1
      if (getch('mid_0','f',mid_0).eq.0)  mid_0=xsou
      if (getch('d_mid','f',d_mid).eq.0)  d_mid=dx_sou
      if (getch('n_h','d',n_h).eq.0)  n_h=1
      if (getch('h_0','f',h_0).eq.0)  h_0=0.
      if (getch('d_h','f',d_h).eq.0)  d_h=dx_sou

      if ((mid_0+(n_mid-1)*d_mid + h_0+(n_h-1)*d_h) 
     *       .gt. (xsou+(nsou-1)*dx_sou))
     *   call erexit('out on the far side')
      if ((mid_0 - h_0+(n_h-1)*d_h) .lt. xsou)
     *   call erexit('out on the near side')

c     }
      end if

C  get information on reflector

      if (auxpar('n1','d',nx_ref,'Refl').eq.0) 
     *   call erexit('n1 must be supplied in Refl')
      if (auxpar('n2','d',n_ref,'Refl').eq.0) 
     *   call erexit('n2 must be supplied in Refl')
C      if (auxpar('o1','f',x0_ref,'Refl').eq.0) 
C     *   call erexit('o1 must be supplied in Refl')
C      if (auxpar('d1','f',dx_ref,'Refl').eq.0) 
C     *   call erexit('d1 must be supplied in Refl')


c check parameters
      if (ndelt .gt. mangle) call erexit('ndelt too large')
      if (maxt .gt. mtim) call erexit('maxt too large')
      if (nsou .gt. msou) call erexit('nsou too large')
      if (nx_ref .gt. mnx_ref) call erexit('nx_ref too large')

      if(cross_time .eq. 1) then
        if(nx_ref .ne. nsou*ndelt) then
					write(0,*) "nx_ref=",nx_ref
					write(0,*) "nsou=",nsou
					write(0,*) "ndelt=",ndelt
          call erexit('nx_ref != nsou*ndelt')
        end if
      end if


C      if (mod(mid_0-xsou,dx_sou) .ne. 0) 
C     *   call erexit('mid_0-xsou not multiple of dx_sou')
C      if (mod(d_mid,dx_sou) .ne. 0) 
C     *   call erexit('d_mid not multiple of dx_sou')
C      if (mod(d_h,dx_sou) .ne. 0) 
C     *   call erexit('d_h not multiple of dx_sou')

      esize4=4
      one=1

      call putch('esize','d',esize)
      call putch('n1','d',2*maxt)
      call putch('o1','f',0.)
      call putch('d1','f',delt)
			if(want_multi .ne. 1) then
      call putch('n2','d',n_h)
      call putch('n3','d',n_mid)
      call putch('n4','d',n_ref)
      call putch('o2','f',h_0)
      call putch('o3','f',mid_0)
      call putch('o4','f',1.)
      call putch('d2','f',d_h)
      call putch('d3','f',d_mid)
      call putch('d4','f',1.)
			end if
      call putch('n_mid','d',n_mid)
      call putch('mid_0','f',mid_0)
      call putch('d_mid','f',d_mid)
      call putch('n_h','d',n_h)
      call putch('h_0','f',h_0)
      call putch('d_h','f',d_h)
      call putch('n_ref','d',n_ref)
      call putch('ref_0','f',1.)
      call putch('d_ref','f',1.)

	
      if(cross_only .ne. 1) then 
c     {
        call auxputch('esize','d',esize,'Raypar_match')
        call auxputch('esize','d',esize4,'Time_match')
        call auxputch('esize','d',esize4,'Time_ref')
        call auxputch('esize','d',esize,'Raymod_match')

        call auxputch('n1','d',2,'Raypar_match')
        call auxputch('n1','d',2*maxt,'Raymod_match')
        call auxputch('o1','f',0.,'Raypar_match')
        call auxputch('o1','f',0.,'Raymod_match')
        call auxputch('d1','f',delt,'Raypar_match')
        call auxputch('d1','f',delt,'Raymod_match')

			if(want_multi .ne. 1) then
        call auxputch('n1','d',n_h,'Time_match')
        call auxputch('n1','d',n_h,'Time_ref')
        call auxputch('n2','d',n_h,'Raymod_match')
        call auxputch('n2','d',n_h,'Raypar_match')

        call auxputch('n2','d',n_mid,'Time_match')
        call auxputch('n2','d',n_mid,'Time_ref')
				call auxputch('n2','d',n_mid,'Pos_match')
        call auxputch('n3','d',n_mid,'Raymod_match')
        call auxputch('n3','d',n_mid,'Raypar_match')

        call auxputch('n3','d',n_ref,'Time_match')
        call auxputch('n3','d',n_ref,'Time_ref')
        call auxputch('n4','d',n_ref,'Raypar_match')
        call auxputch('n4','d',n_ref,'Raymod_match')

        call auxputch('o1','f',h_0,'Time_match')
        call auxputch('o1','f',h_0,'Time_ref')
        call auxputch('o2','f',h_0,'Raymod_match')
        call auxputch('o2','f',h_0,'Raypar_match')

        call auxputch('o2','f',mid_0,'Time_match')
        call auxputch('o2','f',mid_0,'Time_ref')
        call auxputch('o3','f',mid_0,'Raypar_match')
        call auxputch('o3','f',mid_0,'Raymod_match')

        call auxputch('o4','f',one,'Raypar_match')
        call auxputch('o4','f',one,'Raymod_match')
        call auxputch('o3','f',one,'Time_match')
        call auxputch('o3','f',one,'Time_ref')

        call auxputch('d1','f',d_h,'Time_match')
        call auxputch('d1','f',d_h,'Time_ref')
        call auxputch('d2','f',d_h,'Raypar_match')
        call auxputch('d2','f',d_h,'Raymod_match')

        call auxputch('d2','f',d_mid,'Time_match')
        call auxputch('d2','f',d_mid,'Time_ref')
        call auxputch('d3','f',d_mid,'Raypar_match')
        call auxputch('d3','f',d_mid,'Raymod_match')

        call auxputch('d3','f',1.,'Time_match')
        call auxputch('d3','f',1.,'Time_ref')
        call auxputch('d4','f',1.,'Raypar_match')
        call auxputch('d4','f',1.,'Raymod_match')
			else
					call auxputch('n1','d',3,'location')
       end if 



c     }
      else
c     {
        call auxputch('esize','d',esize4,'Time_cross')
        call auxputch('esize','d',esize,'Xv_cross')
        call auxputch('esize','d',esize,'P_cross')

        call auxputch('n1','d',nsou,'Time_cross')
        call auxputch('n2','d',n_ref,'Time_cross')
        call auxputch('n3','d',one,'Time_cross')
        call auxputch('d1','f',1.,'Time_cross')
        call auxputch('d2','f',1.,'Time_cross')
        call auxputch('d3','f',1.,'Time_cross')
        call auxputch('o1','f',0.,'Time_cross')
        call auxputch('o2','f',0.,'Time_cross')
        call auxputch('o3','f',0.,'Time_cross')

        call auxputch('n1','d',nsou,'Xv_cross')
        call auxputch('n2','d',n_ref,'Xv_cross')
        call auxputch('n3','d',one,'Xv_cross')
        call auxputch('d1','f',1.,'Xv_cross')
        call auxputch('d2','f',1.,'Xv_cross')
        call auxputch('d3','f',1.,'Xv_cross')
        call auxputch('o1','f',0.,'Xv_cross')
        call auxputch('o2','f',0.,'Xv_cross')
        call auxputch('o3','f',0.,'Xv_cross')

        call auxputch('n1','d',nsou,'P_cross')
        call auxputch('n2','d',n_ref,'P_cross')
        call auxputch('n3','d',one,'P_cross')
        call auxputch('d1','f',1.,'P_cross')
        call auxputch('d2','f',1.,'P_cross')
        call auxputch('d3','f',1.,'P_cross')
        call auxputch('o1','f',0.,'P_cross')
        call auxputch('o2','f',0.,'P_cross')
        call auxputch('o3','f',0.,'P_cross')

        call auxputch('n1','d',nsou,'Dip_cross')
        call auxputch('n2','d',n_ref,'Dip_cross')
        call auxputch('n3','d',one,'Dip_cross')
        call auxputch('d1','f',1.,'Dip_cross')
        call auxputch('d2','f',1.,'Dip_cross')
        call auxputch('d3','f',1.,'Dip_cross')
        call auxputch('o1','f',0.,'Dip_cross')
        call auxputch('o2','f',0.,'Dip_cross')
        call auxputch('o3','f',0.,'Dip_cross')

				call auxputch('esize','d',8,'Pos_match')
				call auxputch('n1','d',n_h,'Pos_match')
				call auxputch('n3','d',n_ref,'Pos_match')
        call auxputch('o1','f',h_0,'Pos_match')
        call auxputch('o3','f',mid_0,'Pos_match')
        call auxputch('d1','f',d_h,'Pos_match')
        call auxputch('d2','f',d_mid,'Pos_match')
        call auxputch('d3','f',1.,'Pos_match')
      end if
c     }

c     close history file     

        bad=0
      count_match=0
C      call hclose()
C      call auxclose('Raypar_match')


c     precompute information on end point of each ray
 
C      idx_ref=1./dx_ref
      nby_ray=maxt*esize
      nby_sou=nby_ray*ndelt

      dipfd=auxin('Dip_refl')


      do iref=1,n_ref
c     {
			  location(3)=iref

c read reflector information
        if (sreed('Refl',refl,nx_ref*8) .ne. nx_ref*8) then
          call erexit("Problems reading Refl file")
        end if

	write(0,*)'Matching rays for reflector #',iref

        vmin_ref=100000000.
        vmax_ref=-100000000.
        do ix_ref=1,nx_ref
c        {
          if(vmin_ref .gt. refl(2,ix_ref)) then
            vmin_ref=refl(2,ix_ref)
          end if
          if(vmax_ref .lt. refl(2,ix_ref)) then
            vmax_ref=refl(2,ix_ref)
          end if
c        }
        end do 
        write(0,*)'vmin_ref=',vmin_ref
        write(0,*)'vmax_ref=',vmax_ref

        if(dipfd .gt. 0) then
c        {
          write(0,*)'read Dip_refl file'
          write(0,*)'dipfd ',dipfd
          if (sreed('Dip_refl',tandip,nx_ref*4) .ne. nx_ref*4) then
            call erexit("Problems reading Dip_refl file")
          end if
c        }
        else
c        {
          write(0,*)'NOT read Dip_refl file'

c comute cos and sin of 2*dip angle
          do ix_ref=2,nx_ref-1
c          {
            tandip(ix_ref)=
     *        (refl(2,ix_ref+1)-refl(2,ix_ref-1))/
     *        (refl(1,ix_ref+1)-refl(1,ix_ref-1))
          end do 
c          }

          tandip(1) =
     *        (refl(2,2)-refl(2,1))/
     *        (refl(1,2)-refl(1,1))

          tandip(nx_ref) =
     *        (refl(2,nx_ref)-refl(2,nx_ref-1))/
     *        (refl(1,nx_ref)-refl(1,nx_ref-1))
c        }
        end if
        do ix_ref=1,nx_ref
c        {
          dip_ang=2*atan(tandip(ix_ref))
          sinrefl(ix_ref)=sin(dip_ang)
          cosrefl(ix_ref)=cos(dip_ang)
c        }
        end do 

        	call snap("Tandip.h",nx_ref,1,tandip)
       
			if(middle .ne. 1) then
c       {
        do isou=1,nsou
c      {
        do idelt=1,ndelt
c       {
          if(iref .ne. 1) then
            nby_seek=nby_sou*(isou-1) + nby_ray*(idelt-1)
            if (sseek('in',nby_seek,0).ne.nby_seek) then
        call erexit("Problems seeking Input file")
            end if
            if (sseek('Raypar',nby_seek,0).ne.nby_seek) then
        call erexit("Problems seeking Raypar file")
            end if
            if (sseek('Raymod',nby_seek,0).ne.nby_seek) then
        call erexit("Problems seeking Raymod file")
            end if
          end if

          if (sreed('in',rtraj(1,1,idelt),nby_ray).ne.nby_ray) then
      call erexit("Problems reading Input file")
          end if

          if (sreed('Raypar',rpar(1,1,idelt),nby_ray).ne.nby_ray) then
      call erexit("Problems reading Raypar file")
          end if

          if (sreed('Raymod',rmod(1,1,idelt),nby_ray).ne.nby_ray) then
      call erexit("Problems reading Raymod file")
          end if

c				if(isou.eq.1 .and. iref.ne.1) then
c        if (sseek('in',0,0).ne.0) then
c        	call erexit("Problems seeking Input file")
c        end if
c        if (sseek('Raypar',0,0).ne.0) then
c        	call erexit("Problems seeking Raypar file")
c        end if
c        if (sseek('Raymod',0,0).ne.0) then
c           call erexit("Problems seeking Raymod file")
c        end if
c        end if
c
c        if (sreed('in',rtraj(1,1,1),nby_sou).ne.nby_sou) then
c      		call erexit("Problems reading Input file")
c        end if
c
c        if (sreed('Raypar',rpar(1,1,1),nby_sou).ne.nby_sou) then
c          call erexit("Problems reading Raypar file")
c        end if
c
cc        if (sreed('Raymod',rmod(1,1,1),nby_sou).ne.nby_sou) then
c          call erexit("Problems reading Raymod file")
c        end if
c
c      do idelt=1,ndelt


          end_ray=.false.
          crossed(idelt,isou) = .false.
          it=it_min
          do while ((.not.crossed(idelt,isou)).and.
     *              (it.lt.maxt-1).and. 
     *              (.not. end_ray)) 
c         {
            it=it+1
            xray1=rtraj(1,it  ,idelt)
            xray2=rtraj(1,it+1,idelt)
            vray1=rtraj(2,it  ,idelt)
            vray2=rtraj(2,it+1,idelt)

            if((vray1 .eq. vray2).and.(xray1 .eq. xray1)) then
              end_ray=.true.
            end if

            if(cross_time .eq. 1) then
              vray1=(it-1.)*delt
              vray2=(it   )*delt
              ix_ref_min=idelt+(isou-1)*ndelt-1
              ix_ref_max=ix_ref_min+1
            else
              ix_ref_min=0
              ix_ref_max=nx_ref-1
            end if

            Dxa=xray2-xray1
            Dva=vray2-vray1

            if (cross_dir .eq. 1) then
              if((vray1 .lt. vmin_ref).and.(vray2 .lt. vmin_ref)) then
                ray_far_ref=.true.
              else
                ray_far_ref=.false.
              end if
            else
              if((vray1 .gt. vmax_ref).and.(vray2 .gt. vmax_ref)) then
                ray_far_ref=.true.
              else
                ray_far_ref=.false.
              end if
            end if

          ix_ref=ix_ref_min
          do while ((.not.crossed(idelt,isou)).and.
     *              (ix_ref.lt.ix_ref_max).and. 
     *              (.not. end_ray).and. 
     *              (.not. ray_far_ref)) 
c         {
            ix_ref=ix_ref+1

c trick to keep the same code for both cases (more expensive)
            if(cross_time .eq. 1) then
              xref1=xray1-sign(sepsi,Dxa)
              xref2=xray2+sign(sepsi,Dxa)
              vref1=refl(2,ix_ref)
              vref2=vref1
            else
              xref1=refl(1,ix_ref  )
              xref2=refl(1,ix_ref+1)
              vref1=refl(2,ix_ref  )
              vref2=refl(2,ix_ref+1)
            end if

            Dxb=xref2-xref1
            Dvb=vref2-vref1

            Dxab=xref1-xray1
            Dvab=vref1-vray1

            det=Dva*Dxb-Dxa*Dvb

            if(dabs(det) .le. depsi) then
c               {
              if((Dxa .eq. Dxb) .and. (Dva .eq. Dvb) .and.
     *           (Dxab .eq. Dvab)) then
                a=0.d0
                b=0.d0
              call seperr('crossing ray is coincident with reflector')
              else
                write(0,*)'POTENTIAL ERROR'
                write(0,*)'ray is parallel to reflctor'
                write(0,*)'det=',det
                write(0,*)'it=',it
                write(0,*)'ix_ref=',ix_ref
                write(0,*)'isou=',isou
                write(0,*)'idelt=',idelt
                write(0,*)'xray1=',xray1
                write(0,*)'xray2=',xray2
                write(0,*)'xref1=',xref1
                write(0,*)'xref2=',xref2
                write(0,*)'vray1=',vray1
                write(0,*)'vray2=',vray2
                write(0,*)'vref1=',vref1
                write(0,*)'vref2=',vref2
                write(0,*)'Dxab=',Dxab
                write(0,*)'Dvab=',Dvab
C                 call seperr('ray is parallel to reflctor')
              end if
c             }
            else
c             {
              a= (Dvab*Dxb - Dxab*Dvb)/det
              b=-(-Dvab*Dxa + Dxab*Dva)/det
              if((a .ge. -sepsi) .and. (a .le. 1.+sepsi) .and.
     *           (b .ge. -sepsi) .and. (b .le. 1.+sepsi)) then
c                  {
               crossed(idelt,isou) = .true.

               xvcross(1,idelt,isou)=xray1 + (a*Dxa)

               if(cross_time .eq. 1) then
                 xvcross(2,idelt,isou)=rtraj(2,it  ,idelt)*(1.-a) +  
     *                                 rtraj(2,it+1,idelt)*(   a) 
               else
                 xvcross(2,idelt,isou)=vray1 + (a*Dva)
               end if

               tcross(idelt,isou)=((it-1.)+a)*delt

C                write(0,*)'xvcross(1)=',xvcross(1,idelt,isou)
C                write(0,*)'xvcross(2)=',xvcross(2,idelt,isou)
C                write(0,*)'tcross(2)=',tcross(idelt,isou)

               pcross_inc(1,idelt,isou) =
     *             rpar(1,it  ,idelt)*(1.-a) +
     *             rpar(1,it+1,idelt)*(   a) 

               pcross_inc(2,idelt,isou) =
     *             rpar(2,it  ,idelt)*(1.-a) +
     *             rpar(2,it+1,idelt)*(   a) 

c							write(0,*) "pcross",idelt,iref,isou,it,a

               velcross_inc(idelt,isou) =
     *             rmod(1,it  ,idelt)*(1.-a) +
     *             rmod(1,it+1,idelt)*(   a) 

               sigcross_inc(idelt,isou) =
     *             rmod(2,it  ,idelt)*(1.-a) +
     *             rmod(2,it+1,idelt)*(   a) 

               dipcross(idelt,isou) =
     *             -tan(asin((pcross_inc  (1,idelt,isou) *
     *                       velcross_inc(  idelt,isou))))


               if(cross_only .ne. 1) then 
c                {
c note that they have the opposite sign that they should because 
c all rays are going down
                 if(.not. ray_in_tau) then
c                {
                 sin_cross=
     *             sinrefl(ix_ref  )*(1.-b) +
     *             sinrefl(ix_ref+1)*(   b) 
                 cos_cross=
     *             cosrefl(ix_ref  )*(1.-b) +
     *             cosrefl(ix_ref+1)*(   b) 

                 pcross_ref(1,idelt,isou) =
     *              pcross_inc(1,idelt,isou)*cos_cross +
     *              pcross_inc(2,idelt,isou)*sin_cross 

                 pcross_ref(2,idelt,isou) =
     *              pcross_inc(1,idelt,isou)*sin_cross -
     *              pcross_inc(2,idelt,isou)*cos_cross

c                }
                 else
c                {

                 if(dipfd .le. 0) then
c                  {

                   tandip_cross=
     *               tandip(ix_ref  )*(1.-b) +
     *               tandip(ix_ref+1)*(   b) 

C 	ATTENTION VALID ONLY FOR FLAT REFLECTORS
                    dip_ang=2*atan(.5*velcross_inc(idelt,isou) *
     *             (tandip_cross-sigcross_inc(idelt,isou)))

                    sin_cross=sin(dip_ang)
                    cos_cross=cos(dip_ang)

c                  }
                 else
c                  {

c in this case the angles are already in depth and not in tau 

                   sin_cross=
     *               sinrefl(ix_ref  )*(1.-b) +
     *               sinrefl(ix_ref+1)*(   b) 
                   cos_cross=
     *               cosrefl(ix_ref  )*(1.-b) +
     *               cosrefl(ix_ref+1)*(   b) 
c                  }
                 end if

                 pcross_ref(1,idelt,isou) =
     *       pcross_inc(1,idelt,isou)*
     *      (cos_cross -
     *       sin_cross*
     *       .5*sigcross_inc(idelt,isou)*velcross_inc(idelt,isou)) +
     *       pcross_inc(2,idelt,isou)*
     *     (
     *       cos_cross*
     *       2.*sigcross_inc(idelt,isou) + 
     *       sin_cross*
     *      (2./velcross_inc(idelt,isou) -
     *       .5*velcross_inc(idelt,isou)*sigcross_inc(idelt,isou)**2)
     *     )
               pcross_ref(2,idelt,isou) =
     *      +pcross_inc(1,idelt,isou)*
     *      (sin_cross*
     *       .5*velcross_inc(idelt,isou)) +
     *       pcross_inc(2,idelt,isou)*
     *     (
     *      -cos_cross +
     *       sin_cross*
     *       .5*sigcross_inc(idelt,isou)*velcross_inc(idelt,isou)
     *     )
c                }
                 end if
c                }
               end if
c                  }
                end if
c             }
          end if
c          }
           end do
c          }
           end do
        
           if (.not. crossed(idelt,isou)) then
             xvcross(1,idelt,isou)=rtraj(1,1,idelt)
             xvcross(2,idelt,isou)=rtraj(2,1,idelt)
             tcross(idelt,isou)=-1.
             if(cross_time .eq. 1) then
               write(0,*)'isou=',isou
               write(0,*)'idelt=',idelt
               call seperr('Ray never crossed time surface')
             end if
           end if


c       }
        end do
c      }
       end do
c      } Loop for middle
       end if
			do isou=1,nsou
      	if(middle .ne. 1) then
       	 idelt= srite("xvcross.H",xvcross(1,1,isou),ndelt*8)
       	 idelt= srite("tcross.H",tcross(1,isou),ndelt*4)
       	 idelt= srite("pcross_inc.H",pcross_inc(1,1,isou),ndelt*8)
       	 idelt= srite("pcross_ref.H",pcross_ref(1,1,isou),ndelt*8)
       	 idelt= srite("dipcross.H",dipcross(1,isou),ndelt*4)
       	 idelt= srite("velcross_inc.H",velcross_inc(1,isou),ndelt*4)
       	 idelt= srite("sigcross_inc.H",sigcross_inc(1,isou),ndelt*4)
       	 idelt= srite("sigcross_ref.H",sigcross_ref(1,isou),ndelt*4)
        else
       	 idelt= sreed("xvcross.H",xvcross(1,1,isou),ndelt*8)
       	 idelt= sreed("tcross.H",tcross(1,isou),ndelt*4)
       	 idelt= sreed("pcross_inc.H",tcross(1,isou),ndelt*4)
       	 idelt= sreed("pcross_ref.H",tcross(1,isou),ndelt*8)
       	 idelt= sreed("dipcross.H",tcross(1,isou),ndelt*4)
       	 idelt= sreed("velcross_inc.H",velcross_inc(1,isou),ndelt*4)
       	 idelt= sreed("sigcross_inc.H",sigcross_inc(1,isou),ndelt*4)
       	 idelt= sreed("sigcross_ref.H",sigcross_ref(1,isou),ndelt*4)
      	end if
      end do

       if(cross_only .eq. 1) then 
c        {
         do isou=1,nsou
c         {
            if (srite('Time_cross',tcross(1,isou),ndelt*4)
     *                              .ne.  ndelt*4) then
              call erexit("Problems writing Time_cross file")
            end if

            if (srite('Xv_cross',xvcross(1,1,isou),ndelt*8)
     *                              .ne. ndelt*8) then
              call erexit("Problems writing Xv_cross file")
            end if
            if (srite('P_cross',pcross_inc(1,1,isou),ndelt*8)
     *                              .ne. ndelt*8) then
              call erexit("Problems writing P_cross file")
            end if
            if (srite('Dip_cross',dipcross(1,isou),ndelt*4)
     *                              .ne. ndelt*4) then
              call erexit("Problems writing Dip_cross file")
            end if
c         }
         end do
c       }
        else
c         {

	write(0,*)'BEGINNING matching rays'

C	call snap("Xcross",mangle,nsou,xcross)
C	call snap("Tcross",mangle,nsou,tcross)

C	call snap("Vcross",mangle,nsou,vcross)
C	call snap("Pxcross_inc",mangle,nsou,pxcross_inc)
C	call snap("Pvcross_inc",mangle,nsou,pvcross_inc)
C	call snap("Pxcross_ref",mangle,nsou,pxcross_ref)
C	call snap("Pvcross_ref",mangle,nsou,pvcross_ref)

      found_multi=.false.
      do imid=1,n_mid
c     {
			  location(2)=imid
        do ih=1,n_h
c       {
			  	location(1)=ih

          sou_inc=(mid_0 + (imid-1)*d_mid) - (h_0 + (ih-1)*d_h) - xsou
          sou_ref=(mid_0 + (imid-1)*d_mid) + (h_0 + (ih-1)*d_h) - xsou
C          if(abs(mod(sou_inc,dx_sou)) .gt. dx_sou*.001) then
C            write(0,*)'mod(sou_inc,dx_sou)=',mod(sou_inc,dx_sou)
C            write(0,*)'sou_inc=',sou_inc
C            call seperr('sou_inc not on grid')
C          else
C            write(0,*)'sou_inc=',sou_inc
            isou_inc= nint(sou_inc/dx_sou) + 1
C          end if
C          if(abs(mod(sou_ref,dx_sou)) .gt. dx_sou*.001) then
C            write(0,*)'mod(sou_ref,dx_sou)=',mod(sou_ref,dx_sou)
C            write(0,*)'sou_ref=',sou_ref
C            call seperr('sou_ref not on grid')
C          else
C            write(0,*)'sou_ref=',sou_ref
            isou_ref= nint(sou_ref/dx_sou) + 1
C          end if

C          write(0,*)'imid=',imid
C          write(0,*)'ih=',ih
C          write(0,*)'isou_inc=',isou_inc
C          write(0,*)'isou_ref=',isou_ref
  
          idelt_inc=0
          found=.false.
C          do while ((idelt_inc .lt. ndelt-1).and.(.not. found))
          do while ((idelt_inc .lt. ndelt-1))
c         {
            idelt_inc=idelt_inc+1
            if(crossed(idelt_inc  ,isou_inc) .and. 
     *         crossed(idelt_inc+1,isou_inc)) then
c           {
               idelt_ref=0
C               do while ((idelt_ref .lt. ndelt-1).and.(.not. found))
               do while ((idelt_ref .lt. ndelt-1))
c              {
                 idelt_ref=idelt_ref+1
                 if(crossed(idelt_ref  ,isou_ref) .and. 
     *              crossed(idelt_ref+1,isou_ref)) then
c                {

                   xa1=xvcross(1,idelt_inc  ,isou_inc)
                   xa2=xvcross(1,idelt_inc+1,isou_inc)
                   xb1=xvcross(1,idelt_ref  ,isou_ref)
                   xb2=xvcross(1,idelt_ref+1,isou_ref)

                   pa1=-pcross_ref(1,idelt_inc  ,isou_inc)
                   pa2=-pcross_ref(1,idelt_inc+1,isou_inc)
                   pb1=pcross_inc(1,idelt_ref  ,isou_ref)
                   pb2=pcross_inc(1,idelt_ref+1,isou_ref)




                   Dxa=xa2-xa1
                   Dpa=pa2-pa1
                   Dxb=xb2-xb1
                   Dpb=pb2-pb1

                   Dxab=xb1-xa1
                   Dpab=pb1-pa1

                   det=Dpa*Dxb-Dxa*Dpb

                   if(dabs(det) .le. depsi) then
c                  {
                     if((Dxa .eq. Dxb) .and. (Dpa .eq. Dpb) .and.
     *                  (Dxab .eq. Dpab)) then
                       a=0.d0
                       b=0.d0
                       found=.true.
                       call seperr('matching rays are same')
                     else
                      
                       det=depsi
                       write(0,*)'ERROR'
                       write(0,*)'idelt_inc=',idelt_inc
                       write(0,*)'idelt_ref=',idelt_ref
                       write(0,*)'xa1=',xa1
                       write(0,*)'xa2=',xa2
                       write(0,*)'xb1=',xb1
                       write(0,*)'xb2=',xb2
                       write(0,*)'Dxab=',Dxab
                       write(0,*)'Dpab=',Dpab
                       call seperr('matching rays are parrallel')
                     end if
c                  }
                   else
c                  {
                     a= (Dpab*Dxb - Dxab*Dpb)/det
                     b=-(-Dpab*Dxa + Dxab*Dpa)/det
                     if((a .ge. -sepsi) .and. (a .le. 1.+sepsi) .and.
     *                  (b .ge. -sepsi) .and. (b .le. 1.+sepsi)) then
c                    {

C	 	write(0,*)'GOOD det=',det
C	 	write(0,*)'GOOD a=',a
C	 	write(0,*)'GOOD b=',b
C	 	write(0,*)'idelt_inc=',idelt_inc
C	 	write(0,*)'idelt_ref=',idelt_ref
C                       write(0,*)'xa1=',xa1
C                       write(0,*)'xa2=',xa2
C                       write(0,*)'xb1=',xb1
C                       write(0,*)'xb2=',xb2
C                       write(0,*)'pa1=',pa1
C                       write(0,*)'pa2=',pa2
C                       write(0,*)'pb1=',pb1
C                       write(0,*)'pb2=',pb2
C     			write(0,*) "ANGLE CROSSS?",pcross_ref(1,idelt_inc  ,isou_inc)
C     *								,pcross_ref(1,idelt_inc +1,isou_inc)
C     *									,pcross_inc(1,idelt_ref ,isou_ref)
C     *									,pcross_inc(1,idelt_ref +1,isou_ref)
C
         if(found .and. (want_multi .ne. 1)) then
c        {
            found_multi=.true.
c        }
         else
c        {
         if(found) found_multi=.true.
              

C					write(0,*) "xa1 xa2 dxa",xa1,xa2,dxa
C					write(0,*) "pa1 pa2 dpa",pa1,pa2,dpa
C					write(0,*) "xb1 xb2 dxb",xa2,xa2,dxa
C					write(0,*) "pb1 pb2 dpb",pb1,pb2,dpb
C					write(0,*) "Dxab Dpab det",Dxab,Dpab,det
         found=.true.
         count_match=count_match+1

         xref=xa1+ (a*Dxa)
          vref=
     *       xvcross(2,idelt_inc  ,isou_inc)*(1.-a) +
     *       xvcross(2,idelt_inc+1,isou_inc)*(a   ) 
          pxref=
     *       pcross_inc(1,idelt_inc  ,isou_inc)*(1.-a) +
     *       pcross_inc(1,idelt_inc+1,isou_inc)*(a   ) 
          pvref=
     *       pcross_inc(2,idelt_inc  ,isou_inc)*(1.-a) +
     *       pcross_inc(2,idelt_inc+1,isou_inc)*(a   ) 
          velref=
     *       velcross_inc(idelt_inc  ,isou_inc)*(1.-a) +
     *       velcross_inc(idelt_inc+1,isou_inc)*(a   ) 
          sigref=
     *       sigcross_inc(idelt_inc  ,isou_inc)*(1.-a) +
     *       sigcross_inc(idelt_inc+1,isou_inc)*(a   ) 

          tmatch_inc=
     *       tcross(idelt_inc  ,isou_inc)*(1.-a) +
     *       tcross(idelt_inc+1,isou_inc)*(a   ) 

          tmatch_ref=
     *       tcross(idelt_ref  ,isou_ref)*(1.-b) +
     *       tcross(idelt_ref+1,isou_ref)*(b   )  

          tmatch_tot=tmatch_inc+tmatch_ref

C	write(0,*)'xref=',xref
C	write(0,*)'vref=',vref
C	write(0,*)'pxref=',pxref
C	write(0,*)'pvref=',pvref
C	write(0,*)'velref=',velref
C	write(0,*)'tmatch_tot=',tmatch_tot
C	write(0,*)'tcross(idelt_inc)=',
C     *   tcross(idelt_inc,isou_inc)
C	write(0,*)'tcross(idelt_inc+1)=',
C     *   tcross(idelt_inc+1,isou_inc)
C	write(0,*)'tcross(idelt_ref)=',
C     *   tcross(idelt_ref,isou_ref)
C	write(0,*)'tcross(idelt_ref+1)=',
C     *   tcross(idelt_ref+1,isou_ref)
C
C	write(0,*)'xcross(idelt_inc)=',
C     *   xcross(idelt_inc,isou_inc)
C	write(0,*)'xcross(idelt_inc+1)=',
C     *   xcross(idelt_inc+1,isou_inc)
C	write(0,*)'xcross(idelt_ref)=',
C     *   xcross(idelt_ref,isou_ref)
C	write(0,*)'xcross(idelt_ref+1)=',
C     *   xcross(idelt_ref+1,isou_ref)
C
C	write(0,*)'pxcross_ref(idelt_inc)=',
C     *   pxcross_ref(idelt_inc,isou_inc)
C	write(0,*)'pxcross_ref(idelt_inc+1)=',
C     *   pxcross_ref(idelt_inc+1,isou_inc)
C	write(0,*)'pxcross_inc(idelt_inc)=',
C     *   pxcross_inc(idelt_inc,isou_inc)
C	write(0,*)'pxcross_inc(idelt_inc+1)=',
C     *   pxcross_inc(idelt_inc+1,isou_inc)

         nby_seek=nby_sou*(isou_inc-1) + nby_ray*(idelt_inc-1)

         if (sseek('in',nby_seek,0).ne.nby_seek) then
	   call erexit("Problems seeking Input file")
         end if
         if (sreed('in',rtraj_inc(1,1,1),nby_ray).ne.nby_ray) then
	   call erexit("Problems reading Input file")
         end if
         if (sreed('in',rtraj_inc(1,1,2),nby_ray).ne.nby_ray) then
	   call erexit("Problems reading Input file")
         end if

         if (sseek('Raypar',nby_seek,0).ne.nby_seek) then
	   call erexit("Problems seeking Raypar file")
         end if
         if (sreed('Raypar',rpar_inc(1,1,1),nby_ray).ne.nby_ray) then
	   call erexit("Problems reading Raypar file")
         end if
         if (sreed('Raypar',rpar_inc(1,1,2),nby_ray).ne.nby_ray) then
	   call erexit("Problems reading Raypar file")
         end if

         if (sseek('Raymod',nby_seek,0).ne.nby_seek) then
           call erexit("Problems seeking Raymod file")
         end if
         if (sreed('Raymod',rmod_inc(1,1,1),nby_ray).ne.nby_ray) then
           call erexit("Problems reading Raymod file")
         end if
         if (sreed('Raymod',rmod_inc(1,1,2),nby_ray).ne.nby_ray) then
           call erexit("Problems reading Raymod file")
         end if

         if(zmap .eq. 1) then
           if (sseek('Rayzmap',nby_seek,0).ne.nby_seek) then
             call erexit("Problems seeking Rayzmap file")
           end if
           if (sreed('Rayzmap',rtraj_zmap_inc(1,1,1),nby_ray)
     *                                       .ne.nby_ray) then
             call erexit("Problems reading Rayzmap file")
           end if
           if (sreed('Rayzmap',rtraj_zmap_inc(1,1,2),nby_ray)
     *                                       .ne.nby_ray) then
             call erexit("Problems reading Rayzmap file")
           end if
         end if

        nby_seek=nby_sou*(isou_ref-1) + nby_ray*(idelt_ref-1)

         if (sseek('in',nby_seek,0).ne.nby_seek) then
	   call erexit("Problems seeking Input file")
         end if
         if (sreed('in',rtraj_ref(1,1,1),nby_ray).ne.nby_ray) then
	   call erexit("Problems reading Input file")
         end if
         if (sreed('in',rtraj_ref(1,1,2),nby_ray).ne.nby_ray) then
	   call erexit("Problems reading Input file")
         end if

         if (sseek('Raypar',nby_seek,0).ne.nby_seek) then
	   call erexit("Problems seeking Raypar file")
         end if
         if (sreed('Raypar',rpar_ref(1,1,1),nby_ray).ne.nby_ray) then
	   call erexit("Problems reading Raypar file")
         end if
         if (sreed('Raypar',rpar_ref(1,1,2),nby_ray).ne.nby_ray) then
	   call erexit("Problems reading Raypar file")
         end if

         if (sseek('Raymod',nby_seek,0).ne.nby_seek) then
           call erexit("Problems seeking Raymod file")
         end if
         if (sreed('Raymod',rmod_ref(1,1,1),nby_ray).ne.nby_ray) then
           call erexit("Problems reading Raymod file")
         end if
         if (sreed('Raymod',rmod_ref(1,1,2),nby_ray).ne.nby_ray) then
           call erexit("Problems reading Raymod file")
         end if

         if(zmap .eq. 1) then
           if (sseek('Rayzmap',nby_seek,0).ne.nby_seek) then
             call erexit("Problems seeking Rayzmap file")
           end if
           if (sreed('Rayzmap',rtraj_zmap_ref(1,1,1),nby_ray)
     *                                       .ne.nby_ray) then
             call erexit("Problems reading Rayzmap file")
           end if
           if (sreed('Rayzmap',rtraj_zmap_ref(1,1,2),nby_ray)
     *                                       .ne.nby_ray) then
             call erexit("Problems reading Rayzmap file")
           end if
         end if

         nt_inc=(tmatch_inc)/delt+1
         nt_ref=(tmatch_ref)/delt+1

         if((nt_inc+nt_ref +1) .gt. 2*maxt) then
           call seperr('maxt too small') 
         end if

         wt=tmatch_inc/delt-(nt_inc-1.)
C          wt=1.
         xa1= rtraj_inc(1,nt_inc  ,1)*(1.-wt) +
     *        rtraj_inc(1,nt_inc+1,1)*(   wt) 
         xb1= rtraj_inc(1,nt_inc  ,2)*(1.-wt) +
     *        rtraj_inc(1,nt_inc+1,2)*(   wt) 
         w_inc=(xref-xa1)/(xb1-xa1)

         if(zmap .eq. 1) then
           va1= rtraj_zmap_inc(2,nt_inc  ,1)*(1.-wt) +
     *        rtraj_zmap_inc(2,nt_inc+1,1)*(   wt) 
           vb1= rtraj_zmap_inc(2,nt_inc  ,2)*(1.-wt) +
     *        rtraj_zmap_inc(2,nt_inc+1,2)*(   wt) 
           vref=va1+w_inc*(vb1-va1)
         end if

C	write(0,*)'tmatch_inc',tmatch_inc
C	write(0,*)'wt',wt
C	write(0,*)'xa1',xa1
C	write(0,*)'xb1',xb1
C	write(0,*)'va1',pa1
C	write(0,*)'vb1',pb1
C	write(0,*)'xref',xref
C	write(0,*)'idelt_inc',idelt_inc
C	write(0,*)'idelt_ref',idelt_ref
C	write(0,*)'VREF',vref
C
         wt=tmatch_ref/delt-(nt_ref-1.)
         xa1= rtraj_ref(1,nt_ref  ,1)*(1.-wt) +
     *        rtraj_ref(1,nt_ref+1,1)*(   wt) 
         xb1= rtraj_ref(1,nt_ref  ,2)*(1.-wt) +
     *        rtraj_ref(1,nt_ref+1,2)*(   wt) 
         w_ref=(xref-xa1)/(xb1-xa1)

         if(zmap .eq. 1) then
           va1= rtraj_zmap_ref(2,nt_ref  ,1)*(1.-wt) +
     *        rtraj_zmap_ref(2,nt_ref+1,1)*(   wt) 
           vb1= rtraj_zmap_ref(2,nt_ref  ,2)*(1.-wt) +
     *        rtraj_zmap_ref(2,nt_ref+1,2)*(   wt) 
           vref=va1+w_ref*(vb1-va1)
         end if

C	write(0,*)'tmatch_ref',tmatch_ref
C	write(0,*)'wt',wt
C	write(0,*)'xa1',xa1
C	write(0,*)'xb1',xb1
C	write(0,*)'va1',pa1
C	write(0,*)'vb1',pb1
C	write(0,*)'xref',xref
C
C	write(0,*)'VREF',vref
C
C	write(0,*)'nt_inc=',nt_inc
C	write(0,*)'nt_ref=',nt_ref
C	write(0,*)'w_inc=',w_inc
C	write(0,*)'w_ref=',w_ref


         do it=1,nt_inc
c        {
             rtraj_match(1,it) =
     *       rtraj_inc(1,it,1)*(1.-w_inc) +
     *       rtraj_inc(1,it,2)*(   w_inc) 

             if(zmap .ne. 1) then
               rtraj_match(2,it) =
     *         rtraj_inc(2,it,1)*(1.-w_inc) +
     *         rtraj_inc(2,it,2)*(   w_inc) 
             else
               rtraj_match(2,it) =
     *         rtraj_zmap_inc(2,it,1)*(1.-w_inc) +
     *         rtraj_zmap_inc(2,it,2)*(   w_inc) 
             end if

C             rpar_match(1,it) =
C     *       rpar_inc(1,it,1)*(1.-w_inc) +
C     *       rpar_inc(1,it,2)*(   w_inc) 

CC             rpar_match(2,it) =
C     *       rpar_inc(2,it,1)*(1.-w_inc) +
C     *       rpar_inc(2,it,2)*(   w_inc) 

             rmod_match(1,it) =
     *       rmod_inc(1,it,1)*(1.-w_inc) +
     *       rmod_inc(1,it,2)*(   w_inc) 

             if(ray_in_tau) then
               rmod_match(2,it) =
     *         rmod_inc(2,it,1)*(1.-w_inc) +
     *         rmod_inc(2,it,2)*(   w_inc) 
             else
               rmod_match(2,it) = 0.
             end if

C	dt=sqrt(.25*((rtraj_match(2,it)-rtraj_match(2,it-1))-
C     *               rmod_match(2,it)*
C     *               (rtraj_match(1,it)-rtraj_match(1,it-1)))**2 +
C     *               (1./rmod_match(1,it)*
C     *               (rtraj_match(1,it)-rtraj_match(1,it-1)))**2 )
C           write(0,*)'dt=',dt
c        }
         end do

         rtraj_match(1,nt_inc+1) = xref
         rtraj_match(2,nt_inc+1) = vref
C!         rpar_match(1,nt_inc+1) = pxref
C!         rpar_match(2,nt_inc+1) = pvref
         rpar_match(1,1) = pxref
         rpar_match(2,1) = pvref
         rmod_match(1,nt_inc+1) = velref
         if(ray_in_tau) then
           rmod_match(2,nt_inc+1) = sigref
         else
           rmod_match(2,nt_inc+1) = 0.
         end if

        rpar_match(1,2)=-rpar_ref(1,nt_ref,1)*(1.-w_ref)
     *   -rpar_ref(1,nt_ref,2)*w_ref
        rpar_match(2,2)=-rpar_ref(2,nt_ref,1)*(1.-w_ref)
     *  -rpar_ref(2,nt_ref,2)*w_ref
				write(0,*) rpar_match(1,1),rpar_match(1,2),rpar_ref(1,nt_ref,1),w_ref
				write(0,*) rpar_match(2,1),rpar_match(2,2)
         do it=1,nt_ref
c        {
             rtraj_match(1,nt_inc+nt_ref+2-it) =
     *       rtraj_ref(1,it,1)*(1.-w_ref) +
     *       rtraj_ref(1,it,2)*(   w_ref) 

             if(zmap .ne. 1) then
               rtraj_match(2,nt_inc+nt_ref+2-it) =
     *         rtraj_ref(2,it,1)*(1.-w_ref) +
     *         rtraj_ref(2,it,2)*(   w_ref) 
             else
               rtraj_match(2,nt_inc+nt_ref+2-it) =
     *         rtraj_zmap_ref(2,it,1)*(1.-w_ref) +
     *         rtraj_zmap_ref(2,it,2)*(   w_ref) 
             end if

C             rpar_match(1,nt_inc+nt_ref+2-it) =
C     *      -rpar_ref(1,it,1)*(1.-w_ref) -
C     *       rpar_ref(1,it,2)*(   w_ref) 

C             rpar_match(2,nt_inc+nt_ref+2-it) =
C     *      -rpar_ref(2,it,1)*(1.-w_ref) -
C     *       rpar_ref(2,it,2)*(   w_ref) 

             rmod_match(1,nt_inc+nt_ref+2-it) =
     *       rmod_ref(1,it,1)*(1.-w_ref) +
     *       rmod_ref(1,it,2)*(   w_ref) 

             if(ray_in_tau) then
               rmod_match(2,nt_inc+nt_ref+2-it) =
     *         rmod_ref(2,it,1)*(1.-w_ref) +
     *         rmod_ref(2,it,2)*(   w_ref) 
             else
               rmod_match(2,nt_inc+nt_ref+2-it) = 0.
             end if


c        }
         end do

         do it=nt_inc+nt_ref+2,2*maxt
c        {
             rtraj_match(1,it) = rtraj_match (1,nt_inc+nt_ref+1) 
             rtraj_match(2,it) = rtraj_match (2,nt_inc+nt_ref+1)
C             rpar_match(1,it) = rpar_match (1,nt_inc+nt_ref+1)
C             rpar_match(2,it) = rpar_match (2,nt_inc+nt_ref+1)
             rmod_match(1,it) = rmod_match (1,nt_inc+nt_ref+1)
             if(ray_in_tau) then
               rmod_match(2,it) = rmod_match (2,nt_inc+nt_ref+1)
             else
               rmod_match(2,it) = 0.
             end if
c        }
         end do

         if (srite('out',rtraj_match,2*nby_ray).ne.2*nby_ray) then
	   call erexit("Problems writinf Output file")
         end if
         if (srite('Raypar_match',rpar_match,16)
     *                                   .ne.16) then
	   call erexit("Problems writinf Raypar_match file")
         end if
         if (srite('Raymod_match',rmod_match,2*nby_ray)
     *                                   .ne.2*nby_ray) then
	   call erexit("Problems writinf Raymod_match file")
         end if

         if (srite('Time_match',tmatch_tot,esize4)
     *                                  .ne. esize4) then
          call erexit("Problems writing Time_match file")
         end if

         if (srite('Time_ref',tmatch_inc,esize4)
     *                                  .ne. esize4) then
          call erexit("Problems writing Time_ref file")
         end if

				
c        }
         if (srite("location",location,12).ne.12) then
           call erexit("Trouble writing out location")
         end if

         if (srite("Pos_match",rtraj_match(1,nt_inc+1),8).ne.8) then
           call erexit("Trouble writing out Pos_match")
         end if
         			end if
c                    }
                     end if
c                  }
                   end if

c                }
                 end if
c             }
              end do
c           }
            end if
c         }
          end do

         if (.not. found) then
c        {
						bad=bad+1
						write(0,*) "BAD RAY ALERT",bad,'sou_inc=',sou_inc+xsou,'sou_ref=',sou_ref+xsou
					
				 	if (srite('out',rtraj_match,2*nby_ray).ne.2*nby_ray) then
     						call erexit("Problems writinf Output file")
         end if
         if (srite('Raypar_match',rpar_match,16).ne.16) then
           call erexit("Problems writinf Raypar_match file")
         end if
         if (srite('Raymod_match',rmod_match,2*nby_ray).ne.2*nby_ray) then
            call erexit("Problems writinf Raymod_match file")
         end if
          tmatch_tot=-1.
          tmatch_ref=-1.

         if (srite('Time_match',tmatch_tot,esize4)
     *                                  .ne. esize4) then
          call erexit("Problems writing Time_match file")
         end if

         if (srite('Time_ref',tmatch_inc,esize4)
     *                                  .ne. esize4) then
          call erexit("Problems writing Time_ref file")
         end if
         if (srite("location",location,12).ne.12) then
           call erexit("Trouble writing out location")
         end if

         if (srite("Pos_match",rtraj_match(1,nt_inc+1),8).ne.8) then
           call erexit("Trouble writing out Pos_match")
         end if

				end if
         
C       }
        end do
C     }
      end do

      end if 
c         }

c     }
       end do

        if(want_multi .eq. 1) then
          write(0,*)'RECORD MULTIPLE ARRIVALS'
          call putch('n2','d',count_match+bad)
          call auxputch('n2','d',count_match+bad,'Raypar_match')
          call auxputch('n2','d',count_match+bad,'Raymod_match')
          call auxputch('n1','d',count_match+bad,'Time_match')
          call auxputch('n2','d',count_match+bad,'location')
          call auxputch('n1','d',count_match+bad,'Pos_match')
        else
          write(0,*)'FOUND MULTIPLE ARRIVALS BUT NOT RECORDES'
        end if


                 
      end


