|
|
|
|
@ -52,33 +52,6 @@ ellipseInitAE(double a, double e)
|
|
|
|
|
return ellipseInitAB(a, b); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
ellipseGetEccentricity(double a, double p) |
|
|
|
|
{ |
|
|
|
|
return sqrt(1 - p / a); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
bool |
|
|
|
|
ellipseValidate(const EllipseParameters& ep) |
|
|
|
|
{ |
|
|
|
|
// TODO: find out why satellite position gets wonky with orbit
|
|
|
|
|
// eccentricity > 0.995 while passing through true anom = 0.
|
|
|
|
|
// maybe divide by 0, or some floating point error?
|
|
|
|
|
return (ep.a > 0 && |
|
|
|
|
ep.b > 0 && |
|
|
|
|
ep.a >= ep.b && |
|
|
|
|
ep.e >= 0 && |
|
|
|
|
ep.e < 0.995); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
bool |
|
|
|
|
ellipsesEqual(EllipseParameters& e1, EllipseParameters& e2) |
|
|
|
|
{ |
|
|
|
|
return (e1.a == e2.a && |
|
|
|
|
e1.b == e2.b && |
|
|
|
|
e1.e == e2.e); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
Ellipse3D |
|
|
|
|
ellipseInit3D(EllipseParameters ep, uint vert_count) |
|
|
|
|
{ |
|
|
|
|
@ -104,72 +77,6 @@ orbitInit(double a, double e)
|
|
|
|
|
return o; |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetRadialDistance(double e, double p, double true_anom) |
|
|
|
|
{ |
|
|
|
|
return p / (1 + e * cos(true_anom)); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetAngularMomentum(double p, double mu) |
|
|
|
|
{ |
|
|
|
|
return sqrt(mu * p); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetSemiMajorAxis(double epsilon, double mu) |
|
|
|
|
{ |
|
|
|
|
return -1 * mu / (2 * epsilon); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetSemiLatusRectum(double h, double mu) |
|
|
|
|
{ |
|
|
|
|
return pow(h, 2) / mu; |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetAngularMomentumFromStateVectors(double r, double v, double gamma) |
|
|
|
|
{ |
|
|
|
|
return r * v * cos(gamma); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetSpecificEnergy(double a, double mu) |
|
|
|
|
{ |
|
|
|
|
return -1 * mu / (2 * a); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetSpecificEnergyFromStateVectors(double r, double v, double mu) |
|
|
|
|
{ |
|
|
|
|
return pow(v, 2) / 2 - mu / r; |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetVelocity(double epsilon, double mu, double r) |
|
|
|
|
{ |
|
|
|
|
return sqrt(2 * (epsilon + mu / r)); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetFlightPathAngle(double e, double true_anom) |
|
|
|
|
{ |
|
|
|
|
return atan(e * sin(true_anom) / (1 + e * cos(true_anom))); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
orbitGetPeriod(double a, double mu) |
|
|
|
|
{ |
|
|
|
|
return 2 * M_PI / sqrt(mu) * sqrt(pow(a, 3)); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
glm::vec2 |
|
|
|
|
polarToRect(double angle, double r) |
|
|
|
|
{ |
|
|
|
|
return glm::vec2(r * cos(angle), r * sin(angle)); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
void |
|
|
|
|
ellipse3DUpdate(EllipseParameters ep, Ellipse3D& e3d) |
|
|
|
|
{ |
|
|
|
|
@ -189,50 +96,50 @@ ellipse3DUpdate(EllipseParameters ep, Ellipse3D& e3d)
|
|
|
|
|
//
|
|
|
|
|
|
|
|
|
|
// FIXME: variable names 'ecc' should be just 'e'
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
getEccAnomFromTrueAnom(double ecc, double true_anom) |
|
|
|
|
{ |
|
|
|
|
return 2 * atan(sqrt((1 - ecc) / (1 + ecc)) * tan(true_anom / 2)); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
getTrueAnomFromEccAnom(double ecc, double ecc_anom) |
|
|
|
|
{ |
|
|
|
|
return 2 * atan(sqrt((1 + ecc) / (1 - ecc)) * tan(ecc_anom / 2)); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
getMeanAnomFromEccAnom(double ecc_anom, double ecc) |
|
|
|
|
{ |
|
|
|
|
return ecc_anom - ecc * sin(ecc_anom); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
getMeanMotion(double mu, double a) |
|
|
|
|
{ |
|
|
|
|
return sqrt(mu / pow(a, 3)); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
getPropagatedMeanAnom(double mean_anom, double mean_motion, double time_step) |
|
|
|
|
{ |
|
|
|
|
return mean_anom + mean_motion * (time_step); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
getInitialTrialValue(double mean_anom, double ecc) |
|
|
|
|
{ |
|
|
|
|
return mean_anom + ecc * sin(mean_anom) |
|
|
|
|
+ ((pow(ecc, 2) / 2) * sin(2 * mean_anom)); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
getTrialError(double ecc, double test_anom, double mean_anom) |
|
|
|
|
{ |
|
|
|
|
return test_anom - ecc * sin(test_anom) - mean_anom; |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
getNextTrialValue(double err, double ecc, double test_anom, double mean_anom) |
|
|
|
|
{ |
|
|
|
|
// compute derivative of the error function
|
|
|
|
|
@ -288,7 +195,7 @@ getPropagatedTrueAnomaly(TwoBodySystem sys,
|
|
|
|
|
// FIXME: organize into interface/internal functions
|
|
|
|
|
|
|
|
|
|
// NOTE: return an equivalent angle between 0 and 2 pi radians
|
|
|
|
|
double |
|
|
|
|
inline double |
|
|
|
|
orbitClampAngle(double theta) |
|
|
|
|
{ |
|
|
|
|
double revs = floor(theta / (2 * M_PI)); |
|
|
|
|
|