# Model Rocket Altitude predictor...

This is a discussion on Model Rocket Altitude predictor... within the C++ Programming forums, part of the General Programming Boards category; Sorry in advance for the large quantity of code. I have gone over this line by line, and can't find ...

1. ## Model Rocket Altitude predictor...

Sorry in advance for the large quantity of code. I have gone over this line by line, and can't find the issue. I'm not sure if it is in the formulae, or in the code.

This is the issue. This program takes a large variety of variables, from frontal area of the rocket, wind resistance, thrust, burn time, gravity, etc, into account. It performs the calculations, and predicts total altitude of the rocket. Fine, great. The problem is, I get garbage altitude. Like a 1-oz model rocket with a 100-Ns thrust motor.... which is bigger than the rocket, and weighs POUNDS, not ounces. This rocket should almost land on the moon. Total altitude: -44 feet. So I've tried everything I can think of, and am now coming for help. Sorry for the large chunks of code, with few comments, but I will explain anything if asked.

Code:
```#include <iostream>
#include <math.h>

using namespace::std;

float velocity_squared(float velocity);
float r_squared (float r);
float root_q (float thrust, float Mass, float g, float coef_k );
float get_l_value(float Mass, float g, float coef_k, float v_s, float thrust);
float get_l_2(float Mass, float g, float coef_k, float v_s);

const float pi = 3.14;
const float Cd = 0.75;
const float e = 2.71828182845904523536;
const float rho = 1.2;
const float g = 9.8;

float l_value; /* this is the value that is used for the equation yc = [+M / (2*coef_k)]*ln([M*g + k*v^2] / [M*g]);
l_value =(M*g + k*v^2) / M*g)*/
float l_val_2;

float velocity; //v, in the equation
float v_s; //v^2;
float area;
float diam_inch;
float r_s; //radius squared, used for pi*r^2;
float Mass; //in kilograms;
float boost_alt; //altitude at motor burnout
float coast_alt;//altitude gained while coasting;
float total_alt; //boost_alt + coast_alt;
float coef_k; //wind resistance temporary holder;
float impulse;
float thrust;
float burn_time;
float q;
float x;
float x_neg;
float m_neg;

int main()
{
cout<<"Enter the mass of the rocket in ounces: \n";
cin>>Mass;
Mass = Mass /16 / 2.2; //conversion to Kg from ounces
system("CLS");
cout<<"Enter the diameter of the rocket, in inches: \n";
cin>> diam_inch;
system("CLS");
r = diam_inch / 2; //useless now
r_squared(r); //ignore this whole  function, it is no longer used
area = pi* (0.5*(diam_inch / 12) * 0.3048) * (0.5*(diam_inch / 12) * 0.3048);
//area in square meters;
coef_k = 0.5 * rho *Cd * area; /*temporary holding variable to implement wind resistance in other equations*/
cout<<"Enter the Impulse for the motor: \n";
cin>>impulse;
system("CLS");
cout<<"Enter the Thrust for the motor: \n";
cin>>thrust;
system("CLS");
burn_time = impulse / thrust;
q = root_q(thrust,  Mass, g, coef_k);
x = 2*coef_k*q/Mass;
x_neg = x * -1;
velocity = q*(1-pow(e,x_neg*burn_time) / (1+pow(e,x_neg*burn_time)));
m_neg = Mass * -1;
v_s = velocity_squared(velocity);
l_value = get_l_value(Mass, g, coef_k, v_s, thrust);
boost_alt = (m_neg / (2*coef_k)) * log(l_value);
l_val_2 = get_l_2(Mass, g, coef_k, v_s);
coast_alt = (Mass / (2*coef_k)* log(l_val_2));
total_alt = boost_alt + coast_alt;
cout<<"Frontal Area of Rocket is " << area <<"square meters. \n";
cout<<"Total altitude: " << total_alt << "feet. \n";
system("PAUSE");
if (total_alt < 0)
{
cout<<"Try putting the motor in the other way. \n"; /*prompted by the consistent negative altitudes.*/
system("PAUSE");
}
return 0;
}

float velocity_squared(float velocity)
{
v_s = velocity*velocity;
return v_s;
}

float r_squared(float r)
{
r_s = r * r;
return r_s;
}

float root_q (float thrust, float Mass, float g, float coef_k)
{
q = ((thrust - Mass * g )/coef_k);
q = sqrt(q);
return q;
}

float get_l_value(float Mass, float g, float coef_k, float v_s, float thrust)
{
l_value =(thrust - Mass*g + coef_k*v_s) / (Mass*g);
return l_value;
}

float get_l_2(float Mass, float g, float coef_k, float v_s)
{
l_val_2 = (Mass * g + coef_k * v_s) / (Mass * g);
return l_val_2;
}```
Again I apologize for the mass of code, but I can't FIND a problem with the coding or formula.

Anyone have any suggestions that might effectively point the motor the right way, at least in simulation? :P

2. have you tried to debug?

try to simplify formula till you can reproduce the calculations by hand,

replace g=9,8 with g=10
coef_k with 1

remove converting Mass and enter it in Kg

Get to the simplest formula possible. And debug your application step-by-step checking each result on paper.

When you get this stripped formula working - start adding back parts you have removed, adding it one-by-one, checking that your new formula is still working

3. Is l ("l_value") supposed to be small? Because it isn't, and your altitude calculations take a negative number (m_neg will be negative) and multiply it by log(l), so if you want your altitude to be positive l needs to be less than one.

Also "pow(e, whatever)" is better written "exp(whatever)".

4. Originally Posted by tabstop
Is l ("l_value") supposed to be small? Because it isn't, and your altitude calculations take a negative number (m_neg will be negative) and multiply it by log(l), so if you want your altitude to be positive l needs to be less than one.

Also "pow(e, whatever)" is better written "exp(whatever)".
to be perfectly honest... I am not sure. l_value was a step in simplification in the line. I was trying to break it down to the least amount of variables in the actual formula in main. these formulae were taken from a website, they are proven (a friend i dont talk with much anymore used them to make his program, which was accurate to within about 500 feet or so). I realize that m_neg will be negative, it's the negative of mass... i will have to do more research and understand more fully every variable in the equation... something i hoped to do after the program functioned at least semi properly, so that I could see what each part did, when it was functioning properly.

the original formula called for 1-exp(-x*t), where e is Euler's number (a constant), to the power of (-x*t). As I understood it, the best way to express this would be pow(e,-x*t), or perhaps...
int x, t;
x = x * -1;
int xt;
x = x*t;
then inserting that into the exponent, ie, 1-pow(e,x);

But, to learn more, why (and how) would I implement 1-exp(e,-x*t), or is that exactly how I would?

5. Originally Posted by kalor_alros
the original formula called for 1-exp(-x*t)
So that's what you type, "1 - exp(-x*t)" (there's no reason to make a separate variable x_neg, when you can just type -x and you're there). Well, I guess you called t "burn_time".

6. > const float e = 2.71828182845904523536;
The useful part of this constant is shown in blue.
Using doubles would extend that into the red part as well.

Have you considered all the overflow and underflow issues, and all the rounding/truncation issues with your calculations?

Whilst you may have the mathematics correct, it does not necessary mean that what you have is computationally accurate. All floats are approximations, and if you're on the wrong end of it, your results are rubbish.

7. so just using 4-5 places after the decimal point for all floats would make this (possibly) more accurate and functional? I'd rather not have to redeclare all of the variable data types to doubles, since well, fundamentally, i'm lazy :P Not really, but I don't really NEED that level of precision for this.

8. Originally Posted by kalor_alros
so just using 4-5 places after the decimal point for all floats would make this (possibly) more accurate and functional? I'd rather not have to redeclare all of the variable data types to doubles, since well, fundamentally, i'm lazy :P Not really, but I don't really NEED that level of precision for this.
If you don't NEED that level of precision, why did you type all those numbers in then?

And anyway, you don't need the constant e in this program anyway so it doesn't much matter.

9. Probably too lazy to read this then What Every Computer Scientist Should Know About Floating-Point Arithmetic

Step through the code one line at a time using a debugger.
And do the same calculations on paper.

Correct use of floating point is no easy walk in the park, you need to be paying attention to the precision problems all the way. You CAN'T just copy/paste the math directly into C code and expect it to always work for you.

10. I do not think this is a precision problem. It seems more likely a problem in the formulae. Do you happen to be using the equations found here:

Rocket Equations

11. Originally Posted by brewbuck
I do not think this is a precision problem. It seems more likely a problem in the formulae. Do you happen to be using the equations found here:

Rocket Equations
If so, then this:
Code:
`l_value =(thrust - Mass*g + coef_k*v_s) / (Mass*g);`
should be this:
Code:
`l_value =(thrust - Mass*g + coef_k*v_s) / (thrust - Mass*g);`

12. Here's an attempt at simulating rocket flight using Euler's method for numerical integration. I don't trust the direct calculation given by the equations shown here. For one thing, air resistance decreases exponentially with height, which makes a big difference in the drag term. Also, gravity decreases slowly with height and I wanted to account for that. If you really want to simulate rockets to the moon you need to treat both gravity and drag more properly.

The output of the program is in CSV format, with the first column being time since launch, the second rocket height above Earth's surface, the third the rocket velocity relative to Earth.

You can adjust the constants in the block at the beginning of main() to suit the situation. DeltaT might be increased significantly without much loss in precision. I haven't bothered to analyze the error of the method.

Code:
```#include <iostream>

#include <cmath>
#include <algorithm>

void EvaluateTimestep( double K, double DeltaT, double Mr, double Me, double G,
double Ft, double *x )
{
// Leap-frog Euler method -- interleave calculation of DeltaX, DeltaV
// Move the position
x[ 0 ] += DeltaT * x[ 1 ];
// Compute the total force
double F = -G * Me * Mr / ( x[ 0 ] * x[ 0 ] );
double DragSign = ( x[ 1 ] < 0 ) ? 1.0 : -1.0;
F += 0.5 * K * ( x[ 1 ] * x[ 1 ] ) * DragSign;
F += Ft;
// Apply acceleration
x[ 1 ] += DeltaT * F / Mr;
}

int main()
{
double RocketMass = 0.1; // kg
double FuelMass = 0.3; // kg
double BurnRate = 0.083; // kg/s
double ExhaustVelocity = 80.0; // m/s
double Thrust = ExhaustVelocity * BurnRate; // N
double BurnoutTime = FuelMass / BurnRate; // s
double DeltaT = 0.01; // s
double AirDensity = 1.292; // kg/m^3
double AirDensityScale = 5600.0; // m -- height at which density is halved
double CrossSection = 2.5e-5; // m^2
double DragCoefficient = 0.25;
double EarthMass = 5.97e24; // kg
double GravConstant = 6.67e-11; // m^3/kg/s^2;
double EarthRadius = 6.3e6; // m
double InitialHeight = 0.0; // m -- start on the ground

double x[ 2 ]; // x[ 0 ] == position, x[ 1 ] == velocity

x[ 0 ] = EarthRadius + InitialHeight; // Initial position
x[ 1 ] = 0.0; // Initial velocity
double T = 0.0;

std::cout << "Time, Position, Velocity" << std::endl;

// Rocket burn
for( ; T < BurnoutTime; T += DeltaT )
{
double Height = x[ 0 ] - EarthRadius;
double XAirDensity = AirDensity * pow( 0.5, Height / AirDensityScale );
double Drag = XAirDensity * CrossSection * DragCoefficient;
std::cout << T << ", " << Height << ", " << x[ 1 ] << std::endl;
EvaluateTimestep( Drag, DeltaT, RocketMass + FuelMass - BurnRate * T,
EarthMass, GravConstant, Thrust, x );

}

// Freefall
for( ; x[ 0 ] >= EarthRadius; T += DeltaT )
{
double Height = x[ 0 ] - EarthRadius;
double XAirDensity = AirDensity * pow( 0.5, Height / AirDensityScale );
double Drag = XAirDensity * CrossSection * DragCoefficient;
std::cout << T << ", " << Height << ", " << x[ 1 ] << std::endl;
EvaluateTimestep( Drag, DeltaT, RocketMass,
EarthMass, GravConstant, 0.0, x );

}
}```