       subroutine caladp(ctitle,nthtas)
	  
c***********************************************************************
c      CALADP.FOR
c      CALculation of Amplitudes and Delays by Proximity
c      Faro, Mon Jan  6 20:32:08 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 irefl(np),jbdry(np)
       integer*4 iret(51)
       integer*4 i,ifail,ih,imax,j,jj,k,l,m
       integer*4 nthtas,nret,erays
       
       real*8 tbdry(2,np)
       
       real*8 aaddat(6)
       
       real*8 xl(2),yl(2)

       real*8 thetai,ctheta
       real*8 rh,zh,zray,dz,tauray
       real*8 rerays
       
       complex*16 decay(np)
       complex*16 zl(2)
       complex*16 aray

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

       aad = 'aad00000'

c***********************************************************************
c      Too much launching angles? do not proceed with calculations:

       if (nthtas.gt.nangle) then
       write(6,*) 'ADP: # 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,*) 'ADP: # 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***********************************************************************

       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
           
           call seikeq(thetai,imax,irefl,decay,jbdry,tbdry)
           call sdyneq(thetai,imax,irefl,decay,jbdry,tbdry)
          
c          Got the ray? test for proximity to each hydrophone 
c         (yes, this is slow, can you figure out a better way  
c          to do it?):

           do j = 1,nra
	      
	      rh = rarray(j)
	      
	      if ((rh.ge.rmin).and.(rh.lt.rmax)) then
	      
c***********************************************************************
c                 Check if the ray is returning back or not; if not 
c                 we can bracket it without problems, otherwise we 
c                 need to know how many times it passed by the given
c                 array range: 

		 if (iretrn.eq.0) 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)
		    
		    yl(1) = tau( ih )
		    yl(2) = tau(ih+1)

		    call lini1d(xl,yl,rh,tauray,dummy)
		    
		    zl(1) = amp( ih )
		    zl(2) = amp(ih+1)

		    call clin1d(xl,zl,rh,aray,dummy)

		    do jj = 1,nza

		       zh = zarray(jj)
		       dz = abs(zray-zh)

		       if (dz.lt.miss) then
		       
		       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) = tauray 
		       aaddat(4) = realpart( aray )
		       aaddat(5) = imagpart( aray )
		       aaddat(6) = caustc(ih)
		       
		       call dtelrv(outfil,1,6,aad,aaddat)
		       			  
		       end if
		       
		    end do 

                 else
		 
		    call ebrckt(imax,r,rh,nret,iret)
		    
		    do l = 1,nret
		    
		       xl(1) =   r(iret(l)  )
		       xl(2) =   r(iret(l)+1)

		       yl(1) =   z(iret(l)  )
		       yl(2) =   z(iret(l)+1)

		       call lini1d(xl,yl,rh,zray,dummy)

		       yl(1) = tau(iret(l)  )
		       yl(2) = tau(iret(l)+1)

		       call lini1d(xl,yl,rh,tauray,dummy)
		    
		       zl(1) = amp(iret(l)  )
		       zl(2) = amp(iret(l)+1)

		       call clin1d(xl,zl,rh,aray,dummy)
		    
		       do jj = 1,nza

		          zh = zarray(jj)
		          dz = abs(zray-zh)

		          if (dz.lt.miss) then
			  
			  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) = tauray 
		          aaddat(4) = realpart( aray )
		          aaddat(5) = imagpart( aray )
		          aaddat(6) = caustc( iret(l) )
			  
			  call dtelrv(outfil,1,6,aad,aaddat)
					       		       
		          end if

		       end do

		    end do
		 
		 end if

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

	      end if 

	   end do 

           end if

       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
