%****************************************************************
%				REFRIN				(L. Braile  6/21/95)
%
%				PROGRAM REFRIN INTERPRETS SEISMIC REFRACTION DATA BASED ON  
%			THE SLOPE-INTERCEPT METHOD.
%
%	 DATA ENTRY---DATA FILE = refrdat.m (or change file name below)
%					                                     
%			Line 1--- TITLE = ['    '];  
%
%           Line 2--- N = [   ];
%			N IS THE NUMBER OF LAYERS OBSERVED FOR A GIVEN TIME-DISTANCE CURVE.
%
%			Line 3--- X = [   ];
%			X IS THE SPREAD LENGTH.
%				UNITS ARE KM AND KM/SEC, M AND M/SEC, OR FT AND FT/SEC; AND SEC
%				OR KFT/s KFT and seconds.
%
%			Line 4---VA(I),(I=1:N) = [    ]; 
%				VA(I) IS A LIST OF THE VELOCITIES  OBSERVED WITH INCREASING DEPTH.
%
%			Line 5---VB(I),(I=1:N)  = [    ];
%					SAME INPUT AS VA(I) BUT FOR THE REVERSE PROFILE
%
%			Line 6---TAI(I), (I=1:N) = [    ];
%					TAI(I) IS A LIST OF THE TIME INTERCEPTS (SECS.) OF EACH LAYER
%					 THE INTERCEPT OF THE FIRST LAYER IS 0.
%
%			Line 7---TBI(I), (I=1:N) = [    ];
%
%					TBI(I) IS A LIST OF THE INTERCEPTS (SECS.) OF THE REVERSE
%					  PROFILE.  SAME INPUT FORMAT AS TAI(I).
%
%			Line 8---ZMAX=[   ];
%					Maximum Depth for Model Plot axis
%
%***************************************************************** 
%  Read input data [CHANGE FILE NAME TO READ APPROPRIATE DATA FILE]
%
refrdat
%
%*****************************************************************  
%	Print input data
%
format compact		% modify print convention for output
'Reversed refraction profile interpretation by Slope - Intercept Method'
TITLE
'Number of Layers = ',N
'Spread Length = ',X
'  Layer #       VA        VB       TA        TB'
for i=1:N
	AB(i,:) = [i VA(i) VB(i) TAI(i) TBI(i)];
end
AB
%
%  	Check for Reciprocity of input data
%
	for I=2:N
		TBB=TAI(I)+X*(1./VA(I)-1./VB(I)); 
		if TBI(I)<=0
			TBI(I)=TBB;
			else
			TAEND=TAI(I)+X/VA(I);
			TBEND=TBI(I)+X/VB(I);
			ERROR=abs(TAEND/TBEND-1.);
		end
		if (ERROR-0.10) >=0
			'intercept times inconsistent for layer ', I
		end
	end
%	Calculate velocity for first layer
%
		V(1)=(VA(1)+VB(1))*0.5;
%	Calculate angles, velocities, depths for layers 2 to N
%
for M=2:N			% Large Loop, over Number of Layers
		K=1;
		ALPH(1)=asin(V(1)/VB(M));
		BETA(1)=asin(V(1)/VA(M));
	if (M-2)<=0
			A(1)=(ALPH(1)+BETA(1))*0.5;
			W(2)=(ALPH(1)-BETA(1))*0.5;
			V(2)=V(1)/sin(A(1));
	else
			A(1)=ALPH(1)-W(2);
			B(1)=BETA(1)+W(2);
		for j=1:N
			K=K+1;
			VV=V(K)/V(K-1);
			P(K)=asin(VV*sin(A(K-1))); 
			Q(K)=asin(VV*sin(B(K-1))); 
			if (K+1-M) <0
				A(K)=P(K)-W(K+1)+W(K);
				B(K)=Q(K)+W(K+1)-W(K);
				ALPH(K)=A(K)+W(K+1); 
				BETA(K)=B(K)-W(K+1);
			else
				A(K)=(P(K)+Q(K))*0.5;
				B(K)=A(K);
				W(K+1)=W(K)+(P(K)-Q(K))*0.5; 
				ALPH(K)=A(K)+W(K+1); 
				BETA(K)=B(K)-W(K+1); 
				V(K+1)=V(K)/sin(A(K));
			break
			end
		end
	end
			KK=K-1;
			HHA=0;
			HHB=0;
			if (KK) >0
			for  I=1:KK
				HH=cos(ALPH(I))+cos(BETA(I)); 
				HH=HH/V(I);
				HHA=HHA+HH*HA(I);
				HHB=HHB+HH*HB(I);
			end
			end
		R=V(K)/(cos(ALPH(K))+cos(BETA(K)));
		HA(K)=R*(TAI(K+1)-HHA); 
		HB(K)=R*(TBI(K+1)-HHB); 
		DA(1)=HA(1);
		DB(1)=HB(1);
		if (K-1) >0
			DA(K)=DA(K-1)+HA(K);
			DB(K)=DB(K-1)+HB(K);
		end
%
end		% End of Large Loop
%
		for J=2:N
	 		W(J)=W(J)*57.2958+0.001;
		end
		HB(N)=0;
		HA(N)=0;
		DB(N)=0;
		DA(N)=0;
		%
%	Print results
%
[TITLE]
['    Layer    Velocity  Thickness Thickness  Dip (¡)  Depth A   Depth B  ', 
 '   Number    of Layer     at A      at B    at top  (to bottom of Layer)']
for i=1:N
	AC(i,:)=[i V(i) HA(i) HB(i) W(i) DA(i) DB(i)];
end
AC
%[V HA HB W DA DB];
%
%   Plot Model
%
NL=N-1;	% Number of lines
XA=0;
XB=X;
hold on
for i=1:NL
	z1=DA(i);
	z2=DB(i);
	x=[XA XB];
	y=[z1 z2];
	plot(x,y,'r-')
end
x=[0 XB];
y=[0 0];
plot(x,y,'r-')
axis('ij')
axis([0 X 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