API Question Getting orbit altitude for a given time from Orbital Elements

[Topper]

Hello, is there a more or less simple way to get the altitude of an orbit for a given number of seconds from now from the orbital elements from the Orbiter API?

 

[jarmonik]

I haven't tested this but you could start by using following API function
GetElements (OBJHANDLE hRef, ELEMENTS &el, ORBITPARAM *prm=0, double mjd_ref=0, int frame=FRAME_ECL) const

and you would need mna2eca() function from this post: link to a post

inline double limit(double x)
{
    if (x>PI2) return fmod(x,PI2); if (x<0) return PI2-fmod(-x,PI2); return x;
}

double mna = prm->MnA + PI2 * time / prm->T;
double eca = mna2eca(limit(mna), el.e);
double alt = el.a * (1.0 - el.e*cos(eca)) - oapiGetSize(hRef);

 

[asbjos]

For elliptical orbit (eccentricity < 1) and spherical planet, see jarmonik's answer. Note that his limit function is already built into Orbiter as the posangle() function.

For hyperbolic orbit and spherical planet, you have to alter the propagation of MnA (mean anomaly) and EcA (eccentric anomaly):

if (el.e > 1.0)
{
    double meanMotion = sqrt(-oapiGetMass(hRef) * GGRAV / (el.a * el.a * el.a)); // el.a is negative for hyperbolic orbit
    double MnAinXseconds = prm.MnA + meanMotion * time;
    double EcAinXseconds = MnA2EcA(MnAinXseconds, el.e); // note that jarmonik's function takes care of hyperbolic case :-)
  
    double alt = el.a * (1.0 - el.e * cos(EcAinXseconds)) - oapiGetSize(hRef);
}

Taking terrain elevation into account is a bit more difficult, as you have to find the longitude and latitude of the point to find the terrain elevation.
Use the same methods as previously (either elliptic or hyperbolic orbit), but then after finding the eccentric anomaly (EcA):

double TrAinXseconds; // TrA is true anomaly
if (el.e < 1.0) TrAinXseconds = 2.0 * atan(sqrt((1.0 + el.e) / (1.0 - el.e)) * tan(EcAinXseconds / 2.0)); // elliptic
else TrAinXseconds = 2.0 * atan(sqrt((el.e + 1.0) / (el.e - 1.0)) * tanh(EcAinXseconds / 2.0)); // hyperbolic

double latitude = asin(sin(el.i) * sin(el.omegab - el.theta + TrAinXseconds));
double currentLongitudeOfAscencion = el.theta - oapiGetPlanetCurrentRotation(hRef);
double angleFromAscendingNode = atan2(cos(el.i) * sin(el.omegab - el.theta + TrAinXseconds), cos(el.omegab - el.theta + TrAinXseconds));
double longitude = normangle(currentLongitudeOfAscencion - PI2 / oapiGetPlanetPeriod(hRef) * time + angleFromAscendingNode);

double alt = el.a * (1.0 - el.e * cos(EcAinXseconds)) - oapiGetSize(hRef) - oapiSurfaceElevation(hRef, longitude, latitude); // radius - planetRad - mountain height

Footnote: this does not take non-spherical gravity into effect (which you can enable/disable in Orbiter). So for Earth, your position can be up to ~10 kilometers off per orbit in the future you're predicting. For the general terrain this is a negligible error, but for a mountain this can give a quite wrong result.
Footnote2: technically, the calculation of TrA breaks down for orbit with eccentricity exactly equal to 1. But due to how computers store numbers and orbital perturbations, you are extremely unlikely to encounter an exact 1. And if it happens, it will be perturbed away from 1 by the next frame, so it will not be a lasting problem.

 

[Topper]

Thank you! And how can I get dt for an angle or distance(ground track) to travel from now?

 

[jarmonik]

double alt = el.a * (1.0 - el.e*cos(eca)) - oapiGetSize(hRef);

if the orbit is hyperbolic then cosh() should be used instead of cos() because the function is using eccentric anomaly instead of true anomaly.

 

[jarmonik]

Thank you! And how can I get dt for an angle or distance(ground track) to travel from now?

Once again, I haven't tested this.

double mnm = sqrt(oapiGetMass(hRef) * GGRAV / fabs(el.a * el.a * el.a));
double eca = tra2eca(posangle(prm->TrA + angle), el.e);
double dm = posangle(eca2mna(eca, el.e) - prm->MnA);

double dt = dm * mnm; // Travel time for "angle" less than one orbit

// For multiple orbits this probably should be:
double dt = ( (dm * mnm) + floor(angle/PI2) * PI2 * mnm );

 

[Topper]

Ok I testet the following:

double SurfaceSpeedMFD::getAltitudeOverBase()
{
    
    OBJHANDLE hRef = mfdi->v->GetGravityRef();
    ORBITPARAM *prm = new ORBITPARAM;
    ELEMENTS el;

    double angle = mfdi->distanceToBasePad / oapiGetSize(hRef); // mfdi->distanceToBasePad  is ground track distance
    
    mfdi->v->GetElements(hRef, el, prm);

    const double mnm = sqrt(oapiGetMass(hRef) * GGRAV / fabs(el.a * el.a * el.a));    
    double eca = ORBITERTOOLS::tra2eca(posangle(prm->TrA + angle), el.e);
    const double dm = posangle(ORBITERTOOLS::eca2mna(eca, el.e) - prm->MnA);
    const double dt = dm * mnm; // Travel time for "angle" less than one orbit

    sprintf(oapiDebugString(), "Angle: %f  dt: %f", angle toDeg, dt);

    const double mna = prm->MnA + PI2 * dt / prm->T;
    eca = ORBITERTOOLS::mna2eca(limit(mna), el.e);
    double alt = el.a * (1.0 - el.e*cosf(eca)) - oapiGetSize(hRef);
    return alt;
}

which shows a very low value for dt...

 

[jarmonik]

If the mean motion is defined like that then one should actually divide by it (not multiply) in this case

double dt = dm / mnm; // Travel time for "angle" less than one orbit
// For multiple orbits this probably should be:
double dt = ( (dm / mnm) + floor(angle/PI2) * PI2 / mnm );

 

[Topper]

Ha thanks again,

now I'm able to draw the flight path (yellow line) easy (before I used a simple, numerical euler integration) and to calculate importantant altitude information
It's much faster now!

Now I can think about to do the burn to flatten the flight path for final approach automaticly for the baseland autopilot ...
(Actual I used basesynch MFD to get the alt iformation and I did the "deorbit burn" manualy)....

double SurfaceSpeedMFD::getAltitudeAtDistance(double distance)
{
    double dt = 0;
    return getAltitudeAtDistance(distance, &dt);
}

double SurfaceSpeedMFD::getAltitudeAtDistance(double distance, double*dt)
{
    OBJHANDLE hRef = mfdi->v->GetGravityRef();
    ORBITPARAM *prm = new ORBITPARAM;
    ELEMENTS el;
    double angle = distance / oapiGetSize(hRef);
    mfdi->v->GetElements(hRef, el, prm);
    const double mnm = sqrt(oapiGetMass(hRef) * GGRAV / fabs(el.a * el.a * el.a));
    double eca = ORBITERTOOLS::tra2eca(posangle(prm->TrA + angle), el.e);
    const double dm = posangle(ORBITERTOOLS::eca2mna(eca, el.e) - prm->MnA);
    *dt = dm / mnm; // Travel time for "angle" less than one orbit
    const double mna = prm->MnA + PI2 * (*dt) / prm->T;
    eca = ORBITERTOOLS::mna2eca(posangle(mna), el.e);
    return el.a * (1.0 - el.e*cosf(eca)) - oapiGetSize(hRef);;
}

void SurfaceSpeedMFD::drawFlightPath(oapi::Sketchpad *skp, const double shipAltY)
{
    skp->SetPen(stdPenYellow);
    skp->MoveTo(width / 2, shipAltY);
    for (int s = 0; s < mfdi->distanceToBasePad && s/distancePerPixel < width; s += distancePerPixel)
    {
        const double x = (width / 2) + s / distancePerPixel;
        const double y = shipAltY - (getAltitudeAtDistance(s) - oapiGetFocusInterface()->GetAltitude()) / distancePerPixel;
        skp->LineTo(x, y);
    }
    const double altAtMaxElevation = getAltitudeAtDistance(mfdi->distanceToMaxElevation);
    double ttbase;
    const double altAtBase = getAltitudeAtDistance(mfdi->distanceToBasePad, &ttbase);
    sprintf(oapiDebugString(), "ALTITUDE AT MAX ELEVATION: %f      ALTITUDE AT BASE: %f   TIME: %f", altAtMaxElevation - mfdi->maxElevation, altAtBase - mfdi->baseElevation, ttbase);
}

One point still missing is to get the offplane value to the base, so is it also possible to get latitude and longitude for a given time? Then I just need to calculate the distance between those and the base coordinates by this?

double latitude = asin(sin(el.i) * sin(el.omegab - el.theta + TrAinXseconds));
double currentLongitudeOfAscencion = el.theta - oapiGetPlanetCurrentRotation(hRef);
double angleFromAscendingNode = atan2(cos(el.i) * sin(el.omegab - el.theta + TrAinXseconds), cos(el.omegab - el.theta + TrAinXseconds));
double longitude = normangle(currentLongitudeOfAscencion - PI2 / oapiGetPlanetPeriod(hRef) * time + angleFromAscendingNode);

And what are the most efficent points to reduce alt and correct offplane?
I guess for alt it's 180° before base and for correct offplane 90°?

Sorry I had no idea how to use such things like orbital elements for such calculations before and I don't understand what all that stuff means but I try to learn it, the language is also a litlle handicap because most is descriped in english beacuse it's realy special stuff... ;-)...


[edit]
I just noticed that I have this little devation of 0.3s between my calculation and basesynch, I guess it's the rotational speed of the moon, It's a mess...

 

[jarmonik]

One point still missing is to get the offplane value to the base, so is it also possible to get latitude and longitude for a given time? Then I just need to calculate the distance between those and the base coordinates by this?

That is possible but the amount of math would be high. But sooner or later you would probably need it anyway.

I usually prefer vectors over angles and trigonometry. You could use oapiEquToLocal() to get the position vector for the base. I don't recall if oapiEqutoLocal() takes the planet rotation in count or not. If not then you would need to use oapiGetRotationMatrix() and vLoc = mul(mRot, vLoc)

inline double angle(const VECTOR3 &v, const VECTOR3 &h)
{
    double x = dotp(unit(v), unit(h));
    if (x>=1.0)  return 0.0; else if (x<=-1.0) return PI;
    return( acos( x ) );
}

VECTOR3 h = unit(crossp(rpos, rvel));  // Orbit Normal
double off_plane_angle = PI05 - angle(h, vLoc);

Of course, to be more precise you would need to rotate vLoc even more to compensate planet rotation during the time it takes to get there. rotm(Axis, Angle) function would likely do the trick.

And what are the most efficent points to reduce alt and correct offplane?
I guess for alt it's 180° before base and for correct offplane 90°?

Correct, assuming that the initial orbit is near circular.

 

[jarmonik]

There is a pretty good Orbit toolkit in this example project of mine: https://www.orbiter-forum.com/threads/orbits-came-for-the-orbiter.33870/
But it's recommended that you write your own instead just to better understand what you are doing and why.

 

[jarmonik]

You can compute peri-focal reference frame _P, _Q, _H where the _H is often labelled as _W in literature.

double mu = oapiGetMass(hRef) * GGRAV;

VECTOR3 _P  = unit(( ( _pos * (length2(_vel) - mu / length(_pos)) ) - (_vel * dotp(_pos, _vel)) ) / mu);
VECTOR3 _H  = unit(crossp(_pos, _vel));
VECTOR3 _Q  = unit(crossp(_H, _P));

The equation for _P is found for Orbiter.pdf

VECTOR3 PositionByEcA(double eca) const
{
    if (ecc < 1.0) return _P * ((cos(eca)-ecc)*sma) + _Q * (sin(eca)*smi);
    return _P * ((cosh(eca)-ecc)*sma) + _Q * (sinh(eca)*smi);
}

VECTOR3 VelocityByTrA(double tra) const
{
    double up = sqrt(oapiGetMass(hRef) * GGRAV / par);
    return _Q * ((ecc + cos(tra)) * up) - _P * (sin(tra) * up);
}