simgear/ephemeris/star.cxx

   1 /**************************************************************************
   2  * star.cxx
   3  * Written by Durk Talsma. Originally started October 1997, for distribution
   4  * with the FlightGear project. Version 2 was written in August and
   5  * September 1998. This code is based upon algorithms and data kindly
   6  * provided by Mr. Paul Schlyter. (pausch@saaf.se).
   7  *
   8  * This library is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU Library General Public
  10  * License as published by the Free Software Foundation; either
  11  * version 2 of the License, or (at your option) any later version.
  12  *
  13  * This library is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * Library General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU Library General Public
  19  * License along with this library; if not, write to the
  20  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  21  * Boston, MA  02111-1307, USA.
  22  *
  23  * $Id$
  24  **************************************************************************/
  25
  26
  27 #ifdef __BORLANDC__
  28 #  define exception c_exception
  29 #endif
  30 #include <math.h>
  31
  32 #include <simgear/debug/logstream.hxx>
  33
  34 #include "star.hxx"
  35
  36
  37 /*************************************************************************
  38  * Star::Star(double mjd)
  39  * Public constructor for class Star
  40  * Argument: The current time.
  41  * the hard coded orbital elements our sun are passed to
  42  * CelestialBody::CelestialBody();
  43  * note that the word sun is avoided, in order to prevent some compilation
  44  * problems on sun systems
  45  ************************************************************************/
  46 Star::Star(double mjd) :
  47     CelestialBody (0.000000,  0.0000000000,
  48                    0.0000,    0.00000,
  49                    282.9404,  4.7093500E-5,
  50                    1.0000000, 0.000000,
  51                    0.016709,  -1.151E-9,
  52                    356.0470,  0.98560025850, mjd)
  53 {
  54     distance = 0.0;
  55 }
  56
  57 Star::Star() :
  58     CelestialBody (0.000000,  0.0000000000,
  59                    0.0000,    0.00000,
  60                    282.9404,  4.7093500E-5,
  61                    1.0000000, 0.000000,
  62                    0.016709,  -1.151E-9,
  63                    356.0470,  0.98560025850)
  64 {
  65     distance = 0.0;
  66 }
  67
  68 Star::~Star()
  69 {
  70 }
  71
  72
  73 /*************************************************************************
  74  * void Star::updatePosition(double mjd, Star *ourSun)
  75  *
  76  * calculates the current position of our sun.
  77  *************************************************************************/
  78 void Star::updatePosition(double mjd)
  79 {
  80   double
  81     actTime, eccAnom,
  82     xv, yv, v, r,
  83     xe, ecl;
  84
  85   updateOrbElements(mjd);
  86
  87   actTime = sgCalcActTime(mjd);
  88   ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 * actTime); // Angle in Radians
  89   eccAnom = sgCalcEccAnom(M, e);  // Calculate the eccentric Anomaly (also known as solving Kepler's equation)
  90
  91   xv = cos(eccAnom) - e;
  92   yv = sqrt (1.0 - e*e) * sin(eccAnom);
  93   v = atan2 (yv, xv);                   // the sun's true anomaly
  94   distance = r = sqrt (xv*xv + yv*yv);  // and its distance
  95
  96   lonEcl = v + w; // the sun's true longitude
  97   latEcl = 0;
  98
  99   // convert the sun's true longitude to ecliptic rectangular
 100   // geocentric coordinates (xs, ys)
 101   xs = r * cos (lonEcl);
 102   ys = r * sin (lonEcl);
 103
 104   // convert ecliptic coordinates to equatorial rectangular
 105   // geocentric coordinates
 106
 107   xe = xs;
 108   ye = ys * cos (ecl);
 109   ze = ys * sin (ecl);
 110
 111   // And finally, calculate right ascension and declination
 112   rightAscension = atan2 (ye, xe);
 113   declination = atan2 (ze, sqrt (xe*xe + ye*ye));
 114 }