% 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
% Input data (change for other models)
n = 4;  % number of layers
th = [1 2 1 3];  % thicknesses of layers (km or m)
v = [2 3 4 4.5];  % layer velocities (km/s or m/s consistent with thicknesses)
x = [0:0.4:15];  % 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
% 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
hold on
for i = 1:n
    plot(x,ttave(i,:),'+r',x,ttrms(i,:),'ob',x,tt(i,:),'-k')
end
title('Reflections travel times: + = Vave formula; o = Vrms formula; -- = Exact times')
xlabel('Distance (km)')
ylabel('Two Way Travel Time (s)')
hold off