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)

Categories

Find more on Seismology in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!