3d plot of Tridimensional parabolic motion

12 views (last 30 days)
Giovanni Curiazio on 31 Jan 2023
Edited: John D'Errico on 1 Feb 2023
Hi, I am trying to get a 3D plot of a parabolic motion. I am trying to simulate debris from an airplane falling from a certain altitude. Obviously the parabolic motion and the absence of other forces is a simplification.
These are the initial conditions.
h_0=14325 m
V_0= 200*.3048 m/s
Obviously I want my object to reach the ground. I created the following code, but the plot is a little strange. I don't understand where or what is wrong
g=9.8;
z0= 14325;
x0=0;
y0=0;
V_0= 200*0.3048;
gamma=0.78;
V_0x=V_0*cos(gamma);
V_0y=V_0*cos(gamma);
V_0z=V_0*sin(gamma);
t = linspace(0,2*V_0*sin(gamma)/g);
x=x0+V_0x*t;
y=y0+V_0y*t-.5*g*t.^2;
z=z0+V_0z-.5*g*t.^2;
indices=z>=0;
x = x(indices);
y = y(indices);
z = z(indices);
plot3(x,y,z)
hold on
plot3(x0,y0,z0,'ro','MarkerSize',5)
xlabel('Posizione x (m)')
ylabel('Posizione y (m)')
zlabel('Posizione z (m)')
grid on

John D'Errico on 31 Jan 2023
Edited: John D'Errico on 1 Feb 2023
It should look strange. Why?
g=9.8;
z0= 14325;
x0=0;
y0=0;
V_0= 200*0.3048;
gamma=0.78;
V_0x=V_0*cos(gamma);
V_0y=V_0*cos(gamma);
V_0z=V_0*sin(gamma);
t = linspace(0,2*V_0*sin(gamma)/g);
x=x0+V_0x*t;
y=y0+V_0y*t-.5*g*t.^2;
z=z0+V_0z-.5*g*t.^2;
indices=z>=0;
x = x(indices);
y = y(indices);
z = z(indices);
plot3(x,y,z)
hold on
plot3(x0,y0,z0,'ro','MarkerSize',5)
xlabel('Posizione x (m)')
ylabel('Posizione y (m)')
zlabel('Posizione z (m)')
grid on
hold off
Your equations of motion are wrong. Let me explain...
First, you appear to have an angle gamma. Should we presume that gamma is the angle above the horizontal of the initial velocity?
ans = 44.6907
So a roughly 45 degree angle above the horizontal.
And then you split that into a vertical component, AND a horizontal component. So we saw this:
V_0x=V_0*cos(gamma);
V_0y=V_0*cos(gamma);
V_0z=V_0*sin(gamma);
But what you did was set the SAME velocity in both the x and y directions. Essentially, that puts too much velocity into the horizontal component, That is wrong, because now the initial velocity is
norm([V_0x,V_0y,V_0z])
ans = 74.7947
But V_0 is
V_0
V_0 = 60.9600
Do you understand the problem? You needed to split the horizontal velocity into two components properly, NOT match the SAME velocities in both x and y.
Ok, but that was a relatively minor problem. Your initial velocity was wrong. What else is wrong?
x=x0+V_0x*t;
y=y0+V_0y*t-.5*g*t.^2;
z=z0+V_0z-.5*g*t.^2;
Do you see what you did?
Next, you have the z equations wrong. You dropped a t in there.
But worse, why is there an acceleration due to gravity in the y direction?
Finally, is there a reason why you only went out as far in time as you did?
2*V_0*sin(gamma)/g
That is, is that really the amount of time necessary for the projectile to hit the ground, having started up a long way above the ground?
I suppose you can make up equations of motion that will work however you want. But then don't expect them to look as if this projectile is moving in a planetary gravitational field. :)
What would I change to fix it?
First, I would decide how to split up the HORIZONTAL component of velocity into two pieces. So what directino is the projective moving between x and y?
g=9.8;
z0= 14325;
x0=0;
y0=0;
V_0= 200*0.3048;
phi = .5236; % the angle between x and y for the initial components of velocity
gamma=0.78;
V_0x=V_0*cos(gamma)*cos(phi);
V_0y=V_0*cos(gamma)*sin(phi);
V_0z=V_0*sin(gamma);
First, is this the correct overall velocity?
norm([V_0x,V_0y,V_0z])
ans = 60.9600
Now V_0 is consistent.
Next, when will the projectile hit the ground? How long? If the z-equation is correct, then this will do:
troots = roots([-.5*g, V_0z, z0])
troots = 2×1
58.6205 -49.8711
Taking the positive root, we can do this:
t = linspace(0,troots(troots > 0),1000);
x=x0+V_0x*t;
y=y0+V_0y*t; % no acceleration in y due to gravity
z=z0+V_0z*t-.5*g*t.^2; % see the difference!
Did it work?
z(end)
ans = 0
% You don't need to worry about negative z anymore.
% indices=z>=0;
% x = x(indices);
% y = y(indices);
% z = z(indices);
plot3(x,y,z)
hold on
plot3(x0,y0,z0,'ro','MarkerSize',5)
xlabel('Posizione x (m)')
ylabel('Posizione y (m)')
zlabel('Posizione z (m)')
grid on
Now you could add in air resistance too. Make things interesting.
Giovanni Curiazio on 1 Feb 2023
OK, I completely understood the errors. But how do I achieve a final altitude of 0? Thank you

Categories

Find more on Programming in Help Center and File Exchange

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by