]> git.mxchange.org Git - simgear.git/blobdiff - simgear/ephemeris/celestialBody.cxx
Add a new node "float-property" to be used in float comparision in effect predicates
[simgear.git] / simgear / ephemeris / celestialBody.cxx
index a46a34dfca3c416f5efcaf9f923c8fb405c6505f..4e8873925b65ef7290aba629afcf390648d4ac29 100644 (file)
@@ -5,28 +5,25 @@
  * September 1998. This code is based upon algorithms and data kindly 
  * provided by Mr. Paul Schlyter. (pausch@saaf.se). 
  *
- * 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.
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Library General Public
+ * License as published by the Free Software Foundation; either
+ * version 2 of the License, or (at your option) any later version.
  *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * General Public License for more details.
+ * Library General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program; if not, write to the Free Software
- * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/
 
 #include <simgear/debug/logstream.hxx>
 
-#ifdef FG_MATH_EXCEPTION_CLASH
-#  define exception c_exception
-#endif
 #include <math.h>
 
 #include "celestialBody.hxx"
@@ -34,7 +31,7 @@
 
 
 /**************************************************************************
- * void CelestialBody::updatePosition(fgTIME *t, Star *ourSun)
+ * void CelestialBody::updatePosition(double mjd, Star *ourSun)
  *
  * Basically, this member function provides a general interface for 
  * calculating the right ascension and declinaion. This function is 
  * position is calculated an a slightly different manner.  
  *
  * arguments:
- * fgTIME t: provides the current time.
+ * double mjd: provides the modified julian date.
  * Star *ourSun: the sun's position is needed to convert heliocentric 
  *               coordinates into geocentric coordinates.
  *
  * return value: none
  *
  *************************************************************************/
-void CelestialBody::updatePosition(FGTime *t, Star *ourSun)
+void CelestialBody::updatePosition(double mjd, Star *ourSun)
 {
   double eccAnom, v, ecl, actTime, 
     xv, yv, xh, yh, zh, xg, yg, zg, xe, ye, ze;
 
-  updateOrbElements(t);
-  actTime = fgCalcActTime(t);
+  updateOrbElements(mjd);
+  actTime = sgCalcActTime(mjd);
 
   // calcualate the angle bewteen ecliptic and equatorial coordinate system
-  ecl = DEG_TO_RAD * (23.4393 - 3.563E-7 *actTime);
+  ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 *actTime);
   
-  eccAnom = fgCalcEccAnom(M, e);  //calculate the eccentric anomaly
+  eccAnom = sgCalcEccAnom(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
@@ -87,8 +84,8 @@ void CelestialBody::updatePosition(FGTime *t, Star *ourSun)
   ze = yg * sin(ecl) + zg * cos(ecl);
   rightAscension = atan2(ye, xe);
   declination = atan2(ze, sqrt(xe*xe + ye*ye));
-  FG_LOG(FG_GENERAL, FG_INFO, "Planet found at : " 
-        << rightAscension << " (ra), " << declination << " (dec)" );
+  /* SG_LOG(SG_GENERAL, SG_INFO, "Planet found at : " 
+        << rightAscension << " (ra), " << declination << " (dec)" ); */
 
   //calculate some variables specific to calculating the magnitude 
   //of the planet
@@ -106,11 +103,11 @@ void CelestialBody::updatePosition(FGTime *t, Star *ourSun)
       tmp = -1.0;
   }
 
-  FV = RAD_TO_DEG * acos( tmp );
-};
+  FV = SGD_RADIANS_TO_DEGREES * acos( tmp );
+}
 
 /****************************************************************************
- * double CelestialBody::fgCalcEccAnom(double M, double e)
+ * double CelestialBody::sgCalcEccAnom(double M, double e)
  * this private member calculates the eccentric anomaly of a celestial body, 
  * given its mean anomaly and eccentricity.
  * 
@@ -135,7 +132,7 @@ void CelestialBody::updatePosition(FGTime *t, Star *ourSun)
  * the eccentric anomaly
  *
  ****************************************************************************/
-double CelestialBody::fgCalcEccAnom(double M, double e)
+double CelestialBody::sgCalcEccAnom(double M, double e)
 {
   double 
     eccAnom, E0, E1, diff;
@@ -151,28 +148,118 @@ double CelestialBody::fgCalcEccAnom(double M, double e)
          diff = fabs(E0 - E1);
          E0 = E1;
        }
-      while (diff > (DEG_TO_RAD * 0.001));
+      while (diff > (SGD_DEGREES_TO_RADIANS * 0.001));
       return E0;
     }
   return eccAnom;
 }
 
+/*****************************************************************************
+ * inline CelestialBody::CelestialBody
+ * public constructor for a generic celestialBody object.
+ * initializes the 6 primary orbital elements. The elements are:
+ * N: longitude of the ascending node
+ * i: inclination to the ecliptic
+ * w: argument of perihelion
+ * a: semi-major axis, or mean distance from the sun
+ * e: eccenticity
+ * M: mean anomaly
+ * Each orbital element consists of a constant part and a variable part that 
+ * gradually changes over time. 
+ *
+ * Argumetns:
+ * the 13 arguments to the constructor constitute the first, constant 
+ * ([NiwaeM]f) and the second variable ([NiwaeM]s) part of the orbital 
+ * elements. The 13th argument is the current time. Note that the inclination
+ * is written with a capital (If, Is), because 'if' is a reserved word in the 
+ * C/C++ programming language.
+ ***************************************************************************/ 
+CelestialBody::CelestialBody(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);
+}
 
+CelestialBody::CelestialBody(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;
+}
 
+/****************************************************************************
+ * inline void CelestialBody::updateOrbElements(double mjd)
+ * given the current time, this private member calculates the actual 
+ * orbital elements
+ *
+ * Arguments: double mjd: the current modified julian date:
+ *
+ * return value: none
+ ***************************************************************************/
+void CelestialBody::updateOrbElements(double mjd)
+{
+  double actTime = sgCalcActTime(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);
+}
 
+/*****************************************************************************
+ * inline double CelestialBody::sgCalcActTime(double mjd)
+ * this private member function returns the offset in days from the epoch for
+ * wich the orbital elements are calculated (Jan, 1st, 2000).
+ * 
+ * Argument: the current time
+ *
+ * return value: the (fractional) number of days until Jan 1, 2000.
+ ****************************************************************************/
+double CelestialBody::sgCalcActTime(double mjd)
+{
+  return (mjd - 36523.5);
+}
 
+/*****************************************************************************
+ * inline void CelestialBody::getPos(double* ra, double* dec)
+ * gives public access to Right Ascension and declination
+ *
+ ****************************************************************************/
+void CelestialBody::getPos(double* ra, double* dec)
+{
+  *ra  = rightAscension;
+  *dec = declination;
+}
 
-
-
-
-
-
-
-
-
-
-
-
-
+/*****************************************************************************
+ * inline void CelestialBody::getPos(double* ra, double* dec, double* magnitude
+ * gives public acces to the current Right ascension, declination, and 
+ * magnitude
+ ****************************************************************************/
+void CelestialBody::getPos(double* ra, double* dec, double* magn)
+{
+  *ra = rightAscension;
+  *dec = declination;
+  *magn = magnitude;
+}