% this script call tau.m which computes distances, turning points,
% and travel-times for rays and plots the results
%
% you can change the range in freqency and angular degree by modifying
% the first two statements in this script
%
% see ray_codes.ps for details
%
% aaron birch may 22 2002

%
% the frequencies to compute at
%
w=2*pi*linspace(0.0030,0.0030,1);

%
% the angular degrees to compute
%
l=70:5:100;



%
% load the model
%
load model

%
% allocate space for the results
%
time=zeros(length(w),length(l));
delta=zeros(size(time));


%
% loop over frequency, and then angular degree, computing for each one
%
for (i=1:length(w))
  i,
  for (j=1:length(l))
    %
    % call the code to do the actual calculation
    %
    [time(i,j),delta(i,j),lower(i,j),upper(i,j)]=tau(w(i),l(j),r,cs,wcs);
  end;
end;

%
% plot the results
%
clf;c=['r' 'g' 'c' 'b' 'm'];
for (i=1:length(w))
  plot(delta(i,:),time(i,:)/60,c(1+mod(i,5)));hold on;
end;
xlabel('\Delta (deg)');    
ylabel('\tau (minutes)');