function refraction_3layers(v1, v2, v3, z1, z2, FIRST_ARRIVALS_ONLY)

%function refraction_3layers(v1, v2, v3, z1, z2, FIRST_ARRIVALS_ONLY)
%
% Creates travel time plots for a 3-layer system with layer velocities
% v1,  v2 and v3, and layer thicknesses z1 and z2.
% The FIRST_ARRIVALS_ONLY FLAG may be set to 1 to plot only the first arrivals. By 
% default, all arrivals are plotted.


if nargin < 6
    FIRST_ARRIVALS_ONLY = 0;
end

% X-positions for geophones (m, relative to source)
x = [0:5:300];    

% Direct wave 
% Travels along ground surface (at velocity v1)
t1 = x./v1;

% Head wave 
% Refracts along z1-z2 boundary
% Travel time depends on velocities of both layers
% and thickness of layer 1 only (thickness of layer 2 irrelevant)
t2 = (2*z1*sqrt(v2^2-v1^2)/(v1*v2))+x./v2;   %% Note slope should be 1/v2!

%Find critical distance...before this value, no head waves arrive at
%surface.
xcrit = 2*z1*v1/(sqrt(v2^2-v1^2));

if isreal(xcrit)
    a = min(find(x>xcrit));
end

% Head wave for layer 3 
% Refracts along z1-z2 boundary
% Travel time depends on velocities of all 3 layers
% and thicknesses of layers 1 and 2.
t3 = (2*z1*sqrt(v3^2-v1^2)/(v1*v3))+(2*z2*sqrt(v3^2-v2^2)/(v2*v3))+ x./v3;

%Find critical distance...before this value, no head waves arrive at
%surface.
xcrit2 = 2*(((z1*v1)/sqrt(v3^2-v1^2))+((z2*v2)/sqrt(v3^2-v2^2)));
b = min(find(x>xcrit2));

crossover1 = ((2*z1*sqrt(v2^2-v1^2))/(v1*v2))/(1/v1-1/v2);
crossover2 = ((2*z1*sqrt(v3^2-v1^1))/(v1*v3)+(2*z2*sqrt(v3^2-v2^2))/(v3*v2)-(2*z1*sqrt(v2^2-v1^2))/(v1*v2))/(1/v2-1/v3)

if isreal(crossover1)
    if FIRST_ARRIVALS_ONLY
        c = max(find(x <= crossover1))
        d = max(find(x <= crossover2))
        plot(x(1:c),t1(1:c)*1000, '.--')
        hold on
        if isreal(t2)
            plot(x(c:d), t2(c:d)*1000, 'r.--')
        end
        plot(x(d:end), t3(d:end)*1000, 'g.--')
    else
        plot(x,t1*1000, '.--')
        hold on
        if isreal(t2)
            plot(x(a:end), t2(a:end)*1000, 'r.--')
        end
        plot(x(b:end), t3(b:end)*1000, 'g.--')
    end
else
    newcross = ((2*z2*sqrt(v3^2-v2^1))/(v2*v3)+(2*z1*sqrt(v3^2-v1^2))/(v3*v1))/(1/v1-1/v3);
    if FIRST_ARRIVALS_ONLY
        c = max(find(x <= newcross));
        plot(x(1:c),t1(1:c)*1000, '.--')
        hold on
        plot(x(c:end), t3(c:end)*1000, 'g.--')
    else
        plot(x,t1*1000, '.--')
        hold on
        if isreal(t2)
            plot(x(a:end), t2(a:end)*1000, 'r.--')
        end
        plot(x(b:end), t3(b:end)*1000, 'g.--')
    end
end


xlabel('GEOPHONE OFFSET (m)')
ylabel('TIME (ms)')
grid
if isreal(t2)
    legend('DIRECT WAVE', 'FIRST HEAD WAVE', 'SECOND HEAD WAVE', 'Location', 'SouthEast')
else
    legend('DIRECT WAVE', 'SECOND HEAD WAVE', 'Location', 'SouthEast')
end

    title(['z1 = ', num2str(z1), ' m; z2 = ', num2str(z2),...
    ' m; v1 = ', num2str(v1), ' m/s; v2 = ', num2str(v2),...
    ' m/s; v3 = ', num2str(v3), ' m/s'])       
axis ([0 300 0 300])
hold off