From e715b83d119db7ccdf5fb96a104693d4fafea62c Mon Sep 17 00:00:00 2001 From: torcs-ng Date: Thu, 2 Aug 2012 15:35:20 +0000 Subject: [PATCH] add celestialbody omit git-svn-id: https://svn.code.sf.net/p/speed-dreams/code/trunk@4829 30fe4595-0a0c-4342-8851-515496e4dcbd Former-commit-id: 1cf05a7c3bc5a759de6b9d21c05187bc4b5b4497 Former-commit-id: c7283a2e08ea39d86daf475130e417b53aabf458 --- src/libs/ephemeris/celestialbody.cpp | 160 +++++++++++++++++++++++++++ src/libs/ephemeris/celestialbody.h | 95 ++++++++++++++++ 2 files changed, 255 insertions(+) create mode 100644 src/libs/ephemeris/celestialbody.cpp create mode 100644 src/libs/ephemeris/celestialbody.h diff --git a/src/libs/ephemeris/celestialbody.cpp b/src/libs/ephemeris/celestialbody.cpp new file mode 100644 index 000000000..58b4ec85e --- /dev/null +++ b/src/libs/ephemeris/celestialbody.cpp @@ -0,0 +1,160 @@ +/**************************************************************************** + + file : celestialbody.cpp + created : Sun Aug 02 22:54:56 CET 2012 + copyright : (C) 1997 Durk Talsma - (C)2012 Xavier Bertaux + email : bertauxx@yahoo.fr + version : $Id: celestialbody.cpp 4811 2012-07-29 17:16:37Z torcs-ng $ + + ****************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +void ePhCelestialBody::updatePosition(double mjd, ePhStar *ourSun) +{ + double eccAnom, v, ecl, actTime, + xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze; + + updateOrbElements(mjd); + actTime = sdCalcActTime(mjd); + + // calcualate the angle bewteen ecliptic and equatorial coordinate system + ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime); + + eccAnom = sdCalcEccAnom(M, e); //calculate the eccentric anomaly + xv = a * (cos(eccAnom) - e); + yv = a * (sqrt (1.0 - e*e) * sin(eccAnom)); + v = atan2(yv, xv); // the planet's true anomaly + r = sqrt (xv*xv + yv*yv); // the planet's distance + + // calculate the planet's position in 3D space + xh = r * (cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i)); + yh = r * (sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i)); + zh = r * (sin(v+w) * sin(i)); + + // calculate the ecliptic longitude and latitude + xg = xh + ourSun->getxs(); + yg = yh + ourSun->getys(); + zg = zh; + + lonEcl = atan2(yh, xh); + latEcl = atan2(zh, sqrt(xh*xh+yh*yh)); + + xe = xg; + ye = yg * cos(ecl) - zg * sin(ecl); + ze = yg * sin(ecl) + zg * cos(ecl); + rightAscension = atan2(ye, xe); + declination = atan2(ze, sqrt(xe*xe + ye*ye)); + + //calculate some variables specific to calculating the magnitude + //of the planet + R = sqrt (xg*xg + yg*yg + zg*zg); + s = ourSun->getDistance(); + + // It is possible from these calculations for the argument to acos + // to exceed the valid range for acos(). So we do a little extra + // checking. + + double tmp = (r*r + R*R - s*s) / (2*r*R); + if ( tmp > 1.0) + { + tmp = 1.0; + } else if ( tmp < -1.0) + { + tmp = -1.0; + } + + FV = SGD_RADIANS_TO_DEGREES * acos( tmp ); +} + +double ePhCelestialBody::sdCalcEccAnom(double M, double e) +{ + double + eccAnom, E0, E1, diff; + + eccAnom = M + e * sin(M) * (1.0 + e * cos (M)); + // iterate to achieve a greater precision for larger eccentricities + if (e > 0.05) + { + E0 = eccAnom; + do + { + E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e *cos(E0)); + diff = fabs(E0 - E1); + E0 = E1; + } + while (diff > (SGD_DEGREES_TO_RADIANS * 0.001)); + return E0; + } + + return eccAnom; +} + +ePhCelestialBody::ePhCelestialBody(double Nf, double Ns, + double If, double Is, + double wf, double ws, + double af, double as, + double ef, double es, + double Mf, double Ms, double mjd) +{ + NFirst = Nf; NSec = Ns; + iFirst = If; iSec = Is; + wFirst = wf; wSec = ws; + aFirst = af; aSec = as; + eFirst = ef; eSec = es; + MFirst = Mf; MSec = Ms; + updateOrbElements(mjd); +} + +ePhCelestialBody::ePhCelestialBody(double Nf, double Ns, + double If, double Is, + double wf, double ws, + double af, double as, + double ef, double es, + double Mf, double Ms) +{ + NFirst = Nf; NSec = Ns; + iFirst = If; iSec = Is; + wFirst = wf; wSec = ws; + aFirst = af; aSec = as; + eFirst = ef; eSec = es; + MFirst = Mf; MSec = Ms; +} + +void ePhCelestialBody::updateOrbElements(double mjd) +{ + double actTime = sdCalcActTime(mjd); + M = SGD_DEGREES_TO_RADIANS * (MFirst + (MSec * actTime)); + w = SGD_DEGREES_TO_RADIANS * (wFirst + (wSec * actTime)); + N = SGD_DEGREES_TO_RADIANS * (NFirst + (NSec * actTime)); + i = SGD_DEGREES_TO_RADIANS * (iFirst + (iSec * actTime)); + e = eFirst + (eSec * actTime); + a = aFirst + (aSec * actTime); +} + +/* return value: the (fractional) number of days until Jan 1, 2000.*/ +/****************************************************************************/ +double ePhCelestialBody::sdCalcActTime(double mjd) +{ + return (mjd - 36523.5); +} + +void ePhCelestialBody::getPos(double* ra, double* dec) +{ + *ra = rightAscension; + *dec = declination; +} + +void ePhCelestialBody::getPos(double* ra, double* dec, double* magn) +{ + *ra = rightAscension; + *dec = declination; + *magn = magnitude; +} diff --git a/src/libs/ephemeris/celestialbody.h b/src/libs/ephemeris/celestialbody.h new file mode 100644 index 000000000..3ef23733e --- /dev/null +++ b/src/libs/ephemeris/celestialbody.h @@ -0,0 +1,95 @@ +/**************************************************************************** + + file : celestialbody.h + created : Sun Aug 02 22:54:56 CET 2012 + copyright : (C) 2000 Curtis L. Olson - (C)2012 Xavier Bertaux + email : bertauxx@yahoo.fr + version : $Id: celestialbody.h 4811 2012-07-29 17:16:37Z torcs-ng $ + + ****************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + + +#ifndef _CELESTIALBODY_H_ +#define _CELESTIALBODY_H_ + +#include + +class ePhStar; + +class ePhCelestialBody +{ +protected: // make the data protected, in order to give the + // inherited classes direct access to the data + double NFirst; /* longitude of the ascending node first part */ + double NSec; /* longitude of the ascending node second part */ + double iFirst; /* inclination to the ecliptic first part */ + double iSec; /* inclination to the ecliptic second part */ + double wFirst; /* first part of argument of perihelion */ + double wSec; /* second part of argument of perihelion */ + double aFirst; /* semimayor axis first part*/ + double aSec; /* semimayor axis second part */ + double eFirst; /* eccentricity first part */ + double eSec; /* eccentricity second part */ + double MFirst; /* Mean anomaly first part */ + double MSec; /* Mean anomaly second part */ + + double N, i, w, a, e, M; /* the resulting orbital elements, obtained from the former */ + + double rightAscension, declination; + double r, R, s, FV; + double magnitude; + double lonEcl, latEcl; + + double sdCalcEccAnom(double M, double e); + double sdCalcActTime(double mjd); + void updateOrbElements(double mjd); + +public: + ePhCelestialBody(double Nf, double Ns, + double If, double Is, + double wf, double ws, + double af, double as, + double ef, double es, + double Mf, double Ms, double mjd); + ePhCelestialBody(double Nf, double Ns, + double If, double Is, + double wf, double ws, + double af, double as, + double ef, double es, + double Mf, double Ms); + + void getPos(double *ra, double *dec); + void getPos(double *ra, double *dec, double *magnitude); + double getRightAscension(); + double getDeclination(); + double getMagnitude(); + double getLon() const; + double getLat() const; + void updatePosition(double mjd, ePhStar *ourSun); +}; + +inline double ePhCelestialBody::getRightAscension() { return rightAscension; } +inline double ePhCelestialBody::getDeclination() { return declination; } +inline double ePhCelestialBody::getMagnitude() { return magnitude; } + +inline double ePhCelestialBody::getLon() const +{ + return lonEcl; +} + +inline double ePhCelestialBody::getLat() const +{ + return latEcl; +} + +#endif // _CELESTIALBODY_H_ +