function refraction_2layers(v1, v2, z, FIRST_ARRIVALS_ONLY);

% refraction_2layers(v1, v2, z);
%
% Creates travel time plots for a two-layers system with layer velocities
% v1 and v2 and layer 1 thickness of z
% The FIRST_ARRIVALS_ONLY FLAG may be set to 1 to plot only the first arrivals. By 
% default, all arrivals are plotted.

if nargin < 4
    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*z*sqrt(v2^2-v1^2)/(v1*v2))+x./v2;   %% Note slope should be 1/v2!

xcrit = 2*z*v1/(sqrt(v2^2-v1^2));
if isreal(xcrit)
    a = min(find(x>xcrit));
end

crossover = ((2*z*sqrt(v2^2-v1^2))/(v1*v2))/(1/v1-1/v2);

b = max(find(x<= crossover))

if FIRST_ARRIVALS_ONLY
    plot(x(1:b),t1(1:b)*1000, '.--')
    hold on
    if isreal(t2)
        plot(x(b:end), t2(b:end)*1000, 'r.--')
    end
else
    plot(x,t1*1000, '.--')
    hold on
    if isreal(t2)
        plot(x(a:end), t2(a:end)*1000, 'r.--')
    end
end

    xlabel('GEOPHONE OFFSET (m)')
    ylabel('TIME (ms)')
    grid on 
    legend('DIRECT WAVE', 'HEAD WAVE')
    title(['z1 = ', num2str(z), ' m; v1 = ', num2str(v1), ' m/s; v2 = ', num2str(v2), ' m/s'])
    axis ([0 300 0 300])
    hold off