(*        Simple Rocket Example              

The method that you should use for opening your data tables  *)

(* this version reads a plain text file with 2 columns of data *)
tabb=ReadList["/usr/people/faculty/physics/jcrumley/labs/schader/SummerResearch-A8a", {Number, Number}]
(* This reads a spreadsheet file *)
tabbb=Import["/home/jcrumley/labs/schader/A8_sample.ods"];
(* This pulls the condenses the 1x751x2 matrix down to 751 x2 *)
tabbbb=ArrayReshape[tabbb, {751,2}];
(*This reads the first sheet (tab) out of the spreadsheet file *)
tabc=Import["/home/jcrumley/labs/schader/A8_sample.ods", {"Data", 1}];

(*  A simple Table for this example            *)

tab={{0,0},{.2,5},{.4,5},{.6,5},{.8,2},{1,0},{10,0}}

(* A made up set of points governing the thrust curve of the rockets engine*)

T[t_]=Interpolation[tab,InterpolationOrder->1][t]

(* The Mathematica command function that approimates a function that fits the given points*)

(* A plot of the rocket engines thrust in newtons as a funcion of time*)

Plot[T[t],{t,0,10}]

IM[t_]=Integrate[T[t],t]

(* The intergral of the thrust, gives us the total impluse of the rocket engine.  As you can see, it reaches it's max when the engine stops*)

Plot[IM[t],{t,0,5}, PlotRange-> All]

(* A plot of the impulse as a function of time*)

m[t_]=(Subscript[m, 0])-α*IM[t]

(Subscript[m, 0])=.059

(* Differential equation governing rockets motion under both powered and coasting flight, providing Thrust term goes to zero when the rocket engine stops.  a[t_]:=((T[t])/(m[t]))-(((Subscript[C, d])*ρ*s*(v^2))/(2*m[t]))-g*)

(* Rocket surface area*)
s=6.834*(10^(-4))

(* Mass Change Constant = Δm/Total Impulse = Δm/ IM(t ≥ 0) *)
α=((.0055)/(IM[2]))

(* Air Density from ideal gas law*)
ρ=1.223

(* Rocket drag coeficcent*)
(Subscript[C, d])=.9

(* Acceleration due to gravity*)
g=9.807

solution=NDSolve[{y''[t]==((T[t])/(m[t]))-(((Subscript[C, d])*ρ*s*y'[t]*Abs[y'[t]])/(2*m[t]))-g,y[0]==0,y'[0]==0},y,{t,0,10}]

(* Solving the differential equation*)

Plot[Evaluate[y[t]/.First[solution]],{t,0,9},PlotRange->All]

(* A plot of the height of the rocket as a function of time*)

yy[t_]=Evaluate[y[t]/.First[solution]]

v[t_]=(∂_t)yy[t]

(* Taking the derivative of the rockets postion function gives us its velocity function*)

Plot[v[t],{t,0,5}]

(* A plot of the rockets velocity as a function of time*)

(* We could get the rockets acceleration two ways.  One is to take the derivative of velocity, the other, which is below, is to go back to our original equation now that we have time dependent functions of velcity and mass.*)

a[t_]:=((T[t])/(m[t]))-(((Subscript[C, d])*ρ*s*v[t]*Abs[v[t]])/(2*m[t]))-g

Plot[a[t],{t,0,10}]

(* A plot of the rockets acceleration as a function of time*)