% simplerays.m
%
% Simple acoustical ray tracer for piecewise-constant sound speed profile.
%
% Wayne H. Slade, 6-Oct-2005

function simplerays


% sound speed vs depth (e.g. 0-100 is 1530, 100-400 is 1490, ... 800-1200
% is 1540, bottom is at 1200) both bottom and air-water interface are
% assumed to be totally reflecting; depth in meters, velocity in meters per
% second

layer_z = [100   400  800 1200]
layer_c = [1530 1490 1500 1540]

% ray source depths and initial zenith angles (first column is depth in
% meters, second column is zenith angle in degrees, positive clockwise) 
rays0 = [100 80;
        100 90;
        100 100]
    
% specify range, in meters
range = 40e3
    
% BOOM! Trace each ray from time of the explosion...
% each matrix in cell array raypaths contains 5 columns, where each row is
% a step in the ray path; first column is step number, second and third are
% the depth and range for the start of the each step, fourth is the zenith angle
% of each step, and the fifth column is the amount it took for the
% acoustical energy to travel each step
for kr = 1:size(rays0,1)
    raypaths{kr} = traceray(layer_z, layer_c, [0 rays0(kr,:)], range);
end

% plot results-- depth (meters) vs. range (kilometers)
linecolors = ['b' 'r' 'k'];
figure, hold on
for kr = 1:length(raypaths)
    plot(raypaths{kr}(:,2)/1000, raypaths{kr}(:,3), [linecolors(kr) '.-'], 'linewidth', 1.5)
end
set(gca, 'box','on', 'fontsize', 16, 'ydir','reverse')
xlabel('Range [km]', 'fontsize', 16)
ylabel('Depth [m]', 'fontsize', 16)


% Add up travel time for each ray--
for kr = 1:length(raypaths)
    disp(['Travel time, ray ' num2str(kr) ' (' num2str(rays0(kr,1)) 'm,' num2str(rays0(kr,2)) 'deg), t = ' num2str(sum(raypaths{kr}(:,5))) ...
        ' seconds, received at z = ' num2str(num2str(raypaths{kr}(end,3)))])
end


keyboard

return

% layer_z, layer_c - vectors of sound speed vs depth
% ray0 - [source_depth, zenith_angle]
% range - maximum range
function ray = traceray(layer_z, layer_c, ray0, range) 

    ray = [1 ray0 0];
    rayd = []; % debug info for each step
    
    ks = 1;

    while 1
            
        ray_x1 = ray(ks,2);
        ray_z1 = ray(ks,3);
        ray_theta1 = ray(ks,4);
        ray_c1 = 0;
        ray_x2 = 0;
        ray_z2 = 0;
        ray_theta2 = 0;
        ray_c2 = 0;

        % (A) ray is headed towards shallower or air-water interface
        % theta1 (zenith) is < 90
        if (ray_theta1>0) & (ray_theta1<90)
            layeridx = find(layer_z<ray_z1, 1,'last');
            if isempty(layeridx)
                % ray is being REFLECTED at air-sea interface
                ray_z2 = 0;
                ray_theta2 = 180 - ray_theta1;
                ray_x2 = ray_x1 + ray_z1/tand(90-ray_theta1);
                ray_c1 = layer_c(1);
            else
                ray_z2 = layer_z(layeridx); % depth of interface we're hitting
                ray_c1 = layer_c(layeridx+1); % headed up, so c1 at layeridx+1
                ray_c2 = layer_c(layeridx); % headed up, so c2 at layeridx
                ray_x2 = ray_x1 + (ray_z1-ray_z2)/tand(90-ray_theta1);

                if (ray_c2/ray_c1)*sind(ray_theta1) > 1
                    % ray is REFLECTED at density interface
                    ray_theta2 = 180 - ray_theta1;
                else
                    % ray is being REFRACTED at density interface
                    ray_theta2 = asind((ray_c2/ray_c1)*sind(ray_theta1));
                end

            end
        end

        % (B) ray is headed towards deeper interface or bottom
        % theta1 (zenith) is < 90
        % depth at end of ray segment is greater than depth at start of ray
        % segment
        if (ray_theta1>90) & (ray_theta1<180)
            layeridx = find(layer_z(1:(end-1))>ray_z1, 1,'first');
            if isempty(layeridx)
                % ray is being REFLECTED at bottom
                ray_z2 = layer_z(end);
                ray_theta2 = 180 - ray_theta1;
                ray_x2 = ray_x1 + (ray_z2-ray_z1)/tand(ray_theta1-90);
                ray_c1 = layer_c(end);
            else
                ray_z2 = layer_z(layeridx); % depth of interface we're hitting
                ray_c1 = layer_c(layeridx); % headed down, so c1 at layeridx
                ray_c2 = layer_c(layeridx+1); % headed down, so c2 at layeridx+1
                ray_x2 = ray_x1 + (ray_z2-ray_z1)/tand(ray_theta1-90);

                if (ray_c2/ray_c1)*sind(ray_theta1) > 1
                    % ray is being REFLECTED at density interface
                    ray_theta2 = 180 - ray_theta1;
                else
                    % ray is being REFRACTED at density interface
                    ray_theta2 = 180 - asind((ray_c2/ray_c1)*sind(ray_theta1));
                end
            end

        end

        if (ray_theta1==90)
            % ray is horizontal
            ray_x2 = range;
            ray_z2 = ray_z1;
            ray_theta2 = 90;
            ray_c1 = layer_c(find(layer_z(1:(end-1))>ray_z1, 1,'first'))
        end
        
        deltat = sqrt((ray_x2-ray_x1)^2 + (ray_z2-ray_z1)^2)/ray_c1;
        
        % debug info for ray
        rayd = [rayd; ks 0 ray_x1/1000 ray_z1 ray_theta1 0 ray_c1 ray_c2 0 ray_x2/1000 ray_z2 ray_theta2];

        ray(ks,5) = deltat;
        ray(ks+1,:) = [ks+1 ray_x2 ray_z2 ray_theta2 0]; % add new ray

        if ray_x2>=range
            break
        end
        
        ks = ks + 1;
        
    end

    % shorten the last ray step to end at x=range
    if ray(end,2)>range
        delz = (range-ray(end-1,2))*tand(ray(end-1,4)-90)
        ray(end,2) = range;
        ray(end,3) = ray(end-1,3) + delz;
        ray(end,4) = 0;
        ray(end,5) = 0;
        ray_c1 = layer_c(find(layer_z(1:(end-1))>ray_z1, 1,'first'));
        deltat = sqrt((ray(end,2)-ray(end-1,2))^2 + (ray(end,3)-ray(end-1,3))^2)/ray_c1;
        ray(end-1,5) = deltat;
    end
    
return

