       subroutine seikeq(thetax,imax,irefl,decay,jbdry,tbdry)

c**********************************************************************
c      Solve EIKonal EQuations
c      Written by Orlando Camargo Rodriguez
c      orodrig@ualg.pt
c      Universidade do Algarve
c      Physics Department
c      Signal Processing Laboratory
c      Gambelas, Fri Mar  6 10:27:38 WET 2015
c****&******************************************************************

c      Include global (common) variables:

       include 'global.for'
       
c      Define local variables:
       
       integer*4 irefl(np),jbdry(np)
       
       integer*4 i,j
       integer*4 jrefl,jrefr,ibdry,ithtas
       integer*4 noj,nrefl,nrk
       
       integer*4 imax,ikill,idn,iup
       
       real*8 tbdry(2,np)

       real*8 roj(npo),zupj(npo),zdnj(npo)
       
       real*8 yold(4),ynew(4),fold(4),fnew(4)
       
       real*8 taui(2),taur(2),taub(2)
       real*8 la(2),lb(2),li(2)
       real*8 e1(2),es(2),eso(2)
       real*8 normal(2)
       real*8 sigma(2),gradc(2)
       real*8 v(2)
       
       real*8 ri,zi,zatii,zbtyi,ziup,zidn
       real*8 sigmar,sigmaz,sigmai,sri,szi
       real*8 ci,cc,cri,czi,crri,czzi,crzi,cnn
       real*8 dr,dz,dsi,ds4,ds5,dic,dt
       real*8 ap,as,rho1,rho2,cp2,cs2
       real*8 dphase,phasei,phaseii
       real*8 theta,thetax
       real*8 t1,t2,k1,k2,phi1,phi2,beamd,omega
       real*8 za
       real*8 erros
       real*8 adBoW
       real*8 dummy
       
       complex*16 decay(np)

       complex*16 refl,refldc
       complex*16 rk1,rk2

c***********************************************************************

       rho1 = 1.0
       omega = 2*pi*freqx
       
c***********************************************************************
c      Define initial conditions:

       ikill = 0
       iup   = 0
       idn   = 0
 
       srefl = 0
       brefl = 0
       orefl = 0
       jrefl = 0

       iretrn = 0

       nrk    = 0
       
       irefl(1) = jrefl

       refldc = (1.0,0.0)

       decay(1) = refldc 

       phase(1) = 0.0 
       
       i = 1
       
       es(1) = cos( thetax )
       es(2) = sin( thetax )
 
       e1(1) = -es(2)
       e1(2) =  es(1)
 
       call csvals(rx,zx,cx,cc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)

       sigmar = sigmai*es(1)
       sigmaz = sigmai*es(2)
       
          r(1) = rx
          z(1) = zx
       cray(1) = cx
        tau(1) = 0.0
          s(1) = 0.0
         ic(1) = 0.0
          
       yold(1) = rx
       yold(2) = zx
       yold(3) = sigmar
       yold(4) = sigmaz
       
       fold(1) = es(1)
       fold(2) = es(2)
       fold(3) = sri
       fold(4) = szi

       rmin = r(1)
       rmax = r(1)

       beamd = 0.0

c-----------------------------------------------------------------------
c      Trace the ray:

1000   if (( r(i).lt.rbox(2) ).and.( r(i).gt.rbox(1) ).and.
     &     ( ikill.eq.0 )) then
 
          dsi = ds

c         Runge-Kutta integration:

2000       call rkf45(dsi,yold,fold,ynew,fnew,ds4,ds5)

           eso(1) = fold(1)
           eso(2) = fold(2)

            es(1) = fnew(1)
            es(2) = fnew(2)
	    
           erros = 1.0-abs(eso(1)*es(1)+eso(2)*es(2))
	    
	    if (erros.gt.5.0e-8) then
	    
	    nrk = nrk + 1

	     if (nrk.gt.100) then
	     write(6,*) 'Runge-Kutta integration: failure in step'
	     write(6,*) 'convergence, aborting calculations...'
	     stop 
	     end if
	    
	    dsi = dsi/2
	    goto 2000
	    end if

c         Test for boundary intersections: 
          
          ri = ynew(1)
          zi = ynew(2)

          if ((ri.gt.rbox(1)).and.(ri.lt.rbox(2))) then
	  call bdryi(nati,rati,zati,aitype,ri,zatii,v,normal)
	  call bdryi(nbty,rbty,zbty,bitype,ri,zbtyi,v,normal)
          else 
	  ikill = 1
	  end if 

c[====================================================================]
c         Check if the ray is still between the boundaries; if not 
c         find the intersection point and calculate the reflection:
          
          if (((zi.lt.zatii).or.(zi.gt.zbtyi)).and.(ikill.eq.0)) then
           
           la(1) = yold(1)
           la(2) = yold(2)
          
           lb(1) = ynew(1)
           lb(2) = ynew(2)

c======================================================================
c         Ray above surface?
           
          if (zi.lt.zatii) then

             call raybi(nati,rati,zati,aitype,la,lb,li)

             ri = li(1)
             zi = li(2)
 
             call bdryi(nati,rati,zati,aitype,ri,zatii,taub,normal)
             
             ibdry = -1

             srefl = srefl + 1
             jrefl = 1

c            Calculate surface reflection:

             call reflct( normal, es, taur, theta )
	     
c            Check if the ray is digging in beyond the surface:

             call bdryi(nati,rati,zati,aitype,ri+ds*taur(1),zatii,v,v)
	     
	     if ((zi+ds*taur(2)).lt.zatii) then
	     ikill = 1
	     end if

c            Surface reflection:

             if (atype.eq.'A') then
             ikill = 1
             elseif (atype.eq.'E') then
               if (aptype.eq.'H') then
	         rho2 = rhoati(1)
	          cp2 =  cpati(1)
		  cs2 =  csati(1)
                  call bdryr(rho1,rho2,ci,cp2,cs2,apa,asa,theta,refla)
               elseif (aptype.eq.'N') then
                 call bdryi(nati,rati,rhoati,aitype,ri,rho2,v,v)
                 call bdryi(nati,rati, cpati,aitype,ri, cp2,v,v)
                 call bdryi(nati,rati, csati,aitype,ri, cs2,v,v)
                 call bdryi(nati,rati, apati,aitype,ri,  ap,v,v)
                 call bdryi(nati,rati, asati,aitype,ri,  as,v,v)
                 lambda = cp2/freqx
		 call cnvnts(ap,lambda,freqx,atiu,adBoW)
                 apa = adBoW
		 lambda = cs2/freqx
                 call cnvnts(as,lambda,freqx,atiu,adBoW)
                 asa = adBoW
                 call bdryr(rho1,rho2,ci,cp2,cs2,apa,asa,theta,refla)
               else
               end if
	       
c                   Beam displacement calculations:
            
c	       lambda = ci/freqx
c               t1 = theta - pi*5.0e-4
c	       t2 = theta + pi*5.0e-4
c	       k1 = omega/ci*cos( 0.5*pi - t1 )
c	       k2 = omega/ci*cos( 0.5*pi - t2 )
c               call bdryr(rho1,rho2,ci,cp2,cs2,apa,asa,t1,rk1)
c	       call bdryr(rho1,rho2,ci,cp2,cs2,apa,asa,t2,rk2)
c	       phi1 = atan2(imagpart(rk1),realpart(rk1))
c	       phi2 = atan2(imagpart(rk2),realpart(rk2))
c	       beamd = (phi1-phi2)/(k2-k1)
	       
             else
             end if
             
	     refl   = refla
             refldc = refldc*refl

c           Kill the ray if the reflection decay is too small: 

            if ( abs(refldc) .lt. 1.0e-6 ) ikill = 1 

          end if
c======================================================================
c         Ray below bottom?
 
          if (zi.gt.zbtyi) then
 
             call raybi(nbty,rbty,zbty,bitype,la,lb,li)
             
             ri = li(1)
             zi = li(2)

             call bdryi(nbty,rbty,zbty,bitype,ri,zbtyi,taub,normal)

c            Invert the normal at the bottom for reflection:
              
             normal(1) = -normal(1)
             normal(2) = -normal(2)

             ibdry = 1

             brefl = brefl + 1
             jrefl = 1

c            Calculate bottom reflection:
             
             call reflct( normal, es, taur, theta )
	     
c            Check if the ray is digging in beyond the bottom:

             call bdryi(nbty,rbty,zbty,bitype,ri+ds*taur(1),zbtyi,v,v)
	     
	     if ((zi+ds*taur(2)).gt.zbtyi) then
	     ikill = 1
	     end if

c            Bottom reflection:

             if (btype.eq.'A') then
             ikill = 1
             elseif (btype.eq.'E') then
               if (bptype.eq.'H') then	         
	         rho2 = rhobty(1)
	          cp2 =  cpbty(1)
		  cs2 =  csbty(1)
                 call bdryr(rho1,rho2,ci,cp2,cs2,apb,asb,theta,reflb)
               elseif (bptype.eq.'N') then
                 call bdryi(nbty,rbty,rhobty,bitype,ri,rho2,v,v)
                 call bdryi(nbty,rbty, cpbty,bitype,ri, cp2,v,v)
                 call bdryi(nbty,rbty, csbty,bitype,ri, cs2,v,v)
                 call bdryi(nbty,rbty,acpbty,bitype,ri,  ap,v,v)
                 call bdryi(nbty,rbty,acsbty,bitype,ri,  as,v,v)
		 lambda = cp2/freqx
                 call cnvnts(ap,lambda,freqx,btyu,adBoW)
                 apb = adBoW
		 lambda = cs2/freqx
                 call cnvnts(as,lambda,freqx,btyu,adBoW)
                 asb = adBoW
                 call bdryr(rho1,rho2,ci,cp2,cs2,apb,asb,theta,reflb)
               else
               end if
	       
c                   Beam displacement calculations:
            
c	       lambda = ci/freqx
c               t1 = theta - pi*5.0e-4
c	       t2 = theta + pi*5.0e-4
c	       k1 = omega/ci*cos( 0.5*pi - t1 )
c	       k2 = omega/ci*cos( 0.5*pi - t2 )
c               call bdryr(rho1,rho2,ci,cp2,cs2,apb,asb,t1,rk1)
c	       call bdryr(rho1,rho2,ci,cp2,cs2,apb,asb,t2,rk2)
c	       phi1 = atan2(imagpart(rk1),realpart(rk1))
c	       phi2 = atan2(imagpart(rk2),realpart(rk2))
c	       beamd = (phi1-phi2)/(k2-k1)

             else
             end if

             refl   = reflb
             refldc = refldc*refl

c           Kill the ray if the reflection decay is too small: 

            if ( abs(refldc) .lt. 1.0e-6 ) ikill = 1 

          end if
c======================================================================

c         Update marching solution and function:
                    
          ri = li(1)
          zi = li(2)
                     
          call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     &    crri,czzi,crzi)
     
          ynew(1) = ri
          ynew(2) = zi
          ynew(3) = sigmai*taur(1)
          ynew(4) = sigmai*taur(2)

          fnew(1) = taur(1)
          fnew(2) = taur(2)
          fnew(3) = sri
          fnew(4) = szi

          end if
c[====================================================================]

c[oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo]
c         Object reflection:

          if (nobj.gt.0) then

c         For each object detect if the ray is inside the object range: 
  
          do j = 1,nobj
          
          noj = no(j) 

          if ((ri.ge.ro(j,1)).and.(ri.lt.ro(j,noj))) then

           do k = 1,noj
           
              roj(k) =  ro(j,k)
             zdnj(k) = zdn(j,k)
             zupj(k) = zup(j,k)
           
           end do        
           
           if (zdnj(1).ne.zupj(1)) then
           write(6,*) 'Lower and upper faces do not start at the same'
           write(6,*) 'depth! aborting calculations...' 
           stop
           end if 

           call bdryi(noj,roj,zdnj,oitype,ri,zidn,v,normal)
           call bdryi(noj,roj,zupj,oitype,ri,ziup,v,normal)

c          Second point is inside the object?

           if ((ynew(2).ge.zidn).and.(ynew(2).lt.ziup)) then

           la(1) = yold(1)
           la(2) = yold(2)
          
           lb(1) = ynew(1)
           lb(2) = ynew(2)

c=============================================================================>
c            Which face was crossed by the ray: upper or lower?
c
c            Case 1: from left to right, beginning in box & end in box:
             if ((ynew(1).gt.yold(1)).and.(yold(1).ge.roj(1)).and.
     &           (ynew(1).lt.roj(noj))) then

                call bdryi(noj,roj,zdnj,oitype,yold(1),zidn,v,normal)
                call bdryi(noj,roj,zupj,oitype,yold(1),ziup,v,normal)

                if (yold(2).lt.zidn) then
                    call raybi(noj,roj,zdnj,oitype,la,lb,li)
                    idn = 1
                    iup = 0
                else 
                    call raybi(noj,roj,zupj,oitype,la,lb,li)
                    idn = 0
                    iup = 1
                end if

c            Case 2: from left to right, beginning out box & end in box
c            (zdnj(1) = zupj(1)):
              elseif ((ynew(1).gt.yold(1)).and.(yold(1).lt.roj(1)).and.
     &               (ynew(1).lt.roj(noj))) then

               za = lb(2)-(lb(2)-la(2))/(lb(1)-la(1))*(lb(1)-roj(1))
                la(1) = roj(1)
                la(2) = za
                if (za.lt.zupj(1)) then
                    call raybi(noj,roj,zdnj,oitype,la,lb,li)
                    idn = 1
                    iup = 0
                else
                    call raybi(noj,roj,zupj,oitype,la,lb,li)
                    idn = 0
                    iup = 1
                end if

c            Case 3: from right to left, beginning in box & end in box:
             elseif ((ynew(1).lt.yold(1)).and.(yold(1).le.roj(noj)).and.
     &	             (ynew(1).ge.roj(1))) then
	     
                call bdryi(noj,roj,zdnj,oitype,yold(1),zidn,v,normal)
                call bdryi(noj,roj,zupj,oitype,yold(1),ziup,v,normal)

                if (yold(2).lt.zidn) then
                    call raybi(noj,roj,zdnj,oitype,la,lb,li)
                    idn = 1
                    iup = 0
                else 
                    call raybi(noj,roj,zupj,oitype,la,lb,li)
                    idn = 0
                    iup = 1
                end if

c            Case 4: from right to left, beginning out box & end in box
c            (zdnj(1) = zupj(1)):
             elseif ((ynew(1).lt.yold(1)).and.(yold(1).gt.roj(noj)).and.
     &               (ynew(1).ge.roj(1))) then

                za = lb(2)-(lb(2)-la(2))/(lb(1)-la(1))*(lb(1)-roj(noj))
                la(1) = roj(noj)
                la(2) = za
                if (za.lt.zupj(noj)) then
                    call raybi(noj,roj,zdnj,oitype,la,lb,li)
                    idn = 1
                    iup = 0
                else
                    call raybi(noj,roj,zupj,oitype,la,lb,li)
                    idn = 0
                    iup = 1
                end if
	     
             else
             write(6,*) 'Object reflection case: ray beginning'
             write(6,*) 'neither behind or between object box,'
             write(6,*) 'check object coordinates...'
             stop
             end if
c=============================================================================>
             ri = li(1)
             zi = li(2)

c            Face reflection: upper or lower? 

             if ((idn.eq.1).and.(iup.eq.0)) then 
             call bdryi(noj,roj,zdnj,oitype,ri,zidn,taub,normal)
             normal(1) = -normal(1)
             normal(2) = -normal(2)
             ibdry = -1
             elseif ((idn.eq.0).and.(iup.eq.1)) then
             call bdryi(noj,roj,zupj,oitype,ri,ziup,taub,normal)
             ibdry =  1
             else
             write(6,*) 'Object reflection case: ray neither being'
             write(6,*) 'reflected on down or up faces,'
             write(6,*) 'check object coordinates...'
             stop
             end if

             orefl = orefl + 1 
             jrefl = 1

c            Calculate object reflection: 

             call reflct( normal, es, taur, theta )

c            Object reflection => get the reflection coefficient
c            (kill the ray is the object is an absorver):

             if (otype(j:j).eq.'A') then
             refl = (0.0,0.0)
             ikill = 1
             elseif (otype(j:j).eq.'R') then
             refl = (1.0,0.0)
             elseif (otype(j:j).eq.'V') then
             refl = (-1.0,0.0)
             elseif (otype(j:j).eq.'E') then
                 rho2 = orho(j)
                  cp2 =  ocp(j)
                  cs2 =  ocs(j)
                  ap  =  oap(j)
                  as  =  oas(j)
		  lambda = cp2/freqx
                 call cnvnts(ap,lambda,freqx,obju(j:j),adBoW)
                  ap  = adBoW
		  lambda = cs2/freqx
                 call cnvnts(as,lambda,freqx,obju(j:j),adBoW)
                  as  = adBoW
                 call bdryr(rho1,rho2,ci,cp2,cs2,ap,as,theta,refl)
             else
             write(6,*) 'Unknown object type (neither A,E,R or V),'
             write(6,*) 'aborting calculations...'
             stop
             end if

            refldc = refldc*refl

c          Kill the ray if the reflection coefficient is too small: 
	
           if ( abs(refl) .lt. 1.0e-5 ) ikill = 1

c          Update marching solution and function:
          
           es(1) = taur(1)
           es(2) = taur(2)
          
           e1(1) = -es(2)
           e1(2) =  es(1)
          
           call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     &     crri,czzi,crzi)
                    
           ynew(1) = ri
           ynew(2) = zi
           ynew(3) = sigmai*es(1)
           ynew(4) = sigmai*es(2)

           fnew(1) = es(1)
           fnew(2) = es(2)
           fnew(3) = sri
           fnew(4) = szi

           end if
          
          end if 
          
          end do 
          
          end if 
           
c[oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo]

c            Beam displacement correction:

c          ri = ynew(1)
c          zi = ynew(2)
	  
c          call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
c     &    crri,czzi,crzi)
     
c          lambda = ci/freqx
	  
c	  if (beamd.gt.lambda) then
c	   ynew(1) = ynew(1) + beamd
c	    r(i+1) =   r(i)  + beamd
c	    z(i+1) =   z(i)           
c	  tau(i+1) = tau(i)  + beamd/ci
c          cray(i+1) =    ci
c             s(i+1) =   s(i) + beamd
c            ic(i+1) =  ic(i) + beamd*cray(i+1)
c         Option 1: refraction along the boundary, then specular reflection	    
c         irefl(i+1) =     0
c         phase(i+1) = phase(i)
c         Option 2: specular reflection, then refraction along the boundary    
c         irefl(i+1) = jrefl
c         phase(i+1) = phase(i) - atan2(imagpart(refl),realpart(refl))
c	 jrefl = 0
c         jbdry(i+1) = ibdry
c	 if (ibdry.eq.-1) then
c	 call bdryi(nati,rati,zati,aitype,ri,zatii,taub,normal)	
c	 else
c	 call bdryi(nbty,rbty,zbty,bitype,ri,zbtyi,taub,normal)
c	 end if
c         tbdry(1,i+1) = taub(1)
c         tbdry(2,i+1) = taub(2)
c	  i = i + 1
c	  beamd = 0.0
c	  end if

c-----------------------------------------------------------------------
          ri     = ynew(1)
          zi     = ynew(2)

          r(i+1) = ynew(1)
          z(i+1) = ynew(2)
          
          es(1) = fnew(1)
          es(2) = fnew(2)
          
          call csvals(ri,zi,ci,cc,sigmai,cri,czi,sri,szi,
     &    crri,czzi,crzi)

          dr = r(i+1) - r(i)
          dz = z(i+1) - z(i)
          
          dsi = sqrt( dr*dr + dz*dz )
          
           tau(i+1) = tau(i) + dsi/ci
          cray(i+1) =    ci
             s(i+1) =   s(i) + dsi
            ic(i+1) =  ic(i) + dsi*cray(i+1)

         irefl(i+1) = jrefl
         jbdry(i+1) = ibdry

         tbdry(1,i+1) = taub(1)
         tbdry(2,i+1) = taub(2)
          
          if (jrefl.eq.1) then
	  phase(i+1) = phase(i) - atan2(imagpart(refl),realpart(refl))
	  else 
	  phase(i+1) = phase(i)
	  end if 
	  
          jrefl = 0
          ibdry = 0
          nrk   = 0
	   
          taub(1) = 0.0
          taub(2) = 0.0

          decay(i+1) = refldc

          do j = 1,4
       
             yold(j) = ynew(j)
             fold(j) = fnew(j)
         
          end do
        
c-----------------------------------------------------------------------
c         Next one:

          i = i + 1

c-----------------------------------------------------------------------
c         Prevent further calculations if there is no more space in 
c         the memory for the ray coordinates:

          if ( i.gt.np ) then

          write(6,*) 'Ray step too small, number of points in ray'
          write(6,*) 'coordinates > ',np,',' 
          write(6,*) 'aborting calculations...'
          stop

          end if

c-----------------------------------------------------------------------
  
          go to 1000

       end if

       imax = i
       
       nrefl = srefl + brefl + orefl

c-----------------------------------------------------------------------
c      Cut the ray at box exit:

        dr =   r(imax)-  r(imax-1)
        dz =   z(imax)-  z(imax-1)
       dic =  ic(imax)- ic(imax-1)
        dt = tau(imax)-tau(imax-1)

       if (r(imax).gt.rbox(2)) then
            z(imax) =   z(imax-1) + (rbox(2)-r(imax-1))* dz/dr
           ic(imax) =  ic(imax-1) + (rbox(2)-r(imax-1))*dic/dr
	  tau(imax) = tau(imax-1) + (rbox(2)-r(imax-1))* dt/dr
           r(imax) = rbox(2)
       end if
       
       if (r(imax).lt.rbox(1)) then
            z(imax) =   z(imax-1) + (rbox(1)-r(imax-1))* dz/dr
           ic(imax) =  ic(imax-1) + (rbox(1)-r(imax-1))*dic/dr
	  tau(imax) = tau(imax-1) + (rbox(1)-r(imax-1))* dt/dr
           r(imax) = rbox(1)
       end if

c-----------------------------------------------------------------------
c      Search for refraction points (refraction angles are zero!), 
c      rmin, rmax and twisting of rays: 
       
       nrefr = 0

       do i = 2,imax-1
       
          rmin = min( rmin, r(i) )
	  rmax = max( rmax, r(i) )
 
          prod = ( z(i+1)-z(i) )*( z(i)-z(i-1) )

          if ( ( prod.lt.0 ).and.(irefl(i).ne.1) ) then
              
              if (nrefr.gt.np2) then
              write(6,*) 'Refraction points: nrefr > ',np2
              write(6,*) 'aborting calculations...'
              stop
              end if
              
              nrefr = nrefr + 1
 
              rrefr(nrefr) = r(i)
              zrefr(nrefr) = z(i)

          end if

          prod = ( r(i+1)-r(i) )*( r(i)-r(i-1) )

          if ( prod.lt.0 ) then
          iretrn = 1
          end if

          prod = 0.0 

       end do
       
       rmin = min( rmin, r(imax) )
       rmax = max( rmax, r(imax) )

c***********************************************************************
c     Back to main:
c***********************************************************************

       return

       end
