% Ray tracing code for an m profile specified in file ray1.pro
% the first line of ray1.pro has an integer telling the # of m values
% the second line of ray1.pro tells the height step in meters for these 
% values
% following lines of ray1.pro have the m values (not in M-units!) at each 
% height starting with zero height
 
clear
close all

%Initial elevation angle (deg)
theta=0.5;

%Total range and horizontal step (km)
range=100;
step=0.1;

% Source ht (m)
ht=5;

load ray1.pro;
nlay=ray1(1);
dlay=ray1(2);
m(1:nlay)=ray1(3:nlay+2);

% Start program

rad=pi/180;

i1=1;

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

xx(i1)=x;
zz(i1)=z;

while (x <=range);

% Get index of refraction at height z through closest value in profile
  nnum=floor(z/dlay)+1;
  if (nnum > nlay)  
    no=1; 
    display('Warning! Going outside of profile!'); 
  else
    no=m(nnum);
  end;

  x=x+step;
  z=z+step/tan(rad*theta);	

  if (abs(step/tan(rad*theta)) >dlay) 
    display('Smaller range step recommended!');
  end;

% Get index of refraction at new height z through closest value in profile
  nz=no;
  if (z>0)  
    nnum=floor(z/dlay)+1;
    if (nnum > nlay)  
      nz=1; 
      display('Warning! Going outside of profile!'); 
    else
      nz=m(nnum);
    end;
  end;
  if ( (abs(sin(rad*theta)*no/nz)>1) | (z <= 0 ) )

     theta=180.-theta;
     num=z/dlay;

     if (num <= 0) 
       z=-z;
%       display('Intersects! Try another initial angle!');
     elseif (theta < 90 )
       z=2.*dlay*(num+1)-z;
     else
       z=2.*dlay*num-z;
     end
  else
     if (theta > 90.)
       theta=180.-asin(sin(rad*theta)*no/nz)/rad;
     else
       theta=asin(sin(rad*theta)*no/nz)/rad;
     end;
  end;
  i1=i1+1;
  xx(i1)=x;
  zz(i1)=z;
end;
 
x=xx';
z=zz';

clear zz;

for i=1:ray1(1)-2,
zz(i)=(i-1)*ray1(2);
end


figure(1)
subplot(211)
plot(1.e6*(ray1(3:ray1(1))-1),zz);
title('Modified Refractive Index Profile')
xlabel('M(h) - M units')
ylabel('Height (m)')
grid on

subplot(212)
line(x/1000,z/1000)
title('Ray paths')
xlabel('Distance (km)')
ylabel('Height (km)')
grid on