<?xml version="1.0" encoding="utf-8" standalone="no"?>
<!DOCTYPE document PUBLIC "-//CNX//DTD CNXML 0.5 plus MathML//EN" "http://cnx.rice.edu/technology/cnxml/schema/dtd/0.5/cnxml_mathml.dtd">
<document xmlns="http://cnx.rice.edu/cnxml" xmlns:md="http://cnx.rice.edu/mdml/0.4" xmlns:bib="http://bibtexml.sf.net/" xmlns:m="http://www.w3.org/1998/Math/MathML" id="new">
  <name>Programming with M-Files: A Rocket Trajectory Analysis Using For Loops</name>
  <metadata>
  <md:version>1.7</md:version>
  <md:created>2006/01/27 17:11:43 US/Central</md:created>
  <md:revised>2006/09/29 18:17:00.082 GMT-5</md:revised>
  <md:authorlist>
      <md:author id="morrell">
      <md:firstname>Darryl </md:firstname>
      
      <md:surname>Morrell</md:surname>
      <md:email>morrell@asu.edu</md:email>
    </md:author>
  </md:authorlist>

  <md:maintainerlist>
    <md:maintainer id="morrell">
      <md:firstname>Darryl </md:firstname>
      
      <md:surname>Morrell</md:surname>
      <md:email>morrell@asu.edu</md:email>
    </md:maintainer>
  </md:maintainerlist>
  
  <md:keywordlist>
    <md:keyword>Engineering data analysis</md:keyword>
    <md:keyword>for loop</md:keyword>
    <md:keyword>MathScript</md:keyword>
    <md:keyword>MATLAB</md:keyword>
    <md:keyword>Octave</md:keyword>
  </md:keywordlist>

  <md:abstract>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.</md:abstract>
</metadata>
  <content>

<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
<link src="http://en.wikipedia.org/wiki/Motion_graphs_and_derivatives">
 <cite>Motion Graphs and Derivatives</cite></link>
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; <link src="http://dmpeli.math.mcmaster.ca/Matlab/Math1J03/LectureNotes/Lecture4_1.htm">
<cite> Lecture 4.1 </cite> by Dr. Dmitry Pelinovsky at McMaster University </link> contains useful information about this approximation.  We will also approximate integrals using the trapezoidal rule; The Wikipedia article
<link src="http://en.wikipedia.org/wiki/Trapezoidal_rule">
 <cite>Trapezium rule</cite></link>
has an explanation of the trapezoidal rule.
</note>

<section id="SecRocket">
	  <name>Trajectory Analysis of an Experimental Homebuilt Rocket</name>


<para id="ParaRocket1">
On his web page
<link src="http://www.nakka-rocketry.net/ff-2.html">
<cite>Richard Nakka's Experimental Rocketry Web Site: Launch Report - Frostfire Two  Rocket</cite></link>, 
Richard Nakka
provides a very detailed narrative of the test firing of his  Frostfire Two homebuilt rocket  and subsequent data analysis. (<link src="http://www.nakka-rocketry.net/">His site</link> 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.)
</para>
</section>
<section id="Sec2">
	  <name>Computation of Velocity and Acceleration from Altitude Data</name>

<para id="xyzzy">
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 <link src="Altitude.txt">Altitude.txt</link>.
We will use this data to estimate velocity and acceleration of the Frostfire Two rocket during its flight.
</para>

<exercise id="E1">
<problem id="P1">
<name>Get the data</name>
<para id="Pa1">
Download the altitude data set in the file <link src="Altitude.txt">Altitude.txt</link> onto your computer  (right click on <link src="Altitude.txt">this link</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.
</para>
</problem>
</exercise>
<para id="S1">
The following sequence of commands will load the data, create a vector <code>t</code> of time values, create a vector <code>s</code> of altitude values, and plot the altitude as a function of time.
<code type="block">
load Altitude.txt -ascii
t = Altitude(:,1);
s = Altitude(:,2);
plot(t,s)
</code>
The plot should be similar to that in <cnxn target="FigAltPlot"/>.
<figure id="FigAltPlot">
  <media type="application/postscript" src="AltPlot.eps">
    <media type="image/png" src="AltPlot.png"/>
  </media>
  <caption>
    Plot of altitude versus time.
  </caption>
</figure>

</para>

<exercise id="E2">
<problem id="P2">
<name>Forward Differences</name>
<para id="Pa2">
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.
</para>
</problem>
<solution>
<para id="S2">This solution is by Scott Jenne; it computes and plots the velocity:
<code type="block">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
</code>
The plot produced by this code is shown in <cnxn target="FigP2Plot"/>.
<figure id="FigP2Plot">
  <media type="application/postscript" src="P2Plot.eps">
    <media type="image/png" src="P2Plot.png"/>
  </media>
  <caption>
    Plot of velocity computed with the forward difference method versus time.
  </caption>
</figure>

</para>
</solution>
</exercise>

<exercise id="E3">
<problem id="P3">
<name>Backward Differences</name>
<para id="Pa3">
Modify your script from <cnxn target="P2"/> to compute velocity and acceleration using backward differences.  Remember to save  your modified script with a different name than your script from <cnxn target="P2"/>.
</para>
</problem>
<solution>
<para id="S3">This solution by Bryson Hinton:
<code type="block">
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
</code>
The plot produced by this code is shown in <cnxn target="FigP3Plot"/>.
<figure id="FigP3Plot">
  <media type="application/postscript" src="P3Plot.eps">
    <media type="image/png" src="P3Plot.png"/>
  </media>
  <caption>
    Plot of velocity computed with the backward difference method versus time.
  </caption>
</figure>


</para>
</solution>
</exercise>

<exercise id="E4">
<problem id="P4">
<name>Central Differences</name>
<para id="Pa4">
Modify your script from <cnxn target="P2"/> to compute velocity and acceleration using central differences.  Remember to save  your modified script with a different name than your script from <cnxn target="P2"/> and <cnxn target="P3"/>.
</para>
</problem>
<solution>
<para id="S4">This code computes the velocity using the central difference formula.
<code type="block">
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)
</code>
The plot produced by this code is shown in <cnxn target="FigP4Plot"/>.
<figure id="FigP4Plot">
  <media type="application/postscript" src="P4Plot.eps">
    <media type="image/png" src="P4Plot.png"/>
  </media>
  <caption>
    Plot of velocity computed with the central difference method versus time.
  </caption>
</figure>

</para>
</solution>
</exercise>

<para id="Blah">
Compare the velocity and acceleration values computed by the different approximations. What can you say about their accuracy?
</para>

<exercise id="E5">
<problem id="P5">
<name>Can it be done without loops?</name>
<para id="Pa5">
Modify your script from <cnxn target="P2"/> to compute velocity and acceleration without using a for loop.
</para>
</problem>
<solution>
<para id="S5">
This code uses the <code>diff</code> function to compute the difference between adjacient elements of <code>t</code> and <code>s</code>, and the <code>./</code> function to divide each element of the altitude differences with the corresponding element of the time differences:
<code type="block">
load Altitude.txt -ascii
t=Altitude(:,1);
s=Altitude(:,2); 
v=diff(s)./diff(t);
plot(t(1:180),v)
</code>
The plot produced by this code is shown in <cnxn target="FigP5Plot"/>.
<figure id="FigP5Plot">
  <media type="application/postscript" src="P5Plot.eps">
    <media type="image/png" src="P5Plot.png"/>
  </media>
  <caption>
    Plot of velocity computed with the forward difference method versus time.  The values in this plot are the same as in <cnxn target="FigP2Plot"/>.
  </caption>
</figure>

</para>
</solution>
</exercise>
</section>
<section id="Sec3">
	  <name>Computation of Velocity and Altitude from Acceleration  Data</name>
<para id="xx1">
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 <link src="Acceleration.txt">Acceleration.txt</link> onto your computer  (right click on <link src="Acceleration.txt">this link</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.
<code type="block">
load Acceleration.txt -ascii
t = Acceleration(:,1);
a = Acceleration(:,2);
plot(t,a)
</code>
The plot should be similar to that in <cnxn target="FigAccPlot"/>.
<figure id="FigAccPlot">
  <media type="application/postscript" src="AccPlot.eps">
    <media type="image/png" src="AccPlot.png"/>
  </media>
  <caption>
    Plot of acceleration versus time.
  </caption>
</figure>
</para>

<exercise id="E2a">
<problem id="P2a">
<name>Trapezoidal Rule</name>
<para id="Pa2a">
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.
</para>
</problem>
<solution>
<para id="S2a">
This solution is by Jonathan Selby:
<code type="block">
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)
</code>
This code creates the plot in <cnxn target="FigVPlot"/>.
<figure id="FigVPlot">
  <media type="application/postscript" src="VPlot.eps">
    <media type="image/png" src="VPlot.png"/>
  </media>
  <caption>
    Plot of  velocity versus time. The velocity is computed by numerically integrating the  measured acceleration.
  </caption>
</figure>
This code can be easily extended to also compute altitude  while it is computing velocity:
<code type="block">
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)
</code>
This code creates the plot in <cnxn target="FigSPlot"/>.
<figure id="FigSPlot">
  <media type="application/postscript" src="SPlot.eps">
    <media type="image/png" src="SPlot.png"/>
  </media>
  <caption>
    Plot of  altitude versus time.
  </caption>
</figure>
</para>
</solution>
</exercise>

<exercise id="E5a">
<problem id="P5a">
<name>Can it be done without loops?</name>
<para id="Pa5a">
Modify your script from <cnxn target="P2a"/> to compute velocity and altitude without using a for loop.
</para>
</problem>
<solution>
<para id="S5a">This solution by Nicholas Gruman uses the <code>cumtrapz</code> function to compute velocity  with the trapezoidal rule:
<code type="block">
load Acceleration.txt -ascii
t=Acceleration(:,1);
A=Acceleration(:,2);
v=cumtrapz(t,A);
</code>
Altitude could also be computed by adding the following line to the end of the previous code:
<code type="block">
s=cumtrapz(t,v);
</code>
</para>
</solution>
</exercise>

</section>
  </content>
  
</document>
