%This script runs an ode solver on the eikonal equation like described in
%the problem set. The velocity model is set by functions profile.m and c.m.
%Initial angle is given by theta in degrees (0 to 90).
%The solution is integrated until the ray hits the surface, the depth of 50,
%or distance 100 or the time 40, whichever comes first. This can be changed by editing the files.
%This program is not foolproof, mainly because the ode solver can run into
%problems with integration. The ray can also get caught in a waveguide, but escape
%again because of numerical instabillities. Comments are scarce and random,
%most should be self explanatory if some thought is given to it. Plotrays.m
%plots some of the solution
%****** F R A G I L E  ---- H A N D L E   W I T H   C A R E ******

%Written by Keli Karason, karason@mit.edu, 8 Nov. 1999

[zz,cc,theta]=getinput;

%Initial and final times.
to=0;
tfinal=40;

%Initial position in space
xo=0;
zo=0;

theta=theta(:);
%Initial slowness vector from initial angles
pxo=sin(theta/180*pi)/veldep(zo);
pzo=cos(theta/180*pi)/veldep(zo);

%Getting lengths
nth=length(theta);	%number of rays to trace
nr=100;

%initializing variables
t=zeros(100,nth);

x=zeros(100,nth);
z=zeros(100,nth);
px=zeros(100,nth);
pz=zeros(100,nth);


%Running the numerical integration and interpolating to get even number of time steps
options=odeset('Events','on');
for j=1:nth
  [tt,yy]=ode45('yprime',[to tfinal],[xo zo pxo(j) pzo(j)]',options);
  t(:,j)=linspace(0,tt(end),nr)';
  x(:,j)= interp1(tt,yy(:,1),t(:,j),'linear');
  z(:,j)= interp1(tt,yy(:,2),t(:,j),'linear');
  px(:,j)=interp1(tt,yy(:,3),t(:,j),'linear');
  pz(:,j)=interp1(tt,yy(:,4),t(:,j),'linear');
  disp([num2str(j) ' out of ' num2str(nth) ' finished']);
end  

%Getting values for wavefronts
twf=(1:1:20);
xwf=zeros(nth,length(twf));
zwf=zeros(nth,length(twf));
for j=1:nth,
  xwf(j,:)=interp1(t(:,j),x(:,j),twf);
  zwf(j,:)=interp1(t(:,j),z(:,j),twf);
end

plotrays