Clear Filters
Clear Filters

Nested for loop problem

1 view (last 30 days)
JDilla
JDilla on 12 Mar 2015
Edited: Guillaume on 12 Mar 2015
I am modelling seismic rays, (I am very new to Matlab, this is my first attempt at solving a problem), and I have to plot a series of seismic rays at various angles (any number I choose up to 20 degrees) that penetrate through 5 layers, and reflect at the boundary of the 4th and 5th layer. Through each layer, the seismic ray needs to refract. I have managed to achieve this with the code below.
HOWEVER, I need to plot a SERIES of rays, I have only managed one. I wanted a ray to propagate at 1 degree, increments of 5 degrees up to i=20 (the initial angle you see plotted). I want to use a nested for loop to achieve this, but I'm not sure where i'm going. Ultimately, I'm trying to achieve a single plot of several seismic rays with their reflections such as below, at various angles, with the maximum angle being 20.
A plot similar to this picture is what I'm trying to achieve; multiple rays and their reflections. The code for the reflections is
reflectionx=h(end);
plot(2*reflectionx-h,depth)
axis ij
%COMPUTE DISTANCE TRAVELLED BY THE RAY IN EACH LAYER
clear
%DEFINE VARIABLES
z1 = 0; %Depth at y=0
v = [3:7]; %Velocities at each subsequent layer
i = 20; %Takeoff Angle in degrees
%COMPUTE DIAGONAL DISTANCE BY THE RAY IN EACH LAYER
for j = 1:5:i
for n = 1:5; %n = number of layers
v = v(n);
i;
d=z1/cosd(i); %calculate diagonal distance
h(n)=z1*sind(i) %calculate horizontal distance
t(n) = d/v;
v = [3:7];
depth(n)=z1
sini = ((sind(i)*(v+1))/v);
i= asind(sini)
z1=z1+2;
end
end
HorizontalDistance=[cumsum(h)];
disp(HorizontalDistance);
plot(h,depth) %Plot horizontal distance vs depth
hold on
reflectionx=h(end);
plot(2*reflectionx-h,depth) %Plot reflected ray
axis ij %Places origin at upper left corner of plot
  4 Comments
JDilla
JDilla on 12 Mar 2015
Thank you, I will make changes soon!
Guillaume
Guillaume on 12 Mar 2015
Edited: Guillaume on 12 Mar 2015
In addition, on the line
sini = ((sind(i)*(v+1))/v)
Did you really mean to do
(v+1) / v
with v a vector? / (or mrdivive) in this case finds the least-square solution x to x*v = v+1 with v = [3 4 5 6 7].

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!