Plotting Several Data Sets in a For Loop

I am working through a problem that requires me to calculate trajectories of an object based on the angle the object is initially thrown. The end goal of this portion of the problem is to have one graph that shows all of the trajectories on one plot. The angles I am trying to iterate through are 0 to 80 in steps of 10 degrees. However, when I run the script, I am getting one plot where all of the data is the same and all of the data (regardless of the angle) has the same output (as evidenced by the legend on the plot). If I change the angle variable to just one value (i.e. to 10, instead of 0:10:80), the program will successfully output a graph for that angle. Can someone help me sort out what my issue is?

clc
clear
% Constant Variables
m=0.175; %mass, kg
D=0.26; %Diameter,
A=(pi/4)*D^2; %Area, m^2
g=9.81; %m/s^2
rho_air=1.23; %Density of air, kg/m^3
AoA_rad=deg2rad(-4); %Angle of Attack for Drag Calc.

%Position Variables
x_prev=0;
y_prev=1;
x_new=0;
y_new=1;

%Time Variables
t_start=0;
t_step=.01;
t_end=7;
t_total=0;
time=t_start:t_step:t_end;

%Loop Variables/Matrices
Vel=14;
Coord=zeros(length(time),2);
count=1;

%Plot information
figure
hold on
title('Frisbee with Constant Velocity and Varying Angles')
xlabel ('Horizontal Distance (m)')
ylabel('Vertical Distance(m)')
grid on
legend 
for j=0:10:80 %Angles to iterate through
    for i=t_start:t_step:t_end
        alpha_rad=deg2rad(j);
        beta=deg2rad(90-j);
        %Defining and Summing Forces
        Cd=0.08+(2.72*alpha_rad)*(alpha_rad-AoA_rad)^2;
        Cl=0.15+(1.4*alpha_rad);
        Fd=(-Cd*rho_air*A*(Vel^2))/2;
        Fl=(Cl*rho_air*A*(Vel^2))/2;
        Sum_Fx=Fd*sin(beta)-Fl*sin(alpha_rad);
        Sum_Fy=Fd*cos(beta)+Fl*cos(alpha_rad)-(m*g);
        
        %Finding Accelerations and New Velocity
        ax=Sum_Fx/m;
        ay=Sum_Fy/m;
        u0x=Vel*cos(alpha_rad);
        v0y=Vel*sin(alpha_rad);
        u0=u0x+i*ax;
        v0=v0y+i*ay;
        %Finding x and y coordinates
        x_new=x_prev+u0x*i+0.5*ax*(i^2);
        y_new=y_prev+v0y*i+0.5*-g*(i^2);
        Coord(count,1)=x_new;
        Coord(count,2)=y_new;
        x_prev=x_new;
        y_prev=y_new;
        count=count+1;
        Vel=sqrt((u0^2)+(v0^2));
    end
    plot(Coord(:,1),Coord(:,2))
    ylim([0 10])
    xlim([0 30])
end

Answer 1

[William Rose]

@George,

You write:

x_new=x_prev+u0x*i+0.5*ax*(i^2);

which is not correct, given how you define variables. You define i as current time. xprev is updated on each inner loop pass. Therefore you should be using dt, the step size, rather than i=time, for this update step. Your equation for y_new has the same issue.

Inspect your inner loop code to identify calculations that can be pulled to the outer loop - that is, anything that is the same on each inner loop pass.

I would use NxM arrays for X, Y, Vel, etc., where N=length(time) and M=number of angles.

I would change the loop variables to integers and use them as indices for the arrays x(i,j) and y(i,j), which you will fill up as you go through the loops.

Here is a simple example for a ball launched at different angles. Your system is more complicated, so your equations for ax, ay will be different.

dt=.01;                     % time step (s)
time=0:dt:2;                % time is not used in the loops, but you could use it to plot(t,x(:,j)), etc 
angle=0:10:80; 
N=length(time); M=length(angle); 
g=-9.8;                     % y accel m/s^2
x=zeros(N,M); y=zeros(N,M); % allocate arrays for x, y
x(1,:)=zeros(1,M);          % initial distance (m)
y(1,:)=ones(1,M);           % initial height (m)
v0=14;                      % initial veloc (m/s)
for j=1:M
    vx=v0*cos(angle(j)*pi/180);
    vy=v0*sin(angle(j)*pi/180);
    ax=0;       % no x acceleration in my simple model
    ay=g;       % y acceleration
    for i=2:N
        x(i,j)=x(i-1,j)+dt*vx;
        y(i,j)=max(y(i-1,j)+dt*vy,0);  % don't allow y<0
        vx=vx+dt*ax;        % update vx
        vy=vy+dt*ay;        % update vy
    end
    plot(x(:,j),y(:,j))
    hold on
    legstr{j}=num2str(angle(j));
end
legend(legstr); grid on;
xlabel('Distance (m)'); ylabel('Height (m)');

That gives you some ideas.

[George]

Thank you! This got me significantly closer. The plots are now separated, and my code appears to be a lot more put together. However, I am running into an issue now where it appears that the acceleration values aren't working correctly, because all of the lines are straight. It's possible I have a sign error, but I would be surprised since the code appeared to be outputting the correct trajectories previous to this edit. Can you help me find this issue?

clc
clear
% Constant Variables
m=0.175; %mass, kg
D=0.26; %Diameter, 
A=(pi/4)*D^2; %Area, m^2
g=9.81; %m/s^2
rho_air=1.23; %Density of air, kg/m^3
AoA_rad=deg2rad(-4); %Angle of Attack for Drag Calc. 
%Time Variables
t_step=.01;
time=0:t_step:8;
angle=0:10:80;
time_length=length(time);
angle_length=length(angle);
x_coord=zeros(time_length,angle_length);
y_coord=zeros(time_length,angle_length);
y_coord(1,:)=ones(1,angle_length);           % initial height (m)
%Loop Variables/Matrices
Vel=14;
%Plot information
figure
for j=1:angle_length %Angles to iterate through
    alpha_rad=deg2rad(j);
    beta=deg2rad(90-j);
    %Defining and Summing Forces
    Cd=0.08+(2.72*alpha_rad)*(alpha_rad-AoA_rad)^2;
    Cl=0.15+(1.4*alpha_rad);
    Fd=(-Cd*rho_air*A*(Vel^2))/2;
    Fl=(Cl*rho_air*A*(Vel^2))/2;
    Sum_Fx=Fd*sin(beta)-Fl*sin(alpha_rad);
    Sum_Fy=Fd*cos(beta)+Fl*cos(alpha_rad)-(m*g);
    
    %Finding Accelerations and New Velocity
    ax=Sum_Fx/m;
    ay=Sum_Fy/m;
    u0x=Vel*cos(alpha_rad);
    v0y=Vel*sin(alpha_rad);
    for i=2:time_length
        %Finding x and y coordinates
        x_coord(i,j)=x_coord(i-1,j)+u0x*t_step+0.5*ax*(t_step^2);
        y_coord(i,j)=y_coord(i-1,j)+v0y*t_step+0.5*-g*(t_step^2);
        u0=u0x+t_step*ax;
        v0=v0y+t_step*ay;
        Vel=sqrt((u0^2)+(v0^2));
    end 
    plot(x_coord(:,j),y_coord(:,j))
    hold on 
    ylim([0 10])
    xlim([0 30])
end
title('Frisbee with Constant Velocity and Varying Angles')
xlabel ('Horizontal Distance (m)')
ylabel('Vertical Distance(m)')
grid on
legend  

[William Rose]

@George,

[edit: fix spelling errors]

I am not sure what is wrong in your code. I think I have a hard time finding the problem because the equations in the code do not look very much like the equations of motion which one would write out on paper. For this reason I recommend writing code in which the code resembles the equations on paper.

The equations on paper are (tell me if I'm wrong):

lift:

drag:

Lift is perpendicular to the velocity vector; drag is anti-parallel to the velocity vector. The instantaneous velocity is , i.e. ,. Then the x and y forces are:

In your model, C_d and C_L are constants, right? (They could depend on the angle of attack, which is likely to change, but that complication does not appear to be included in your model - and that is fine.)

[Voss]

@George: This was good advice from @William Rose:

"Inspect your inner loop code to identify calculations that can be pulled to the outer loop - that is, anything that is the same on each inner loop pass."

And I see that you followed it, but I think you pulled a little too much out of the inner loop. Specifically, notice that Vel is recalculated on the last line of the inner loop, so anything that depends on Vel must remain in the inner loop.

See below for the code with that change. (I also disabled setting xlim and ylim, so I can see the entire span of the curves.)

clc
clear
% Constant Variables
m=0.175; %mass, kg
D=0.26; %Diameter, 
A=(pi/4)*D^2; %Area, m^2
g=9.81; %m/s^2
rho_air=1.23; %Density of air, kg/m^3
AoA_rad=deg2rad(-4); %Angle of Attack for Drag Calc. 
%Time Variables
t_step=.01;
time=0:t_step:8;
angle=0:10:80;
time_length=length(time);
angle_length=length(angle);
x_coord=zeros(time_length,angle_length);
y_coord=zeros(time_length,angle_length);
y_coord(1,:)=ones(1,angle_length);           % initial height (m)
%Loop Variables/Matrices
Vel=14;
%Plot information
figure
for j=1:angle_length %Angles to iterate through
    alpha_rad=deg2rad(j);
    beta=deg2rad(90-j);
    %Defining and Summing Forces
    Cd=0.08+(2.72*alpha_rad)*(alpha_rad-AoA_rad)^2;
    Cl=0.15+(1.4*alpha_rad);
    
    for i=2:time_length
        Fd=(-Cd*rho_air*A*(Vel^2))/2;
        Fl=(Cl*rho_air*A*(Vel^2))/2;
        Sum_Fx=Fd*sin(beta)-Fl*sin(alpha_rad);
        Sum_Fy=Fd*cos(beta)+Fl*cos(alpha_rad)-(m*g);
        
        %Finding Accelerations and New Velocity
        ax=Sum_Fx/m;
        ay=Sum_Fy/m;
        u0x=Vel*cos(alpha_rad);
        v0y=Vel*sin(alpha_rad);
        
        %Finding x and y coordinates
        x_coord(i,j)=x_coord(i-1,j)+u0x*t_step+0.5*ax*(t_step^2);
        y_coord(i,j)=y_coord(i-1,j)+v0y*t_step+0.5*-g*(t_step^2);
        u0=u0x+t_step*ax;
        v0=v0y+t_step*ay;
        Vel=sqrt((u0^2)+(v0^2));
    end
    
    plot(x_coord(:,j),y_coord(:,j))
    hold on 
%     ylim([0 10])
%     xlim([0 30])
end
title('Frisbee with Constant Velocity and Varying Angles')
xlabel ('Horizontal Distance (m)')
ylabel('Vertical Distance(m)')
grid on
legend  

[William Rose]

@George,

I get the images below (right plot is zoomed in version of left plot):

The high-launch angle trajectories look weird. The frisbee is launched with a forward velocity component, but it travels backward for a little while. Why? Because the lift force (which is perpendicular to velocity) is directed mainly in the negative x direction.

[George]

@William Rose and @Voss

Thank you both. @William Rose what is the difference between the two codes? I understand why @Voss made the change, but I can't seem to get a change that replicates what you got.

@Voss - I am also concerned that the objects seem to take a while to come back down. Does that also seem odd to you?

[William Rose]

@George, @Voss, my code is below. It is my attempt to implement the equaitons as I wrote them in an earlier comment. I think these are the equaitons that you are trying to solve. And I did it with a first order Euler method, which is more or less what you were doing - although your method is partly second order, because you incoporate acceleration in the update for x and y.

I don;t know why @Voss's modification of your code gives different answers than mine. If you really want to know, you could allocate arrays for Fx(i,j) and Fy(i,j) in each script and use them to save the forces. Then you can plot the forces in the two scripts for the different angles. This would likely provide insight.

[William Rose]

@George,

Your script incldues the lines

u0x=Vel*cos(alpha_rad);
v0y=Vel*sin(alpha_rad);

After that, you use u0x and v0y to update x and y. In other words, the equations above are the x, y components of instantaneous velocity. I think the equations are incorrect, because they use alpha_rad, which is the initial launch angle. Alpha_rad does not get updated as the trajectory changes. It should. Or you can do what I did, which is to update vx and vy, based on the accelerations ax and ay.

vx=vx+dt*Fx/m;                  % update vx
vy=vy+dt*Fy/m;                  % update vy

Use vx and vy to update the positions:

x(i,j)=x(i-1,j)+dt*vx;          % update x
y(i,j)=y(i-1,j)+dt*vy;          % update y

That may explain the difference.