simple example applying perturbed effect on orbit propogation?

Question

I need to apply perturbation effect to orbit propagation. I am programmer not a physicist. I’ve already wrote program to propagate not perturbation orbit and show satellite ground track. So I know how to find state vector (RV – position, velocity) from orbital elements also state vector in PQW and IJK frames. But with those formulas I hardly understand when and where apply perturbation effects. At this moment I understand:

Adrag = -1/2*p*(Cd*A)/m * sqr(v)*Iv
Where:
p – atmospheric density;
Cd – drag coefficient;
A – crosssectional area of satellite perpendicular to the velocity vector;
m – mass of the satellite;
v – velocity of the satellite relative to the atmosphere;
Iv – unit vector in the direction of the satellite’s velocity.

This formula give me an atmospheric drag acceleration vector. Then I must apply acceleration vector to position vector R and I will get perturbed satellite position. After I need to apply acceleration vector to some formulas to change orbital elements. Necessary formulas I found in “Fundamentals of astrodynamics and applications” by David A. Vallado, paragraph “8.3.2 Gaussian VOP (Nonconservative and Conservative Effects)”. Also there is algorithm in “8.7.1 Application: Perturbed Two-body Propogation” but as I understand I cant use them because I don’t know both mean motion derivatives, I find out that they only can be known in TLE and can’t be calculated with orbital elements. But I don’t sure if they are what I need.

Can someone write simple example applying atmospheric drag (or another perturbed effect) to unperturbed orbit? Or if its hard, can you advise me a forum where I can ask such question and get answer?

Answer 1

First you need input data. Input can be presented as orbit elements (semimajor axis, eccentricity, inclination, true anomaly, longitude of ascending node, perigee), or state vector (position and velocity). If you choose orbit elements they need to be translated to state vector in ECI coordinates (find orbit to ECI).

Second choose integrator for propagation (three types: analytic, semi-analytic, numerical). If you choose numerical propagator as I did. Now choose integrator (a lot of them) but with order higher than Euler, because it doesn’t provide necessary accuracy. For simple example and understanding I choose Runge Kutta 4 (RK4).

Third, choose perturbed effects. The most simple of them is gravity (find 2-Body problem) satellite_position * -MU / Magnitude(satellite_position)^3 Gravity effect used in all perturbed orbit calculations where MU – Earth gravity = 398600.4415. In most cases when you see “only-drag”, “only-solar radiation pressure” etc, it’s mean that also include gravity, so technically they are not “only”.

When choosing perturbed effect you must decide how to calculate it, for example even gravity has different models, or in magnetic Earth field model. Atmospheric drag can be calculated in two ways (I find only two):

1) Calculation depending on time like described in “Montenbruck, E. Gill; Satellite Orbits - Models, Methods, and Applications”. And also choose model like Harris-Priester

2) Calculation with already calculated table of atmosphere density values, like did these gentleman https://github.com/komrad36 in one of his programs. For table values you also can choose any “already tables” you like: “U.S. Standard Atmosphere”, or Russian “ГОСТ 4401-81” etc.

I choose table version. My implementation in c++, main loop:

QDateTime dt = QDateTime (QDate(2018,05,11),QTime(14,30,30));
double JD = GetJulDate(dt);
double Mjd = MJD(dt);
double h = 0.5 // integration step
for (double tCur = tBegin; tCur < tEnd; tCur += tStep){
        // for table version
        RK4 (h, sat);
        // for time depend version
    //RK4 (h, sat, Mjd);
}

For integration step h is better to choose 0.5. Mjd – modified Julian date. sat – is object which contains satellite data (in my case only position, velocity, size and mass). Runge kutta 4 with time t, if you don’t use time just send only 1 parameter in function “acceleration”:

void RK4(const double& h, Satellite& sat, double &t){
    Satellite k1,k2,k3,k4;
    double MJDstep = (0.5 * h)/Tm;
    Satellite yy;
    yy = sat;

    k1 = acceleration (t, yy) * h;
    yy.loc = sat.loc + .5 * k1.loc;
    yy.vel = sat.vel + .5 * k1.vel;

    k2 = acceleration (t + MJDstep, yy) * h;
    yy.loc = sat.loc + h * .5 * k2.loc;
    yy.vel = sat.vel + h * .5 * k2.vel;

    k3 = acceleration (t + MJDstep, yy) * h;
    yy.loc = sat.loc + h * k3.loc;
    yy.vel = sat.vel + h * k3.vel;

    k4 = acceleration (t + h/Tm, yy) * h;
    sat.loc += h/6 * (k1.loc + 2 * k2.loc + 2 * k3.loc + k4.loc);
    sat.vel += h/6 * (k1.vel + 2 * k2.vel + 2 * k3.vel + k4.vel);
}

Line MJDstep = (0.5 * h)/Tm is prestep calculation where Tm = 1440. This is for translate time in seconds to modified Julian date. Or any other time unit of your chose. And acceleration function:

Satellite acceleration (double const &t, Satellite &sat) {
    Satellite f;
    f.loc = sat.vel;
    double p = Magnitude(sat.loc);
    f.vel = sat.loc * -MU /(p*p*p);

    double Cd = 2.2;
    f.vel += AccelDrag(sat, Cd);
    // for time version
    // f.vel += AccelDrag2(sat, Cd, t);
    return f;
}

the last but not the least (and obvious) use same units everywhere. If you choose meters translate length units to meters. If seconds for time translate all time to seconds.