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
This commit is contained in:
torcs-ng
parent
08e2b6a1a6
commit
e715b83d11
2 changed files with 255 additions and 0 deletions
160
src/libs/ephemeris/celestialbody.cpp
Normal file
160
src/libs/ephemeris/celestialbody.cpp
Normal file
|
|
@ -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;
|
||||
}
|
||||
95
src/libs/ephemeris/celestialbody.h
Normal file
95
src/libs/ephemeris/celestialbody.h
Normal file
|
|
@ -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 <constants.h>
|
||||
|
||||
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_
|
||||
|
||||
Loading…
Reference in a new issue