% program txrefl.m
% txrefl calculates approximate reflection travel times versus distance
% for a stack of n plane, homogeneous, isotropic laters using the
% t-squared, x-squared equation and both everage velocities and rms
% velocities.  txrefl also calculates "exact" travel times by tracing
% rays through the model with appropriate Snell's law calculations.  The
% ray tracing is performed iteratively for each distance until the
% calculated distance is within a tolerance (TEST) of the requested
% distance using function rfltim.
% L. Braile, 2/7/03
clear all
% Input data (change for other models)
n = 5;  % number of layers
%change the first number below to 2.0 for part "b".
th = [0.5 1.0 1.0 2.0 3.0];  % thicknesses of layers (km or m)
v = [1.5 2 3 4 5];  % layer velocities (km/s or m/s consistent with thicknesses)
x = [0:0.4:12];  % distances for t-x (km or m consistent with thicknesses)
% Begin calculations
nx = max(size(x));  % number of distances
t = th./v;  % one way vertical travel times in layers
% Plot vel vs depth
d=cumsum(th);	
% depth to top of each layer
d=[0 d];		
for i=1:n
	z(2*i-1)=d(i);
	z(2*i)=d(i+1);
	vel(2*i-1)=v(i);
	vel(2*i)=v(i);
end
plot(vel,z,'k-')
axis('ij')
axis([0 6 0 9])
xlabel('Velocity (km/s)')
ylabel('Depth (km)')
pause;
% calculate average velocities
for i = 1:n
    sum1 = sum(v(1:i).*t(1:i));  % numerator
    sum2 = sum(t(1:i));  % denominator
    Vave(i) = sum1/sum2;
end
% calculate rms velocities
for i = 1:n
    sum1 = sum(v(1:i).*v(1:i).*t(1:i));  % numerator
    sum2 = sum(t(1:i));  % denominator
    Vrms(i) = sqrt(sum1/sum2);
end
% calculate vertical incidence travel times (two way)
for i = 1:n
    t0(i) = 2*sum(t(1:i));
end
% calculate travel times for reflections for layers 1 to n for all
% distances using both average velocity and rms velocity formulas
ttave = zeros(n,nx);
ttrms = zeros(n,nx);
tt = zeros(n,nx);
for i = 1:n
    ttave(i,:) = sqrt(t0(i)*t0(i) + x.*x/(Vave(i)*Vave(i)));
    ttrms(i,:) = sqrt(t0(i)*t0(i) + x.*x/(Vrms(i)*Vrms(i)));
end
% calculate exact travel times using rfltim.m
dp = zeros(n,1);
time = zeros(n);
for j = 1:nx
    dist = x(j);
    time = rfltim(th,v,dp,n+1,dist,0.001);
    tt(:,j) = time';
end
% plot results
figure
axis('ij')
axis ([0 10 0 8])
title('Reflections travel times: + = Vave formula; o = Vrms formula; -- = Exact times')
xlabel('Distance (km)')
ylabel('Two Way Travel Time (s)')
hold on
for i = 1:n
    plot(x,ttave(i,:),'+r',x,ttrms(i,:),'ob',x,tt(i,:),'-k')
    pause;
end


hold off

    