       subroutine calcpr(ctitle,nthtas)
	  
c***********************************************************************
c      CALCPR.FOR
c      CALculation of Coherent Acoustic Pressure
c      Gambelas, Tue Jan  7 11:43:17 WET 2014
c      Written by Orlando Camargo Rodriguez
c****&******************************************************************
c      Note: 
c      readin.for will stop calculations if the array type is incorrect.
c      Therefore, if the array is no horizontal, vertical or linear then 
c      it should be rectangular.
c****&******************************************************************

c      Include global (common) variables:

       include 'global.for'

c      Define local variables:

       character(len=128) header       
       character(len=60 ) ctitle
       character(len=8  ) v1
       
       integer*4 irefl(np),jbdry(np)
       integer*4 iret(51)
       integer*4 i,ih,imax,j,jj,k,nthtas,nret

       real*8 tbdry(2,np)
       real*8 p1d(2,nhyd)
       real*8 es(2),e1(2),deltar(2)
       
       real*8 thetai,ctheta
       real*8 rh,zh,n,nxn,dzdr,qray
       real*8 sray,dr,dz,q0,w

       complex*16 decay(np),paux(npress)
       
       complex*16 aray,press
       
c***********************************************************************

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

       header = 'Signal Processing LABoratory: TRACEO
     & acoustic ray tracing program - Coherent PRessure'
     
      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='cpr.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***********************************************************************

       omega  = 2.0*pi*freqx

        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)
	   
	   q0 = cx/( pi*dtheta/180.0 )
	   
c          Now that the ray has been calculated let's determine the
c          ray influence at each point of the array:

           if (artype.eq.'HRY') then
               
	       zh = zarray(1)
               
	       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)

		    call grayp(imax,irefl,ih,q0,rh,dzdr,tauray,zray,
     &		    aray,qray,w)

 		    call gpress(ih,rh,zh,tauray,zray,dzdr,aray,w,press)

 		    pres1d(j) = pres1d(j) + press

		  else

		    call ebrckt(imax,r,rh,nret,iret)
		    
		    do jj = 1,nret

		       call grayp(imax,irefl,iret(jj),q0,rh,dzdr,
     &                 tauray,zray,aray,qray,w)
 		       call gpress(iret(jj),rh,zh,tauray,zray,dzdr,
     &                 aray,w,press)

 		       pres1d(j) = pres1d(j) + press

 		    end do

		  end if
c***********************************************************************

		  end if

               end do

	   elseif (artype.eq.'VRY')  then

	       rh = rarray(1)
	       
	       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)
	
		 do j = 1,nza
	
 		    zh = zarray(j)

		    call grayp(imax,irefl,ih,q0,rh,dzdr,tauray,zray,
     &		    aray,qray,w)
 		    call gpress(ih,rh,zh,tauray,zray,dzdr,aray,w,press)

 		    pres1d(j) = pres1d(j) + press

 		 end do
		 
		 else
		  
		 call ebrckt(imax,r,rh,nret,iret)
		    
		 do jj = 1,nret
 		    
		 do k = 1,nza

 		    zh = zarray(k)
	
		    call grayp(imax,irefl,iret(jj),q0,rh,dzdr,
     &              tauray,zray,aray,qray,w)
 		    call gpress(iret(jj),rh,zh,tauray,zray,dzdr,
     &              aray,w,press)

 		    pres1d(k) = pres1d(k) + press

 		 end do
		    
		 end do

		 end if

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

	       end if

           elseif (artype.eq.'LRY')  then
	   
               do j = 1,nra
           
	          rh = rarray(j)
                  zh = zarray(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)

		    call grayp(imax,irefl,ih,q0,rh,dzdr,tauray,zray,
     &		    aray,qray,w)
 		    call gpress(ih,rh,zh,tauray,zray,dzdr,aray,w,press)

 		    pres1d(j) = pres1d(j) + press
		  
		  else
		  
		    call ebrckt(imax,r,rh,nret,iret)
		    
		    do jj = 1,nret
 		    
		    do k = 1,nza

 		       zh = zarray(k)
	
		       call grayp(imax,irefl,iret(jj),q0,rh,dzdr,
     &                 tauray,zray,aray,qray,w)
 		       call gpress(iret(jj),rh,zh,tauray,zray,dzdr,
     &                 aray,w,press)

 		       pres1d(j) = pres1d(j) + press

 		    end do
		    
		    end do

		  end if
		  
c***********************************************************************
                  end if
		  
	       end do

	   else
               
	       do j = 1,nra
	       
                  rh = rarray(j)
		  
c                 Check first if the array range is inside the min 
c                 and max ranges of the ray:

		  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)
	
 		    do k = 1,nza

 		       zh = zarray(k)
	
		       call grayp(imax,irefl,ih,q0,rh,dzdr,tauray,zray,
     &		       aray,qray,w)
 		       call gpress(ih,rh,zh,tauray,zray,dzdr,aray,w,
     &		       press)

 		       pres2d(k,j) = pres2d(k,j) + press

 		    end do
		  
		  else
		  
		    call ebrckt(imax,r,rh,nret,iret)
		    
		    do jj = 1,nret
 		    
		    do k = 1,nza

 		       zh = zarray(k)
	
		       call grayp(imax,irefl,iret(jj),q0,rh,dzdr,
     &                 tauray,zray,aray,qray,w)
 		       call gpress(iret(jj),rh,zh,tauray,zray,dzdr,
     &                 aray,w,press)

 		       pres2d(k,j) = pres2d(k,j) + press

 		    end do
		    
		    end do

		  end if
		  
c***********************************************************************
		  
		  end if

               end do

           end if

           end if

       end do

c***********************************************************************
c      Write the data to the output file:
c***********************************************************************

        v1 = 'pressure'
	
        if (artype.ne.'RRY') then
	
	    jj = max(nra,nza)
	    
	    call dtelcm(outfil,jj,1,v1,pres1d)
		
	else

c         Yes, this is not elegant, but it works!!!!!

            jj = 0

	    do j = 1,nra
	    
	       do i = 1,nza
	       
	       jj = jj + 1
	       
	       paux(jj) = pres2d(i,j)
	       
	       end do
	    
	    end do
	    
            call dtelcm(outfil,nza,nra,v1,paux)
	    
	end if

c***********************************************************************
c      The Happy End:
c***********************************************************************

       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
