File Exchange

image thumbnail

raytrace

version 1.6.0.0 (5.63 KB) by Val Schmidt
A function to model acoustic refraction due to variations in sound speed.

12 Downloads

Updated 12 Feb 2014

View Version History

View License

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 ray-trace.

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.

Cite As

Val Schmidt (2021). raytrace (https://www.mathworks.com/matlabcentral/fileexchange/26253-raytrace), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (18)

Sebastian

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

alex

Hi do this work for doubly curved lens??

mark gray

mark gray

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

Walter Allensworth

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 tab-delimited column format.

Hope this helps.

Walter Allensworth

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!

Val Schmidt

Thanks Yoann! And Thanks for the code post. Always interesting to see another way to do it!

Yoann Ladroit

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:end-1))./(z(2:end)-z(1:end-1));
g_down=[g_down g_down(end)];

g_up=-(c(1:end-1)-c(2:end))./(z(1:end-1)-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(idx-1,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(1-xi^2*c(idx)^2))/(1+sqrt(1-xi^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(i-1));
t_ray(i+1)=t_ray(i)+(t_ray(i)-t_ray(i-1));
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(1-xi^2*c(idx)^2)-sqrt(1-xi^2*c(idx_plus)^2)));
t_ray(i+1)=t_ray(i)+1/(g(idx))*log(c(idx_plus)/c(idx)*(1+sqrt(1-xi^2*c(idx)^2))/(1+sqrt(1-xi^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(1-xi^2*c(idx)^2));
t_ray(i+1)=t_ray(i)+2/(g(idx))*log(1/(c(idx)*abs(xi))*(1+sqrt(1-xi^2*c(idx)^2)));
z_ray(i+1)=z(nanmax(idx-1,1));

if beta_ray(i+1)>0
g=g_down;
else
g=g_up;
end
end
end
end

end


end

end

Val Schmidt

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).

Haree V

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.

kubas121

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.

Chen zhifei

As a beginner, I got a lot from your work. thank you very much.

anastasia

Getting back to say that it worked!!! Thank a lot!!!!

anastasia

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.

anastasia

Geoff

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?

Saurabh Singh

Very very useful. Thanks a lot, Val. :)

Prakash Manandhar

Thanks, this is fantastic.

MATLAB Release Compatibility
Created with R2009a
Compatible with any release
Platform Compatibility
Windows macOS Linux
Acknowledgements

Inspired: Water sound speed calculator

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!