In ocean acoustics a common task is to model the propagation of sound through the water column. The speed of sound varies with temperature, salinity and pressure (depth). To a first approximation, the ocean is horizontally stratified and these parameters only vary with depth. Variations in sound speed in layers of the water column cause a "ray" of acoustic energy to be refracted (to bend) as the ray passes through it.
When given a starting location, launch angle, vertical sound speed profile and total travel time, the attached function, raytrace.m, models the propagation of sound though a water column. It handles both reflection from the surface of the water and the bottom (the deepest depth given in the sound speed profile) as well as caustics (places in which the ray is bent back in the direction from which it came). Raytrace automatically plots the provided sound speed profile and the results of the raytrace.
Although raytrace.m is designed for ocean acoustics, the mathematics for optical refraction is the same and where the assumptions apply raytrace.m will produce reasonable results.
Val Schmidt (2021). raytrace (https://www.mathworks.com/matlabcentral/fileexchange/26253raytrace), MATLAB Central File Exchange. Retrieved .
Inspired: Water sound speed calculator
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!Create scripts with code, output, and formatted text in a single executable document.
Hello Val,
i would really like to try out your code witht he provided example, but it keeps telling me
that my maximum recursion limit is reached. If i increase the limit, the programm crashes and
by reducing the maximum time and/or the bounces it still gives me this error.
Could you give me a hint on where my mistake lies.
Regards,
Basti
Hi do this work for doubly curved lens??
Hello All,
As a Matlab novice I am having trouble getting this to work for both my own collected SSP data and the example Val has kindly provided, and feel like I am not understanding a basic element of the process.
The error message I receive is:
 Not enough input arguments.
Error in ETX_raytrace_setup (line 136)
if zz(1) ~= 0
As I've said this occurs for both the example and my own data. I have tried to set up my variables as follows:
zz = 1.5:.25:12.5
cc = identical array size of corresponding SSP values at each depth
theta0 = 70
xo = 0
zo = 1.5
tt = 0.0243
Any help or guidance to correct usage would be much appreciated.
Many Thanks,
Mark
Great program!
A tip when using this code...
I had problems while using a real SVP from a CZ environment. Results were sporadic.
I had to upsample the SVP to 1 meter resolution, and then it worked fine.
Example:
load svp.txt
zsvp=svp(:,1)';
csvp=svp(:,2)';
zz = 0:1:max(zsvp);
cc=interp1(zsvp,csvp,zz,'linear');
svp.txt is ascii text of depth,speed pairs in tabdelimited column format.
Hope this helps.
Another tip when using this code...
The code will take a vector of angles and plot the rays for all of them at the same time.
Example:
instead of theta0=0;
use
theta0=[40:5:5];
Hope this helps!
Thanks Yoann! And Thanks for the code post. Always interesting to see another way to do it!
Hi Val,
Thanks for that, it works fine for me. Gives very similar results to another way of modelling it, which is always reassuring!
function [r_ray,t_ray,z_ray,beta_ray]=compute_acoustic_ray_path(z,c,r_ray,t_ray,z_ray,beta_ray,T)
[z,ia,ib] = unique(z);
if size(z,1)==1
c=c(ia);
else
c=c(ib);
end
[z,idx_sort]=sort(z);
c=c(idx_sort);
g_down=(c(2:end)c(1:end1))./(z(2:end)z(1:end1));
g_down=[g_down g_down(end)];
g_up=(c(1:end1)c(2:end))./(z(1:end1)z(2:end));
g_up=[g_up(1) g_up];
if beta_ray>0
g=g_down;
else
g=g_up;
end
xi=(cos(beta_ray(1))/c(1));
while t_ray(end)<=T
%t_ray(end)
[~,idx]=nanmin(abs(z_ray(end)z));
i=length(r_ray);
if beta_ray(i)<0
idx_plus=nanmax(idx1,1);
elseif beta_ray(i)>0
idx_plus=nanmin(idx+1,length(z));
else
idx_plus=idx;
end
if idx==idx_plus&&idx~=1
return;
end
if xi==0 %case where the ray is propagating vertically
beta_ray(i+1)=beta_ray(i);
r_ray(i+1)=r_ray(idx);
t_ray(i+1)=t_ray(i)+1/(g(idx))*log(c(idx_plus)/c(idx)*(1+sqrt(1xi^2*c(idx)^2))/(1+sqrt(1xi^2*c(idx_plus)^2)));
z_ray(i+1)=z(idx_plus);
else
if idx_plus==1&&idx==1 %case where it hits the surface
if g(idx)==0%case where velocity stays the same
beta_ray(i+1)=beta_ray(i);
r_ray(i+1)=r_ray(i)+(z(2)z(1)/tan(beta_ray(i+1)));
t_ray(i+1)=t_ray(i)+((z(2)z(1))/c(1));
z_ray(i+1)=z(1);
if beta_ray(i+1)>0
g=g_down;
else
g=g_up;
end
else
beta_ray(i+1)=pi/2(beta_ray(i)pi/2);
r_ray(i+1)=r_ray(i)+(r_ray(i)r_ray(i1));
t_ray(i+1)=t_ray(i)+(t_ray(i)t_ray(i1));
z_ray(i+1)=z(1);
if beta_ray(i+1)>0
g=g_down;
else
g=g_up;
end
end
else
if g(idx)==0%case where velocity stays the same
beta_ray(i+1)=beta_ray(i);
r_ray(i+1)=r_ray(i)+((z(idx_plus)z(idx))/tan(beta_ray(i+1)));
t_ray(i+1)=t_ray(i)+((z(idx_plus)z(idx))/c(idx));
z_ray(i+1)=z(idx_plus);
else
if (xi*c(idx_plus))^2<1
beta_ray(i+1)=sign(beta_ray(i))*acos(c(idx_plus)*xi);
r_ray(i+1)=r_ray(i)+(1/(xi*g(idx))*(sqrt(1xi^2*c(idx)^2)sqrt(1xi^2*c(idx_plus)^2)));
t_ray(i+1)=t_ray(i)+1/(g(idx))*log(c(idx_plus)/c(idx)*(1+sqrt(1xi^2*c(idx)^2))/(1+sqrt(1xi^2*c(idx_plus)^2)));
z_ray(i+1)=z(idx_plus);
else %case where the ray got to the critical angle
beta_ray(i+1)=beta_ray(i);
xi=(cos(beta_ray(i+1))/c(idx));
r_ray(i+1)=r_ray(i)+2/(xi*g(idx))*(sqrt(1xi^2*c(idx)^2));
t_ray(i+1)=t_ray(i)+2/(g(idx))*log(1/(c(idx)*abs(xi))*(1+sqrt(1xi^2*c(idx)^2)));
z_ray(i+1)=z(nanmax(idx1,1));
if beta_ray(i+1)>0
g=g_down;
else
g=g_up;
end
end
end
end
end
end
end
If you're having trouble start with the example you'll find commented out in the code and reproduced here:
zz = 0:1:5000;
cc = 1520 + [zz*.05];
cc(751:end) = cc(750) + (zz(751:end)zz(750))*.014;
% Conduct the raytrace and plot result
[x z t d] = raytrace(0,0,5,120,zz,cc,true);
Compare this to what you've provided in your own work. Are you specifying the same number of arguments? Are there the same number depths as sound speed values in your profile? etc.
If you get an error "Warning, not enough bounces..." this means that you've specified a really long travel time for the water depth and the code is predicting many reflections (more than 10). This may reflections almost certainly reflects a mistake on your part and you should shorten your travel time to something reasonable. However, if it is really what you want to do, you can edit the code to increase the number of bounces (this is a variable).
I get this message when I run the code:
Error using raytrace (line 136)
Not enough input arguments.
I am new to matlab. Please tell a solution.
Thank you.
This is absolutely amazing. It's not only great implementation but also the best explanation for this subject. If only it contained transmission loss it'd be even more great.
As a beginner, I got a lot from your work. thank you very much.
Getting back to say that it worked!!! Thank a lot!!!!
I have trouble making this work.
I am new to matlab and don't know what i am doing wrong. i am making a raytrace.m that has my speed profile but when i run the function this error appears:
function [xxf zzf ttf ddf] = raytrace(xo,zo, theta0,tt,zz,cc,plotflag)

Error: Function definitions are not permitted in this context.
I have had trouble making this work with a fairly standard 13 point ssp:
Warning: Not enough bounces specified for time requested.
> In raytrace at 261
I'm calling it repeatedly for launch angles 2 to 2 in steps of 0.5. I changed diff(z{:}) to diff(z) and the same for c because it didn't seem to be giving the right result. dt doesn't seem to be changing from the value zero so does this mean the time isn't advancing?
Very very useful. Thanks a lot, Val. :)
Thanks, this is fantastic.