diff --git a/src/orbits.cpp b/src/orbits.cpp index 9ee1fee..b5e0560 100644 --- a/src/orbits.cpp +++ b/src/orbits.cpp @@ -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)); diff --git a/src/orbits.h b/src/orbits.h index c9ed76c..83c41c4 100644 --- a/src/orbits.h +++ b/src/orbits.h @@ -82,14 +82,32 @@ ellipseInitAB(double a, double b); EllipseParameters ellipseInitAE(double a, double e); -double -ellipseGetEccentricity(double a, double p); - -bool -ellipseValidate(const EllipseParameters& ep); +inline double +ellipseGetEccentricity(double a, double p) +{ + return sqrt(1 - p / a); +} -bool -ellipsesEqual(EllipseParameters& e1, EllipseParameters& e2); +inline 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); +} + +inline bool +ellipsesEqual(EllipseParameters& e1, EllipseParameters& e2) +{ + return (e1.a == e2.a && + e1.b == e2.b && + e1.e == e2.e); +} // NOTE: create vertices for a 3d ellipse // NOTE: all vertices are in the x/y plane with z = 0 @@ -102,14 +120,72 @@ ellipse3DUpdate(EllipseParameters ep, Ellipse3D& e3d); OrbitalElements orbitInit(double a, double e); -double -orbitGetVelocity(double epsilon, double mu, double r); +inline double +orbitGetAngularMomentum(double p, double mu) +{ + return sqrt(mu * p); +} -double -orbitGetFlightPathAngle(double e, double true_anom); +inline double +orbitGetSemiMajorAxis(double epsilon, double mu) +{ + return -1 * mu / (2 * epsilon); +} -double -orbitGetPeriod(double a, double mu); +inline double +orbitGetSemiLatusRectum(double h, double mu) +{ + return pow(h, 2) / mu; +} + +inline double +orbitGetAngularMomentumFromStateVectors(double r, double v, double gamma) +{ + return r * v * cos(gamma); +} + +inline double +orbitGetSpecificEnergy(double a, double mu) +{ + return -1 * mu / (2 * a); +} + +inline double +orbitGetSpecificEnergyFromStateVectors(double r, double v, double mu) +{ + return pow(v, 2) / 2 - mu / r; +} + +inline double +orbitGetVelocity(double epsilon, double mu, double r) +{ + return sqrt(2 * (epsilon + mu / r)); +} + +inline double +orbitGetFlightPathAngle(double e, double true_anom) +{ + return atan(e * sin(true_anom) / (1 + e * cos(true_anom))); +} + +inline double +orbitGetPeriod(double a, double mu) +{ + return 2 * M_PI / sqrt(mu) * sqrt(pow(a, 3)); +} + +// NOTE: aka) the trajectory equation (eq. 2.45) +inline double +orbitGetRadialDistance(double e, double p, double true_anom) +{ + return p / (1 + e * cos(true_anom)); +} + +inline glm::vec2 +polarToRect(double angle, double r) +{ + return glm::vec2(r * cos(angle), r * sin(angle)); +} void orbitUpdate(OrbitalElements& o, double a, double e); @@ -153,30 +229,5 @@ getPropagatedTrueAnomaly(TwoBodySystem sys, double initial_anom, double time_step); // NOTE: in seconds -// NOTE: aka) the trajectory equation (eq. 2.45) -double -orbitGetRadialDistance(double e, double p, double true_anom); - -double -orbitGetAngularMomentum(double p, double mu); - -double -orbitGetSemiMajorAxis(double epsilon, double mu); - -double -orbitGetSemiLatusRectum(double h, double mu); - -double -orbitGetAngularMomentumFromStateVectors(double r, double v, double gamma); - -double -orbitGetSpecificEnergy(double a, double mu); - -double -orbitGetSpecificEnergyFromStateVectors(double r, double v, double mu); - -glm::vec2 -polarToRect(double true_anom, double r); - double orbitGetTimeOfFlight(TwoBodySystem sys, double theta_begin, double theta_end);