%  program txcurv.m
%  txcurv calculates reflection and refraction (head wave) travel times
%  for an n-layered half space of plane, homogeneous, dipping layers.
%  The last layer veocity corresponds to an infinite half space.
% L. Braile, 2/10/03.  Modified from FORTRAN program TXCURV.
%
% Input data (change for other models)
n = 5;  % number of layers
% last layer thickness should be 999 to simulate half space
th = [1 2 1 3 999];  % thicknesses of layers (km or m)
v = [2 3 4 4.5 6];  % layer velocities (km/s or m/s consistent with thicknesses)
dp = [-.5 0 -0.5 -0.5 0];  % dip angle (degrees) is negative for updip and positive 
%  for downdip.  Dip angles refer to the bottom of the layer. dp(n) = 0.
mindist = 0;  % minimum distance for t-x calculation
maxdist = 15;  % maximum distance for t-x calculation
dx = 0.4;  % distance increment
redvel = 6;  % reducing velocity, if redvel = 0, use actual travel times
% end of input
%
x = [mindist:dx:maxdist];  % distances for t-x (km or m consistent with thicknesses)
% Begin calculations
%
% calculate refraction (head wave) travel times
nx = max(size(x));  % number of distances
time = zeros(n,1);
ttrefr = zeros(n,nx);
for i = 1:nx
    dist = x(i);
    time = rfrtim(th,v,dp,n,dist);
    ttrefr(:,i) = time';
end
%
% calculate reflection travel times
% calculate exact travel times using rfltim.m
time = zeros(n-1);
ttrefl = zeros(n-1,nx);
for i = 1:nx
    dist = x(i);
    time = rfltim(th,v,dp,n,dist,0.001);
    ttrefl(:,i) = time';
end
%  apply reducing velocity if redvel > 0
if redvel > 0
    for i = 1:nx
        ttrefr(:,i) = ttrefr(:,i) - x(i)/redvel;
        ttrefl(:,i) = ttrefl(:,i) - x(i)/redvel;
    end
end
% plot  t-x curve results
figure
hold on
for i = 1:n
    plot(x,ttrefr(i,:),'-k') % plot direct wave and n-1 head wave arrivals
end
for i = 1:n-1
    plot(x,ttrefl(i,:),'-k') % plot n-1 reflections from bottom of first n-1 layers
end
title('Refraction (Head Wave) and Reflection Travel Times')
xlabel('Distance (km)')
if redvel == 0
    ylabel('Time (s)')
    else
    ylabel('Reduced Time (s)')
end
hold off
pause
figure
% calculate depths at the "A" end (left side of model) and
% depths and thicknesses at the "B" end (right side of model
DA = zeros(n,1);
DB = DA;
for i = 1:n-1
    DA(i) = sum(th(1:i));
end
for i = 1:n-1
    DB(i) = DA(i) + tan(dp(i));
end
thB(1) = DB(1);
if n > 2
    for i = 2:n-1
        thB(i) = DB(i) - DB(i-1);
    end
end
% find maximum depth
ZMAXA = max(DA);
ZMAXB = max(DB);
ZMAX = max([ZMAXA,ZMAXB]);
ZMAX = ZMAX + 0.1*ZMAX;
% plot model
NL=n-1;	% Number of lines
XA=0;
XB=maxdist;
hold on
for i=1:NL
	z1=DA(i);
	z2=DB(i);
	x2=[XA XB];
	y=[z1 z2];
	plot(x2,y,'r-')
end
x1=[0 XB];
y=[0 0];
plot(x1,y,'r-')
axis('ij')
axis([0 maxdist 0 ZMAX])
% label velocities
DA(n)=ZMAX;
DB(n)=ZMAX;
damin=0;
dbmin=0;
for i=1:n
	xave=XB/2;
	zave=(damin+dbmin+DA(i)+DB(i))/4;
	label=num2str(v(i));
	text(xave,zave,label)
	damin=DA(i);
	dbmin=DB(i);
end
hold off