%----------------------------------------------------------------------
% linray
% program to draw rays for linear gradient models
% Coded by Fred Duennebier
%
% velocity= vzero+k * z
% Adapted to make interactive by G. Ito 3/03
%----------------------------------------------------------------------
clear all;
fig='figure(1)';
eval(fig); clf;
%figure(2); clf;
xmax=60; ymax=20;	%min max bounds in distance to plot	
dv=500;         %tick increment on velocity plot
vmax=4000;		%maximum velocity to plot
tmax=0.04;		%maximum time to plot
iother=1;		%plot every iother ray in raypath plot

%----------------------------------------------------------------------
%parameters that can be adjusted if needed
%----------------------------------------------------------------------
startangle=0.1;			%SPECIFY: min. take-off angle from vertical
endangle=89;  			%SPECIFY: max. take-off angle from vertical
pointsper=30;  			% SPECIFY points to plot per ray segment
numberofrays=200;  		% SPECIFY: number of rays

imodify=3;
%----------------------------------------------------------------------
%Loop to play with velocities
%----------------------------------------------------------------------
%startangle=90-max_takeoff_angle;  endangle=90-min_takeoff_angle; 

while (imodify~=2)

%----------------------------------------------------------------------
%Define velocity model
%----------------------------------------------------------------------
vmod;

%----------------------------------------------------------------------
%Calculate ray paths
%----------------------------------------------------------------------

calc_raypath;
iray=[1:iother:numberofrays];

%----------------------------------------------------------------------
%plot raypaths and T vs X
%----------------------------------------------------------------------
rc=['  b'; 'r-.'; '  g'; 'k-.'; '  b'; 'r-.'; '  g'; 'k-.'];
subplot(222); hold off;
iray=[1:iother:numberofrays];
ww=w(:,iray); uu=u(:,iray);
mww=max(ww);
clear jcolor
for k=1:numlayer;
 irc=min([length(rc) k]);
 jj=find(mww>zees(k)&mww<zees(k+1));
 plot(uu(:,jj)-rayxbegin,ww(:,jj),rc(irc,:)); hold on;
end
%plot(u,w);	
xlay=[0 xmax]'*ones([1 numlayer+1]); 
ylay=[1 1]'*zees;
hold on;	
plot(xlay,ylay,'--k')
xlabel ('Shot-Receiver Distance (m)')
ylabel ('Depth (m)')
axis ([0 xmax 0 ymax]);
grid on
set(gca, 'ydir', 'reverse')
subplot(121)
orient landscape


%figure(2); clf;
subplot(224); hold off;
plot(dist-rayxbegin,raytime,'k-'); hold on;
rc=['bo-'; 'ro-'; 'go-'; 'ko-'; 'bo-'; 'ro-'; 'go-'; 'ko-'];
rc=['b.'; 'r.'; 'g.'; 'k.'; 'b.'; 'r.'; 'g.'; 'k.'];
for k=1:numlayer;
 irc=min([length(rc) k]);
 jj=find(mww>zees(k)&mww<zees(k+1));
 if (sum(jj)>0)
   plot(dist(jj)-rayxbegin,raytime(jj),rc(irc,:)); 
 end
end
grid on;
xlabel ('Shot-Receiver Distance (m)')
ylabel ('Travel Time (s)')
axis([0 xmax 0 tmax]);
grid on;


iloop=1;
disp('Enter:')
disp( 'return= modify velocity model');
disp( '1=      revert to previous velocity model and pick new ones');
disp(['2=      EXIT and write velocities to output file']);
imodify=input('--->');
if (sum(size(imodify))==0); imodify=0; end;
end

outfile=input('Enter name of file to write velocities: >','s');
out=[zees; vees];
fid=fopen([outfile '.vmod'],'w');
fprintf(fid,'%10.3f %10.3f \n',out);
fclose(fid);