Trajectory Code

Today, let’s cover the basics of writing a trajectory code. There are many different levels of fidelity for an analysis, but today we will start at the basic 1-DOF code, just worrying about altitude. Also, we will assume constant drag coefficient, thrust, and gravity. If you want a more detailed code, the first step is to add thrust and drag coefficients that vary with altitude. The next step would be to add a 2nd degree of freedom with downrange, as well as altitude. At this point, if you are doing just basic sizing, you could stop. Then, a 3-DOF model is made by adding angle of attack with vehicle center of gravity and center of pressure. Then, add wind into the 3D model as well as thrust vector response. This is the most complicated that I have done and it is sufficient for sizing gimbals and fins in a variety of wind loads. From here, you would add the next 3-DOF and then probably add a Monte Carlo analysis on top of that.

If you are just using normal hobby rockets, you could get away with using a code like RASAero and not have to write your own. But, let’s be honest, if you are reading this blog, you probably like doing the code for the fun of it. I still recommend using a well tested code like RASAero to ballpark the first few answers your analysis gives you.

The basics of the code is Newton’s Second Law that the Sum of Forces = Mass x Acceleration and the knowledge that velocity is the integral of acceleration and position is the integral of velocity. So we start with a known position, velocity, and mass and solve for all of the forces (drag and thrust) to find the acceleration. Then we go to the next step and find a new velocity as a function of acceleration and the old velocity, new position as a function of old position, old velocity, and acceleration, new mass as a function of old mass and mass flow rate, and then start the cycle over again.

X_new = X_old + V_old*dt + 1/2 * A_old  * dt^2

V_new = V_old + A_old * dt

A_new = Force_new / Mass_new

For the forces, we set Thrust = Thrust when thrusting, and Thrust = 0 otherwise. Drag is more complicated as it is a function of speed and air density. For air density, we use the old standby: the 1976 Standard Atmosphere. Annoyingly, density does not curve fit well, so I curve fit pressure and temperature (they are both exponential), then calculated density. And gravity is always 9.8 m/s pointing to the earth

Force_new = Thrust – Drag – Gravity

Drag = 1/2 * Density * Velocity * abs(Velocity) * Cd * Area

Density = Pressure * Molar Mass / Temperature / Gas Constant

The velocity is multiplied by the absolute value of velocity instead of just squaring to retain the proper sign; otherwise, drag acts as thrust when the vehicle is coming back down to earth.

So you pretty much just take all of these equations and add them together into one integrated iterator and you are good to go. I have put an example for Earendel done in Openoffice Calc below. It has an OK accuracy, achieving 86 km altitude and a max speed of 1200 m/s vs. a more accurate 2D model with variable thrust and drag and finer stepping that achieves 1299 m/s max speed and an altitude of 106 km. Not too shabby for a simple code whipped up in an hour.

Simple Trajectory

Thanks to Professor Anderson at Purdue for teaching me this in undergrad at Purdue.