You are using an out of date browser. It may not display this or other websites correctly.

You should upgrade or use an alternative browser.

You should upgrade or use an alternative browser.

- Thread starter Topper
- Start date

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

and you would need

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:

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);
}
```

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:

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

if the orbit is

Last edited:

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 );
```

- Joined
- Mar 28, 2008

- Messages
- 653

- Reaction score
- 10

- Points
- 33

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

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 );
```

- 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)....

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?

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

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);
}
```

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:

That is possible but the amount of math would be high. But sooner or later you would probably need it anyway.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?

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.

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:

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.

But it's recommended that you write your own instead just to better understand what you are doing and why.

Last edited:

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);
}
```

- Replies
- 27

- Views
- 394

- Replies
- 12

- Views
- 547

- Replies
- 1

- Views
- 131