SO3Engine/SO3Utils/SO3Astronomy.cpp

Go to the documentation of this file.

/*
-----------------------------------------------------------------------------
This source file is part of OpenSpace3D
For the latest info, see http://www.openspace3d.com

Copyright (c) 2010 I-maginer

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser 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 MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA, or go to
http://www.gnu.org/copyleft/lesser.txt

You may alternatively use this source under the terms of a specific version of
the OpenSpace3D Unrestricted License provided you have obtained such a license from
I-maginer.
-----------------------------------------------------------------------------
*/

/*
Taken from Caelum (http://www.ogre3d.org/wiki/index.php/Caelum)
Planetary position computation from http://stjarnhimlen.se/comp/ppcomp.html#5
Sunrise/sunset & moonrise/moonset computation from http://stjarnhimlen.se/comp/riset.html
There's also usefull informations in "A physically-based night sky model" concerning the moon position calculation.
*/

#include "SO3Astronomy.h"

namespace SO3
{

const double SAstronomy::PI     = 3.1415926535897932384626433832795029L;
const double SAstronomy::J2000  = 2451545.0;

double SAstronomy::RadToDeg(double value)
{
  return value * 180 / PI;
}

double SAstronomy::DegToRad(double value)
{
  return value * PI / 180;
}

double SAstronomy::SinDeg(double x)
{
  return std::sin(DegToRad(x));
}

double SAstronomy::CosDeg(double x)
{
  return std::cos(DegToRad(x));
}

double SAstronomy::Atan2Deg(double y, double x)
{
  return RadToDeg(std::atan2(y, x));
}

double SAstronomy::NormalizeDegrees(double value)
{
  value = fmod(value, 360);
  if (value < double(0))
  {
    value += double(360);
  }
  return value;
}

void SAstronomy::ConvertEclipticToEquatorialRad(double lon, double lat, double& rasc, double& decl)
{
  double ecl = SAstronomy::DegToRad(23.439281);
  
  double x = cos(lon) * cos(lat);
        double y = cos(ecl) * sin(lon) * cos(lat) - sin(ecl) * sin(lat);
        double z = sin(ecl) * sin(lon) * cos(lat) + cos(ecl) * sin(lat);

  double r = sqrt(x * x + y * y);
  rasc = atan2(y, x);
  decl = atan2(z, r);
} 

void SAstronomy::ConvertRectangularToSpherical(double x, double y, double z, double& rasc, double& decl, double& dist)
{
  dist = sqrt (x * x + y * y + z * z);
  rasc = Atan2Deg(y, x);
  decl = Atan2Deg(z, sqrt (x * x + y * y));
}

void SAstronomy::ConvertSphericalToRectangular(double rasc, double decl, double dist, double& x, double& y, double& z)
{
  x = dist * CosDeg(rasc) * CosDeg(decl);
  y = dist * SinDeg(rasc) * CosDeg(decl);
  z = dist * SinDeg(decl);
}

void SAstronomy::ConvertEquatorialToHorizontal(double jday, double longitude, double latitude, double rasc, double decl, double& azimuth, double& altitude)
{
  double d = jday - 2451543.5;
  double w = double(282.9404 + 4.70935E-5 * d);
  double M = double(356.0470 + 0.9856002585 * d);

  // Sun's mean longitude
  double L = w + M;

  // Universal time of day in degrees.
  double UT = double(fmod(d, 1) * 360);
  double hourAngle = longitude + L + double(180) + UT - rasc;

  double x = CosDeg(hourAngle) * CosDeg(decl);
  double y = SinDeg(hourAngle) * CosDeg(decl);
  double z = SinDeg(decl);

  double xhor = x * SinDeg(latitude) - z * CosDeg(latitude);
  double yhor = y;
  double zhor = x * CosDeg(latitude) + z * SinDeg(latitude);

  azimuth = Atan2Deg(yhor, xhor) + double(180);
  altitude = Atan2Deg (zhor, sqrt (xhor * xhor + yhor * yhor));
}

void SAstronomy::GetEclipticSunMeanLongitude(double jday, double& sunlon)
{
  // 2451544.5 == Astronomy::getJulianDayFromGregorianDateTime(2000, 1, 1, 0, 0, 0));
  // 2451543.5 == Astronomy::getJulianDayFromGregorianDateTime(1999, 12, 31, 0, 0, 0));
  double d = jday - 2451543.5;

  // Sun's Orbital elements:
  // argument of perihelion
  double w = double(282.9404 + 4.70935E-5 * d);
  // eccentricity (0=circle, 0-1=ellipse, 1=parabola)
  double e = 0.016709 - 1.151E-9 * d;
  // mean anomaly (0 at perihelion; increases uniformly with time)
  double M = double(356.0470 + 0.9856002585 * d);
  // Obliquity of the ecliptic.
  //double oblecl = double(23.4393 - 3.563E-7 * d);

  // Eccentric anomaly
  double E = M + RadToDeg(e * SinDeg(M) * (1 + e * CosDeg(M)));

  // Sun's Distance(R) and true longitude(L)
  double xv = CosDeg(E) - e;
  double yv = SinDeg(E) * sqrt(1 - e * e);
  sunlon = Atan2Deg(yv, xv) + w;
}

void SAstronomy::GetEclipticSunPosition(double jday, double& lambda, double& beta)
{
  double sunlon, sunlat = 0;
  GetEclipticSunMeanLongitude(jday, sunlon);

  lambda = DegToRad(sunlon);
        beta = DegToRad(sunlat);
}

void SAstronomy::GetEquatorialSunPosition(double jday, double& sunRightAscension, double& sunDeclinaison)
{
  double lambda, beta;
  GetEclipticSunPosition(jday, lambda, beta);
  ConvertEclipticToEquatorialRad(lambda, beta, sunRightAscension, sunDeclinaison);
        sunRightAscension = RadToDeg(sunRightAscension);
        sunDeclinaison = RadToDeg(sunDeclinaison);
}

void SAstronomy::GetHorizontalSunPosition(double jday, double longitude, double latitude, double& azimuth, double& altitude)
{
  // Horizontal spherical.
  double sunRightAscension, sunDeclinaison;
  GetEquatorialSunPosition(jday, sunRightAscension, sunDeclinaison);
  ConvertEquatorialToHorizontal(jday, longitude, latitude, sunRightAscension, sunDeclinaison, azimuth, altitude);
}

void SAstronomy::GetHorizontalSunPosition(double jday, Ogre::Degree longitude, Ogre::Degree latitude, Ogre::Degree& azimuth, Ogre::Degree& altitude)
{
  double az, al;
  GetHorizontalSunPosition(jday, longitude.valueDegrees(), latitude.valueDegrees(), az, al);
  azimuth = Ogre::Degree(static_cast<float>(az));
  altitude = Ogre::Degree(static_cast<float>(al));
}

void SAstronomy::GetEclipticMoonPositionRad(double jday, double& lon, double& lat)
{
  // Julian centuries since January 1, 2000
        double T = (jday - 2451545.0L) / 36525.0L; 
        double lprim = 3.8104L + 8399.7091L * T;
        double mprim = 2.3554L + 8328.6911L * T;
        double m = 6.2300L + 648.3019L * T;
        double d = 5.1985L + 7771.3772L * T;
        double f = 1.6280L + 8433.4663L * T;
        lon = lprim
              + 0.1098L * sin(mprim)
        + 0.0222L * sin(2.0L * d - mprim)
        + 0.0115L * sin(2.0L * d)
                                + 0.0037L * sin(2.0L * mprim)
        - 0.0032L * sin(m)
        - 0.0020L * sin(2.0L * f)
        + 0.0010L * sin(2.0L * d - 2.0L * mprim)
                                + 0.0010L * sin(2.0L * d - m - mprim)
        + 0.0009L * sin(2.0L * d + mprim)
        + 0.0008L * sin(2.0L * d - m)
                                + 0.0007L * sin(mprim - m)
        - 0.0006L * sin(d)
        - 0.0005L * sin(m + mprim);
                
  lat = + 0.0895L * sin(f)
        + 0.0049L * sin(mprim + f)
        + 0.0048L * sin(mprim - f)
        + 0.0030L * sin(2.0L * d - f)
        + 0.0010L * sin(2.0L * d + f - mprim)
        + 0.0008  * sin(2.0L * d - f - mprim)
        + 0.0006L * sin(2.0L * d + f);
}

void SAstronomy::GetHorizontalMoonPosition(double jday, double longitude, double latitude, double& azimuth, double& altitude)
{
  // Ecliptic spherical
  double lonecl, latecl;
  SAstronomy::GetEclipticMoonPositionRad(jday, lonecl, latecl);

  // Equatorial spherical
  double rasc, decl; 
  SAstronomy::ConvertEclipticToEquatorialRad(lonecl, latecl, rasc, decl);
  
  // Radians to degrees (all angles are in radians up to this point)
  rasc = RadToDeg(rasc);
  decl = RadToDeg(decl);

  // Equatorial to horizontal
  SAstronomy::ConvertEquatorialToHorizontal(jday, longitude, latitude, rasc, decl, azimuth, altitude);
}

void SAstronomy::GetHorizontalMoonPosition(double jday, Ogre::Degree longitude, Ogre::Degree latitude, Ogre::Degree& azimuth, Ogre::Degree& altitude)
{
  double az, al;
  GetHorizontalMoonPosition(jday, longitude.valueDegrees(), latitude.valueDegrees(), az, al);
  azimuth = Ogre::Degree(static_cast<float>(az));
  altitude = Ogre::Degree(static_cast<float>(al));
}

void SAstronomy::GetMoonPhase(double jday, float& moonPhase)
{
  // Calculates julian days since January 22, 2008 13:36 (full moon) and divides by the time between lunations (synodic month)
  ScopedHighPrecissionFloatSwitch precissionSwitch;
  double T = (jday - 2454488.0665L) / 29.531026L;

  T = fabs(fmod(T, 1));
  moonPhase = static_cast<float>(-fabs(-4 * T + 2) + 2);
}

int SAstronomy::GetJulianDayFromGregorianDate(int year, int month, int day)
{
  // Formulas from http://en.wikipedia.org/wiki/Julian_day
  // These are all integer divisions, but I'm not sure it works
  // correctly for negative values.
  int a = (14 - month) / 12;
  int y = year + 4800 - a;
  int m = month + 12 * a - 3;
  return day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045;
}

double SAstronomy::GetJulianDayFromGregorianDateTime(int year, int month, int day, int hour, int minute, double second)
{
  ScopedHighPrecissionFloatSwitch precissionSwitch;
  int jdn = GetJulianDayFromGregorianDate(year, month, day);
  
  // These are NOT integer divisions.
  double jd = jdn + (hour - 12) / 24.0 + minute / 1440.0 + second / 86400.0;
  return jd;
}

double SAstronomy::GetJulianDayFromGregorianDateTime(int year, int month, int day, double secondsFromMidnight)
{
  int jdn = GetJulianDayFromGregorianDate(year, month, day);
  double jd = jdn + secondsFromMidnight / 86400.0 - 0.5;
  return jd;
}

void SAstronomy::GetGregorianDateFromJulianDay(int julianDay, int& year, int& month, int& day)
{
  // From http://en.wikipedia.org/wiki/Julian_day
  int J = julianDay;
  int j = J + 32044;
  int g = j / 146097;
  int dg = j % 146097;
  int c = (dg / 36524 + 1) * 3 / 4;
  int dc = dg - c * 36524;
  int b = dc / 1461;
  int db = dc % 1461;
  int a = (db / 365 + 1) * 3 / 4;
  int da = db - a * 365;
  int y = g * 400 + c * 100 + b * 4 + a;
  int m = (da * 5 + 308) / 153 - 2;
  int d = da - (m + 4) * 153 / 5 + 122;
  year = y - 4800 + (m + 2) / 12;
  month = (m + 2) % 12 + 1;
  day = d + 1;
}

void SAstronomy::GetGregorianDateTimeFromJulianDay(double julianDay, int& year, int& month, int& day, int& hour, int& minute, double& second)
{
  // Integer julian days are at noon.
  // static_cast<int)(floor( is more precise than Ogre::Math::IFloor.
  // Yes, it does matter.
  int ijd = static_cast<int>(floor(julianDay + 0.5));
  GetGregorianDateFromJulianDay(ijd, year, month, day);

  double s = (julianDay + 0.5 - ijd) * 86400.0;
  hour = static_cast<int>(floor(s / 3600));
  s -= hour * 3600;
  minute = static_cast<int>(floor(s / 60));
  s -= minute * 60;
  second = s;
}

void SAstronomy::GetGregorianDateFromJulianDay(double julianDay, int& year, int& month, int& day)
{
  int hour;
  int minute;
  double second;
  GetGregorianDateTimeFromJulianDay(julianDay, year, month, day, hour, minute, second);
}

#if (OGRE_PLATFORM == OGRE_PLATFORM_WIN32) && (OGRE_COMPILER == OGRE_COMPILER_MSVC)
int SAstronomy::EnterHighPrecissionFloatingPointMode()
{
  int oldMode = ::_controlfp (0, 0);
  ::_controlfp (_PC_64, _MCW_PC);
  return oldMode;
}

void SAstronomy::RestoreFloatingPointMode(int oldMode)
{
  ::_controlfp (oldMode, _MCW_PC);
}
#else
int SAstronomy::EnterHighPrecissionFloatingPointMode()
{
  // Meaningless
  return 0xC0FFEE;
}

void SAstronomy::RestoreFloatingPointMode(int oldMode)
{
  // Useless check.
  assert(oldMode == 0xC0FFEE);
}
#endif

}

---