clear
close all

% Ray tracing in a linear refractive index atmosphere

%Initial elevation angle (deg)
theta=0.5;

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

% Source ht (m)
ht=5;

%Initial M value
m0=350;

% Effective earth radius multiplier
k=1.33;

% Start program

rad=pi/180;
dm=(1/k)*m0/6370e3;

i1=1;

x=0;
z=ht;
theta=90-theta;
range=range*1000;
step=step*1000;
dlay=1.e-6;

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

while (x <=range),
  no=m0+z*dm;
  x=x+step;
  z=z+step/tan(rad*theta);	
  nz=m0+z*dm;

  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;
 
ray1(1)=200;
ray1(2)=0.2;
z=0;
i1=3;
for i=1:200
  ray1(i1)=m0+z*dm;
  z=z+0.1;
  i1=i1+1;
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