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

Topper

Addon Developer
Addon Developer
Donator
Joined
Mar 28, 2008
Messages
653
Reaction score
10
Points
33
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

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,390
Reaction score
461
Points
83
Website
users.kymp.net
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

C:
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);
 
Last edited:

asbjos

tuanibrO
Addon Developer
Joined
Jun 22, 2011
Messages
696
Reaction score
254
Points
78
Location
This place called "home".
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):
C++:
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):
C++:
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.
 
Last edited:

Topper

Addon Developer
Addon Developer
Donator
Joined
Mar 28, 2008
Messages
653
Reaction score
10
Points
33
Thank you! And how can I get dt for an angle or distance(ground track) to travel from now?
 

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,390
Reaction score
461
Points
83
Website
users.kymp.net
C:
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.
 
Last edited:

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,390
Reaction score
461
Points
83
Website
users.kymp.net
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.

C:
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

Addon Developer
Addon Developer
Donator
Joined
Mar 28, 2008
Messages
653
Reaction score
10
Points
33
Ok I testet the following:

Code:
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...

1625087377551.png
 

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,390
Reaction score
461
Points
83
Website
users.kymp.net
If the mean motion is defined like that then one should actually divide by it (not multiply) in this case

C:
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

Addon Developer
Addon Developer
Donator
Joined
Mar 28, 2008
Messages
653
Reaction score
10
Points
33
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)....

Code:
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);
}

1625173977441.png

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?
Code:
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...
 
Last edited:

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,390
Reaction score
461
Points
83
Website
users.kymp.net
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)

C:
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 ) );
}

C:
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.
 
Last edited:

jarmonik

Well-known member
Orbiter Contributor
Addon Developer
Beta Tester
Joined
Mar 28, 2008
Messages
2,390
Reaction score
461
Points
83
Website
users.kymp.net
You can compute peri-focal reference frame _P, _Q, _H where the _H is often labelled as _W in literature.

C:
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


Code:
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);
}

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