Yet another propagation problem (atmospheric drag mostly)

[cristiapi]

If you don't need accuracy, your

result = sqrt(x*x + y*y + z*z);
// R3 - Earth radius = 6378.137
height = result - R3;

is good.

 

[Raccoon]

so if my calculating formula for height is good. Why with this results height is rising?

 

[cristiapi]

Because the satellite is in a point of the orbit where it goes from the perigee to apogee, the orbit is highly elliptical and the atmospheric drag is very small. In the long run, the orbit will decay (see my graph).

 

[Raccoon]

I understand what you mean. With only "g" height changes. But with atmospheric drag its only rising, how long must be "run"? Your graph shows only decreasing, it is solo atmospheric drag?

 

[cristiapi]

I understand what you mean. With only "g" height changes. But with atmospheric drag its only rising, how long must be "run"?

Not sure to understand.

Your graph shows only decreasing, it is solo atmospheric drag?

Yes.

 

[Raccoon]

With RK4 implementation of atmospheric drag will be like this?

void acc(Satellite &y, Satellite &f, double Mjd_TT,const QMatrix3x3 &T) {
    f.loc = y.vel;
    double Cd = 2.2;
    f.vel = AccelDrag(Mjd_TT, y.loc, y.vel, T, y.size, y.mass, Cd);
}

void RK4step(const double& h, Satellite& y, const double &Mjd_TT, const QMatrix3x3 &T) {
    Satellite k1,k2,k3,k4;
 
    Satellite yy;
    acc(y , k1, Mjd_TT, T); //for(int i= 0; i < 2; i++) yy[i]= y[i] + stp * .5 * k1[i];
    yy.loc = y.loc + h * .5 * k1.loc;
    yy.vel = y.vel + h * .5 * k1.vel;
    acc(yy, k2, Mjd_TT, T);//for(int i= 0; i < 2; i++) yy[i]= y[i] + stp * .5 * k2[i];
    yy.loc = y.loc + h * .5 * k2.loc;
    yy.vel = y.vel + h * .5 * k2.vel;
    acc(yy, k3, Mjd_TT, T);//for(int i= 0; i < 2; i++) yy[i]= y[i] + stp * k3[i];
    yy.loc = y.loc + h * .5 * k3.loc;
    yy.vel = y.vel + h * .5 * k3.vel;
    acc(yy, k4, Mjd_TT, T);//for(int i= 0; i < 2; i++) y[i] += stp/6 * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]);
    y.loc += h/6 * (k1.loc + 2 * k2.loc + 2 * k3.loc + k4.loc);
    y.vel += h/6 * (k1.vel + 2 * k2.vel + 2 * k3.vel + k4.vel);
}

and main loop is:

Satellite sat(QVector3D(7000000.0, 0.0, 500000.0), QVector3D(0.0, 9000.0, 0.0),           62.5, 8000.0);
    QVector< QVector3D > satellite;

    double Mjd_UTC;
    double Mjd_TT;
    QVector <double> height;
    double h = 10.0;
    QDateTime dt = QDateTime(QDate(2002,4,24),QTime(0,0,0,0));
    for (int j = 0.0; j < 5000; j++) {
        Mjd_UTC = MJD(dt);
        Mjd_TT  = Mjd_UTC + (32.0 + 32.184)/86400.0;

        QMatrix3x3 T = NutationMatrix(Mjd_TT) * PreccessionMatrix(MJD_J2000, Mjd_TT);
        RK4step(h, sat, Mjd_TT, T);
        dt = dt.addSecs(h);

        // R3 - Earth radius  = 6378.137
        height.push_back(Magnitude2(sat.loc) - R3*1000.0);
        satellite.push_back(sat.loc);
        satellite.push_back(sat.vel);
    }

 

[cristiapi]

I see an error here:

    acc(yy, k3, Mjd_TT, T);//for(int i= 0; i < 2; i++) yy[i]= y[i] + stp * k3[i];
    yy.loc = y.loc + h * .5 * k3.loc;
    yy.vel = y.vel + h * .5 * k3.vel;

where do you see ".5"? Anyway, my sample code is not good for you because you need the full version with the time:
https://en.wikipedia.org/wiki/Runge–Kutta_methods#The_Runge–Kutta_method
notice f(tn + h/2,...).
You could try to implement that code, but you should first understand what that code does.

 

[Raccoon]

I fix error.
Solo atmospheric drag acceleration:

Satellite acc(double const &t, Satellite &y) {
    Satellite f;
    f.loc = y.vel;
    double Cd = 2.2;
    QMatrix3x3 T = NutationMatrix(t) * PreccessionMatrix(MJD_J2000, t);
    f.vel = AccelDrag(t, y.loc, y.vel, T, y, y.size, y.mass, Cd);
    return f;
}

Runge kutta with time will be:

void RK4step(const double& h, Satellite& y, double &t){
    Satellite k1,k2,k3,k4, acceleration;
Satellite yy;
    k1 = acc(t, y) * h;    //for(int i= 0; i < 2; i++) yy[i]= y[i] + stp * .5 * k1[i];
    yy.loc = (y.loc + .5 * k1.loc);
    yy.vel = (y.vel + .5 * k1.vel);
    yy.mass = y.mass;
    yy.size = y.size;
    k2 = acc(t + 0.5 * h, yy) * h;    //for(int i= 0; i < 2; i++) yy[i]= y[i] + stp * .5 * k2[i];
    yy.loc = y.loc + h * .5 * k2.loc;
    yy.vel = y.vel + h * .5 * k2.vel;
    k3 = acc(t + 0.5 * h, yy) * h;    //for(int i= 0; i < 2; i++) yy[i]= y[i] + stp * k3[i];
    yy.loc = y.loc + h * k3.loc;
    yy.vel = y.vel + h * k3.vel;
    k4 = acc(t + h, yy) * h;     //for(int i= 0; i < 2; i++) y[i] += stp/6 * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]);
    y.loc += h/6 * (k1.loc + 2 * k2.loc + 2 * k3.loc + k4.loc);
    y.vel += h/6 * (k1.vel + 2 * k2.vel + 2 * k3.vel + k4.vel);
}

Main loop:

    double Mjd_UTC;
    double Mjd_TT;
    QVector <double> height;
    double h = 10.0;
    QDateTime dt = QDateTime(QDate(2002,4,24),QTime(0,0,0,0));
    for (int j = 0.0; j < 100; j++) {
        Mjd_UTC = MJD(dt);
        Mjd_TT  = Mjd_UTC + (32.0 + 32.184)/86400.0;

        RK4step(h, sat, Mjd_TT);
        dt = dt.addSecs(h);

        // R3 - Earth radius  = 6378.137
        height.push_back(Magnitude2(sat.loc) - R3*1000.0);
        satellite.push_back(sat.loc);
        satellite.push_back(sat.vel);
    }

Do I need send to RK modificated julian date (Mjd) or just send seconds starting, for example, from 0?
What atmospheric model do you use?
If its not hard can you post your implementation of solo atmospherid drag acceleration? I guess i dont calculate it right.

 

[cristiapi]

Your are still mixing code because you don't understand what RK4 does. Stop that random coding and take your time to understand RK4 without asking someone else to debug your random code.

 

[Raccoon]

Hm, RK find value between to point (interpolating) to get close to real result. Higher rank then higher accuracy. And my code not random, i guess you dont understand what I ask. I dont ask anyone to debug my code, i want explanation how to calculate solo perturbed effects.

 

[cristiapi]

It's really hard for me to understand what you write.
You may be interested in this implementation of RK4:
http://gafferongames.com/game-physics/integration-basics/