program ray

c Ray tracing code for an m profile specified in file ray1.pro
c the first line of ray1.pro has an integer telling the # of m values
c the second line of ray1.pro tells the height step in meters for these 
c values
c following lines of ray1.pro have the m values (not in M-units!) at each 
c height starting with zero height
c Code outputs to files x.out and z.out, use rayplot.m in Matlab to
c plot the rays
 
	implicit none

	integer mmax,i,nlay,junk,num
	parameter (mmax=2000)
	real*8 m(mmax),dlay,theta
	real*8 pi,rad,x,z,ht,range,step,n,no

	common m,dlay,nlay

	pi=4.*atan(1.d0)
	rad=pi/180.

	write(*,*)'Initial elevation angle? (deg)'
        read(*,*)theta
	write(*,*)'range, horiz. step? (km)'
	read(*,*)range,step
	write(*,*)'Source ht? (m)'
	read(*,*)ht

	open(unit=8,file='ray1.pro',status='old')
	open(unit=9,file='x.out',status='unknown')
	open(unit=10,file='z.out',status='unknown')

	read(8,*)nlay
        if (nlay.gt.mmax) then
          write(*,*)'Need bigger value of mmax!'
	  stop
        endif
	read(8,*)dlay
	do i=1,nlay
	  read(8,*)m(i)
        enddo

	x=0
	z=ht
	theta=90-theta
	range=range*1000
	step=step*1000

	write(9,*)x
	write(10,*)z

	do while (x.lt.range)

	  no=n(z)

 	  x=x+step
          z=z+step/dtan(rad*theta)	

	  if (dabs(step/dtan(rad*theta)).gt.dlay) 
     .	    write(*,*)'Try a smaller range step!'

	  if ((dabs(dsin(rad*theta)*no/n(z)).gt.1).or.(z.le.0.)) then
	    theta=180.-theta
	    num=z/dlay
	    if (num.le.0) then
	      z=-z
	    elseif (theta.lt.90.) then
              z=2.*dlay*(num+1)-z
	    else
              z=2.*dlay*num-z
	    endif
	  else
	    if (theta.gt.90.) then
	      theta=180.-dasin(dsin(rad*theta)*no/n(z))/rad
	    else
	      theta=dasin(dsin(rad*theta)*no/n(z))/rad
	    endif
	  endif

   	  write(9,*)x
	  write(10,*)z
	  write(11,*)theta
	enddo

	end

c------------------------------------------

	function n(z)

	integer mmax,nlay,num
        parameter (mmax=2000)
	real*8 m(mmax),dlay,n,z

	common m,dlay,nlay
	
	num=int(z/dlay)+1

	if (num.gt.nlay) then
	  n=1.
	  write(*,*)'Warning! Going above specified profile!'
	else
	  n=m(num)
	endif

	return
	end