       subroutine calari(ctitle,nthtas)

c***********************************************************************
c      CALARI.FOR
c      CALculation of All Ray Information
c      Faro, Mon Jan  6 14:41:38 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  ) ray,rinfo,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 i,imax,j,k,nthtas
       
       real*8 raydat(np,6)
       real*8  tbdry(2,np)
       real*8   refr(2,np2)
       real*8 rayinf(5)
       
       real*8 thetai,ctheta

       complex*16 decay(np)

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

       ray   = 'ray00000'
       rinfo = 'inf00000'       

c***********************************************************************
       
       header = 'Signal Processing LABoratory: TRACEO
     & acoustic ray tracing program - All Ray Information'
     
      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='ari.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     Write the angles:

      v1 = 'ray_elev'
	   
      call dtelrv(outfil,1,nthtas,v1,thetas)

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

        do i = 1,nthtas
       
c.......................................................................
c          Yes, this is ugly, but it works! 

	   if (i.lt.10) then
           j = 4
	   write (ci, '(I1)'), i
	      ray(  8:8) = ci       
	      rinfo(8:8) = ci
	   elseif ((i.ge.10).and.(i.lt.100)) then
	   j = 5
	   write (cii, '(I2)'), i
	      ray(  7:8) = cii
	      rinfo(7:8) = cii
	   elseif ((i.ge.100).and.(i.lt.1000))  then
	   j = 6
	   write (ciii, '(I3)'), i
	      ray(  6:8) = ciii
	      rinfo(6:8) = ciii
	   elseif ((i.ge.1000).and.(i.lt.10000))  then
	   j = 7
	   write (civ, '(I4)'), i
	      ray(  5:8) = civ
	      rinfo(5:8) = civ
	   else
	   j = 8
	   write (cv, '(I5)'), i
	      ray(  4:8) = cv
	      rinfo(4:8) = cv
	   end if 
c.......................................................................

           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	    Instead of saving each parameter independently let's 
c           put the information in matrices (4th row => real 
c           part of the amplitude, 5th row => complex part)

	    do k = 1,imax
	    raydat(k,1) =   r(k)
	    raydat(k,2) =   z(k)
	    raydat(k,3) = tau(k)
	    raydat(k,4) = realpart( amp(k) )
	    raydat(k,5) = imagpart( amp(k) )
	    raydat(k,6) = caustc(k)
	    end do
	    
	    call dtelm6(outfil,6,imax,ray,raydat)
	    	     
	    rayinf(1) = iretrn 
	    rayinf(2) =  srefl 
	    rayinf(3) =  brefl 
	    rayinf(4) =  orefl 
	    rayinf(5) =  nrefr
	    
            call dtelrv(outfil,1,5,rinfo,rayinf)
	    	    
	   end if

       end do
       
c***********************************************************************
c      Close the OUTFIL:

       close(outfil)
    
c***********************************************************************
c     Back to main:
c***********************************************************************

       return

       end
