How to create nested for loops that has two changing variables to output one variable?

This is the code I am trying to run, and I am unsure how to write nested for loops to run over a certain time period and length.

Context: I am trying to find what the function will look like after the time ends and how r and hz evolve. Ultimately I will want to plot the r and hz together.

%Intial Conditions
rho=(3.6*30*24*60*60); %3.6 months
u=10^9;
u0=2*10^-5; 
v=1*10^6;
r0=2000; %Intial radius
h0=180; %Intial height
V=(3/4)*pi()*(r0^2)*h0; %volume
v0=1*10^6;
g=1.31;
tau=((3/4)^5)*(((pi()^3)*v0*(r0^8))/(g*V^3));
%% Theta 
for t=1:172800; %Using the time for 2 days
    theta(t)=rho*(1-exp(-t/rho));
end

%% Solving function  
for r=1:3000;
    for ta=1:theta;
        hz(r,ta)=(((4*V)/(3*pi()*r0^2)).*(1./((1+(ta./tau)).^1/4)).*(1-(((r.^2)./(r0^2)).*(1./(1+((ta./tau).^1/4))))).^(1/3));
    end
end

[Walter Roberson]

Close but consider

hz(r,ta)=(((4*V)/(3*pi()*r0^2)).*(1./((1+(theta(ta)./tau)).^1/4)).*(1-(((r.^2)./(r0^2)).*(1./(1+((theta(ta)./tau).^1/4))))).^(1/3));

Also, are you sure that you want to have an initial radius but then vary r in steps of 1 from 1 to 3000 ? Consider creating a vector of r values, perhaps with linspace(), and indexing that vector inside the loop.

[Alyssa Mills]

I'm trying to find height as a function of radius. Here is the equation I am trying to use where I eventually want to test different times (t) and radii (r).

[Walter Roberson]

Your code already tests different angles (ta) and different radii (r ) . If you want to test different tau as well you would add another loop and probably add another index on the output.

Answer 1

[Walter Roberson]

rho=(3.6*30*24*60*60); %3.6 months
u=10^9;
u0=2*10^-5; 
v=1*10^6;
r0=2000; %Intial radius
h0=180; %Intial height
V=(3/4)*pi()*(r0^2)*h0; %volume
v0=1*10^6;
g=1.31;
tau=((3/4)^5)*(((pi()^3)*v0*(r0^8))/(g*V^3));
%% Theta 
for t=1:172800; %Using the time for 2 days
    theta(t)=rho*(1-exp(-t/rho));
end
%% Solving function
hz = zeros(3000, length(theta));
for r=1:3000;
    for ta=1:length(theta);
        hz(r,ta)=(((4*V)/(3*pi()*r0^2)).*(1./((1+(theta(ta)./tau)).^1/4)).*(1-(((r.^2)./(r0^2)).*(1./(1+((theta(ta)./tau).^1/4))))).^(1/3));
    end
end

This is rather slow and you should try to vectorize it.

[Alyssa Mills]

Thank you. Although for some reason when I plot the hz over r, I don't get the expected curve that is associated with the equation as shown in the attachment. (Not sure if I need to resubmit as a new question..first time using the forum).

[Walter Roberson]

Those look to me to be four different tau values, but one fixed angle (value not shown.)

[Alyssa Mills]

The tau values are the same for each curve, but the theta's (angles) change because they are dependent on time. I made the curve for t=2 days, but I can no longer remember how, and how I use the code to plot the correct plot. When I plot the hz you shared, I get straight lines. The r^2/r0^2 dictacte the curvature ( I plotted this separately to find what component of the equation accounts for the curvature).

[Walter Roberson]

In the graph, the line t = 3.6 months, how does the t there correspond to any variable in your equation for hz ? And how does any of that relate to the t=2 days curve ?

[Alyssa Mills]

The t=3.6 months in the graph is still the time it takes the height/thickness to get to that end radius of 3000. The rho value is a constant incase you were seeing that 3.6 months which is used in the theta calculation that is then fed into the height/thickness (hz) equation.

[Walter Roberson]

To confirm, what is labeled as t in your graph, such as t = 3.6 months, is the rho variable in your code? So the plot is over 4 different rho values?

[Alyssa Mills]

The t's labeled in the graph are the t variable in the code which is time. So the plot is over 4 different time values which translate to 4 theta values as the time is used in the theta function. Rho is a constant that you can ignore for simplicity (it just helps calculate my theta value).

[Walter Roberson]

tvals = [2, 27, 3.6*365/12, 1.5*365] * 24 * 60 * 60;    %seconds
rvals = 1 : 3000;
rho=(3.6*30*24*60*60); %3.6 months
u=10^9;
u0=2*10^-5; 
v=1*10^6;
r0=2000; %Intial radius
h0=180; %Intial height
V=(3/4)*pi()*(r0^2)*h0; %volume
v0=1*10^6;
g=1.31;
tau=((3/4)^5)*(((pi()^3)*v0*(r0^8))/(g*V^3));
%% Theta 
for tidx = 1 : length(tvals)
    theta(tidx)=rho*(1-exp(-tvals(tidx)/rho));
end
%% Solving function
hz = zeros(3000, length(theta));
for ridx = 1 : length(rvals)
    r = rvals(ridx);
    for tidx = 1:length(theta)
        ta = theta(tidx);
        hz(r,tidx)=(((4*V)/(3*pi()*r0^2)).*(1./((1+(ta./tau)).^1/4)).*(1-(((r.^2)./(r0^2)).*(1./(1+((ta./tau).^1/4))))).^(1/3));
    end
end

subplot(1,2,1)
plot(rvals, real(hz))
legend(sprintfc('%g', tvals/(24*60*60)))
subplot(1,2,2)
plot(rvals, imag(hz))
legend(sprintfc('%g', tvals/(24*60*60)))

[Alyssa Mills]

Thank you so much for the help! You helped vectorize which I am sort of familiar with but couldn't get it exactly down. I'll try to play with this to make the curves match exactly to the previous plot. I had gotten the exact curve of t=2 days previously (using the conditions I have above), but I sadly didn't save that script.