%% ==========================================================================
%% GEOL335 Lab 3 by .... (write your name, NSID)
%% ... write date ....
%% ==========================================================================

%% Use comments like this to explain the purposes of various steps, 
%% meanings of variables, and ideas of approaches
%% ------------------------------------------------------------------

%% in Matlab, you can extract data columns from Excel file like shown below
%% (uncomment the following lines):

%data = readtable("lab3_table.xlsx")

%data        % look at the structure 'data'

	% Note that from the given table, Matlab will have NaN (not-a-number)
	% values loaded from the second line of the spreadsheet.
	% these spurious values can be skipped as shown below:
	
%x = data.Distance(2:end);
%z = data.Elevation(2:end);

%tA = data.ShotA(2:end);
%tB = data.ShotB(2:end);
%tC = data.ShotC(2:end);
%tD = data.ShotD(2:end);

%% in Octave (with package 'io' loaded), numerical data from same file can be loaded
%%  

data = xlsread("lab3_table.xlsx");	% this 'data' is a simple matrix

x = data(:,1);
z = data(:,2);

tA = data(:,3);
tB = data(:,4);
tC = data(:,5);
tD = data(:,6);

% specify X locations of shots

xA = 0;       
xB = 84;
xC = 180;
xD = 276;

% extract reciprocal times for pairs of shots
% check in the printouts that the reciprocal times are equal

trAD = tA(find(x==xD))  % time from A at the position of shot D
trDA = tD(find(x==xA))  % vice versa. This time should be the same (reciprocal time)

% ... add more shot pairs ...

%% Task 1: plotting measured times
%% ============================================================================

	% create a figure with two subplots showing travel times and elevation
	
figure(1);	
clf

subplot(2,1,1);
hold on

plot(x,tA,'-+r')
plot(x,tB,'-ob')
plot(x,tC,'-+g')
plot(x,tD,'-ok')

title('Trave-time data')
ylabel('Travel time [s]')
%xlabel('x [m]')

elevation_subplot = subplot(2,1,2);	% note that you can save the 'handle' of the subplot
					% and use it later in task 7

plot(x,z)
title('Elevation')
xlabel('x [m]')
ylabel('z [m]')

%% Task 2: Extract forward and backward (reversed) headwave segments
%% =========================================================================

	% it is convenient to use boolean arrays (consisting of values 0 and 1) 
	% for selecting the various subsets of data
	% values = 1 mean the location is selected and 0 means it is not selected
	
hwfA = x > 8;	% boolean vector for selecting forward headwave from shot A
hwbB = x < 65;	% headwave backward from B
hwfB = x > 100;	% headwave forward from B

hwbC = x < 150;  % headwave backward from C
...

hwbD = x < 255;

...

	% direct waves from shot C
	
dwfC = x >= xC & x < 225;	% forward
%dwbC = x <= xC & ...;		% similarly define the backward branch

ind_dwfC = find(dwfC);		% corresponding row indices in data table
%ind_dwbC = ... 

offset_dwfC = x(ind_dwfC) - xC;	% signed source-receiver offsets for these waves

% ..... add similar for backward direct-wave branch ...

%% Task 3: Measure V1 from shot C 
%% =========================================================================

	% specify the range of direct waves from shot C
	
ind_dwC = find( x > 150 & x < 225 );	% check this! 


					
	% plot direct wave times given by dwfC and dwbC 
	% versus ABSOLUTE VALUEs of source-receiver OFFSETS:
	
figure(3)
clf

hold on

%plot(abs(x(ind_dwfC)),tC...)	% complete this line and add line for backward wave!
				% comment whether they are close?

title("Direct-wave travel times")
xlabel('abs(offset) [m]')
ylabel('travel time [s]')

	% obtain V1 by least-squares fitting
	% in the following line, see how 'polyfit' works 
	% and modify it to include BOTH the forward and backward branches.
	% you can combine branches using indices again: ind_C = [ ind_dwC, ind_dfC ]
	
p1 = polyfit(offset_dwfC,tC(ind_dwfC),1); % least-squares fitting by a linear function

p1	% This should be a 2-element array. 
	% The first element is the slowness, the second is the intercept
	
V1 = 1.0/p1(1)		% the desired velocity

%% Task 4: Extracting reversed head wave segments
%% =========================================================================

	% Use function find(...) to obtain indices of the overlappng segments 
	% of forward and backward headwaves
	
ind_AD = find(hwfA & hwbD);	% overlap of forward hw A and backward hw D
ind_BD = find(hwfB & hwbD);	% B and D
ind_AC = find(hwfA & hwbC);	% A and C

% ... add other pairs ...

	% plot all segments on one plot
	% Evaluate the coverage?
	
figure(4);

clf
hold on

plot(x(ind_AD),tA(ind_AD),'-r')
plot(x(ind_AD),tD(ind_AD),'-g')

plot(x(ind_BD),tB(ind_BD),'-m')
plot(x(ind_BD),tD(ind_BD),'-b')

plot(x(ind_AC),tA(ind_AC),'-k')
plot(x(ind_AC),tC(ind_AC),'-k')

title("Selected head-wave travel times")
xlabel('x [m]')
ylabel('Time [s]')

%% ... other pairs?

%% Tasks 5 and 6: Plus and Minus travel times
%% =========================================================================

	% for the segments in Task 4, 
	% calculate delay times (half of the PLUS function minus reciprocal time)
	% and the MINUS time functions
	
delay_AD = 0.5*(tA(ind_AD)+tD(ind_AD) - trAD);
tminus_AD = tA(ind_AD)-tD(ind_AD);

% ... add pairs of times for other shots ...

	% plot subplots of delay times and the minus-time analysis functions
	
figure(41)      % plus-minus times in subplots
clf

subplot(2,1,1)
hold on

plot(x(ind_AD),delay_AD,'-r')

% ... add other curves ...

ylabel('Delay time [s]')
xlabel('x [m]')
set(gca,'ydir','reverse')	% this reverses the delay-time axis downward, so that
				% it resembles depth

subplot(2,1,2)
hold on

plot(x(ind_AD),tminus_AD,'-r')

% ... add other ...

ylabel('Minus time [s]')
xlabel('x [m]')

	% in the MINUS curve, select a 2-3 intervals of X
	% and measure the slope p2 and velocity V2 at each of these points,
	% similar to Task 3 above:
	
V21 = 1000.0;	% I've put these values for testing, but you replace them with yur results
V22 = 1200.0;

x1 = 75.0; x2 = 200.0;	% middles of intervals at which V21 and V22 are measured

	% from the above values, interpolate V2 at every point in x:
	
V2 = interp1([x1,x2],[V21,V22],x,'linear','extrap');

%% Task 7: Depth and elevation of the refractor
%% =========================================================================

	% calculate the critical angle at each point of the profile.
	% Note that 'dot' operations are needed for calculations with each element
	% of the arrays:
	
i_crit = asin(V1*V2.^(-1));	% angles in radians

	% calculate depth to the refractor for each delay-time curve above:
	
depth_AD = delay_AD .* V1 ./ cos(i_crit(ind_AD));

% do the same for other delay curves ...

	% plot the inverted refractor positions on the elevations plot in figure(1) 

figure(1);
subplot(elevation_subplot);
hold on

plot(x(ind_AD),z(ind_AD)-depth_AD,'-r')

% .... add estimates from other pairs of headwave travel-time curves ...