Programming with M-Files: A Rocket Trajectory Analysis Using For Loops
Summary: This is an example of using m-files to solve a simple engineering data analysis problem in which velocity and acceleration are computed from altitude data from a home-built rocket test launch.
Note:
This example requires an understanding of the relationships between position, velocity, and acceleration of an object moving in a straight line. The Wikipedia article
Motion Graphs and Derivatives
has a clear explanation of these relationships, as well as a discussion of average and instantaneous velocity and acceleration and the role derivatives play in these relationships. Also, in this example, we will approximate derivatives with forward, backward, and central differences;
Lecture 4.1 by Dr. Dmitry Pelinovsky at McMaster University contains useful information about this approximation. We will also approximate integrals using the trapezoidal rule; The Wikipedia article
Trapezium rule
has an explanation of the trapezoidal rule.
Trajectory Analysis of an Experimental Homebuilt Rocket
On his web page
Richard Nakka's Experimental Rocketry Web Site: Launch Report - Frostfire Two Rocket,
Richard Nakka
provides a very detailed narrative of the test firing of his Frostfire Two homebuilt rocket and subsequent data analysis. (His site provides many detailed accounts of tests of rockets and rocket motors. Some rocket launches were not as successful as the Frostfire Two launch; his site provides very interesting post-flight analysis of all launches.)
Computation of Velocity and Acceleration from Altitude Data
In this section, we will use m-files to analyze the altitude data extracted from the plot "Altitude and Acceleration Data from R-DAS" on Richard Nakk's web page. This data is in the file Altitude.txt.
We will use this data to estimate velocity and acceleration of the Frostfire Two rocket during its flight.
Problem 1:
Get the data
Download the altitude data set in the file Altitude.txt onto your computer (right click on this link). The file is formatted as two columns: the first column is time in seconds, and the second column is altitude in feet. Load the data and plot the altitude as a function of time.
The following sequence of commands will load the data, create a vector
t of time values, create a vector s of altitude values, and plot the altitude as a function of time.
load Altitude.txt -ascii t = Altitude(:,1); s = Altitude(:,2); plot(t,s)The plot should be similar to that in Figure 1.
Figure 1:
Plot of altitude versus time.
|
Problem 2:
Forward Differences
Write a script that uses a for loop to compute velocity and acceleration from the altitude data using forward differences. Your script should also plot the computed velocity and acceleration as function of time.
Solution 2
This solution is by Scott Jenne; it computes and plots the velocity:
load Altitude.txt -ascii
t=Altitude(:,1);
s=Altitude(:,2);
for n=1:180;
v=((s(n+1))-s(n))/((t(n+1))-t(n))
hold on
plot(t(n),v,'o')
end
The plot produced by this code is shown in Figure 2.
Figure 2:
Plot of velocity computed with the forward difference method versus time.
|
Problem 3:
Backward Differences
Solution 3
This solution by Bryson Hinton:
load altitude.txt -ascii
t=altitude(:,1);
s=altitude(:,2);
hold on
for x=2:181
v(x)=(s(x)-s(x-1))/(t(x)-t(x-1));
plot(t(x),v(x),'b.')
end
The plot produced by this code is shown in Figure 3.
Figure 3:
Plot of velocity computed with the backward difference method versus time.
|
Problem 4:
Central Differences
Solution 4
This code computes the velocity using the central difference formula.
load Altitude.txt -ascii
t=Altitude(:,1);
s=Altitude(:,2);
for n=2:180
v(n-1)=(s(n+1)-s(n-1))/(t(n+1)-t(n-1));
end
plot(t(2:180),v)
The plot produced by this code is shown in Figure 4.
Figure 4:
Plot of velocity computed with the central difference method versus time.
|
Compare the velocity and acceleration values computed by the different approximations. What can you say about their accuracy?
Problem 5:
Can it be done without loops?
Modify your script from 2 to compute velocity and acceleration without using a for loop.
Solution 5
This code uses the
diff function to compute the difference between adjacient elements of t and s, and the ./ function to divide each element of the altitude differences with the corresponding element of the time differences:
load Altitude.txt -ascii t=Altitude(:,1); s=Altitude(:,2); v=diff(s)./diff(t); plot(t(1:180),v)The plot produced by this code is shown in Figure 5.
Figure 5:
Plot of velocity computed with the forward difference method versus time. The values in this plot are the same as in Figure 2.
|
Computation of Velocity and Altitude from Acceleration Data
In this section, we will use m-files to analyze the acceleration data extracted from the plot "Altitude and Acceleration Data from R-DAS" on Richard Nakk's web page. Download the acceleration data set in the file Acceleration.txt onto your computer (right click on this link). The first column is time in seconds, and the second column is acceleration in g's. The following commands load the data and plot the acceleration as a function of time.
load Acceleration.txt -ascii t = Acceleration(:,1); a = Acceleration(:,2); plot(t,a)The plot should be similar to that in Figure 6.
Figure 6:
Plot of acceleration versus time.
|
Problem 6:
Trapezoidal Rule
Write a script that uses a for loop to compute velocity and altitude from the acceleration data using the trapezoidal rule. Your script should also plot the computed velocity and altitude as function of time.
Solution 6
This solution is by Jonathan Selby:
This code can be easily extended to also compute altitude while it is computing velocity:
load Acceleration.txt -ascii
t=Acceleration (:,1);
a=Acceleration (:,2);
v(1)=0;
for n=1:181
v(n+1)=(t(n+1)-t(n))*(a(n+1)+a(n))/2+v(n);
end
plot(t,v)
This code creates the plot in Figure 7.
Figure 7:
Plot of velocity versus time. The velocity is computed by numerically integrating the measured acceleration.
|
load Acceleration.txt -ascii
t=Acceleration (:,1);
a=Acceleration (:,2);
v(1)=0; % Initial velocity
s(1)=0; % Initial altitude
for n=1:181
v(n+1)=(t(n+1)-t(n))*(a(n+1)+a(n))/2+v(n);
s(n+1)=(t(n+1)-t(n))*(v(n+1)+v(n))/2+s(n);
end
plot(t,s)
This code creates the plot in Figure 8.
Figure 8:
Plot of altitude versus time.
|
Problem 7:
Can it be done without loops?
Modify your script from 6 to compute velocity and altitude without using a for loop.
Solution 7
This solution by Nicholas Gruman uses the
cumtrapz function to compute velocity with the trapezoidal rule:
load Acceleration.txt -ascii t=Acceleration(:,1); A=Acceleration(:,2); v=cumtrapz(t,A);Altitude could also be computed by adding the following line to the end of the previous code:
s=cumtrapz(t,v);
More about this content: Metadata | Version History | Cite This Content
This work is licensed by Darryl Morrell. See the Creative Commons License about permission to reuse this material.
Last edited by Darryl Morrell on Sep 29, 2006 6:10 pm GMT-5.