clear all
prem=load('prem.txt');
p=input('Input ray parameter: ');
z=prem(:,1);
Vp=prem(:,4);
nprem=size(z,1);
Cx=1/p;
np=0;
for i=1:nprem-1
    u(i)=2/(Vp(i)+Vp(i+1));
    if Vp(i) <= Cx && Vp(i+1) > Cx && np == 0
        np=i;
        zp=z(i)+(Cx-Vp(i))*(z(i+1)-z(i))/(Vp(i+1)-Vp(i));
        u(i)=2/(Vp(i)+Cx);
    end
end
x=0;
t=0;
for i=1:np-1
    x=x+2*p*(z(i+1)-z(i))/sqrt(u(i)^2-p^2);
    t=t+2*u(i)^2*(z(i+1)-z(i))/sqrt(u(i)^2-p^2);
end
x=x+2*p*(zp-z(np))/sqrt(u(np)^2-p^2)
t=t+2*u(np)^2*(zp-z(np))/sqrt(u(np)^2-p^2)
