program ray

c Ray tracing code for a linear m profile versus height, specified by k 
c Code outputs to files x.out and z.out, use rayplot.m in Matlab to
c plot the rays
 
	implicit none

	integer i,nlay,junk,num
	real*8 k,pi,rad,dlay,theta
	real*8 x,z,ht,range,step,n,no,m0,dm

	common m0,dm

	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
c	write(*,*)'Initial M value, deriv (/m) ?'
c	read(*,*)m0,dm
	write(*,*)'Initial M value, k ?'
	read(*,*)m0,k

	dm=(1/k)*m0/6370e3

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

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

	dlay=1.e-6

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

	do while (x.lt.range)

	  no=n(z)

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

c	  if (dabs(step/dtan(rad*theta)).gt.dlay) 
c     .	    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
	      write(*,*)'Intersects! Try another initial angle!'
c	      write(9,*)'Intersects!'
c	      write(10,*)'Intersects!'
c	      stop
	    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

	open(unit=12,file='ray1.pro',status='unknown')
	write(12,*)200
	write(12,*)0.2
	z=0.
	do i=1,200
	  write(12,*)m0+z*dm
	  z=z+0.1
	enddo

	end

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

	function n(z)

	integer nlay,num
	real*8 m0,dm,n,z

	common m0,dm
	
	n=m0+z*dm

	return
	end