Thread: Model Rocket Altitude predictor...

  1. #1
    Registered User
    Join Date
    Aug 2009
    Posts
    19

    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; //radius
    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. #2
    Hurry Slowly vart's Avatar
    Join Date
    Oct 2006
    Location
    Rishon LeZion, Israel
    Posts
    6,788
    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
    All problems in computer science can be solved by another level of indirection,
    except for the problem of too many layers of indirection.
    – David J. Wheeler

  3. #3
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,336
    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. #4
    Registered User
    Join Date
    Aug 2009
    Posts
    19
    Quote Originally Posted by tabstop View Post
    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. #5
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,336
    Quote Originally Posted by kalor_alros View Post
    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. #6
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    39,660
    > 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.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.

  7. #7
    Registered User
    Join Date
    Aug 2009
    Posts
    19
    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. #8
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,336
    Quote Originally Posted by kalor_alros View Post
    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. #9
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    39,660
    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.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.

  10. #10
    Officially An Architect brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,396
    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
    Code:
    //try
    //{
    	if (a) do { f( b); } while(1);
    	else   do { f(!b); } while(1);
    //}

  11. #11
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,336
    Quote Originally Posted by brewbuck View Post
    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. #12
    Officially An Architect brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,396
    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 );
    	
        }
    }
    Code:
    //try
    //{
    	if (a) do { f( b); } while(1);
    	else   do { f(!b); } while(1);
    //}

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. questions on multiple thread programming
    By lehe in forum C Programming
    Replies: 11
    Last Post: 03-27-2009, 07:44 AM
  2. Replies: 14
    Last Post: 06-28-2006, 01:58 AM
  3. Model not showing up on screen...
    By Shamino in forum Game Programming
    Replies: 14
    Last Post: 03-09-2006, 08:00 PM
  4. multiple file loading. so fruturated! help!
    By psychopath in forum Game Programming
    Replies: 5
    Last Post: 05-09-2005, 05:13 PM
  5. Help with file reading/dynamic memory allocation
    By Quasar in forum C++ Programming
    Replies: 4
    Last Post: 05-17-2004, 03:36 PM