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

c**********************************************************************
c      Solve DYNamic EQuations
c      (ray tracing subroutine for cylindrical symmetry)
c      Written by Orlando Camargo Rodriguez
c      orodrig@ualg.pt
c      Universidade do Algarve
c      Physics Department
c      Signal Processing Laboratory
c      Faro, Sun Jun 15 13:49:31 WEST 2014
c****&******************************************************************

c      Include global (common) variables:

       include 'global.for'
       
c      Define local variables:
       
       integer*4 irefl(np),jbdry(np)
       
       integer*4 i,imax
       integer*4 ibdry
       
       real*8 tbdry(2,np)
       
       real*8 es(2),normal(2),taub(2)
       real*8 gradc(2),ngradc(2),dgradc(2)
       real*8 sigma(2),sigman(2)
       
       real*8 thetax,ci,cii,cxc,cnn
       real*8 cri,czi,crri,czzi,crzi
       real*8 ri,zi,si,sigmai,zatii,zbtyi
       real*8 cnj,csj,rm,rn
       real*8 dr,dz,dsi
       real*8 drds,dzds
       real*8 alpha
       real*8 prod
       real*8 mynan,mnsone
       
       complex*16 decay(np)
       complex*16 ap_aq

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

       lambda = cx/freqx
       
c      Let me know if you can figure out a better way to define a NaN:

       mnsone = -1.0 
       mynan  = sqrt( mnsone )

c***********************************************************************
       
c      Define initial conditions:

       p(1) = 1.0
       q(1) = 0.0

       caustc(1) = 0.0 
       
       ibdry = 1

c      Get Thorpe attenuation in dB/m: 

       call thorpe(freqx,alpha)

c-----------------------------------------------------------------------
c      Solve the Dynamic Equations:

       do i = 1,imax-1
       
       ri = r(i+1)
       zi = z(i+1)

       call csvals(ri,zi,cii,cxc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)
       
       ngradc(1) = cri
       ngradc(2) = czi
       
       ri = r(i)
       zi = z(i)

       call csvals(ri,zi,ci,cxc,sigmai,cri,czi,sri,szi,
     & crri,czzi,crzi)
       
       gradc(1) = cri
       gradc(2) = czi
       
       dgradc(1) = ngradc(1) - gradc(1)
       dgradc(2) = ngradc(2) - gradc(2)
       
       dr = r(i+1) - r(i)
       dz = z(i+1) - z(i)
       
       dsi = sqrt( dr*dr + dz*dz )
       
          es(1) = dr/dsi
          es(2) = dz/dsi

       sigma(1) = sigmai*es(1)
       sigma(2) = sigmai*es(2)
       
       drdn = -es(2)
       dzdn =  es(1)
       
       cnn = ( drdn**2 )*crri + 2.0*drdn*dzdn*crzi + ( dzdn**2 )*czzi

       ci = cray(i)
       
       cxc = ci*ci
       
       if (irefl(i+1).ne.1) then
       
          p(i+1) = p(i) - q(i)*(cnn/cxc)*dsi
          q(i+1) = q(i) + p(i)*   ci    *dsi

c         This is the Dragon: Bellhop's refraction correction:

	  sigman(1) = -sigma(2)
	  sigman(2) =  sigma(1)
	  
	  call dot( dgradc, sigman, cnj )
	  call dot( dgradc, sigma , csj )
	  
	  if (sigma(2).ne.0.0) then 
	   
	  rm  = sigma(1)/sigma(2)
	  rn  = -( rm * ( 2.0*cnj - rm*csj )/cii )
	  
	  p(i+1) = p(i) + q(i)*rn
          
	  end if 
	  
       elseif (irefl(i+1).eq.1) then

          ibdry = jbdry(i+1)

          taub(1) = tbdry(1,i+1)
          taub(2) = tbdry(2,i+1)

	  call rnrefl(ibdry,sigma,taub,gradc,ci,rn)

          p(i+1) = p(i) + q(i)*rn
          q(i+1) = q(i)
	  
       else

       write(6,*) 'Solving dynamic equations:'
       write(6,*) 'irefl neither equal 1 or 0!'
       write(6,*) 'aborting calculations...'
       stop

       end if
       
       prod = q(i)*q(i+1)
       
       if ((prod.le.0.0).and.(q(i).ne.0.0)) then
       caustc(i+1) = caustc(i) + pi/2.0
       else
       caustc(i+1) = caustc(i)
       end if 

       end do

c-----------------------------------------------------------------------
c      Amplitude calculation:
       
       do i = 2,imax

       ap_aq = cray(1)*cos(thetax)*cray(i)/( ic(i)*q(i) )

       amp(i) = sqrt( ap_aq )*decay(i)*exp( -alpha*s(i) )

       end do 

       amp(1) = amp(2)

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

       return

       end
