% ******************************************************
% Rays
% ******************************************************

% The equations we're solving are:
% r'    = c rho
% z'    = c zeta
% rho'  = -c_r / c2
% zeta' = -c_z / c2

send = 100000;		% arclength for rays

ntheta   = 31;		% number of rays
theta = pi / 180 * linspace( -14.0, 14.0, ntheta );

zs = 1000.0;		% source depth
c0 = 1501.4;		% sound speed at source depth

for ith = 1:ntheta	% loop over take-off angle

   % ray initial condition:
   x0 = [ 0.0   zs   cos( theta( ith ) ) / c0   sin( theta( ith ) ) / c0 ];

   % now solve the DE to trace the ray
   [ s, x ] = ode45( 'rayf', 0.0, send, x0 );

   plot( x( : , 1 ), x( : , 2 ) );
   hold on;	% hold the old rays on screen when plotting new rays
end
hold off;

% label the plot

xlabel( 'Range (m)' )
ylabel( 'Depth (m)' )
view( 0, -90 );	% flip plot so that z-axis is pointing down