       subroutine caladr(ctitle,nthtas)

c***********************************************************************
c      CALADR.FOR
c      CALculation of Amplitudes and Delays by Regula falsi
c      Faro, Tue Jan  7 00:30:12 WET 2014
c      Written by Orlando Camargo Rodriguez
c****&******************************************************************

c      Include global (common) variables:

       include 'global.for'

c      Define local variables:
       
       character(len=128) header       
       character(len=60 ) ctitle
       character(len=8  ) aad,v1
       character(len=5  ) cv
       character(len=4  ) civ
       character(len=3  ) ciii
       character(len=2  ) cii
       character(len=1  ) ci
   
       integer*4 inan(nangle,nhyd2)
       integer*4 irefl(np),jbdry(np)
       integer*4 i,ifail,ih,imax,il,ir,j,jj,k,l,m
       integer*4 nfail,neig,nrays,nthtas,ntrial,erays
       
       real*8 tbdry(2,np)
       
       real*8 mz(nangle,nhyd2)
       
       real*8 f(nangle),thtas(nangle)
       real*8 thetal(np2),thetar(np2)
       real*8 aaddat(6)
       real*8 xl(2),yl(2)
       
       real*8 thetai,ctheta
       real*8 rh,zh,zray
       real*8 thtae,thtal,thtar,thta0
       real*8 fl,fr,f0,prod
       real*8 mynan,mnsone
       real*8 rerays
       
       complex*16 decay(np)

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

       aad = 'aad00000'

c***********************************************************************
c      Let me know if you can figure out a better way to define a NaN:

       mnsone = -1.0 
       mynan  = sqrt( mnsone )
       
c***********************************************************************
c      Too much launching angles? do not proceed with calculations:

       if (nthtas.gt.nangle) then
       write(6,*) 'ADR: # of launching angles > ',nangle
       write(6,*) 'aborting calculations...'
       stop 
       end if
       
c***********************************************************************
c      Too much hydrophones along range? do not proceed with calculations:

       if (nra.gt.nhyd2) then
       write(6,*) 'ADR: # of  receivers along range > ',nhyd2
       write(6,*) 'aborting calculations...'
       stop 
       end if
       
c***********************************************************************

       erays = 0 
       
       header = 'Signal Processing LABoratory: TRACEO
     & acoustic ray tracing program - amplitude and delay calculations'
     
      header(117:124) = char(ichar(' '))
      header(125:126) = char(0) // char(1)
       
      write(header(127:128),'(a)') int(19785,kind=2)
      
      open(unit=outfil,file='aad.mat',form='unformatted',
     & access='stream')
     
c     Open the file and write the header:

      do i = 1, 16
         j = 8*(i-1) + 1
         write(outfil) header(j:j+7)
      end do

c***********************************************************************
c      Solve the EIKonal and the Dynamic sets of EQuations:
c***********************************************************************

       nrays = 0 

       do i = 1,nthtas

           thetai = -thetas(i)*pi/180.0

           ctheta = abs( cos( thetai ) ) 

c          Trace a ray as long as it is neither 90 or -90:
           
           if (ctheta.gt.1.0e-7) then
           
	   nrays = nrays + 1 
	   
	   thtas(nrays) = thetai 
	   
           call seikeq(thetai,imax,irefl,decay,jbdry,tbdry)
           call sdyneq(thetai,imax,irefl,decay,jbdry,tbdry)
	   
	   if (iretrn.eq.1) then
	   write(6,*) 'Eigenray search by Regula Falsi:'
	   write(6,*) 'a returning ray detected for ',thetas(i)
	   write(6,*) 'aborting calculations...'
	   stop
	   end if 
          
c          Got the ray? then fill the matrix of depths: 

           do j = 1,nra
	      
	      rh = rarray(j)
	      
	      if ((rh.ge.rmin).and.(rh.lt.rmax)) then
	      
	      call brcket(imax,r,rh,ih)
	      
	      xl(1) =   r( ih )
	      xl(2) =   r(ih+1)
	      
	      yl(1) =   z( ih )
	      yl(2) =   z(ih+1)
	      
	      call lini1d(xl,yl,rh,zray,dummy)
	      
	        mz(i,j) = zray
	      inan(i,j) = 0
	      
	      else
 
	        mz(i,j) = mynan
	      inan(i,j) = 1
	      
	      end if 
	      
	   end do 
	   
           end if
	   
       end do

c***********************************************************************
c      Angle loop finished? proceed searching eigenrays at each point 
c      of the array:
c***********************************************************************

       do i = 1,nra
       
          rh = rarray(i)
       
          do j = 1,nza 
	  
	     zh = zarray(j)

c            At each range calculate the difference between 
c            the hydrophone and ray depths:

	     do k = 1,nrays
	     
	        f(k) = zh - mz(k,i) 
		
	     end do

c            Determine the number of eigenrays by looking at 
c            signal variations (or zero values) of the difference:

	     neig = 0 
	     
	     do k = 1,nrays-1
	
	     	fl = f( k )
	     	fr = f(k+1)
	     	prod = fl*fr

                il = inan( k ,i)
		ir = inan(k+1,i)
	     	
		if ((il.eq.0).and.(ir.eq.0)) then 
		
		   if ((fl.eq.0.0).and.(fr.ne.0.0)) then
	
		   neig = neig + 1
		   thetal(neig) = thtas( k )
		   thetar(neig) = thtas(k+1)
	
		   elseif ((fr.eq.0.0).and.(fl.ne.0.0)) then

		   neig = neig + 1
		   thetal(neig) = thtas( k )
		   thetar(neig) = thtas(k+1)
	
		   elseif (prod.lt.0.0) then
	
		   neig = neig + 1
		   thetal(neig) = thtas( k )
		   thetar(neig) = thtas(k+1)
	
		   else
	
		   end if
	        
		end if 
		
		if (neig.gt.np2) then
		write(6,*) 'ADR: # of arrivals at (',rh,',',zh,')'
		write(6,*) '> ',np2
                write(6,*) 'aborting calculations...'
                stop
		end if
		
	     end do
	     
c            Time to find eigenrays; either we are lucky or 
c            we need to apply regula falsi:

             nfail = 0 

c            If neig = 0 the loop is iddle:
	     
             do l = 1,neig
	     
	        rbox(2) = rh
		 
		thtal = thetal(l)
		thtar = thetar(l)
		
		call seikeq(thtal,imax,irefl,decay,jbdry,tbdry)
		
		fl = z(imax) - zh
		
		call seikeq(thtar,imax,irefl,decay,jbdry,tbdry)
		
		fr = z(imax) - zh
		
		if (abs(fl).le.miss) then
		
		    thtae = thtal
		    ifail = 0
		    
		elseif (abs(fr).le.miss) then
		
		    thtae = thtar
		    ifail = 0
		    
		else
		 
		    ntrial = 0

1000                thta0 = thtar - fr*( thtal - thtar )/( fl - fr )

                    if (ntrial.gt.21) then
		    ifail = 1
		    nfail = nfail + 1
		    write(prtfil,*) rh,zh
	            write(prtfil,*) 'Eigenray search failure,'
	            write(prtfil,*) 'skipping to next case...'
		    end if

                    call seikeq(thta0,imax,irefl,decay,jbdry,tbdry)

		    f0 = z(imax) - zh
		    
		    if (abs(f0).lt.miss) then

		        ifail = 0
		        thtae = thta0

		    elseif (ifail.ne.1) then

		        prod = fl*f0
			ntrial = ntrial + 1

			if ( prod.lt.0.0 ) then
			thtar = thta0
			fr    = f0
			else
			thtal = thta0
			fl    = f0
			end if

			goto 1000
			
		    else

		    end if
		    
		end if 
		
		if (ifail.ne.1) then
		
		   call seikeq(thtae,imax,irefl,decay,jbdry,tbdry)
		   call sdyneq(thtae,imax,irefl,decay,jbdry,tbdry)
                   
	           erays = erays + 1
			  
c#######################################################################
		  if (erays.lt.10) then
		  write (ci, '(I1)'), erays
		  aad(8:8) = ci
		  elseif ((erays.ge.10).and.(erays.lt.100)) then
		  write (cii, '(I2)'), erays
		  aad(7:8) = cii
		  elseif ((erays.ge.100).and.(erays.lt.1000))  then
		  write (ciii, '(I3)'), erays
		  aad(6:8) = ciii
		  elseif ((erays.ge.1000).and.(erays.lt.10000))  then
		  write (civ, '(I4)'), erays
		  aad(5:8) = civ
		  else
		  write (cv, '(I5)'), erays
		  aad(4:8) = cv
		  end if
c#######################################################################
		          
		  aaddat(1) = rh
		  aaddat(2) = zh
		  aaddat(3) = tau(imax)
		  aaddat(4) = realpart( amp(imax) )
		  aaddat(5) = imagpart( amp(imax) )
		  aaddat(6) = caustc(imax)
		  
		  call dtelrv(outfil,1,6,aad,aaddat)
		   
	        end if
		
	     end do 

	  end do  
       
       end do 
       
c***********************************************************************
c      The Happy End:
c***********************************************************************

       rerays = 1.0*erays 
       
       v1 = 'neigrays'
	   
       call dtelrv(outfil,1,1,v1,rerays)
       
       v1 = 'rarray2d'
	   
       call dtelrv(outfil,1,nra,v1,rarray)

       v1 = 'zarray2d'
	   
       call dtelrv(outfil,1,nza,v1,zarray)
       
c***********************************************************************
c      Close the OUTFIL:

       close(outfil)

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

       return

       end
