function [dx,dt,irtr] = layerxt(p,h,utop,ubot)

// After Peter Shearer's fortran routine LAYERXT
//

// inputs :
// p = ray parameter
// h = layer thickness
// utop= slowness at top of layer
// ubot= slowness at bottom of layer
// output :
// dx = change in distance
// dt = change in travel time
// irtr = return code
// = -1, zero thickness layer
// = 0, ray turned above layer
// = 1, ray passed through layer
// = 2, ray turned in layer, 1 leg counted in dt and dx.
//
// method:
// This code assumes a linear velocity gradient between depth points and
// calculates the one way travel time and distance increments using
// equations (4.21) and (4.22) on page 41 of Introduction to Seismology by Shearer.
// When a layer has a constant velocity, it uses equations (4.18) and (4.19) on page 40.



if (p >= utop) then // ray turned above layer
    dx = 0.;
    dt = 0.;
    irtr = 0;
    return;
elseif (h == 0.) then // zero thickness layer
    dx = 0.;
    dt = 0.;
    irtr = -1;
    return;
end

u1 = utop;
u2 = ubot;
v1 = 1.0/u1;
v2 = 1.0/u2;
b = (v2-v1)/h; // the slope of the velocity gradient

eta1 = sqrt(u1*u1 - p*p);

if (b == 0.0) then // constant velocity layer
    dx = (h*p)/eta1;
    dt = (h*u1*u1)/eta1;
    irtr = 1;
    return
end

x1 = eta1/(u1*b*p);
tau1 = (log((u1 + eta1)/p) - eta1/u1) / b;

if (p >= ubot) then
    dx = x1;
    dtau = tau1;
    dt = dtau + p*dx;
    irtr = 2;
    return
end

irtr = 1;

eta2 = sqrt(u2*u2 - p*p);
x2 = eta2/(u2*b*p);
tau2=(log(( u2 + eta2)/p) - eta2/u2) / b;
dx = x1 - x2;
dtau = tau1 - tau2;
dt = dtau+p*dx;

return

endfunction