Modeling velocity and height of rocket c++

Question

Hi We are supposed to model the height and velocity of a rocket in c++ for our final project. Having user input the total flight time and delta time for the points during flight that they wish to measure. The following is the code I have written for this project. The velocity is supposed to start positive and after 60 seconds at which point there is no fuel left and thus no thrust the velocity should start becoming negative. However both my height and velocity are coming out as negative from the start and reaching negative infinite by the end.

#include <iostream>
using namespace std;
int main()
{
    float *v;
    float *h;
    float a;

    double mass=0.0, thrust, time, dt;
    double g = 32.2;
    double K = 0.008;

    cout << "enter time";
    cin >> time;
    cout << "enter dt";
    cin >> dt;

    a = (time/dt);
    v = new float[a];
    h = new float[a];

    v[0] = 0;
    h[0] = 0;
    float tt = 0;

    // for loop to calculate velocity and time
    for(int i = 0; i <= (time/dt) ; i++)
    {
        tt = dt + tt;

        if( tt <= 60)
        {
            mass = (3000-(40*tt)/g);
            thrust = 7000;
        }
        if ( tt > 60)
        {
            mass = 3000/g;
            thrust = 0;
        }

        // these are the formulas for velocity and height position our prof gave us
        v[i+1] = v[i] - (K/mass)*v[i]*v[i-1] * dt + (thrust/mass - g)*dt;
        h[i+1] = v[i+1] * dt + h[i];
    }

    // for loop to output
    for(int i = 0; i <= (time/dt); i++)
    {
        cout << i << " - " << "Velocity:" << v[i+1] << " Position:" << h[i+1] <<endl;
    }
    return 0;
}

sample output:

enter time120
enter dt.01
0 - Velocity:-0.298667 Position:-0.00298667
1 - Velocity:-0.597333 Position:-0.00896
2 - Velocity:-0.895999 Position:-0.01792
3 - Velocity:-1.19467 Position:-0.0298666
4 - Velocity:-1.49333 Position:-0.0448
5 - Velocity:-1.792 Position:-0.0627199
6 - Velocity:-2.09066 Position:-0.0836266
7 - Velocity:-2.38933 Position:-0.10752

<...i left out a lot of numbers in the middle to not make this post too long...>

11994 - Velocity:-inf Position:-inf
11995 - Velocity:-inf Position:-inf
11996 - Velocity:-inf Position:-inf
11997 - Velocity:-inf Position:-inf
11998 - Velocity:-inf Position:-inf
11999 - Velocity:-inf Position:-inf
12000 - Velocity:-inf Position:-inf
Program ended with exit code: 0

I have compared with my friends who are getting good results and we can not determine a difference between their code and mine. I have the rest of my program complete and working fine I just cannot figure out why my calculations are wrong

Answer 1

Thank you guys for your help I was able to solve it.

    #include <iostream>

    using namespace std;

   int main()
  {
  double *v;
double *h;
double g = 32.2;
double K = .008;
double mass;
double t;
double dt;
double tt = 0;
double thrust;

cout << "t \n";
cin >> t;
cout << "dt \n";
cin >> dt;

double a = t/dt;
v = new double[a];
h = new double[a];
v[0] = 0;
h[0] = 0;
for (int i = 0; i <= a; i++)
{

    tt = tt+dt;

    if( tt == 0)
    {
        thrust = 7000;
        mass = (3000 -40*dt)/g;
        v[i+1] = v[i] + ((thrust/mass)-g)*dt;
    }
    else if( tt > 0 && tt < 60)
    {
        thrust = 7000;
        mass = (3000 - 40 *tt)/g;
        v[i+1] = v[i] - ((K/mass)*v[i]*v[i-1] * dt) + ((thrust/mass) - g)*dt;
    }
    else if (tt > 60)
    {
        thrust = 0;
        mass = 600/g;
        v[i+1] = v[i] - ((K/mass)*v[i]*v[i-1] * dt) + ((thrust/mass) - g)*dt;
    }
    h[i+1] = v[i+1] * dt + h[i];

}



cout << " end results \n";
for(int i = 0; i <= a; i++)
{
    cout << i << " v - " << v[i] << " h - " << h[i] <<endl;

}

return 0;
}

New results:

t 120 dt .01 end results 0 v - 0 h - 0 1 v - 0.429434 h - 0.00429434 2 v - 0.858967 h - 0.012884 3 v - 1.2886 h - 0.02577 4 v - 1.71833 h - 0.0429534 5 v - 2.14817 h - 0.064435 6 v - 2.5781 h - 0.090216 7 v - 3.00813 h - 0.120297

You can see below at 60s where the velocity changes due to no more thrust 5997 - Velocity:890.361 5998 - Velocity:890.392 5999 - Velocity:890.422 6000 - Velocity:886.697 6001 - Velocity:882.985 6002 - Velocity:879.302

Answer 2

Ignoring the out-of-bounds access to v[-1] when i is zero, there is something wrong with your thrust, mass, or g.

thrust is 7000, mass is 3000 at time = 0. That means thrust/mass is just over 2. With g=32 (really? you are doing rocketry calculations in imperial units?), that means the rocket never has enough thrust to counter gravity, and just sits on the pad.

---

Edit: That would be reality. Because this is a fairly simple simulation, and doesn't include a "pad", in the model the rocket starts free-falling to the centre of the earth.

Answer 3

I am not 100% convinced about the forumla, I dont understand why it contains a v[i] and a v[i-1] term. Anyhow, even if correct, in the first iteration (i==0) you are accesing out of bounds of the velocity array: v[i-1]. That is undefined behaviour.

To fix this, either review the formula, does it really contain a v[i-1] term? ..or start the iteration at i=1 (and initialize v[0] and v[1]).

Answer 4

You are using v[i-1] but i starts out at 0, therefore this calculation is going to use whatever happens to be at v[-1]. I suggest you initialize i to 1 (and then check all the uses of i to ensure that the correct array elements will be used).