%---------------------------------------------------------------
%					t2x2			(L. Braile 1/20/95, revised 2/24/01)
%---------------------------------------------------------------
% t2x2 uses the t squared - x squared method to calculate a 
% 1-D velocity model from t and x reflection data.
%---------------------------------------------------------------
% Inputs are title, nrefl, number of relection x-t pairs
% for each reflection and the x-t data in matrices.
% Inputs are contained in a m-files (sample shown below)
% which is called by the program to input the data.
%---------------------------------------------------------------
% Sample input file txdat.m:
%
%	txtitle=['test t2-x2 input'];
%	nrefl=[4];
%	pts=[5 5 4 6];
%	maxn=max(pts);
%	x1=zeros(nrefl,maxn);
%	x1(1,1:6)=	[1 	2 	3 	4 	5 	0];
%	t1(1,1:6)=	[.2 .25 .3 .4 .55 	0];
%	x1(2,1:6)= [1 	2	3	4	5	0];
%	t1(2,1:6)=	[.4 .42 .44	.47	.52 0];
%	x1(3,1:6)=	[1 	2 	3 	4 	0 	0];
%	t1(3,1:6)=	[.6 .61 .623 .64 	0 	0];
%	x1(4,1:6)= [1 	2	3	4	5	6];
%	t1(4,1:6)=	[.8 .81 .82 .832 .845 .86];
%	x=x1';
%	t=t1';
%---------------------------------------------------------------
%
% The matrices t and x must be full (zeros are used to fill
% out the matrices for missing x-t pairs) and are of dimension
% (max number of data points for a reflector, nrefl).  The x-t 
% data are in the columns of the x and t matrices.
%
% Input the data: title, nrefl, pts, x, t
%
%--------------------------------------------------------------- 
txdat       % CHANGE HERE FOR DIFFERENT INPUT FILE.
%---------------------------------------------------------------
% For each reflection, use the x-t data to calculate t squared
% and x squared, the rms velocity from the slope of the least
% squares line fit to the t^2 - x^2 data and the two-way travel
% time from the intecept time.  Then, calculate the interval
% velocity from the Dix equation.  Plot data and least squares
% travel time fits and print out the data, fit, and derived model.
%
%---------------------------------------------------------------
for i=1:nrefl
	npts=pts(i);
	for j=1:npts;
		x2(j)=x(j,i)^2;
		t2(j)=t(j,i)^2;
	end
	mode=0;
	sigmay=zeros(npts,1);
%---------------------------------------------------------------
% Calculate least squares fit to t^2 - x^2 data for this reflection
	[a,sigmaa,b,sigmab,r]=linfit(x2,t2,sigmay,npts,mode);
	t0(i)=sqrt(a);
	sl(i)=b;
	if i==1
		v=sqrt(1/sl(i));
		h=v*t0(i)/2;
	else
		v=sqrt((1/(t0(i)-t0(i-1)))*(t0(i)/sl(i)-t0(i-1)/sl(i-1)));
		h=h + v*(t0(i)-t0(i-1))/2;
	end
	vrms=sqrt(1/sl(i));
	hsave(i)=h;
	asave(i)=t0(i);
	t0std(i)=sigmaa^.5;
	vrmssv(i)=vrms;
	vrmsstd(i)=sigmab^.5;
	vsave(i)=v;
end
%---------------------------------------------------------------
% Print output table
txtitle
['x data']
x
['t data']
t
layer=1:nrefl;
['     layer  depth(km)   t0(s)  std t0(s)  vrms(km/s) std vrms  v layer']
A=zeros(nrefl,7);
for i=1:nrefl;
A(i,:)=[layer(i),hsave(i),asave(i),t0std(i),vrmssv(i),vrmsstd(i),vsave(i)]; 
end
A
plot(x,t,'or')
xlabel(' DISTANCE (km)')
ylabel('TIME (s)')
title('t^2-x^2 ANALYSIS, o = DATA, -- = CALCULATED')
hold on
for i=1:nrefl
	a2=asave(i)^2;
	b=sl(i);
	xt=x(1:pts(i),i);
	z=feval('txline',xt,a2,b);
	plot(xt,z,'k')
end
hold off
pause
% Plot velocity - depth models
% Layer velocities
zl = [0];
vl = [];
for i=1:nrefl
   zl = [zl hsave(i) hsave(i)];
   vl = [vl vsave(i) vsave(i)];
end
vl = [vl vsave(nrefl)];
% Vrms velocities
zr = [0 hsave];
vr = [vrmssv(1) vrmssv];
% plot data
plot(zl,vl,'-r')
hold on
plot(zr,vr,'ob',zr,vr,'--b')
xlabel(' Depth (km)')
ylabel('Velocity (km/s)')
hold off
% end of t2x2 
%---------------------------------------------------------------