## working with two variables

on 23 Nov 2011

### Sven (view profile)

working with two variables i'm not sure if this is possible, but is it possible to make a 3D plot using an equation based on two variables?
here's my script: -
r_dir=5 %Distance of direct sound path (m)
r_a=0.15 %Radius of rotating drum (m)
d_dir=-5 %Position of centre of rotating drum (on axis) (m)
l_dir=0 %Position of centre of rotating drum off axis/90 degrees from the axis - in the same plane as the rotating drum width ways
t=[0.05] %Time (s)
r_d=5 %Direct distance between receiver and centre of rotating drum (m)
f_r=5 %Frequency/Rotations per second of rotating drum (Hz)
w_r=2*pi*f_r %Angular frequency of the rotating drum
l_f=r_a*sin(w_r*t) %Position of opening on the rotating drum - on the l axis (90 degrees to the axis, so 'off' axis) (m)
d_f=r_a*cos(w_r*t) %Position of opening on the rotating drum - on the d axis (on axis) (m)
r_pd=sqrt((l_f-l_dir).^2+(d_f-d_dir).^2) %Path difference between direct sound path distance and the diffracted sound path distance (m)
r_f=r_a+r_pd %Diffracted sound path distance (m)
c=350 %Speed of sound at 28 degrees celcius with 90% relative humidity
t_y=(r_f-r_dir)/c %Time delay between the arrival of the direct sound and the diffracted sound (sec)
f=20:5:20000 %Frequency of sound from 20 Hz to 20 kHz
w=2*pi*f; %Angular frequency of the sound
AR=r_dir./r_f %Amplitued Ratio - Pressure amplitude of the direct sound divided by the pressure amplitude of the diffracted sound (same as diffracted sound path distance divided by direct sound path distance in metres)
MF=sqrt(AR.^2+1+AR*2*cos(t_y*w))
plot(f,MF)
I picked t=0.05 sec as a constant in the above example, but I'd like to have it so that is a variable too - hence leading to some kind of graph with three dimensions. It would then be a 3D graph (surface) - one axis for the result, one for frequency of sound and one for time.

### Sven (view profile)

on 23 Nov 2011

It's not the most efficient code, but I've kept things similar to your original input for clarity. Think now of "MF" as a grid or results, with one pixel for each combination of time t, and frequency f:
f = 20:5:20000; %Frequency of sound from 20 Hz to 20 kHz
w=2*pi*f; %Angular frequency of the sound
t = 0:0.001:0.1; % Set of t locations
MF = zeros(length(t),length(f));
r_dir=5; %Distance of direct sound path (m)
r_a=0.15; %Radius of rotating drum (m)
d_dir=-5; %Position of centre of rotating drum (on axis) (m)
l_dir=0; %Position of centre of rotating drum off axis/90 degrees from the axis - in the same plane as the rotating drum width ways
c=350; %Speed of sound at 28 degrees celcius with 90% relative humidity
r_d=5; %Direct distance between receiver and centre of rotating drum (m)
f_r=5; %Frequency/Rotations per second of rotating drum (Hz)
w_r=2*pi*f_r; %Angular frequency of the rotating drum
for i = 1:length(t)
l_f=r_a*sin(w_r*t(i)); %Position of opening on the rotating drum - on the l axis (90 degrees to the axis, so 'off' axis) (m)
d_f=r_a*cos(w_r*t(i)); %Position of opening on the rotating drum - on the d axis (on axis) (m)
r_pd=sqrt((l_f-l_dir).^2+(d_f-d_dir).^2); %Path difference between direct sound path distance and the diffracted sound path distance (m)
r_f=r_a+r_pd; %Diffracted sound path distance (m)
t_y=(r_f-r_dir)/c; %Time delay between the arrival of the direct sound and the diffracted sound (sec)
AR=r_dir./r_f; %Amplitued Ratio - Pressure amplitude of the direct sound divided by the pressure amplitude of the diffracted sound (same as diffracted sound path distance divided by direct sound path distance in metres)
MF(i,:) = sqrt(AR.^2+1+AR*2*cos(t_y*w));
end
figure
subplot(2,1,1), imagesc(f,t,MF), colorbar, xlabel f, ylabel t
subplot(2,1,2), surf(f,t,MF,'EdgeColor','none'), xlabel f, ylabel t, zlabel result

Tom

### Tom (view profile)

on 23 Nov 2011
Thanks a lot for this. I would be really grateful if you would be willing to help me a little further by just helping me to understand a couple of the coding steps that are new to me. These are: -
row 4 - MF = zeros(length(t),length(f));
row 13 & 21 - 'for' and 'end'
row 13 & 20 - i = 1:length(t) & MF(i,:) = ...
and then the different plot types (I just watched a tutorial video on subplots so I get that bit - it's just the bits after the subplot bit).
If you know some good tutorial videos on these concepts - or just some good sticky posts - that would be perfect.
Many thanks
Tom
Sven

### Sven (view profile)

on 23 Nov 2011
Sure, think of it this way... in your original code, your output (MF) was a 1-by-M array of results, where N was the number of frequencies (in f) that you wanted to sample at. In that original code, time (t) was just a single scalar number (ie, 1-by-1). What you then said you wanted to do was repeat your output for *many* different measures of "t". So I also defined "t" as a vector. If you think about how your final output (MF) would look *now*, you'll see that it will be like repeating your original 1-by-M array of results N times (where N is the number of time samples in "t").
The best way to *store* all those results would be in a N-by-M matrix. Each row gives all the frequency results for a given single value of time.
My line 4 command: "MF = zeros(length(t),length(f));"
It simply initialises an empty matrix full of zeros to be the exact size (M-by-N) that we want our final output to be. This is a good MATLAB coding practise when your results will be a big matrix. It's called "preallocation".
Lines 13 and 21. They define a "for loop". Think of looping over every single value for "t" and basically doing your original script (which calculated your MF output for "1" value of t)
Line 13: The "colon" operator. Type "1:20" into the command line, you'll see it makes a vector of all the numbers from 1 to 20.
I think you'll find that the keywords I covered ("for loop", "colon operator", "preallocation" (which is a little more advanced)) are covered well in the "Getting Started" section of the doc. It sounds like you're on the right track though. Hope that's helpful.