exec('vels.sci');

xf=20; z0=0.6;

xset('window', 0);
for i=-10:2.5:10                                                  
  theta=i*%pi/180; A=cos(theta)/depths(z0);                    
  deff('yprim=acous(t, y)', 'yprim=[y(2); -depthsp(y(1))/(A^2 * depths(y(1))^3)]')
  z=ode([z0; tan(theta)], 0, 0:1:60, acous);               
  plot2d(0:1:60, -z(1,:))                                     
end        

function zf=ray(theta)
theta=theta*%pi/180; A=cos(theta)/depths(z0);
deff('yprim=acous(t, y)', 'yprim=[y(2); -depthsp(y(1))/(A^2 * depths(y(1))^3)]');
z=ode([z0; tan(theta)], 0, xf, acous);
zf=z(1,$)-0.9;
endfunction  

xset('window',1)
fplot2d((-10:.2:10), ray) 
[sol, v,info]= fsolve(3,ray)