#include <functional>#include <algorithm>#include <cmath>#include <cassert>#include <celmath/mathlib.h>#include <celmath/solve.h>#include "astro.h"#include "orbit.h"#include "body.h"Include dependency graph for orbit.cpp:
Go to the source code of this file.
Typedefs | |
| typedef pair< double, double > | Solution |
Functions | |
| static EllipticalOrbit * | StateVectorToOrbit (const Point3d &position, const Vec3d &v, double mass, double t) |
|
|
|
|
||||||||||||||||||||
|
Definition at line 252 of file orbit.cpp. References cube(), astro::kilometersToAU(), Vector3< T >::length(), PI, radToDeg(), sin(), sqrt(), square(), Vector3< T >::x, Vector3< T >::y, and Vector3< T >::z. Referenced by MixedOrbit::MixedOrbit(). 00256 {
00257 Vec3d R = position - Point3d(0.0, 0.0, 0.0);
00258 Vec3d L = R ^ v;
00259 double magR = R.length();
00260 double magL = L.length();
00261 double magV = v.length();
00262 L *= (1.0 / magL);
00263
00264 Vec3d W = L ^ (R / magR);
00265
00266 double G = astro::G * 1e-9; // convert from meters to kilometers
00267 double GM = G * mass;
00268
00269 // Compute the semimajor axis
00270 double a = 1.0 / (2.0 / magR - square(magV) / GM);
00271
00272 // Compute the eccentricity
00273 double p = square(magL) / GM;
00274 double q = R * v;
00275 double ex = 1.0 - magR / a;
00276 double ey = q / sqrt(a * GM);
00277 double e = sqrt(ex * ex + ey * ey);
00278
00279 // Compute the mean anomaly
00280 double E = atan2(ey, ex);
00281 double M = E - e * sin(E);
00282
00283 // Compute the inclination
00284 double cosi = L * Vec3d(0, 1.0, 0);
00285 double i = 0.0;
00286 if (cosi < 1.0)
00287 i = acos(cosi);
00288
00289 // Compute the longitude of ascending node
00290 double Om = atan2(L.x, L.z);
00291
00292 // Compute the argument of pericenter
00293 Vec3d U = R / magR;
00294 double s_nu = (v * U) * sqrt(p / GM);
00295 double c_nu = (v * W) * sqrt(p / GM) - 1;
00296 s_nu /= e;
00297 c_nu /= e;
00298 Vec3d P = U * c_nu - W * s_nu;
00299 Vec3d Q = U * s_nu + W * c_nu;
00300 double om = atan2(P.y, Q.y);
00301
00302 // Compute the period
00303 double T = 2 * PI * sqrt(cube(a) / GM);
00304 T = T / 86400.0; // Convert from seconds to days
00305
00306 #if 0
00307 cout << "a: " << astro::kilometersToAU(a) << '\n';
00308 cout << "e: " << e << '\n';
00309 cout << "i: " << radToDeg(i) << '\n';
00310 cout << "Om: " << radToDeg(Om) << '\n';
00311 cout << "om: " << radToDeg(om) << '\n';
00312 cout << "M: " << radToDeg(M) << '\n';
00313 cout << "T: " << T << '\n';
00314 #endif
00315
00316 return new EllipticalOrbit(a * (1 - e), e, i, Om, om, M, T, t);
00317 }
|
1.4.1