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);
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);
}
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
C:double alt = el.a * (1.0 - el.e*cos(eca)) - oapiGetSize(hRef);
Thank you! And how can I get dt for an angle or distance(ground track) to travel from now?
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 );
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;
}
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 );
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);
}
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);
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?
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);
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°?
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));
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);
}