[prev in list] [next in list] [prev in thread] [next in thread]
List: kde-commits
Subject: [marble] src/plugins/render/satellites: [SoCiS 2013] Port satellites plugin to astro library
From: Rene Kuettner <rene () bitkanal ! net>
Date: 2013-12-18 16:02:39
Message-ID: E1VtJa3-0005cR-9f () scm ! kde ! org
[Download RAW message or body]
Git commit 9fb7c441ddc08eb9c78c545e048a6fe3c453a624 by Rene Kuettner, on behalf of \
Marek Hakala. Committed on 18/12/2013 at 15:57.
Pushed by renek into branch 'master'.
[SoCiS 2013] Port satellites plugin to astro library
This patch adapts the satellites plugin to the core astronomy
library (libastro). It also removes the legacy version of the
library from the plugin's directory.
REVIEW: 114268
M +2 -7 src/plugins/render/satellites/CMakeLists.txt
M +1 -1 src/plugins/render/satellites/SatellitesMSCItem.h
M +1 -1 src/plugins/render/satellites/SatellitesModel.cpp
D +0 -9 src/plugins/render/satellites/mex/README
D +0 -2278 src/plugins/render/satellites/mex/astrolib.cpp
D +0 -146 src/plugins/render/satellites/mex/astrolib.h
D +0 -669 src/plugins/render/satellites/mex/attlib.cpp
D +0 -109 src/plugins/render/satellites/mex/attlib.h
D +0 -727 src/plugins/render/satellites/mex/planetarySats.cpp
D +0 -104 src/plugins/render/satellites/mex/planetarySats.h
http://commits.kde.org/marble/9fb7c441ddc08eb9c78c545e048a6fe3c453a624
diff --git a/src/plugins/render/satellites/CMakeLists.txt \
b/src/plugins/render/satellites/CMakeLists.txt index 03650a4..b1b611e 100644
--- a/src/plugins/render/satellites/CMakeLists.txt
+++ b/src/plugins/render/satellites/CMakeLists.txt
@@ -1,6 +1,7 @@
PROJECT(SatellitesPlugin)
INCLUDE_DIRECTORIES(
+ ${CMAKE_SOURCE_DIR}/src/lib/astro
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_BINARY_DIR}
${QT_INCLUDE_DIR}
@@ -14,11 +15,6 @@ set( sgp4_SRCS
sgp4/sgp4unit.cpp
sgp4/sgp4io.cpp )
-set( mex_SRCS
- mex/astrolib.cpp
- mex/attlib.cpp
- mex/planetarySats.cpp )
-
set( satellites_SRCS
TrackerPluginModel.cpp
TrackerPluginItem.cpp
@@ -39,6 +35,5 @@ qt_add_resources( satellites_SRCS satellites.qrc )
marble_add_plugin( SatellitesPlugin
${satellites_SRCS}
- ${sgp4_SRCS}
- ${mex_SRCS} )
+ ${sgp4_SRCS} )
diff --git a/src/plugins/render/satellites/SatellitesMSCItem.h \
b/src/plugins/render/satellites/SatellitesMSCItem.h index f1fff67..e1eddfd 100644
--- a/src/plugins/render/satellites/SatellitesMSCItem.h
+++ b/src/plugins/render/satellites/SatellitesMSCItem.h
@@ -16,7 +16,7 @@
#include <QString>
#include <QDateTime>
-#include "mex/planetarySats.h"
+#include <planetarySats.h>
class QColor;
diff --git a/src/plugins/render/satellites/SatellitesModel.cpp \
b/src/plugins/render/satellites/SatellitesModel.cpp index 5d73653..91ced5d 100644
--- a/src/plugins/render/satellites/SatellitesModel.cpp
+++ b/src/plugins/render/satellites/SatellitesModel.cpp
@@ -19,7 +19,7 @@
#include "GeoDataStyle.h"
#include "sgp4/sgp4io.h"
-#include "mex/planetarySats.h"
+#include <planetarySats.h>
#include <locale.h>
diff --git a/src/plugins/render/satellites/mex/README \
b/src/plugins/render/satellites/mex/README deleted file mode 100644
index 2903c53..0000000
--- a/src/plugins/render/satellites/mex/README
+++ /dev/null
@@ -1,9 +0,0 @@
-
-This code allows to calculate spacecraft orbits/positions around other
-planets and has been contributed by Gerhard Holtkamp. See license
-headers in the source code files for detailed license information.
-
- astrolib - Standard astronomy routines
- attlib - Routines to handle matrix and vector operations in 3D
-planetarySats - Spacecraft orbit/position calculation routines
-
diff --git a/src/plugins/render/satellites/mex/astrolib.cpp \
b/src/plugins/render/satellites/mex/astrolib.cpp deleted file mode 100644
index a1d5cbf..0000000
--- a/src/plugins/render/satellites/mex/astrolib.cpp
+++ /dev/null
@@ -1,2278 +0,0 @@
-//
-// This file is part of the Marble Virtual Globe.
-//
-// This program is free software licensed under the GNU LGPL. You can
-// find a copy of this license in LICENSE.txt in the top directory of
-// the source code.
-//
-// Copyright 2012 Gerhard HOLTKAMP
-//
-
-/* =========================================================================
- Procedures needed for standard astronomy programs.
- The majority of procedures in this unit are taken from Montenbruck,
- Pfleger "Astronomie mit dem Personal Computer", Springer Verlag, 1989
- as well as from the "Explanatory Supplement to the Astronomical Almanac"
- University Science Books, Mill Valley, California, 1992
- and modified correspondingly.
-
- Linux C++ version
- Copyright : Gerhard HOLTKAMP 11-JAN-2012
- ========================================================================= */
-
-#include <cmath>
-using namespace std;
-
-#include "attlib.h"
-#include "astrolib.h"
-
-double frac (double f)
- { return fmod(f,1.0); }
-
-/*--------------- Function ddd ------------------------------------------*/
-
-double ddd (int d, int m, double s)
- {
- /*
- Conversion from Degrees, Minutes, Seconds into decimal degrees
-
- INPUT :
- d : degrees
- m : minutes
- s : seconds (and fractions)
-
- OUTPUT :
- dd : decimal degrees
- */
- double dd;
-
- if ( (d < 0) || (m < 0) || (s < 0)) dd = -1.0;
- else dd = 1.0;
- dd = dd * (fabs(double(d)) + fabs(double(m))/60.0 + fabs(s)/3600.0);
-
- return dd;
- }
-
-/*--------------- Function dms ------------------------------------------*/
-
-void dms (double dd, int &d, int &m, double &s)
- /*
- Conversion from decimal degrees into Degrees, Minutes, Seconds
-
- INPUT :
- dd : decimal degrees
-
- OUTPUT :
- d : degrees
- m : minutes
- s : seconds (and fractions)
- */
- {
- double d1;
-
- d1 = fabs(dd);
- // d = int(floor(d1));
- d = int(d1);
- d1 = (d1 - d) * 60.0;
- // m = int(floor(d1));
- m = int(d1);
- s = (d1 - m) * 60.0;
- if (dd < 0)
- {
- if (d != 0) d = -d;
- else
- {
- if (m != 0) m = -m;
- else s = -s;
- };
- };
- }
-
- /*------------ FUNCTION mjd ----------------------------------------------*/
-
-double mjd (int day, int month, int year, double hour)
-/*
- Modified Julian Date ( MJD = Julian Date - 2400000.5)
- valid for every date
- Julian Calendar up to 4-OCT-1582,
- Gregorian Calendar from 15-OCT-1582.
- */
- {
- double a;
- long int b, c;
-
- a = 10000.0 * year + 100.0 * month + day;
- if (month <= 2)
- {
- month = month + 12;
- year = year - 1;
- };
- if (a <= 15821004.1)
- {
- b = ((year+4716)/4) - 1181;
- if (year < -4716)
- {
- c = year + 4717;
- c = -c;
- c = c / 4;
- c = -c;
- b = c - 1182;
- };
- }
- else b = (year/400) - (year/100) + (year/4);
- // { b = -2 + floor((year+4716)/4) - 1179;}
- // else {b = floor(year/400) - floor(year/100) + floor(year/4);};
-
- a = 365.0 * year - 679004.0;
- a = a + b + int(30.6001 * (month + 1)) + day + hour / 24.0;
-
- return a;
- }
-
-/*---------------- Function julcent ---------------------------------------*/
-
-double julcent (double mjuld)
- {
- /*
- Julian Centuries since J2000.0 of Modified Julian Date mjuld.
- */
- return (mjuld - 51544.5) / 36525.0;
- }
-
-/*---------------- Function caldat -----------------------------------------*/
-
-void caldat (double mjd, int &day, int &month, int &year, double &hour)
- /*
- Calculate the calendar date from the Modified Julian Date
-
- INPUT :
- mjd : Modified Julian Date (Julian Date - 2400000.5)
-
- OUTPUT :
- day, month, year : corresponding date
- hour : corresponding hours of the above date
- */
- {
- long int b, c, d, e, f, jd0;
-
- jd0 = long(mjd + 2400001.0);
- if (jd0 < 2299161) c = jd0 + 1524; /* Julian calendar */
- else
- { /* Gregorian calendar */
- b = long (( jd0 - 1867216.25) / 36524.25);
- c = jd0 + b - (b/4) + 1525;
- };
-
- if (mjd < -2400001.0) // special case for year < -4712
- {
- if (mjd == floor(mjd)) jd0 = jd0 + 1;
- c = long((-jd0 - 0.1)/ 365.25);
- c = c + 1;
- year = -4712 - c;
- d = c / 4;
- d = c * 365 + d; // number of days from JAN 1, -4712
- f = d + jd0; // day of the year
- if ((c % 4) == 0) e = 61;
- else e = 60;
- if (f == 0)
- {
- year = year - 1;
- month = 12;
- day = 31;
- f = 500; // set as a marker
- };
- if (f < e)
- {
- if (f < 32)
- {
- month = 1;
- day = f;
- }
- else
- {
- month = 2;
- day = f - 31;
- };
- }
- else
- {
- if (f < 500)
- {
- f = f - e;
- month = long((f + 123.0) / 30.6001);
- day = f - long(month * 30.6001) + 123;
- month = month - 1;
- };
- };
- }
- else // normal case
- {
- d = long ((c - 122.1) / 365.25);
- e = 365 * d + (d/4);
- f = long ((c - e) / 30.6001);
- day = c - e - long(30.6001 * f);
- month = f - 1 - 12 * (f / 14);
- year = d - 4715 - ((7 + month) / 10);
- };
-
- hour = 24.0 * (mjd - floor(mjd));
- }
-
-/*---------------- Function DefTdUt --------------------------------------*/
-double DefTdUt (int yi)
- {
- /*
- Get a suitable default for value of TDT - UT in year yi in seconds.
- */
-
- const int td[9] = {55,55,56,56,57,58,58,59,60};
- double t, result;
- int yy, yr;
-
- yr = yi;
-
- if (yr > 1899)
- {
- if (yr >= 2000) // this corrects for observations until DEC 2006
- {
- if (yr > 2006) yr -= 1;
- if (yr > 2006) yr -= 1; // no leap second at the end of 2007
- if (yr > 2007) yr -= 1; // no leap second at the end of 2009
- if (yr > 2007) yr -= 1; // no leap second at the end of 2010
- yr -= 6;
- if (yr < 1999) yr = 1999;
- };
- t = yr - 2000;
- t = t / 100.0;
- result = (((-339.84*t - 516.12)*t -160.22)*t + 92.23)*t + 71.28;
- if (yr > 2013)
- {
- t = yr - 2013;
- t = t / 100.0;
- result = (27.5*t + 75.0)*t + 73.0;
- }
- else if (yr > 1994)
- {
- result -= 6.28;
- }
- else if (yr > 1985)
- {
- yy = yr - 1986;
- result = td[yy];
- };
- };
-
- if (yr < 1900)
- {
- t = yr - 1800;
- t = t / 100;
- if (yr > 1649)
- {
- t = t - 0.19;
- result = 5.156 + 13.3066 * t * t;
- if (yr > 1864) result = -0.6*(yr - 1865) + 6;
- if (yr > 1884) result = -0.2*(yr - 1885) - 6;
- }
- else
- {
- if (yr > 947) result = 25.5 * t*t;
- else result = 1360 + (44.3*t + 320.0) * t;
- };
- };
-
-// round to full seconds
-
- if (result < 0)
- {
- t = -result;
- t += 0.5;
- result = - floor(t);
- }
- else result = floor(result + 0.5);
-
- result += 0.184; // because of TT = TAI + 32.184;
-
- return result;
- }
-
-/*---------------- Function lsidtim -------------------------------------*/
-
-double lsidtim (double jd, double lambda, double ep2)
- {
- /* Calculate the Apparent Local Mean Sidereal Time (in decimal hours) for
- the Modified Julian Date jd and geographic longitude lambda (in degrees)
- and with the correction (AST - MST) eps (given in seconds)
- */
- double lmst;
- double t, ut, gmst;
- int mjd0;
-
- mjd0 = int(jd);
- ut = (jd - mjd0)*24.0;
- t = (mjd0 - 51544.5) / 36525.0;
- gmst = 6.697374558 + 1.0027379093*ut
- + (8640184.812866 + (0.093104 - 6.2e-6*t)*t)*t/3600.0;
- lmst = 24.0 * frac((gmst + lambda/15.0) / 24.0);
- lmst = lmst + ep2 / 3600.0;
-
- return lmst;
- }
-
-/*---------------- Function eps -----------------------------------------*/
-
-double eps (double t) // obliquity of ecliptic
- {
- /*
- Obliquety of ecliptic eps (in radians)
- at time t (in Julian Centuries from J2000.0)
- */
- double tp;
-
- tp = 23.43929111 - (46.815+(0.00059-0.001813*t)*t)*t/3600.0;
- tp = 1.74532925199e-2 * tp;
-
- return tp;
- }
-
-/*---------------- Function eclequ -----------------------------------------*/
-
-Vec3 eclequ (double t, Vec3& r1) // ecliptic -> equatorial
- {
- /*
- Convert position vector r1 from ecliptic into equatorial coordinates
- at t (in Julian Centuries from J2000.0 )
- */
-
- Mat3 m;
- Vec3 r2;
-
- m = xrot (-eps(t));
- r2 = mxvct (m, r1);
- return r2;
- }
-
-/*---------------- Function equecl -----------------------------------------*/
-
-Vec3 equecl (double t, Vec3& r1) // equatorial -> ecliptic
- {
- /*
- Convert position vector r1 from equatorial into ecliptic coordinates
- at t (in Julian Centuries from J2000.0)
- */
-
- Mat3 m;
- Vec3 r2;
-
- m = xrot (eps(t));
- r2 = mxvct (m, r1);
- return r2;
- }
-
-/*---------------- Function pmatecl -----------------------------------------*/
-
-Mat3 pmatecl (double t1, double t2) // ecl. precession
- {
- /*
- Calculate ecliptic precession matrix a from t1 to t2
- (times in Julian Centuries from J2000.0)
- */
-
- const double secrad = 4.8481368111e-6; // arcsec -> radians
-
- double ppi, pii, pa, dt;
- Mat3 m1, m2, a;
-
- dt = t2 - t1;
- ppi = 174.876383889 + (((3289.4789+0.60622*t1)*t1) +
- ((-869.8089-0.50491*t1) + 0.03536*dt)*dt) / 3600.0;
- ppi = ppi * 1.74532925199e-2;
- pii = ((47.0029-(0.06603-0.000598*t1)*t1) +
- ((-0.03302+0.000598*t1)+0.00006*dt)*dt)*dt * secrad;
- pa = ((5029.0966+(2.22226-0.000042*t1)*t1) +
- ((1.11113-0.000042*t1)-0.000006*dt)*dt)*dt * secrad;
-
- pa = ppi + pa;
- m1 = zrot (-pa);
- m2 = xrot (pii);
- m1 = m1 * m2;
- m2 = zrot (ppi);
- a = m1 * m2;
-
- return a;
- }
-
-/*---------------- Function pmatequ --------------------------------------*/
-
-Mat3 pmatequ (double t1, double t2) // equ. precession
- {
- /*
- Calculate equatorial precession matrix a from t1 to t2
- (times in Julian Centuries from J2000.0)
- */
- const double secrad = 4.8481368111e-6; // arcsec -> radians
-
- double dt, zeta, z, theta;
- Mat3 m1, m2, a;
-
- dt = t2 - t1;
- zeta = ((2306.2181+(1.39656-0.000139*t1)*t1) +
- ((0.30188-0.000345*t1)+0.017998*dt)*dt)*dt * secrad;
- z = zeta + ((0.7928+0.000411*t1)+0.000205*dt)*dt*dt * secrad;
- theta = ((2004.3109-(0.8533+0.000217*t1)*t1) -
- ((0.42665+0.000217*t1)+0.041833*dt)*dt)*dt * secrad;
-
- m1 = zrot (-z);
- m2 = yrot (theta);
- m1 = m1 * m2;
- m2 = zrot (-zeta);
- a = m1 * m2;
-
- return a;
- }
-
-/*---------------- Function nutmat --------------------------------------*/
-
-Mat3 nutmat (double t, double& ep2, bool hpr)
- {
- /*
- Calculate nutation matrix a from mean to true equatorial coordinates
- at time t (in Julian Centuries from J2000.0)
- Also calculates the correction ep2 for apparent sidereal time in sec
-
- if hpr is true a high precision is used, otherwise a low precision
- (only the first 50 terms of the nutation theory are used)
- */
- const int ntb1 = 15;
- const int tb1[ntb1][5] =
- {
- { 0, 0, 0, 0, 1}, // 1
- { 0, 0, 0, 0, 2}, // 2
- { 0, 0, 2,-2, 2}, // 9
- { 0, 1, 0, 0, 0}, // 10
- { 0, 1, 2,-2, 2}, // 11
- { 0,-1, 2,-2, 2}, // 12
- { 0, 0, 2,-2, 1}, // 13
- { 0, 2, 0, 0, 0}, // 16
- { 0, 2, 2,-2, 2}, // 18
- { 0, 0, 2, 0, 2}, // 31
- { 1, 0, 0, 0, 0}, // 32
- { 0, 0, 2, 0, 1}, // 33
- { 1, 0, 2, 0, 2}, // 34
- { 1, 0, 0, 0, 1}, // 38
- { -1, 0, 0, 0, 1}, // 39
- };
-
- const int ntb2 = 35;
- const int tb2[ntb2][5] =
- {
- { -2, 0, 2, 0, 1}, // 3
- { 2, 0,-2, 0, 0}, // 4
- { -2, 0, 2, 0, 2}, // 5
- { 1,-1, 0,-1, 0}, // 6
- { 0,-2, 2,-2, 1}, // 7
- { 2, 0,-2, 0, 1}, // 8
- { 2, 0, 0,-2, 0}, // 14
- { 0, 0, 2,-2, 0}, // 15
- { 0, 1, 0, 0, 1}, // 17
- { 0,-1, 0, 0, 1}, // 19
- { -2, 0, 0, 2, 1}, // 20
- { 0,-1, 2,-2, 1}, // 21
- { 2, 0, 0,-2, 1}, // 22
- { 0, 1, 2,-2, 1}, // 23
- { 1, 0, 0,-1, 0}, // 24
- { 2, 1, 0,-2, 0}, // 25
- { 0, 0,-2, 2, 1}, // 26
- { 0, 1,-2, 2, 0}, // 27
- { 0, 1, 0, 0, 2}, // 28
- { -1, 0, 0, 1, 1}, // 29
- { 0, 1, 2,-2, 0}, // 30
- { 1, 0, 0,-2, 0}, // 35
- { -1, 0, 2, 0, 2}, // 36
- { 0, 0, 0, 2, 0}, // 37
- { -1, 0, 2, 2, 2}, // 40
- { 1, 0, 2, 0, 1}, // 41
- { 0, 0, 2, 2, 2}, // 42
- { 2, 0, 0, 0, 0}, // 43
- { 1, 0, 2,-2, 2}, // 44
- { 2, 0, 2, 0, 2}, // 45
- { 0, 0, 2, 0, 0}, // 46
- { -1, 0, 2, 0, 1}, // 47
- { -1, 0, 0, 2, 1}, // 48
- { 1, 0, 0,-2, 1}, // 49
- { -1, 0, 2, 2, 1}, // 50
- };
-
- const double tb3[ntb1][4] =
- {
- {-171996.0,-174.2, 92025.0, 8.9 }, // 1
- { 2062.0, 0.2, -895.0, 0.5 }, // 2
- { -13187.0, -1.6, 5736.0, -3.1 }, // 9
- { 1426.0, -3.4, 54.0, -0.1 }, // 10
- { -517.0, 1.2, 224.0, -0.6 }, // 11
- { 217.0, -0.5, -95.0, 0.3 }, // 12
- { 129.0, 0.1, -70.0, 0.0 }, // 13
- { 17.0, -0.1, 0.0, 0.0 }, // 16
- { -16.0, 0.1, 7.0, 0.0 }, // 18
- { -2274.0, -0.2, 977.0, -0.5 }, // 31
- { 712.0, 0.1, -7.0, 0.0 }, // 32
- { -386.0, -0.4, 200.0, 0.0 }, // 33
- { -301.0, 0.0, 129.0, -0.1 }, // 34
- { 63.0, 0.1, -33.0, 0.0 }, // 38
- { -58.0, -0.1, 32.0, 0.0 }, // 39
- };
-
- const double tb4[ntb2][2] =
- {
- { 46.0, -24.0}, // 3
- { 11.0, 0.0 }, // 4
- { -3.0, 1.0 }, // 5
- { -3.0, 0.0 }, // 6
- { -2.0, 1.0 }, // 7
- { 1.0, 0.0 }, // 8
- { 48.0, 1.0 }, // 14
- {-22.0, 0.0 }, // 15
- {-15.0, 9.0 }, // 17
- {-12.0, 6.0 }, // 19
- { -6.0, 3.0 }, // 20
- { -5.0, 3.0 }, // 21
- { 4.0, -2.0 }, // 22
- { 4.0, -2.0 }, // 23
- { -4.0, 0.0 }, // 24
- { 1.0, 0.0 }, // 25
- { 1.0, 0.0 }, // 26
- { -1.0, 0.0 }, // 27
- { 1.0, 0.0 }, // 28
- { 1.0, 0.0 }, // 29
- { -1.0, 0.0 }, // 30
- {-158.0,-1.0 }, // 35
- {123.0,-53.0 }, // 36
- { 63.0, -2.0 }, // 37
- {-59.0, 26.0 }, // 40
- {-51.0, 27.0 }, // 41
- {-38.0, 16.0 }, // 42
- { 29.0, -1.0 }, // 43
- { 29.0,-12.0 }, // 44
- {-31.0, 13.0 }, // 45
- { 26.0, -1.0 }, // 46
- { 21.0,-10.0 }, // 47
- { 16.0, -8.0 }, // 48
- {-13.0, 7.0 }, // 49
- {-10.0, 5.0 }, // 50
- };
-
- const double secrad = 4.8481368111e-6; // arcsec -> radians
- const double p2 = 2.0 * M_PI;
-
- double ls, lm, d, f, n, dpsi, deps, ep0;
- int j;
- Mat3 m1, m2, a;
-
- if (hpr)
- {
- lm = 2.355548393544 +
- (8328.691422883903 + (0.000151795164 + 0.000000310281*t)*t)*t;
- ls = 6.240035939326 +
- (628.301956024185 + (-0.000002797375 - 0.000000058178*t)*t)*t;
- f = 1.627901933972 +
- (8433.466158318464 + (-0.000064271750 + 0.000000053330*t)*t)*t;
- d = 5.198469513580 +
- (7771.377146170650 + (-0.000033408511 + 0.000000092115*t)*t)*t;
- n = 2.182438624361 +
- (-33.757045933754 + (0.000036142860 + 0.000000038785*t)*t)*t;
-
- lm = fmod(lm,p2);
- ls = fmod(ls,p2);
- f = fmod(f,p2);
- d = fmod(d,p2);
- n = fmod(n,p2);
-
- dpsi = 0.0;
- deps = 0.0;
- for(j=0; j<ntb1; ++j)
- {
- ep0 = tb1[j][0]*lm + tb1[j][1]*ls + tb1[j][2]*f + tb1[j][3]*d + tb1[j][4]*n;
- dpsi = dpsi + (tb3[j][0]+tb3[j][1]*t) * sin(ep0);
- deps = deps + (tb3[j][2]+tb3[j][3]*t) * cos(ep0);
- };
- for(j=0; j<ntb2; ++j)
- {
- ep0 = tb2[j][0]*lm + tb2[j][1]*ls + tb2[j][2]*f + tb2[j][3]*d + tb2[j][4]*n;
- dpsi = dpsi + tb4[j][0] * sin(ep0);
- deps = deps + tb4[j][1] * cos(ep0);
- };
- dpsi = 1.0e-4 * dpsi * secrad;
- deps = 1.0e-4 * deps * secrad;
- }
-
- else // low precision
- {
- ls = p2 * frac (0.993133+99.997306*t); // mean anomaly sun
- d = p2 * frac (0.827362+1236.853087*t); // diff long. moon-sun
- f = p2 * frac (0.259089+1342.227826*t); // dist. node
- n = p2 * frac (0.347346 - 5.372447*t); // long. node
-
- dpsi = (-17.2*sin(n) - 1.319*sin(2*(f-d+n)) - 0.227*sin(2*(f+n))
- + 0.206*sin(2*n) + 0.143*sin(ls)) * secrad;
- deps = (+9.203*cos(n) + 0.574*cos(2*(f-d+n)) + 0.098*cos(2*(f+n))
- -0.09*cos(2*n) ) * secrad;
- };
-
- ep0 = eps (t);
- ep2 = ep0 + deps;
- m1 = xrot (ep0);
- m2 = zrot (-dpsi);
- m1 *= m2;
- m2 = xrot (-ep2);
- a = m2 * m1;
- ep2 = dpsi * cos (ep2);
- ep2 *= 13750.9870831; // convert radians into time-seconds
-
- return a;
- }
-
-/*---------------- Function nutecl --------------------------------------*/
-
-Mat3 nutecl (double t, double& ep2) // nutation matrix (ecliptic)
- {
- /*
- Calculate nutation matrix a from mean to true ecliptic coordinates
- at time t (in Julian Centuries from J2000.0)
- Also calculates the correction ep2 for apparent sidereal time in sec
- */
-
- const double secrad = 4.8481368111e-6; // arcsec -> radians
- const double p2 = 2.0 * M_PI;
-
- double ls, d, f, n, dpsi, deps, ep0;
- Mat3 m1, m2, a;
-
- ls = p2 * frac (0.993133+99.997306*t); // mean anomaly sun
- d = p2 * frac (0.827362+1236.853087*t); // diff long. moon-sun
- f = p2 * frac (0.259089+1342.227826*t); // dist. node
- n = p2 * frac (0.347346 - 5.372447*t); // long. node
-
- dpsi = (-17.2*sin(n) - 1.319*sin(2*(f-d+n)) - 0.227*sin(2*(f+n))
- + 0.206*sin(2*n) + 0.143*sin(ls)) * secrad;
-
- deps = (+9.203*cos(n) + 0.574*cos(2*(f-d+n)) + 0.098*cos(2*(f+n))
- -0.09*cos(2*n) ) * secrad;
-
- ep0 = eps (t);
- ep2 = ep0 + deps;
- m1 = xrot (-deps);
- m2 = zrot (-dpsi);
- a = m1 * m2;
- ep2 = dpsi * cos (ep2);
- ep2 *= 13750.9870831; // convert radians into time-seconds
-
- return a;
- }
-
-/*---------------- Function pmatequ --------------------------------------*/
-
-Mat3 PoleMx (double xp, double yp)
- {
- /* Returns Polar Motion matrix.
- xp and yp are the coordinates of the Celestial Ephemeris Pole
- with respect to the IERS Reference Pole in arcsec as published
- in the IERS Bulletin B.
- */
- const double arctrd = M_PI / (180.0*3600.0);
- Mat3 res;
- double xr, yr;
-
- xr = xp * arctrd;
- yr = yp * arctrd;
- res.assign(1.0,0.0,xr,0.0,1.0,-yr,-xr,yr,1.0);
-
- return res;
- }
-
-/*---------------- Function aberrat --------------------------------------*/
-
-Vec3 aberrat (double t, Vec3& ve) // aberration
- {
- /*
- Correct position vector ve for aberration into new position va
- at time t (in Julian Centuries from J2000.0)
- */
- const double p2 = 2.0 * M_PI;
-
- double l, cl, d0;
- Vec3 va;
-
- d0 = abs(ve);
- l = p2 * frac(0.27908+100.00214*t);
- cl = cos(l)*d0;
- va[0] = ve[0] - 9.934e-5 * sin(l) * d0;
- va[1] = ve[1] + 9.125e-5 * cl;
- va[2] = ve[2] + 3.927e-5 * cl;
-
- return va;
- }
-
-/*--------------------- Function GeoPos ----------------------------------*/
-
-Vec3 GeoPos (double jd, double ep2, double lat, double lng, double ht)
- {
- /* Return the geocentric vector (in the equatorial system) of the
- geographic position given by the latitude lat and longitude lng
- (in radians) and height ht (in m) at the MJD-time jd (UT).
- ep2 : correction for apparent sidereal time in sec
-
- The length unit of the vector is in terms of the equatorial
- Earth radius (6378.137 km)
- */
- double const e2 = 6.69438499959e-3;
- double np, h, sp;
-
- Vec3 r;
-
- sp = sin(lat);
- h = ht / 6378.137e3;
- np = 1.0 / (sqrt(1.0 - e2*sp*sp));
- r[2] = ((1.0 - e2)*np + h)*sp;
-
- sp = (np + h) * cos(lat);
- np = lsidtim(jd, (lng*180.0/M_PI), ep2) * M_PI/12.0;
-
- r[0] = sp * cos(np);
- r[1] = sp * sin(np);
-
- return r;
- }
-
-Vec3 GeoPos (double jd, double ep2, double lat, double lng, double ht,
- double xp, double yp)
- {
- /* Return the geocentric vector (in the equatorial system) of the
- geographic position given by the latitude lat and longitude lng
- (in radians) and height ht (in m) at the MJD-time jd (UT).
-
- ep2 : correction for apparent sidereal time in sec
- xp, yp: coordinates of polar motion in arcsec.
-
- The length unit of the vector is in terms of the equatorial
- Earth radius (6378.137 km)
- */
- double const e2 = 6.69438499959e-3;
- double np, h, sp;
- Mat3 prx;
- Vec3 r;
-
- sp = sin(lat);
- h = ht / 6378.137e3;
- np = 1.0 / (sqrt(1.0 - e2*sp*sp));
- r[2] = ((1.0 - e2)*np + h)*sp;
- sp = (np + h) * cos(lat);
- r[0] = sp * cos(lng);
- r[1] = sp * sin(lng);
-
- // correct for polar motion
- if ((xp != 0) || (yp != 0))
- {
- prx = mxtrn(PoleMx(xp, yp));
- r = mxvct(prx, r);
- };
-
- // convert into equatorial
- np = lsidtim(jd, 0.0, ep2) * M_PI/12.0;
- prx = zrot(-np);
- r = mxvct(prx, r);
-
- return r;
- }
-
-/*--------------------- Function EquHor ----------------------------------*/
-
-Vec3 EquHor (double jd, double ep2, double lat, double lng, Vec3 r)
- {
- /* convert vector r from the equatorial system into the horizontal system.
- jd = MJD-time (UT)
- ep2 : correction for apparent sidereal time in sec
- lat, lng : geographic latitude and longitude (in radians)
- */
- double lst;
- Vec3 s;
- Mat3 mx;
-
- lst = lsidtim(jd, (lng*180.0/M_PI), ep2) * M_PI/12.0;
- mx = zrot(lst);
- s = mxvct(mx, r);
- mx = yrot(M_PI/2.0 - lat);
- s = mxvct(mx, s);
-
- return s;
- }
-
-/*--------------------- Function HorEqu ----------------------------------*/
-
-
-Vec3 HorEqu (double jd, double ep2, double lat, double lng, Vec3 r)
- {
- /* convert vector r from the horizontal system into the equatorial system.
- jd = MJD-time (UT)
- ep2 : correction for apparent sidereal time in sec
- lat, lng : geographic latitude and longitude (in radians)
- */
- double lst;
- Vec3 s;
- Mat3 mx;
-
- mx = yrot(lat - M_PI/2.0);
- s = mxvct(mx, r);
- lst = -lsidtim(jd, (lng*180.0/M_PI), ep2) * M_PI/12.0;
- mx = zrot(lst);
- s = mxvct(mx, s);
-
- return s;
- }
-
-/*--------------------- Function AppPos ----------------------------------*/
-
-void AppPos (double jd, double ep2, double lat, double lng, double ht,
- int solsys, Vec3 r, double& azim, double& elev, double& dist)
- {
- /* get apparent position in the horizontal system
- jd = MJD-time (UT)
- ep2 : correction for apparent sidereal time in sec
- lat, lng : geographic latitude and longitude (in radians)
- ht : height above normal in meters.
- solsys : = 1 if object is in solar system and parallax has to be
- taken into account, 0 otherwise.
- r = vector of celestial object. The unit of length of this vector
- has to be in terms of the equatorial Earth radius (6378.14 km)
- if solsys = 1, otherwise it's arbitrary.
- azim : azimuth in radians (0 is to the North).
- elev : elevation in radians
- dist : distance (if solsys = 1; otherwise abs(r))
- */
- Vec3 s;
-
- // correct for topocentric position (parallax)
- if (solsys) s = r - GeoPos(jd, ep2, lat, lng, ht);
- else s = r;
-
- s = EquHor(jd, ep2, lat, lng, s);
- s = carpol(s);
- dist = s[0];
- elev = s[2];
- azim = M_PI - s[1];
- }
-
-/*--------------------- Function AppRADec ----------------------------------*/
-
-void AppRADec (double jd, double ep2, double lat, double lng,
- double azim, double elev, double& ra, double& dec)
- {
- /* get apparent position in the horizontal system
- jd = MJD-time (UT)
- ep2 : correction for apparent sidereal time in sec
- lat, lng : geographic latitude and longitude (in radians)
- azim : azimuth in radians (0 is to the North).
- elev : elevation in radians
- ra : Right Ascension (in radians)
- dec : Declination (in radians)
- */
- Vec3 s;
-
- s[0] = 1.0;
- s[1] = M_PI - azim;
- s[2] = elev;
- s = polcar(s);
- s = HorEqu(jd, ep2, lat, lng, s);
- s = carpol(s);
- dec = s[2];
- ra = s[1];
- }
-
-/*------------------------- Refract ---------------------------------*/
-
-double Refract (double h, double p, double t)
- {
- /* Calculate atmospheric refraction with low precision
- h = height (in radians) of object
- p = presure (in millibars)
- t = temperature (in degrees Celsius)
-
- RETURN: refraction angle in radians
- */
- double const raddeg = 180.0 / M_PI;
- double r;
-
- r = h * raddeg;
- r = (r + 7.31 / (r + 4.4)) / raddeg;
- r = (0.28*p/(t+273.0)) * 0.0167 / tan(r);
- r = r / raddeg;
-
- return r;
- }
-
-/*---------------- Function eccanom --------------------------------------*/
-
-double eccanom (double man, double ecc)
- {
- /*
- Solve Kepler equation for eccentric anomaly of elliptical orbits.
- man : Mean Anomaly (in radians)
- ecc : eccentricity
- */
- const double p2 = 2.0*M_PI;
- const double eps = 1E-11;
- const int maxit = 15;
-
- double m, e, f;
- int i, mi;
-
- m=man/p2;
- mi = int(m);
- m=p2*(m-mi);
- if (m < 0) m=m+p2;
- if (ecc < 0.8) e=m;
- else e=M_PI;
- f = e - ecc*sin(e) - m;
- i=0;
-
- while ((fabs(f) > eps) && (i < maxit))
- {
- e = e - f / (1.0 - ecc*cos(e));
- f = e - ecc*sin(e) - m;
- i = i + 1;
- }
-
- return e;
- }
-
-
-/*---------------- Function hypanom --------------------------------------*/
-
-double hypanom (double mh, double ecc)
- {
- /*
- Solve Kepler equation for eccentric anomaly of hyperbolic orbits.
- mh : Mean Anomaly
- ecc : eccentricity
- */
- const double eps = 1E-11;
- const int maxit = 15;
-
- double h, f;
- int i;
-
- h = log (2.0*fabs(mh)/ecc+1.8);
- if (mh < 0.0) h = -h;
- f = ecc * sinh(h) - h - mh;
- i = 0;
-
- while ((fabs(f) > eps*(1.0+fabs(h+mh))) && (i < maxit))
- {
- h = h - f / (ecc*cosh(h) - 1.0);
- f = ecc*sinh(h) - h - mh;
- i = i + 1;
- }
-
- return h;
- }
-
-
-/*------------------ Function ellip ------------------------------------*/
-
-void ellip (double gm, double t0, double t, double a, double ecc,
- double m0, Vec3& r1, Vec3& v1)
- {
- /*
- Calculate position r1 and velocity v1 of for elliptic orbit at time t.
- gm : Gravitational Constant
- t0 : epoch
- a : semim-ajor axis
- ecc : eccentricity
- m0 : Mean Anomaly at epoch (in radians).
- The units must be consistent with each other.
- */
- double m, e, fac, k, c, s;
-
- if (fabs(a) < 1e-60) a = 1e-60;
- k = gm / a;
- if (k >= 0) k = sqrt(k);
- else k = 0; // just in case
-
- m = k * (t - t0) / a + m0;
- e = eccanom(m, ecc);
- fac = sqrt(1.0 - ecc*ecc);
- c = cos(e);
- s = sin(e);
- r1.assign (a*(c - ecc), a*fac*s, 0.0);
- m = 1.0 - ecc*c;
- v1.assign (-k*s/m, k*fac*c/m, 0.0);
- }
-
-/*---------------- Function hyperb --------------------------------------*/
-
-void hyperb (double gm, double t0, double t, double a, double ecc,
- Vec3& r1, Vec3& v1)
- {
- /*
- Calculate position r1 and velocity v1 of for hyperbolic orbit at time t.
- gm : Gravitational Constant
- t0 : time of perihelion passage
- a : semi-major axis
- ecc : eccentricity
- The units must be consistent with each other.
- */
- double mh, h, fac, k, c, s;
-
- a = fabs(a);
- if (a < 1e-60) a = 1e-60;
- k = gm / a;
- if (k >= 0) k = sqrt(k);
- else k = 0; // just in case
-
- mh = k * (t - t0) / a;
- h = hypanom (mh, ecc);
- fac = sqrt(ecc*ecc-1.0);
- c = cosh(h);
- s = sinh(h);
- r1.assign (a*(ecc-c), a*fac*s, 0.0);
- mh = ecc*c - 1.0;
- v1.assign (-k*s/mh, k*fac*c/mh, 0.0);
- }
-
-/*---------------- Function parab --------------------------------------*/
-
- void stumpff (double e2, double& c1, double& c2, double& c3)
- {
- /*
- Calculation of Stumpff functions c1=sin(e)/c,
- c2=(1-cos(e))/e^2 and c3=(e-sin(e))/e^3 for the
- argument e2=e^2 (e: eccentric anomaly)
- */
- const double eps=1E-12;
- double n, add;
-
- c1=0.0; c2=0.0; c3=0.0; add=1.0; n=1.0;
- do
- {
- c1=c1+add; add=add/(2.0*n);
- c2=c2+add; add=add/(2.0*n+1);
- c3=c3+add; add=-e2*add;
- n=n+1.0;
- } while (abs(add)>eps);
- }
-
-void parab (double gm, double t0, double t, double q, double ecc,
- Vec3& r1, Vec3& v1)
- {
- /*
- Calculate position r1 and velocity v1 of for parabolic
- (or near parabolic) orbit at time t using Stumpff's method.
- gm : Gravitational Constant
- t0 : time of perihelion passage
- q : perihelion distance
- ecc : eccentricity
- The units must be consistent with each other.
- */
- const double eps = 1E-9;
- const int maxit = 15;
-
- double e2, e20, fac, c1, c2, c3, k, tau, a, u, u2;
- double x, y;
- int i, fle;
-
- ecc = fabs(ecc);
- e2=0.0; fac=0.5*ecc; i=0;
-
- q = fabs(q);
- if (q < 1e-40) q = 1e-40;
- k = gm / (q*(1+ecc));
-
- if (k >= 0) k = sqrt(k);
- else k = 0; // just in case
-
- tau = gm / (q*q*q);
- if (tau >= 0) tau= 1.5*sqrt(tau)*(t-t0);
- else tau = 0;
-
- do
- {
- i = i + 1;
- e20 = e2;
- if (fac < 0.0) a = -1.0;
- else a = tau * sqrt(fac);
- a = sqrt(a*a+1.0)+a;
- if (a > 0.0) a = exp(log(a)/3.0);
- if (a == 0.0) u = 0.0;
- else u = a - 1.0/a;
- u2 = u*u;
- if (fac != 0) e2 = u2*(1.0-ecc)/fac;
- else e2 = 1;
- stumpff (e2, c1, c2, c3);
- fac = 3.0*ecc*c3;
- fle = (abs(e2-e20) < eps) || (i > maxit);
- } while (!fle);
-
- if (fac != 0)
- {
- tau = q*(1.0+u2*c2*ecc/fac);
- x = q*(1.0-u2*c2/fac);
- y = (1.0+ecc)/fac;
- if (y >= 0) y = q*sqrt(y)*u*c1;
- else y = 0;
- r1.assign (x, y, 0.0);
- v1.assign (-k*y/tau, k*(x/tau+ecc), 0.0);
- }
- else // just in case
- {
- r1.assign (0, 0, 0);
- v1.assign (0, 0, 0);
- };
- }
-
-/*---------------- Function kepler --------------------------------------*/
-
-void kepler (double gm, double t0, double t, double m0, double a, double ecc,
- double ran, double aper, double inc, Vec3& r1, Vec3& v1)
- {
- /*
- Calculate position r1 and velocity v1 of body at time t as an
- undisturbed 2-body Kepler problem.
-
- gm : Gravitational Constant
- t0 : time of perihelion passage (hyperbolic or near-parabolic)
- epoch of m0 for elliptical orbits.
- m0 : Mean Anomaly at epoch for elliptical orbits in degrees in the
- range 0 to 360. If m0 < 0 the calculations will be done as
- near-parabolic. Set m0 = 0 for hyperbolic orbits.
- a : semi-major axis for elliptical (positive) or hyperbolic orbits
- (negative). If m0 < 0, a signifies perihelion distance q instead
- ecc : eccentricity
- ran : Right Ascension of Ascending Node (in degrees)
- aper : Argument of Perihelion (in degrees)
- inc : Inclination (in degrees)
-
- The units must be consistent with each other.
-
- NOTE : If ecc = 1 the orbit will always be calculated as a parabolic one.
- If ecc < 1 (elliptical orbits) two possibilities exist:
- If m0 >= 0 the calculation will be done in the classical
- way with m0 signifying the mean anomaly at time t0 and
- a the semi-major axis.
- If m0 < 0 the orbit will be calculated as a near-parabolic
- orbit with a signifying the perihelion distance at time t0.
- (This latter method will be useful for high eccentricity
- comet orbits for which typically q is given instead of a)
- If ecc > 1 (hyperbolic orbits) two possibilities exist:
- If m0 >= 0 (the value is ignored for hyperbolic orbits)
- a classical hyperbolic calculation is performed with a
- signifying the semi-major axis (negative for hyperbolas).
- If m0 < 0 the orbit will be calculated as a near-parabolic
- orbit with a signifying the perihelion distance. (This can
- be useful for comet orbits which would rarely have an
- eccentricity >> 1)
- In either case t0 signifies the time of perihelion passage.
- */
- const double dgrd = M_PI / 180.0;
- enum {ell, par, hyp} kepc;
- Mat3 m1,m2;
-
- kepc = ell; // just to keep the compiler happy
-
- // convert into radians
- m0 *= dgrd;
- ran *= dgrd;
- aper *= dgrd;
- inc *= dgrd;
-
- // calculate the position and velocity within the orbit plane
- if (ecc == 1.0) kepc = par;
- if (ecc < 1.0)
- {
- if (m0 < 0.0) kepc = par;
- else kepc = ell;
- }
- if (ecc > 1.0)
- {
- if (m0 < 0.0) kepc = par;
- else kepc = hyp;
- }
-
- switch (kepc)
- {
- case ell : ellip (gm, t0, t, a, ecc, m0, r1, v1); break;
- case par : parab (gm, t0, t, a, ecc, r1, v1); break;
- case hyp : hyperb (gm, t0, t, a, ecc, r1, v1); break;
- }
-
- // convert into reference plane
- m1 = zrot (-aper);
- m2 = xrot (-inc);
- m1 *= m2;
- m2 = zrot (-ran);
- m2 = m2 * m1;
-
- r1 = mxvct (m2, r1);
- v1 = mxvct (m2, v1);
- }
-
-/*---------------- Function oscelm --------------------------------------*/
-
-void oscelm (double gm, double t, Vec3& r1, Vec3& v1,
- double& t0, double& m0, double& a, double& ecc,
- double& ran, double& aper, double& inc)
- {
- /*
- Get the osculating Kepler elements at epoch t from position r1 and
- velocity v1 of body.
-
- gm : Gravitational Constant. If set negative, its absolute value will
- be taken as gm and the perihelion distance will be returned
- instead of the semimajor axis.
- t0 : time of perihelion passage
- m0 : Mean Anomaly at epoch for elliptical orbits in degrees in the
- range 0 to 360. Set m0 = 0 for hyperbolic orbits.
- m0 will be set to -1 if perihelion distance is to be returned
- instead of a.
- a : semi-major axis for elliptical (positive) or hyperbolic orbits
- (negative). If m0 < 0, a signifies perihelion distance q instead.
- If gm is negative, the perihelion distance q will be returned
- ecc : eccentricity
- ran : Right Ascension of Ascending Node (in degrees)
- aper : Argument of Perihelion (in degrees)
- inc : Inclination (in degrees)
-
- The units must be consistent with each other.
- */
- const double rddg = 180.0/M_PI;
- const double p2 = 2.0*M_PI;
-
- Vec3 c;
- double cabs, u, cu, su, p, r;
- int pflg; // 1, if perihelion distance to be returned
-
- pflg = 0;
- if (gm < 0.0)
- {
- gm = -gm;
- pflg = 1;
- };
-
- if (gm < 1e-60) gm = 1e-60;
- c = r1 * v1;
- cabs = abs(c);
- if (fabs(cabs) < 1e-40) cabs = 1e-40;
- ran = atan20 (c[0], -c[1]);
- inc = c[2]/cabs;
- if (fabs(inc) <= 1.0) inc = acos (inc);
- else inc = 0;
-
- r = abs(r1);
- if (fabs(r) < 1e-40) r = 1e-40;
- su = sin(inc);
- if (su != 0.0) su = r1[2]/su;
- cu = r1[0]*cos(ran)+r1[1]*sin(ran);
- u = atan20 (su, cu); // argument of latitude
-
- a = abs(v1);
- a = 2.0/r - a*a/gm;
- if (fabs(a) < 1.0E-30) ecc = 1.0; // parabolic
- else
- {
- a = 1.0/a; // semimajor axis
- ecc = 0;
- };
-
- p = cabs*cabs/gm; // semilatum rectum
-
- if (ecc == 1.0)
- {
- p = p / 2.0; // perihelion distance
- a = 2*p;
- }
- else
- {
- ecc = 1.0 - p/a;
- if (ecc >= 0) ecc = sqrt(ecc);
- else ecc = 0;
- p = p / (1.0 + ecc); // perihelion distance
- }
-
- if (fabs(a) > 1e-60) cu = (1.0 - r/a);
- else cu = 0; // just in case
- su = dot(r1,v1) / sqrt(fabs(a)*gm);
- if (ecc < 1.0)
- {
- m0 = atan20(su,cu); // eccentric anomaly
- su = sin(m0);
- cu = cos(m0);
- aper = 1.0-ecc*ecc;
- if (aper >= 0) aper = atan20(sqrt(aper)*su,(cu-ecc)); // true anomaly
- m0 = m0 - ecc*su; // mean anomaly
- }
- else if (ecc > 1.0)
- {
- su = su / ecc;
- m0 = su + sqrt(su*su + 1.0);
- if (m0 >= 0) m0 = log(m0); // = asinh (su); hyperbolic anomaly
- aper = (ecc+1.0)/(ecc-1.0);
- if (aper >= 0) aper = 2.0 * atan(sqrt(aper)*tanh(m0/2.0));
- m0 = ecc*su-m0;
- }
-
- if (ecc != 1.0)
- {
- aper = u - aper; // argument of perihelion
- u = fabs(a);
- t0 = u/gm;
- if (t0 >= 0) t0 = t - u*sqrt(t0)*m0; // time of perihelion transit
- else t0 = t;
- }
- else // parabolic
- {
- pflg = 1;
-
- aper = 2.0 * atan(su); // true anomaly
- aper = u - aper; // argument of perihelion
- t0 = 2.0*p*p*p/gm;
- if (t0 >= 0) t0 = t - sqrt(t0) * (su*su*su/3.0 + su);
- else t0 = t;
- }
-
- if (m0 < 0.0) m0 = m0 + p2;
- if (ran < 0.0) ran = ran + p2;
- if (aper < 0.0) aper = aper + p2;
- m0 = rddg * m0;
- ran = rddg * ran;
- aper = rddg * aper;
- inc = rddg * inc;
-
- if (ecc > 1.0) m0 = 0.0;
- if (pflg)
- {
- a = p;
- m0 = -1.0;
- }
- }
-
-/*-------------------- QuickSun --------------------------------------*/
-
-Vec3 QuickSun (double t) // low precision position of the Sun at time t
- {
- /* Low precision position of the Sun at time t given in
- Julian centuries since J2000.
- Valid only between 1950 and 2050.
-
- Returns the position vector (in A.U.) in ecliptic coordinates
-
- */
- double n, g, l;
- Vec3 rs;
-
- n = 36525.0 * t;
- l = 280.460 + 0.9856474*n; // mean longitude
- g = 357.528 + 0.9856003*n; // mean anomaly
- l = 2.0*M_PI * frac(l/360.0);
- g = 2.0*M_PI * frac(g/360.0);
- l = l + (1.915*sin(g) + 0.020*sin(2.0*g)) * 1.74532925199e-2; // ecl.long.
- g = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2.0*g); // radius
- rs[0] = g * cos(l);
- rs[1] = g * sin(l);
- rs[2] = 0;
-
- return rs;
- }
-
-/*-------------------- Class Sun200 --------------------------------------*/
- /*
- Ecliptic coordinates (in A.U.) and velocity (in A.U./day)
- of the sun for Equinox of Date given in Julian centuries since J2000.
- ======================================================================
- */
-Sun200::Sun200 ()
- { }
-
-Vec3 Sun200::position (double t) // position of the Sun at time t
- {
- Vec3 rs, vs;
-
- state (t, rs, vs);
-
- return rs;
- }
-
-void Sun200::state (double t, Vec3& rs, Vec3& vs)
- {
- /* State vector rs (position) and vs (velocity) of the Sun in
- ecliptic of date coordinates at time t (in Julian Centruries
- since J2000).
- */
- const double p2 = 2.0 * M_PI;
-
- double l, b, r;
- int i;
-
- tt = t;
-
- dl = 0.0; dr = 0.0; db = 0.0;
- m2 = p2 * frac(0.1387306 + 162.5485917*t);
- m3 = p2 * frac(0.9931266 + 99.9973604*t);
- m4 = p2 * frac(0.0543250 + 53.1666028*t);
- m5 = p2 * frac(0.0551750 + 8.4293972*t);
- m6 = p2 * frac(0.8816500 + 3.3938722*t);
- d = p2 * frac(0.8274 + 1236.8531*t);
- a= p2 * frac(0.3749 + 1325.5524*t);
- uu = p2 * frac(0.2591 + 1342.2278*t);
-
- c3[1] = 1.0; s3[1] = 0.0;
- c3[2] = cos(m3); s3[2] = sin(m3);
- c3[0] = c3[2]; s3[0] = -s3[2];
-
- for (i=3; i<9; ++i)
- addthe(c3[i-1],s3[i-1],c3[2],s3[2],c3[i],s3[i]);
- pertven(); pertmar(); pertjup(); pertsat(); pertmoo();
-
- dl = dl + 6.4 * sin(p2*(0.6983 + 0.0561*t))
- + 1.87 * sin(p2*(0.5764 + 0.4174*t))
- + 0.27 * sin(p2*(0.4189 + 0.3306*t))
- + 0.20 * sin(p2*(0.3581 + 2.4814*t));
- l = p2 * frac(0.7859453 + m3/p2 +((6191.2+1.1*t)*t+dl)/1296.0E3);
- r = 1.0001398 - 0.0000007*t + dr*1E-6;
- b = db * 4.8481368111E-6;
-
- cl= cos(l); sl=sin(l); cb=cos(b); sb=sin(b);
- rs[0]=r*cl*cb; rs[1]=r*sl*cb; rs[2]=r*sb;
-
- /* velocity calculation
- gms=2.9591202986e-4 AU^3/d^2
- e=0.0167086
- sqrt(gms/a)=1.7202085e-2 AU/d
- sqrt(gms/a)*sqrt(1-e^2)=1.71996836e-2 AU/d
- nu=M + 2e*sin(M) = M + 0.0334172 * sin(M) */
-
- uu = m3 + 0.0334172*sin(m3); // eccentric Anomaly E
- d = cos(uu);
- uu = sin(uu);
- a = 1.0 - 0.0167086*d;
- vs[0] = -1.7202085e-2*uu/a; // velocity in orbit plane
- vs[1] = 1.71996836e-2*d/a;
- uu = atan2 (0.9998604*uu, (d-0.0167086)); // true anomaly
- d = cos(uu);
- uu = sin(uu);
- dr = d*vs[0]+uu*vs[1];
- dl = (d*vs[1]-uu*vs[0]) / r;
-
- vs[0] = dr*cl*cb - dl*r*sl*cb;
- vs[1] = dr*sl*cb + dl*r*cl*cb;
- vs[2] = dr*sb;
- }
-
-void Sun200::addthe (double c1, double s1, double c2, double s2,
- double& cc, double& ss)
- {
- cc=c1*c2-s1*s2;
- ss=s1*c2+c1*s2;
- }
-
-void Sun200::term (int i1, int i, int it, double dlc, double dls, double drc,
- double drs, double dbc, double dbs)
- {
- if (it == 0) addthe (c3[i1+1],s3[i1+1],c[i+8],s[i+8],u,v);
- else
- {
- u=u*tt;
- v=v*tt;
- }
- dl = dl + dlc*u + dls*v;
- dr = dr + drc*u + drs*v;
- db = db + dbc*u + dbs*v;
- }
-
-void Sun200::pertven() // Kepler terms and perturbations by Venus
- {
- int i;
-
- c[8]=1.0; s[8]=0.0; c[7]=cos(m2); s[7]=-sin(m2);
- for (i=7; i>2; i--)
- addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
-
- term (1, 0, 0, -0.22, 6892.76, -16707.37, -0.54, 0.0, 0.0);
- term (1, 0, 1, -0.06, -17.35, 42.04, -0.15, 0.0, 0.0);
- term (1, 0, 2, -0.01, -0.05, 0.13, -0.02, 0.0, 0.0);
- term (2, 0, 0, 0.0, 71.98, -139.57, 0.0, 0.0, 0.0);
- term (2, 0, 1, 0.0, -0.36, 0.7, 0.0, 0.0, 0.0);
- term (3, 0, 0, 0.0, 1.04, -1.75, 0.0, 0.0, 0.0);
- term (0, -1, 0, 0.03, -0.07, -0.16, -0.07, 0.02, -0.02);
- term (1, -1, 0, 2.35, -4.23, -4.75, -2.64, 0.0, 0.0);
- term (1, -2, 0, -0.1, 0.06, 0.12, 0.2, 0.02, 0.0);
- term (2, -1, 0, -0.06, -0.03, 0.2, -0.01, 0.01, -0.09);
- term (2, -2, 0, -4.7, 2.9, 8.28, 13.42, 0.01, -0.01);
- term (3, -2, 0, 1.8, -1.74, -1.44, -1.57, 0.04, -0.06);
- term (3, -3, 0, -0.67, 0.03, 0.11, 2.43, 0.01, 0.0);
- term (4, -2, 0, 0.03, -0.03, 0.1, 0.09, 0.01, -0.01);
- term (4, -3, 0, 1.51, -0.4, -0.88, -3.36, 0.18, -0.1);
- term (4, -4, 0, -0.19, -0.09, -0.38, 0.77, 0.0, 0.0);
- term (5, -3, 0, 0.76, -0.68, 0.3, 0.37, 0.01, 0.0);
- term (5, -4, 0, -0.14, -0.04, -0.11, 0.43, -0.03, 0.0);
- term (5, -5, 0, -0.05, -0.07, -0.31, 0.21, 0.0, 0.0);
- term (6, -4, 0, 0.15, -0.04, -0.06, -0.21, 0.01, 0.0);
- term (6, -5, 0, -0.03, -0.03, -0.09, 0.09, -0.01, 0.0);
- term (6, -6, 0, 0.0, -0.04, -0.18, 0.02, 0.0, 0.0);
- term (7, -5, 0, -0.12, -0.03, -0.08, 0.31, -0.02, -0.01);
- }
-
-void Sun200::pertmar() // Kepler terms and perturbations by Mars
- {
- int i;
-
- c[7] = cos(m4); s[7] = -sin(m4);
- for (i=7; i>0; i--)
- addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
- term (1, -1, 0, -0.22, 0.17, -0.21, -0.27, 0.0, 0.0);
- term (1, -2, 0, -1.66, 0.62, 0.16, 0.28, 0.0, 0.0);
- term (2, -2, 0, 1.96, 0.57, -1.32, 4.55, 0.0, 0.01);
- term (2, -3, 0, 0.4, 0.15, -0.17, 0.46, 0.0, 0.0);
- term (2, -4, 0, 0.53, 0.26, 0.09, -0.22, 0.0, 0.0);
- term (3, -3, 0, 0.05, 0.12, -0.35, 0.15, 0.0, 0.0);
- term (3, -4, 0, -0.13, -0.48, 1.06, -0.29, 0.01, 0.0);
- term (3, -5, 0, -0.04, -0.2, 0.2, -0.04, 0.0, 0.0);
- term (4, -4, 0, 0.0, -0.03, 0.1, 0.04, 0.0, 0.0);
- term (4, -5, 0, 0.05, -0.07, 0.2, 0.14, 0.0, 00);
- term (4, -6, 0, -0.1, 0.11, -0.23, -0.22, 0.0, 0.0);
- term (5, -7, 0, -0.05, 0.0, 0.01, -0.14, 0.0, 0.0);
- term (5, -8, 0, 0.05, 0.01, -0.02, 0.1, 0.0, 0.0);
- }
-
-void Sun200::pertjup() // Kepler terms and perturbations by Jupiter
- {
- int i;
-
- c[7] = cos(m5); s[7] = -sin(m5);
- for (i=7; i>4; i--)
- addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
- term (1, -1, 0, 0.01, 0.07, 0.18, -0.02, 0.0, -0.02);
- term (0, -1, 0, -0.31, 2.58, 0.52, 0.34, 0.02, 0.0);
- term (1, -1, 0, -7.21, -0.06, 0.13, -16.27, 0.0, -0.02);
- term (1, -2, 0, -0.54, -1.52, 3.09, -1.12, 0.01, -0.17);
- term (1, -3, 0, -0.03, -0.21, 0.38, -0.06, 0.0, -0.02);
- term (2, -1, 0, -0.16, 0.05, -0.18, -0.31, 0.01, 0.0);
- term (2, -2, 0, 0.14, -2.73, 9.23, 0.48, 0.0, 0.0);
- term (2, -3, 0, 0.07, -0.55, 1.83, 0.25, 0.01, 0.0);
- term (2, -4, 0, 0.02, -0.08, 0.25, 0.06, 0.0, 0.0);
- term (3, -2, 0, 0.01, -0.07, 0.16, 0.04, 0.0, 0.0);
- term (3, -3, 0, -0.16, -0.03, 0.08, -0.64, 0.0, 0.0);
- term (3, -4, 0, -0.04, -0.01, 0.03, -0.17, 0.0, 0.0);
- }
-
-void Sun200::pertsat() // Kepler terms and perturbations by Saturn
- {
- c[7] = cos(m6); s[7] = -sin(m6);
- addthe(c[7],s[7],c[7],s[7],c[6],s[6]);
- term (0, -1, 0, 0.0, 0.32, 0.01, 0.0, 0.0, 0.0);
- term (1, -1, 0, -0.08, -0.41, 0.97, -0.18, 0.0, -0.01);
- term (1, -2, 0, 0.04, 0.1, -0.23, 0.1, 0.0, 0.0);
- term (2, -2, 0, 0.04, 0.1, -0.35, 0.13, 0.0, 0.0);
- }
-
-void Sun200::pertmoo() // corrections for Earth-Moon center of gravity
- {
- dl = dl + 6.45*sin(d) - 0.42*sin(d-a) + 0.18*sin(d+a) + 0.17*sin(d-m3) - \
0.06*sin(d+m3);
- dr = dr + 30.76*cos(d) - 3.06*cos(d-a) + 0.85*cos(d+a) - 0.58*cos(d+m3) + \
0.57*cos(d-m3);
- db = db + 0.576*sin(uu);
- }
-
-/*--------------------- Class moon200 ----------------------------------*/
-
- /*
- Position vector (in Earth radii) of the Moon referred to Ecliptic
- for Equinox of Date.
- t is the time in Julian centuries since J2000.
- */
-
-Moon200::Moon200 ()
- { }
-
-Vec3 Moon200::position (double t) // position of the Moon at time t
- {
- Vec3 rm;
-
- const double arc = 206264.81; // 3600*180/pi = ''/r
- const double pi2 = M_PI * 2.0;
-
- double lambda, beta, r, fac;
-
- minit(t);
- solar1(); solar2(); solar3(); solarn(n); planetary(t);
-
- lambda = pi2 * frac ((l0+dlam/arc) / pi2);
- if (lambda < 0) lambda = lambda + pi2;
- s = f + ds / arc;
-
- fac = 1.000002708 + 139.978 * dgam;
- beta = (fac*(18518.511+1.189+gam1c)*sin(s)-6.24*sin(3*s)+n);
- beta = beta * 4.8481368111e-6;
- sinpi = sinpi * 0.999953253;
- r = arc / sinpi;
-
- rm.assign (r, lambda, beta);
- rm = polcar(rm);
-
- return rm;
- }
-
-
-void Moon200::addthe (double c1, double s1, double c2, double s2,
- double& c, double& s)
- {
- c=c1*c2-s1*s2;
- s=s1*c2+c1*s2;
- }
-
-double Moon200::sinus (double phi)
- {
- /* sin(phi) in units of 1r = 360ø */
- return sin (2.0*M_PI * frac(phi));
- }
-
-void Moon200::long_periodic (double t)
- {
- /* long periodic changes of mean elements */
-
- double s1, s2, s3, s4, s5, s6, s7;
-
- s1 = sinus (0.19833 + 0.05611*t);
- s2 = sinus (0.27869 + 0.04508*t);
- s3 = sinus (0.16827 - 0.36903*t);
- s4 = sinus (0.34734 - 5.37261*t);
- s5 = sinus (0.10498 - 5.37899*t);
- s6 = sinus (0.42681 - 0.41855*t);
- s7 = sinus (0.14943 - 5.37511*t);
- dl0 = 0.84*s1 + 0.31*s2 + 14.27*s3 + 7.26*s4 + 0.28*s5 + 0.24*s6;
- dl = 2.94*s1 + 0.31*s2 + 14.27*s3 + 9.34*s4 + 1.12*s5 + 0.83*s6;
- dls = -6.4*s1 - 1.89*s6;
- df = 0.21*s1+0.31*s2+14.27*s3-88.7*s4-15.3*s5+0.24*s6-1.86*s7;
- dd = dl0 - dls;
- dgam = -3332E-9 * sinus (0.59734 - 5.37261*t)
- -539E-9 * sinus (0.35498 - 5.37899*t)
- -64E-9 * sinus (0.39943 - 5.37511*t);
- }
-
-void Moon200::minit(double t)
- {
- /* calculate mean elements l (moon), F (node distance)
- l' (sun) and D (moon's elongation) */
-
- const double arc = 206264.81; // 3600*180/pi = ''/r
- const double pi2 = M_PI * 2.0;
-
- int i, j, max;
- double t2, arg, fac;
-
- max = 0; // just to keep the compiler happy
- arg = 0;
- fac = 0;
-
- t2 = t * t;
- dlam = 0; ds = 0; gam1c = 0; sinpi = 3422.7;
- long_periodic (t);
- l0 = pi2*frac(0.60643382+1336.85522467*t-0.00000313*t2) + dl0/arc;
- l = pi2*frac(0.37489701+1325.55240982*t+0.00002565*t2) + dl/arc;
- ls = pi2*frac(0.99312619+99.99735956*t-0.00000044*t2) + dls/arc;
- f = pi2*frac(0.25909118+1342.22782980*t-0.00000892*t2) + df/arc;
- d = pi2*frac(0.82736186+1236.85308708*t-0.00000397*t2) + dd/arc;
- for (i=0; i<4; ++i)
- {
- switch (i)
- {
- case 0: arg=l; max=4; fac=1.000002208; break;
- case 1: arg=ls; max=3; fac=0.997504612-0.002495388*t; break;
- case 2: arg=f; max=4; fac=1.000002708+139.978*dgam; break;
- case 3: arg=d; max=6; fac=1.0; break;
- }
- co[6][i] = 1.0; co[7][i] = cos (arg) * fac;
- si[6][i] = 0.0; si[7][i] = sin (arg) * fac;
- for (j = 2; j <= max; ++j)
- addthe(co[j+5][i],si[j+5][i],co[7][i],si[7][i],co[j+6][i],si[j+6][i]);
- for (j = 1; j <= max; ++j)
- {
- co[6-j][i]=co[j+6][i];
- si[6-j][i]=-si[j+6][i];
- }
- }
- }
-
-void Moon200::term (int p, int q, int r, int s, double& x, double& y)
- {
- // calculate x=cos(p*arg1+q*arg2...); y=sin(p*arg1+q*arg2...)
- int i[4];
- int k;
-
- i[0]=p; i[1]=q; i[2]=r; i[3]=s; x=1.0; y=0.0;
- for (k=0; k<4; k++)
- if (i[k]!=0) addthe(x,y,co[i[k]+6][k],si[i[k]+6][k],x,y);
- }
-
-void Moon200::addsol(double coeffl, double coeffs, double coeffg,
- double coeffp, int p, int q, int r, int s)
- {
- double x, y;
- term (p,q,r,s,x,y);
- dlam = dlam + coeffl*y; ds = ds + coeffs*y;
- gam1c = gam1c + coeffg*x; sinpi = sinpi + coeffp*x;
- }
-
-void Moon200::solar1()
- {
- addsol ( 13.902, 14.06, -0.001, 0.2607, 0, 0, 0, 4);
- addsol ( 0.403, -4.01, 0.394, 0.0023, 0, 0, 0, 3);
- addsol ( 2369.912, 2373.36, 0.601, 28.2333, 0, 0, 0, 2);
- addsol ( -125.154, -112.79, -0.725, -0.9781, 0, 0, 0, 1);
- addsol ( 1.979, 6.98, -0.445, 0.0433, 1, 0, 0, 4);
- addsol (191.953, 192.72, 0.029, 3.0861, 1, 0, 0, 2);
- addsol (-8.466, -13.51, 0.455, -0.1093, 1, 0, 0, 1);
- addsol (22639.5, 22609.07, 0.079, 186.5398, 1, 0, 0, 0);
- addsol (18.609, 3.59, -0.094, 0.0118, 1, 0, 0, -1);
- addsol (-4586.465, -4578.13, -0.077, 34.3117, 1, 0, 0, -2);
- addsol (3.215, 5.44, 0.192, -0.0386, 1, 0, 0, -3);
- addsol (-38.428, -38.64, 0.001, 0.6008, 1, 0, 0, -4);
- addsol (-0.393, -1.43, -0.092, 0.0086, 1, 0, 0, -6);
- addsol (-0.289, -1.59, 0.123, -0.0053, 0, 1, 0, 4);
- addsol (-24.420, -25.1, 0.04, -0.3, 0, 1, 0, 2);
- addsol (18.023, 17.93, 0.007, 0.1494, 0, 1, 0, 1);
- addsol (-668.146, -126.98, -1.302, -0.3997, 0, 1, 0, 0);
- addsol (0.56, 0.32, -0.001, -0.0037, 0, 1, 0, -1);
- addsol (-165.145, -165.06, 0.054, 1.9178, 0, 1, 0, -2);
- addsol (-1.877, -6.46, -0.416, 0.0339, 0, 1, 0, -4);
- addsol (0.213, 1.02, -0.074, 0.0054, 2, 0, 0, 4);
- addsol (14.387, 14.78, -0.017, 0.2833, 2, 0, 0, 2);
- addsol (-0.586, -1.2, 0.054, -0.01, 2, 0, 0, 1);
- addsol (769.016, 767.96, 0.107, 10.1657, 2, 0, 0, 0);
- addsol (1.75, 2.01, -0.018, 0.0155, 2, 0, 0, -1);
- addsol (-211.656, -152.53, 5.679, -0.3039, 2, 0, 0, -2);
- addsol (1.225, 0.91, -0.03, -0.0088, 2, 0, 0, -3);
- addsol (-30.773, -34.07, -0.308, 0.3722, 2, 0, 0, -4);
- addsol (-0.57, -1.4, -0.074, 0.0109, 2, 0, 0, -6);
- addsol (-2.921, -11.75, 0.787, -0.0484, 1, 1, 0, 2);
- addsol (1.267, 1.52, -0.022, 0.0164, 1, 1, 0, 1);
- addsol (-109.673, -115.18, 0.461, -0.949, 1, 1, 0, 0);
- addsol (-205.962, -182.36, 2.056, 1.4437, 1, 1, 0, -2);
- addsol (0.233, 0.36, 0.012, -0.0025, 1, 1, 0, -3);
- addsol (-4.391, -9.66, -0.471, 0.0673, 1, 1, 0, -4);
- }
-
-void Moon200::solar2()
- {
- addsol (0.283, 1.53, -0.111, 0.006, 1, -1, 0, 4);
- addsol (14.577, 31.7, -1.54, 0.2302, 1, -1, 0, 2);
- addsol (147.687, 138.76, 0.679, 1.1528, 1, -1, 0, 0);
- addsol (-1.089, 0.55, 0.021, 0.0, 1, -1, 0, -1);
- addsol (28.475, 23.59, -0.443, -0.2257, 1, -1, 0, -2);
- addsol (-0.276, -0.38, -0.006, -0.0036, 1, -1, 0, -3);
- addsol (0.636, 2.27, 0.146, -0.0102, 1, -1, 0, -4);
- addsol (-0.189, -1.68, 0.131, -0.0028, 0, 2, 0, 2);
- addsol (-7.486, -0.66, -0.037, -0.0086, 0, 2, 0, 0);
- addsol (-8.096, -16.35, -0.74, 0.0918, 0, 2, 0, -2);
- addsol (-5.741, -0.04, 0.0, -0.0009, 0, 0, 2, 2);
- addsol (0.255, 0.0, 0.0, 0.0, 0, 0, 2, 1);
- addsol (-411.608, -0.2, 0.0, -0.0124, 0, 0, 2, 0);
- addsol (0.584, 0.84, 0.0, 0.0071, 0, 0, 2, -1);
- addsol (-55.173, -52.14, 0.0, -0.1052, 0, 0, 2, -2);
- addsol (0.254, 0.25, 0.0, -0.0017, 0, 0, 2, -3);
- addsol (0.025, -1.67, 0.0, 0.0031, 0, 0, 2, -4);
- addsol (1.06, 2.96, -0.166, 0.0243, 3, 0, 0, 2);
- addsol (36.124, 50.64, -1.3, 0.6215, 3, 0, 0, 0);
- addsol (-13.193, -16.4, 0.258, -0.1187, 3, 0, 0, -2);
- addsol (-1.187, -0.74, 0.042, 0.0074, 3, 0, 0, -4);
- addsol (-0.293, -0.31, -0.002, 0.0046, 3, 0, 0, -6);
- addsol (-0.29, -1.45, 0.116, -0.0051, 2, 1, 0, 2);
- addsol (-7.649, -10.56, 0.259, -0.1038, 2, 1, 0, 0);
- addsol (-8.627, -7.59, 0.078, -0.0192, 2, 1, 0, -2);
- addsol (-2.74, -2.54, 0.022, 0.0324, 2, 1, 0, -4);
- addsol (1.181, 3.32, -0.212, 0.0213, 2, -1, 0, 2);
- addsol (9.703, 11.67, -0.151, 0.1268, 2, -1, 0, 0);
- addsol (-0.352, -0.37, 0.001, -0.0028, 2, -1, 0, -1);
- addsol (-2.494, -1.17, -0.003, -0.0017, 2, -1, 0, -2);
- addsol (0.36, 0.2, -0.012, -0.0043, 2, -1, 0, -4);
- addsol (-1.167, -1.25, 0.008, -0.0106, 1, 2, 0, 0);
- addsol (-7.412, -6.12, 0.117, 0.0484, 1, 2, 0, -2);
- addsol (-0.311, -0.65, -0.032, 0.0044, 1, 2, 0, -4);
- addsol (0.757, 1.82, -0.105, 0.0112, 1, -2, 0, 2);
- addsol (2.58, 2.32, 0.027, 0.0196, 1, -2, 0, 0);
- addsol (2.533, 2.4, -0.014, -0.0212, 1, -2, 0, -2);
- addsol (-0.344, -0.57, -0.025, 0.0036, 0, 3, 0, -2);
- addsol (-0.992, -0.02, 0.0, 0.0, 1, 0, 2, 2);
- addsol (-45.099, -0.02, 0.0, -0.0010, 1, 0, 2, 0);
- addsol (-0.179, -9.52, 0.0, -0.0833, 1, 0, 2, -2);
- addsol (-0.301, -0.33, 0.0, 0.0014, 1, 0, 2, -4);
- addsol (-6.382, -3.37, 0.0, -0.0481, 1, 0, -2, 2);
- addsol (39.528, 85.13, 0.0, -0.7136, 1, 0, -2, 0);
- addsol (9.366, 0.71, 0.0, -0.0112, 1, 0, -2, -2);
- addsol (0.202, 0.02, 0.0, 0.0, 1, 0, -2, -4);
- }
-
-void Moon200::solar3()
- {
- addsol (0.415, 0.1, 0.0, 0.0013, 0, 1, 2, 0);
- addsol (-2.152, -2.26, 0.0, -0.0066, 0, 1, 2, -2);
- addsol (-1.44, -1.3, 0.0, 0.0014, 0, 1, -2, 2);
- addsol (0.384, -0.04, 0.0, 0.0, 0, 1, -2, -2);
- addsol (1.938, 3.6, -0.145, 0.0401, 4, 0, 0, 0);
- addsol (-0.952, -1.58, 0.052, -0.0130, 4, 0, 0, -2);
- addsol (-0.551, -0.94, 0.032, -0.0097, 3, 1, 0, 0);
- addsol (-0.482, -0.57, 0.005, -0.0045, 3, 1, 0, -2);
- addsol (0.681, 0.96, -0.026, 0.0115, 3, -1, 0, 0);
- addsol (-0.297, -0.27, 0.002, -0.0009, 2, 2, 0, -2);
- addsol (0.254, 0.21, -0.003, 0.0, 2, -2, 0, -2);
- addsol (-0.25, -0.22, 0.004, 0.0014, 1, 3, 0, -2);
- addsol (-3.996, 0.0, 0.0, 0.0004, 2, 0, 2, 0);
- addsol (0.557, -0.75, 0.0, -0.009, 2, 0, 2, -2);
- addsol (-0.459, -0.38, 0.0, -0.0053, 2, 0, -2, 2);
- addsol (-1.298, 0.74, 0.0, 0.0004, 2, 0, -2, 0);
- addsol (0.538, 1.14, 0.0, -0.0141, 2, 0, -2, -2);
- addsol (0.263, 0.02, 0.0, 0.0, 1, 1, 2, 0);
- addsol (0.426, 0.07, 0.0, -0.0006, 1, 1, -2, -2);
- addsol (-0.304, 0.03, 0.0, 0.0003, 1, -1, 2, 0);
- addsol (-0.372, -0.19, 0.0, -0.0027, 1, -1, -2, 2);
- addsol (0.418, 0.0, 0.0, 0.0, 0, 0, 4, 0);
- addsol (-0.330, -0.04, 0.0, 0.0, 3, 0, 2, 0);
- }
-
-void Moon200::addn (double coeffn, int p, int q, int r, int s,
- double& n, double&x, double& y)
- {
- term (p,q,r,s,x,y);
- n=n+coeffn*y;
- }
-
-void Moon200::solarn (double& n)
- {
- // perturbation N of ecliptic latitude
- double x, y;
-
- n = 0.0;
- addn (-526.069, 0, 0, 1, -2, n, x, y);
- addn (-3.352, 0, 0, 1, -4, n, x, y);
- addn (44.297, 1, 0, 1, -2, n, x, y);
- addn (-6.0, 1, 0, 1, -4, n, x, y);
- addn (20.599, -1, 0, 1, 0, n, x, y);
- addn (-30.598, -1, 0, 1, -2, n, x, y);
- addn (-24.649, -2, 0, 1, 0, n, x, y);
- addn (-2.0, -2, 0, 1, -2, n, x, y);
- addn (-22.571, 0, 1, 1, -2, n, x, y);
- addn (10.985, 0, -1, 1, -2, n, x, y);
- }
-
-void Moon200::planetary (double t)
- {
- // perturbations of the ecliptic longitude by Venus and Jupiter
-
- dlam = dlam
- + 0.82*sinus(0.7736 -62.5512*t)+0.31*sinus(0.0466-125.1025*t)
- + 0.35*sinus(0.5785-25.1042*t)+0.66*sinus(0.4591+1335.8075*t)
- + 0.64*sinus(0.313-91.568*t)+1.14*sinus(0.148+1331.2898*t)
- + 0.21*sinus(0.5918+1056.5859*t)+0.44*sinus(0.5784+1322.8595*t)
- + 0.24*sinus(0.2275-5.7374*t)+0.28*sinus(0.2965+2.6929*t)
- + 0.33*sinus(0.3132+6.3368*t);
- }
-
-/*-------------------- Class Eclipse --------------------------------------*/
-/*
- Calculalations for Solar and Lunar Eclipses
- */
-
-Eclipse::Eclipse()
- {
- // some assignments just in case
- rs.assign(1,0,0);
- rm.assign(1,0,0);
- eshadow.assign (1,0,0);
- rint.assign (1,0,0);
- d_umbra = 0;
- ep2 = 0;
- }
-
-int Eclipse::solar (double jd, double tdut, double& phi, double& lamda)
- {
- /* Calculates various items about a solar eclipse.
- INPUT:
- jd = Modified Julian Date (UT)
- tdut = TDT - UT in sec
-
- OUTPUT:
- phi, lamda: Geographic latitude and longitude (in radians)
- of center of shadow if central eclipse
-
- RETURN:
- 0 : No eclipse
- 1 : Partial eclipse
- 2 : Non-central annular eclipse
- 3 : Non-central total eclipse
- 4 : Annular eclipse
- 5 : Total eclipse
- */
-
- const double flat = 0.996633; // flatting of the Earth
- const double ds = 218.245445; // diameter of Sun in Earth radii
- const double dm = 0.544986; // diameter of Moon in Earth radii
- double s0, s, dlt, r2, r0;
- int phase;
- Vec3 ve;
-
- // get the apparent equatorial coordinates of the sun and the moon
- equ_sun_moon(jd, tdut);
-
- rs[2]/=flat; // adjust for flatting of the Earth
- rm[2]/=flat;
- rint.assign ();
- lamda = 0;
- phi = 0;
-
- // intersect shadow axis with Earth
- eshadow = rm - rs;
- eshadow = vnorm(eshadow); // direction vector of shadow
-
- s0 = - dot(rm, eshadow); // distance Moon - fundamental plane
- r2 = dot (rm,rm);
- dlt = s0*s0 + 1.0 - r2;
- r0 = 1.0 - dlt;
- if (r0 > 0) r0 = sqrt (r0);
- else r0 = 0; // distance center of Earth - shadow axis
-
- r2 = abs(rs - rm);
- d_umbra = (ds - dm) * s0 / r2 - dm;// diameter of umbra at fundamental plane
- d_penumbra = (ds + dm) * s0 / r2 + dm;
-
- // get phase of eclipse
- if (r0 < 1.0)
- {
- if (dlt > 0) dlt = sqrt(dlt);
- else dlt = 0;
- s = s0 - dlt; // distance Moon - fundamental plane
- d_umbra = (ds - dm) * s / r2 - dm; // diameter of umbra at surface
- rint = rm + s * eshadow;
- rint[2] *= flat; // vector to intersection
- ve = carpol(rint);
- lamda = ve[1] - lsidtim(jd,0,ep2)*0.261799387799; // geographic coordinates
- if (lamda > M_PI) lamda -= 2.0*M_PI;
- if (lamda < (-M_PI)) lamda += 2.0*M_PI;
- phi = sqrt(rint[0]*rint[0] + rint[1]*rint[1])*0.993305615;
- phi = atan2(rint[2],phi);
-
- if (d_umbra > 0) phase = 4; // central annular eclipse
- else phase = 5; // central total eclipse
- }
- else
- {
- if (r0 < (1.0 + 0.5 * fabs(d_umbra)))
- {
- if (d_umbra > 0) phase = 2; // non-central annular eclipse
- else phase = 3; // non-central total eclipse
- }
- else
- {
- if (r0 < (1.0 + 0.5*d_penumbra)) phase = 1; // partial eclipse
- else phase = 0; // no eclipse
- }
- }
-
- rs[2]*=flat; // restore from flatting of the Earth
- rm[2]*=flat;
-
- return phase;
- }
-
-void Eclipse::maxpos (double jd, double tdut, double& phi, double& lamda)
- {
- /* Calculates the geographic position of the place of maximum eclipse
- at the time jd for a non-central eclipse. (NOTE that in case of
- a central eclipse the maximum position will be calculated by
- function solar. Do not use this function in that case as the
- result would be incorrect! No check is being done in this routine
- whether an eclipse is visible at all. Use maxpos only in case of a
- confirmed partial or non-central eclipse.
-
- INPUT:
- jd = Modified Julian Date (UT)
- tdut = TDT - UT in sec
-
- OUTPUT:
- phi, lamda: Geographic latitude and longitude (in radians)
- of maximum eclipse at the time
- */
-
- const double flat = 0.996633; // flatting of the Earth
- double s0;
- Vec3 ve;
-
- // get the apparent equatorial coordinates of the sun and the moon
- equ_sun_moon(jd, tdut);
-
- rs[2]/=flat; // adjust for flatting of the Earth
- rm[2]/=flat;
- rint.assign ();
- lamda = 0;
- phi = 0;
-
- // intersect shadow axis with Earth
- eshadow = rm - rs;
- eshadow = vnorm(eshadow); // direction vector of shadow
-
- s0 = - dot(rm, eshadow); // distance Moon - fundamental plane
- rint = rm + s0 * eshadow;
- rint = vnorm(rint); // normalize to 1 Earth radius
- rint[2] *= flat; // vector to position of maximum eclipse
- ve = carpol(rint);
-
- lamda = ve[1] - lsidtim(jd,0,ep2)*0.261799387799; // geographic coordinates
- if (lamda > M_PI) lamda -= 2.0*M_PI;
- if (lamda < (-M_PI)) lamda += 2.0*M_PI;
- phi = sqrt(rint[0]*rint[0] + rint[1]*rint[1])*0.993305615;
- phi = atan2(rint[2],phi);
-
- rs[2]*=flat; // restore from flatting of the Earth
- rm[2]*=flat;
- }
-
-void Eclipse::penumd (double jd, double tdut, Vec3& vrm, Vec3& ves,
- double& dpn, double& pang)
- {
- /* Calculates various items needed for finding out the northern and
- southern border of the penumbra during a solar eclipse.
-
- INPUT:
- jd = Modified Julian Date (UT)
- tdut = TDT - UT in sec
-
- OUTPUT:
- vrm = position vector of Moon adjusted for flattening
- ves = unit vector pointing into the direction of the center of the shadow
- dpn = diameter of the penumbra at the fundamental plane
- pang = angle of penumbral half-cone in radians
- */
-
- const double flat = 0.996633; // flatting of the Earth
- const double ds = 218.245445; // diameter of Sun in Earth radii
- const double dm = 0.544986; // diameter of Moon in Earth radii
- double s0, r2;
-
- // get the apparent equatorial coordinates of the sun and the moon
- equ_sun_moon(jd, tdut);
- rs[2]/=flat; // adjust for flatting of the Earth
- rm[2]/=flat;
-
- // intersect shadow axis with Earth
- eshadow = rm - rs;
- pang = abs(eshadow);
- eshadow = vnorm(eshadow); // direction vector of shadow
- ves = eshadow;
- vrm = rm;
-
- s0 = - dot(rm, eshadow); // distance Moon - fundamental plane
-
- r2 = abs(rs - rm);
- dpn = (ds + dm) * s0 / r2 + dm;
-
- // penumbral angle
- pang = asin((dm + ds) / (2.0 * pang));
-
- rs[2]*=flat; // restore from flatting of the Earth
- rm[2]*=flat;
- }
-
-void Eclipse::umbra (double jd, double tdut, Vec3& vrm, Vec3& ves,
- double& dpn, double& pang)
- {
- /* Calculates various items needed for finding out the northern and
- southern border of the umbra during a solar eclipse.
-
- INPUT:
- jd = Modified Julian Date (UT)
- tdut = TDT - UT in sec
-
- OUTPUT:
- vrm = position vector of Moon adjusted for flattening
- ves = unit vector pointing into the direction of the center of the shadow
- dpn = diameter of the umbra at the fundamental plane
- pang = angle of umbral half-cone in radians
- */
-
- const double flat = 0.996633; // flatting of the Earth
- const double ds = 218.245445; // diameter of Sun in Earth radii
- const double dm = 0.544986; // diameter of Moon in Earth radii
- double s0, r2;
-
- // get the apparent equatorial coordinates of the sun and the moon
- equ_sun_moon(jd, tdut);
- rs[2]/=flat; // adjust for flatting of the Earth
- rm[2]/=flat;
-
- // intersect shadow axis with Earth
- eshadow = rm - rs;
- pang = abs(eshadow);
- eshadow = vnorm(eshadow); // direction vector of shadow
- ves = eshadow;
- vrm = rm;
-
- s0 = - dot(rm, eshadow); // distance Moon - fundamental plane
-
- r2 = abs(rs - rm);
- dpn = (ds - dm) * s0 / r2 - dm;
-
- // umbral angle
- pang = asin((ds - dm) / (2.0 * pang));
-
- rs[2]*=flat; // restore from flatting of the Earth
- rm[2]*=flat;
- }
-
-void Eclipse::equ_sun_moon(double jd, double tdut)
- {
- /* Get the equatorial coordinates of the Sun and the Moon and
- store them in rs and rm.
- Also store the time t in Julian Centuries and the correction
- ep2 for the Apparent Sidereal Time.
- jd = Modified Julian Date (UT)
- tdut = TDT - UT in sec
- */
- double ae = 23454.77992; // 149597870.0/6378.14 = 1AE -> Earth Radii
- Mat3 mx;
-
- t = julcent (jd) + tdut / 3.15576e9; // =(86400.0 * 36525.0);
- rs = sun.position(t);
- rm = moon.position(t);
- rs = eclequ(t,rs);
- rm = eclequ(t,rm);
-
- // correct Moon coordinates for center of figure
- mx = zrot(-2.4240684e-6); // +0.5" in longitude
- rm = mxvct(mx,rm);
- mx = yrot(-1.2120342e-6); // -0.25" in latitude
- rm = mxvct(mx,rm);
- mx = nutmat(t,ep2); // nutation
- rs = mxvct(mx,rs); // apparent coordinates
- rs = aberrat(t,rs);
- rs *= ae;
- rm = mxvct(mx,rm);
- }
-
-double Eclipse::duration (double jd, double tdut, double& width)
- {
- /* Get the duration of a central eclipse at the center point
- for MJD jd and TDT-UT tdut in seconds.
- Also return the width of the umbra in km.
- A call to solar with the respective jd must have been done
- prior to calling this function !!!
- */
- const double omega = 4.3755e-3; // radial velocity of Earth in rad/min.
-
- double dur, lm, pa, umbold;
- Vec3 rold, eold, rsold, rmold;
- Mat3 mx;
-
- // save old values for jd
- rold = rint;
- eold = eshadow;
- umbold = d_umbra;
- rsold = rs;
- rmold = rm;
-
- dur = 0.1; // 0.1 min
- if (solar(jd+dur/1440.0,tdut, lm, pa) > 3)
- {
- mx = zrot(dur*omega);
- rint = mxvct(mx,rint);
- rint -= rold;
- pa = dot (rint, eold);
- lm = dot (rint, rint) - pa*pa;
- if (lm > 0) lm = sqrt(lm);
- else lm = 0;
- if (lm > 0) dur = fabs(umbold) / lm * dur * 60.0;
- else dur = 0;
- }
-
- else dur = -1;
-
- // restore old values for jd
- d_umbra = umbold;
- eshadow = eold;
- eold = rold*rint; // direction perpendicular of apparent shadow movement
- rint = rold;
- rs = rsold;
- rm = rmold;
-
- // get width of umbra at center location
- rold = vnorm (eold);
- pa = dot(rold,eshadow);
- if (pa > 1.0) pa = 1.0;
- if (pa < -1.0) pa = -1.0;
- pa = fabs(sin(acos(pa)));
- if (pa < 0.00001) pa = 0.00001;
- width = d_umbra / pa * 6378.14; // allow for projection
- width = fabs(width);
-
- return dur;
- }
-
-Vec3 Eclipse::GetRSun () // get Earth - Sun vector in Earth radii
- {
- return rs;
- }
-
-Vec3 Eclipse::GetRMoon () // get Earth - Moon vector in Earth radii
- {
- return rm;
- }
-
-double Eclipse::GetEp2 () // get the ep2 value
- {
- return ep2;
- }
-
-int Eclipse::lunar (double jd, double tdut)
- {
- /* check whether lunar eclipse is in progress at
- Modified Julian Date jd (UT).
- tdut = TDT - UT in sec
-
- RETURN:
- 0 : No eclipse
- 1 : Partial Penumbral Eclipse
- 2 : Total Penumbral Eclipse
- 3 : Partial Umbral Eclipse
- 4 : Total Umbral Eclipse
- */
-
- const double dm = 0.544986; // diameter of Moon in Earth radii
- const double ds = 218.245445; // diameter of Sun in Earth radii
-
- double umbra, penumbra;
- double r2, s0, sep;
- int phase;
- Vec3 v1, v2;
-
- // get position of Sun and Moon
- equ_sun_moon(jd, tdut);
-
- // get radius of umbra and penumbra
- r2 = abs(rs);
- s0 = abs (rm);
- umbra = 1.02*fabs((ds - 2.0) * s0 / r2 - 2.0)*0.5; // radius of umbra
- penumbra = 1.02*fabs((ds + 2.0) * s0 / r2 + 2.0)*0.5;//radius of penumbra
- /* (the factor 1.02 allows for enlargment of shadow due to
- Earth's atmosphere) */
-
- // get angular separation of center of shadow and Moon
- r2 = abs(rm);
- sep = dot(rs,rm)/(abs(rs)*r2);
- if (fabs(sep) > 1.0) sep = 1.0;
- sep = acos(sep); // in radians
- sep = fabs(tan(sep)*r2); // distance of Moon and shadow in Earth radii
-
- // Now check the kind of eclipse
- if (sep < (umbra - dm/2.0)) phase = 4;
- else if (sep < (umbra + dm/2.0)) phase = 3;
- else if (sep < (penumbra - dm/2.0)) phase = 2;
- else if (sep < (penumbra + dm/2.0)) phase = 1;
- else phase = 0;
-
- return phase;
- }
-
-
diff --git a/src/plugins/render/satellites/mex/astrolib.h \
b/src/plugins/render/satellites/mex/astrolib.h deleted file mode 100644
index 8a878e8..0000000
--- a/src/plugins/render/satellites/mex/astrolib.h
+++ /dev/null
@@ -1,146 +0,0 @@
-//
-// This file is part of the Marble Virtual Globe.
-//
-// This program is free software licensed under the GNU LGPL. You can
-// find a copy of this license in LICENSE.txt in the top directory of
-// the source code.
-//
-// Copyright 2012 Gerhard HOLTKAMP
-//
-
-#if !defined(__astrolib_h)
-#define __astrolib_h
-
-#include "attlib.h"
-
-double ddd (int d, int m, double s); // deg, min, sec -> decimal degrees
-void dms (double dd, int &d, int &m, double &s); // dec deg -> deg, min, sec
-double mjd (int day, int month, int year, double hour); // modified Julian date
-double julcent (double mjuld); // Julian centuries since 2000.0
-void caldat (double mjd, int &day, int &month, int &year, double &hour);
-double DefTdUt (int yr); // default value for TDT - UT in seconds
-double lsidtim (double jd, double lambda, double ep2); // Sidereal Time
-
-double eps (double t); // obliquity of ecliptic
-Vec3 eclequ (double t, Vec3& r1); // ecliptic -> equatorial
-Vec3 equecl (double t, Vec3& r1); // equatorial -> ecliptic
-Mat3 pmatecl (double t1, double t2); // ecl. precession
-Mat3 pmatequ (double t1, double t2); // equ. precession
-Mat3 nutmat (double t, double& ep2, bool hipr = false); // nutation (equatorial)
-Mat3 nutecl (double t, double& ep2); // nutation matrix (ecliptic)
-Mat3 PoleMx (double xp, double yp); // Polar motion matrix
-Vec3 aberrat (double t, Vec3& ve); // aberration
-
-Vec3 GeoPos (double jd, double ep2, double lat, double lng, double ht);
-Vec3 GeoPos (double jd, double ep2, double lat, double lng, double ht,
- double xp, double yp);
-Vec3 EquHor (double jd, double ep2, double lat, double lng, Vec3 r);
-Vec3 HorEqu (double jd, double ep2, double lat, double lng, Vec3 r);
-void AppPos (double jd, double ep2, double lat, double lng, double ht,
- int solsys, Vec3 r, double& azim, double& elev, double& dist);
-void AppRADec (double jd, double ep2, double lat, double lng,
- double azim, double elev, double& ra, double& dec);
-double Refract (double h, double p = 1015.0, double t = 15.0); // refraction
-
-double eccanom (double man, double ecc); // eccentric anomaly
-double hypanom (double mh, double ecc); // hyperbolic anomaly
-void ellip (double gm, double t0, double t, double a, double ecc,
- double m0, Vec3& r1, Vec3& v1); // elliptic state vector
-void hyperb (double gm, double t0, double t, double a, double ecc,
- Vec3& r1, Vec3& v1); // hyperbolic state vector
-void parab (double gm, double t0, double t, double q, double ecc,
- Vec3& r1, Vec3& v1); // elliptic state vector
-void kepler (double gm, double t0, double t, double m0, double a, double ecc,
- double ran, double aper, double inc, Vec3& r1, Vec3& v1);
-void oscelm (double gm, double t, Vec3& r1, Vec3& v1,
- double& t0, double& m0, double& a, double& ecc,
- double& ran, double& aper, double& inc);
-
-Vec3 QuickSun (double t); // low precision position of the Sun at time t
-
-class Sun200 // Calculating the Sun in epoch J2000.0 coordinates
- {
- public:
- Sun200();
- Vec3 position (double t); // position of the Sun
- void state (double t, Vec3& rs, Vec3& vs); // State vector of the Sun
-
- private:
- double c3[9], s3[9];
- double c[9], s[9];
- double m2, m3, m4, m5, m6;
- double d, a, uu, tt;
- double cl, sl, cb, sb;
- double u, v, dl, dr, db;
- void addthe (double c1, double s1, double c2, double s2,
- double& cc, double& ss);
- void term (int i1, int i, int it, double dlc, double dls, double drc,
- double drs, double dbc, double dbs);
- void pertven();
- void pertmar();
- void pertjup();
- void pertsat();
- void pertmoo();
- };
-
-class Moon200 // Calculating the position of the Moon in J2000.0
- {
- public:
- Moon200();
- Vec3 position (double t); // position of the Moon
-
- private:
- double dgam, dlam, n, gam1c, sinpi;
- double l0, l, ls, f, d, s;
- double dl0, dl, dls, df, dd, ds;
- double co[13][4];
- double si[13][4];
- void addthe (double c1, double s1, double c2, double s2,
- double& c, double& s);
- double sinus (double phi);
- void long_periodic (double t);
- void minit(double t);
- void term (int p, int q, int r, int s, double& x, double& y);
- void addsol(double coeffl, double coeffs, double coeffg,
- double coeffp, int p, int q, int r, int s);
- void solar1();
- void solar2();
- void solar3();
- void addn (double coeffn, int p, int q, int r, int s,
- double& n, double&x, double& y);
- void solarn (double& n);
- void planetary (double t);
- };
-
-class Eclipse // Eclipse Calculations
- {
- public:
- Eclipse();
- int solar (double jd, double tdut, double& phi, double& lamda);
- void maxpos (double jd, double tdut, double& phi, double& lamda);
- void penumd (double jd, double tdut, Vec3& vrm, Vec3& ves,
- double& dpn, double& pang);
- void umbra (double jd, double tdut, Vec3& vrm, Vec3& ves,
- double& dpn, double& pang);
- double duration (double jd, double tdut, double& width);
- Vec3 GetRSun (); // get Earth - Sun vector in Earth radii
- Vec3 GetRMoon (); // get Earth - Moon vector in Earth radii
- double GetEp2 (); // get the ep2 value
- int lunar (double jd, double tdut);
-
- private:
- Sun200 sun;
- Moon200 moon;
- Vec3 rs, rm; // position of the Sun and the Moon
- Vec3 eshadow; // unit vector in direction of shadow
- Vec3 rint; // intersection shadow axis - Earth surface
- double t; // time in Julian Centuries
- double ep2; // correction for Apparent Sideral Time
- double d_umbra; // diameter of umbra in Earth Radii
- double d_penumbra; // diameter of penumbra in Earth Radii
- void equ_sun_moon(double jd, double tdut);
- };
-
-#endif // __astrolib_h sentry.
-
-
diff --git a/src/plugins/render/satellites/mex/attlib.cpp \
b/src/plugins/render/satellites/mex/attlib.cpp deleted file mode 100644
index 0264bd7..0000000
--- a/src/plugins/render/satellites/mex/attlib.cpp
+++ /dev/null
@@ -1,669 +0,0 @@
-//
-// This file is part of the Marble Virtual Globe.
-//
-// This program is free software licensed under the GNU LGPL. You can
-// find a copy of this license in LICENSE.txt in the top directory of
-// the source code.
-//
-// Copyright 2012 Gerhard HOLTKAMP
-//
-
-/***********************************************************************
- 3-Dim Vector and Matrix Definitions and Operations
-
- Linux C++ version
-
- Author: Gerhard HOLTKAMP 15-APR-2010
- ***********************************************************************/
-
-//#include <iostream>
-#include <cmath>
-using namespace std;
-
-#include "attlib.h"
-
-/*********************************************************************/
-double atan20 (double y, double x)
-{
- /* redefine atan2 so that it does'nt crash when both x and y are 0 */
- double result;
-
- if ((x == 0) && (y == 0)) result = 0;
- else result = atan2 (y, x);
-
- return result;
-}
-
-
- Vec3::Vec3(double x, double y, double z)
- {assign(x, y, z);}
-
- Vec3::Vec3 (const Vec3& c)
- {
- for (int j = 0; j<3; ++j) v[j] = c.v[j];
- }
-
-
- void Vec3::assign (double x, double y, double z)
- {
- v[0]=x; v[1]=y; v[2]=z;
- }
-
- double& Vec3::operator [] (unsigned index)
- {
- return (index < 3) ? * (v + index) : *v;
- }
-
- Vec3& Vec3::operator = (const Vec3& c)
- {
- for (int j = 0; j<3; ++j) v[j] = c.v[j];
- return *this;
- }
-
- Vec3& Vec3::operator += (const Vec3& c)
- {
- for (int j = 0; j<3; ++j) v[j] = v[j] + c.v[j];
- return *this;
- }
-
- Vec3& Vec3::operator -= (const Vec3& c)
- {
- for (int j = 0; j<3; ++j) v[j] = v[j] - c.v[j];
- return *this;
- }
-
- Vec3& Vec3::operator *= (const Vec3& c) // cross product
- {
- Vec3 result;
- result.v[0] = v[1] * c.v[2] - v[2] * c.v[1];
- result.v[1] = v[2] * c.v[0] - v[0] * c.v[2];
- result.v[2] = v[0] * c.v[1] - v[1] * c.v[0];
- for (int j = 0; j<3; ++j) v[j] = result.v[j];
- return *this;
- }
-
- Vec3& Vec3::operator *= (double r)
- {
- for (int j = 0; j<3; ++j) v[j] = r * v[j];
- return *this;
- }
-
- Vec3& Vec3::operator /= (double r)
- {
- double q;
-
- if (r < 1E-100) q = 0.0;
- else q = 1.0 / r;
- for (int j = 0; j<3; ++j) v[j] = q * v[j];
- return *this;
- }
-
- double abs(const Vec3& c) // absolute value
- {
- double result=0;
-
- for (int j=0; j<3; ++j) result += c.v[j]*c.v[j];
- return sqrt(result);
- }
-
- double dot (const Vec3& c1, const Vec3& c2) // dot product
- {
- double result=0;
-
- for (int j = 0; j<3; ++j) result += c1.v[j]*c2.v[j];
- return result;
- }
-
- Vec3 operator + (const Vec3& c1, const Vec3& c2)
- {
- Vec3 result;
-
- for (int j = 0; j<3; ++j) result.v[j] = c1.v[j] + c2.v[j];
- return result;
- }
-
- Vec3 operator - (const Vec3& c1, const Vec3& c2)
- {
- Vec3 result;
-
- for (int j = 0; j<3; ++j) result.v[j] = c1.v[j] - c2.v[j];
-
- return result;
- }
-
- Vec3 operator * (double r, const Vec3& c1)
- {
- Vec3 result;
-
- for (int j = 0; j<3; ++j) result.v[j] = r * c1.v[j];
- return result;
- }
-
- Vec3 operator * (const Vec3& c1, double r)
- {
- Vec3 result;
-
- for (int j = 0; j<3; ++j) result.v[j] = c1.v[j] * r;
- return result;
- }
-
- Vec3 operator * (const Vec3& c1, const Vec3& c2) // cross product
- {
- Vec3 result;
-
- result.v[0] = c1.v[1] * c2.v[2] - c1.v[2] * c2.v[1];
- result.v[1] = c1.v[2] * c2.v[0] - c1.v[0] * c2.v[2];
- result.v[2] = c1.v[0] * c2.v[1] - c1.v[1] * c2.v[0];
- return result;
- }
-
- Vec3 operator / (const Vec3& c1, double r)
- {
- double q;
- Vec3 result;
-
- if (r < 1E-100) q = 0.0;
- else q = 1.0 / r;
- for (int j = 0; j<3; ++j) result.v[j] = q * c1.v[j];
- return result;
- }
-
- Vec3 vnorm (const Vec3& c) // norm vector
- {
- int j;
- double q=0;
- Vec3 result;
-
- for (j=0; j<3; ++j) q += c.v[j]*c.v[j];
- q = sqrt(q);
- if (q < 1E-100) q = 0.0;
- else q = 1.0 / q;
- for (j=0; j<3; ++j) result.v[j] = c.v[j]*q;
- return result;
- }
-
- Vec3 carpol (const Vec3& c) // Cartesian -> Polar
- /* Convert vector from cartesian coordinates c = [x, y, z]
- into polar coordinates v = [rho, phi, theta]
- where rho is the radius, phi the "azimuth" and theta the
- "declination" (equatorial plane = 0, pole = 90ø)
- */
- {
- double x, y, z, r;
- Vec3 result;
-
- x = c.v[0];
- y = c.v[1];
- z = c.v[2];
- r = x*x + y*y;
- result.v[0] = sqrt (r + z*z);
- result.v[1] = atan20 (y, x);
- if (result.v[1] < 0) result.v[1] += 2.0*M_PI;
- r = sqrt (r);
- result.v[2] = atan20 (z, r);
-
- return result;
- }
-
- Vec3 polcar (const Vec3& c) // Polar -> Cartesian
- /* Convert vector from polar coordinates c = [rho, phi, theta]
- where rho is the radius, phi the "azimuth" and theta the
- "declination" (equatorial plane = 0, pole = 90ø)
- into vector v = [x, y, z] of cartesian coordinates.
- */
- {
- double r;
- Vec3 result;
-
- r = c.v[0] * cos(c.v[2]);
- result.v[0] = r * cos(c.v[1]);
- result.v[1] = r * sin(c.v[1]);
- result.v[2] = c.v[0] * sin(c.v[2]);
-
- return result;
- }
-
- ostream& operator << (ostream& os, const Vec3& c)
- {
- os << "[" << c.v[0] << "," << c.v[1] << "," << c.v[2] << "]";
- return os;
- }
-
-/*******************************************************************/
-
-Mat3::Mat3 (double x)
- {
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = x;
- }
-
-Mat3::Mat3 (const Mat3& c)
- {
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = c.m[i][j];
- }
-
-void Mat3::assign (double x11, double x12, double x13,
- double x21, double x22, double x23,
- double x31, double x32, double x33)
- {
- m[0][0] = x11;
- m[0][1] = x12;
- m[0][2] = x13;
- m[1][0] = x21;
- m[1][1] = x22;
- m[1][2] = x23;
- m[2][0] = x31;
- m[2][1] = x32;
- m[2][2] = x33;
- }
-
-void Mat3::assign (double x[3][3])
- {
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = x[i][j];
- }
-
-void Mat3::PutMij (double x, int i, int j)
- {
- if ((i>0) && (i<4) && (j>0) && (j<4)) m[i-1][j-1] = x;
- }
-
-double Mat3::GetMij (int i, int j)
- {
- double result;
-
- if ((i>0) && (i<4) && (j>0) && (j<4)) result = m[i-1][j-1];
- else result = 0;
- return result;
- }
-
-
-Mat3& Mat3::operator = (const Mat3& c)
- {
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = c.m[i][j];
-
- return (*this);
- }
-
-Mat3& Mat3::operator += (const Mat3& c)
- {
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = m[i][j] + c.m[i][j];
- return *this;
- }
-
-Mat3& Mat3::operator -= (const Mat3& c)
- {
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = m[i][j] - c.m[i][j];
- return *this;
- }
-
-Mat3& Mat3::operator *= (const Mat3& c)
- {
- int i, j, k;
- double r;
- Mat3 b;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j)
- {
- r = 0;
- for (k=0; k<3; ++k) r = r + c.m[i][k] * m[k][j];
- b.m[i][j] = r;
- }
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = b.m[i][j];
-
- return *this;
- }
-
-Mat3& Mat3::operator *= (double r)
- {
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = r*m[i][j];
- return *this;
- }
-
-Mat3& Mat3::operator /= (double r)
- {
- int i, j;
- double q;
-
- if (r < 1E-100) q = 0.0;
- else q = 1.0 / r;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) m[i][j] = m[i][j] * q;
- return *this;
- }
-
-Mat3 mxcon (double r) // constant matrix with all elements of value r
- {
- Mat3 result;
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) result.m[i][j] = r;
-
- return result;
- }
-
-Mat3 mxidn () // identity matrix
- {
- Mat3 result;
- int i;
-
- result = mxcon (0);
- for (i=0; i<3; ++i) result.m[i][i] = 1.0;
-
- return result;
- }
-
-Mat3 mxtrn (const Mat3& m1) // returns transposed of matrix m1
- {
- Mat3 result;
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) result.m[i][j] = m1.m[j][i];
-
- return result;
- }
-
-double mxdet (const Mat3& c) // returns determinant of matrix m1
- {
- double result;
-
- result = c.m[0][0]*c.m[1][1]*c.m[2][2] + c.m[0][1]*c.m[1][2]*c.m[2][0]
- + c.m[0][2]*c.m[1][0]*c.m[2][1] - c.m[0][2]*c.m[1][1]*c.m[2][0]
- - c.m[0][0]*c.m[1][2]*c.m[2][1] - c.m[0][1]*c.m[1][0]*c.m[2][2];
-
- return result;
- }
-
-Mat3 operator + (const Mat3& c1, const Mat3& c2)
- {
- Mat3 result;
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) result.m[i][j] = c1.m[i][j] + c2.m[i][j];
- return result;
- }
-
-Mat3 operator - (const Mat3& c1, const Mat3& c2)
- {
- Mat3 result;
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) result.m[i][j] = c1.m[i][j] - c2.m[i][j];
- return result;
- }
-
-Mat3 operator * (double r, const Mat3& c1)
- {
- Mat3 result;
- int i, j;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) result.m[i][j] = c1.m[i][j] * r;
- return result;
- }
-
-
-Mat3 operator * (const Mat3& c1, double r)
- {
- int i, j;
- Mat3 result;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) result.m[i][j] = c1.m[i][j] * r;
- return result;
- }
-
-Mat3 operator * (const Mat3& c1, const Mat3& c2)
- {
- int i, j, k;
- double r;
- Mat3 result;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j)
- {
- r = 0;
- for (k=0; k<3; k++) r = r + c1.m[i][k] * c2.m[k][j];
- result.m[i][j] = r;
- };
-
- return result;
- }
-
-Mat3 operator / (const Mat3& c1, double r)
- {
- Mat3 result;
- int i, j;
- double q;
-
- if (r < 1E-100) q = 0.0;
- else q = 1.0 / r;
-
- for (i=0; i<3; ++i)
- for (j=0; j<3; ++j) result.m[i][j] = c1.m[i][j] * q;
- return result;
- }
-
-Vec3 mxvct (const Mat3& m1, Vec3& v1)
- {
- int i, j;
- double r;
- Vec3 result;
-
- for (i=0; i<3; ++i)
- {
- r = 0;
- for (j=0; j<3; ++j) r = r + m1.m[i][j] * v1[j];
- result[i] = r;
- };
- return result;
- }
-
-Mat3 xrot (double a)
- {
- // Rotation matrix around x-axis with angle a (in radians).
- double c, s;
- Mat3 m1;
-
- c = cos (a);
- s = sin (a);
- m1.m[0][0] = 1.0;
- m1.m[0][1] = 0.0;
- m1.m[0][2] = 0.0;
- m1.m[1][0] = 0.0;
- m1.m[1][1] = c;
- m1.m[1][2] = s;
- m1.m[2][0] = 0.0;
- m1.m[2][1] = -s;
- m1.m[2][2] = c;
-
- return m1;
- }
-
-Mat3 yrot (double a)
- {
- // Rotation matrix around y-axis with angle a (in radians).
- double c, s;
- Mat3 m1;
-
- c = cos (a);
- s = sin (a);
- m1.m[0][0] = c;
- m1.m[0][1] = 0.0;
- m1.m[0][2] = -s;
- m1.m[1][0] = 0.0;
- m1.m[1][1] = 1.0;
- m1.m[1][2] = 0.0;
- m1.m[2][0] = s;
- m1.m[2][1] = 0.0;
- m1.m[2][2] = c;
-
- return m1;
- }
-
-Mat3 zrot (double a)
- {
- // Rotation matrix around z-axis with angle a (in radians).
- double c, s;
- Mat3 m1;
-
- c = cos (a);
- s = sin (a);
- m1.m[0][0] = c;
- m1.m[0][1] = s;
- m1.m[0][2] = 0.0;
- m1.m[1][0] = -s;
- m1.m[1][1] = c;
- m1.m[1][2] = 0.0;
- m1.m[2][0] = 0.0;
- m1.m[2][1] = 0.0;
- m1.m[2][2] = 1.0;
-
- return m1;
- }
-
-Mat3 csmx (double p, double y, double r)
- {
- // Form cosine matrix m from pitch p, yaw y and roll r angles.
- // The angles are to be given in radians.
- // Roll is rotation about the x-axis, pitch about y, yaw about z.
- //
- Mat3 pitch, yaw, roll, result;
-
- pitch = yrot (p);
- yaw = zrot (y);
- roll = xrot (r);
- result = yaw * pitch;
- result *= roll;
-
- return result;
- }
-
-void gpyr (const Mat3& m1, double& p, double& y, double& r)
- {
- // Get pitch p, yaw y and roll r angles (in radians) from
- // cosine matrix m1.
-
- y = asin (m1.m[0][1]);
- r = atan20 (-m1.m[2][1], m1.m[1][1]);
- p = atan20 (-m1.m[0][2], m1.m[0][0]);
- }
-
-void vcpy (Vec3& v, double& p, double& y)
- {
- // Convert direction given by cartesion vector v into a corresponding
- // pitch (p) / yaw (y) sequence. The angles are in radians.
-
- p = atan20 (-v[2], v[0]);
- y = atan20 (v[1], sqrt (v[0]*v[0] + v[2]*v[2]));
- }
-
-void vcrp (Vec3& v, double& p, double& r)
- {
- // Convert direction given by cartesian vector v into a corresponding
- // roll (r) / pitch (p) sequence. The angles are in radians.
-
- r = atan20 (v[1], -v[2]);
- p = M_PI / 2.0 - atan20 (v[0], sqrt (v[1]*v[1] + v[2]*v[2]));
- }
-
-void mxevc (const Mat3& m, double& a, Vec3& v)
- {
- // Get eigenvalue a and eigenvector v from cosine matrix m.
- // The eigenangle a is in radians.
-
- double ri, rj, rk, q1, q2, q3, q4;
-
- // using the 3-1-3 rotation matrix, first get the corresponding
- // angles ri, rj, rk of rotation about the x-, y-, z-axis.
- ri = atan20 (m.m[2][0], -m.m[2][1]);
- rj = 0.5 * acos (m.m[2][2]);
- rk = atan20 (m.m[0][2], m.m[1][2]);
-
- // now convert to quaternions
- q4 = sin (rj);
- q1 = q4 * cos (0.5 * (ri - rk));
- q2 = q4 * sin (0.5 * (ri - rk));
- q4 = cos (rj);
- q3 = q4 * sin (0.5 * (ri + rk));
- q4 = q4 * cos (0.5 * (ri + rk));
-
- // now get the eigen angle and eigen vector
- v.assign(q1, q2, q3);
- ri = abs (v);
- if (ri == 0)
- {
- // treat singularity for 3-1-3 matrix
- v.assign (0.0, 0.0, 1.0);
- q4 = 0.5 * sqrt (1.0 + m.m[0][0] + m.m[1][1] + m.m[2][2]);
- }
- else v /= ri;
-
- a = 2.0 * acos (q4);
- }
-
-Mat3 mxrox (double& a, Vec3& v)
- {
- // Convert eigenvalue a (eigen angle in radians) and eigenvector v
- // into a corresponding cosine rotation matrix m.
-
- double q1, q2, q3, q4, q12, q22, q32, q42;
- Mat3 result;
-
- // calculate quaternions and their squares
- q4 = sin ( 0.5 * a);
- q1 = v[0] * q4;
- q2 = v[1] * q4;
- q3 = v[2] * q4;
- q4 = cos (0.5 * a);
- q12 = q1 * q1;
- q22 = q2 * q2;
- q32 = q3 * q3;
- q42 = q4 * q4;
-
- // now get the matrix elements
- result.assign ((q12 - q22 - q32 + q42), (2.0 * (q1*q2 + q3*q4)),
- (2.0 * (q1*q3 - q2*q4)), (2.0 * (q1*q2 - q3*q4)),
- (-q12 + q22 - q32 + q42),(2.0 * (q2*q3 + q1*q4)),
- (2.0 * (q1*q3 + q2*q4)), (2.0 * (q2*q3 - q1*q4)),
- (-q12 - q22 + q32 + q42));
-
- return result;
- }
-
-
-ostream& operator << (ostream& os, const Mat3& c)
- {
- os << "[" << c.m[0][0] << "," << c.m[0][1] << "," << c.m[0][2] << "]" << endl
- << "[" << c.m[1][0] << "," << c.m[1][1] << "," << c.m[1][2] << "]" << endl
- << "[" << c.m[2][0] << "," << c.m[2][1] << "," << c.m[2][2] << "]" << endl;
- return os;
- }
-
diff --git a/src/plugins/render/satellites/mex/attlib.h \
b/src/plugins/render/satellites/mex/attlib.h deleted file mode 100644
index 7e209cf..0000000
--- a/src/plugins/render/satellites/mex/attlib.h
+++ /dev/null
@@ -1,109 +0,0 @@
-//
-// This file is part of the Marble Virtual Globe.
-//
-// This program is free software licensed under the GNU LGPL. You can
-// find a copy of this license in LICENSE.txt in the top directory of
-// the source code.
-//
-// Copyright 2012 Gerhard HOLTKAMP
-//
-
-#if !defined(__attlib_h)
-#define __attlib_h
-
-/***********************************************************************
- 3-Dim Vector and Matrix Definitions and Operations
-
- License: GNU LGPL Version 2+
-
- Author: Gerhard HOLTKAMP 14-JAN-2012
- ***********************************************************************/
-
-#include <iostream>
-
-double atan20 (double y, double x);
-
-class Vec3
-{
- private:
- double v[3];
-
- public:
-
- friend class Mat3;
-
- explicit Vec3(double x=0, double y=0, double z=0);
- Vec3 (const Vec3& c);
- void assign (double x=0, double y=0, double z=0);
- double& operator [] (unsigned index);
- Vec3& operator = (const Vec3& c);
- Vec3& operator += (const Vec3& c);
- Vec3& operator -= (const Vec3& c);
- Vec3& operator *= (const Vec3& c); // cross product
- Vec3& operator *= (double r);
- Vec3& operator /= (double r);
- friend double abs(const Vec3& c); // absolute value
- friend double dot (const Vec3& c1, const Vec3& c2); // dot product
- friend Vec3 operator + (const Vec3& c1, const Vec3& c2);
- friend Vec3 operator - (const Vec3& c1, const Vec3& c2);
- friend Vec3 operator * (double r, const Vec3& c1);
- friend Vec3 operator * (const Vec3& c1, double r);
- friend Vec3 operator * (const Vec3& c1, const Vec3& c2); // cross product
- friend Vec3 operator / (const Vec3& c1, double r);
- friend Vec3 vnorm(const Vec3& c); // norm vector
- friend Vec3 carpol (const Vec3& c); // Cartesian -> Polar
- friend Vec3 polcar (const Vec3& c); // Polar -> Cartesian
- friend std::ostream& operator << (std::ostream& os, const Vec3& c);
-};
-
-/********************************************************************/
-
-// class Mat3: public Vec3
-class Mat3
-{
- public:
- double m[3][3];
-
- explicit Mat3(double x=0);
- Mat3 (const Mat3& c);
- void assign (double x11, double x12, double x13, double x21, double x22,
- double x23, double x31, double x32, double x33);
- void assign (double x[3][3]); // assign matrix
- void PutMij (double x, int i, int j); // put single matrix element
- double GetMij (int i, int j); // get single matrix element
- Mat3& operator = (const Mat3& c);
- Mat3& operator += (const Mat3& c);
- Mat3& operator -= (const Mat3& c);
- Mat3& operator *= (const Mat3& c);
- Mat3& operator *= (double r);
- Mat3& operator /= (double r);
- friend Mat3 mxtrn (const Mat3& m1); // transposed matrix
- friend double mxdet (const Mat3& c); // determinant
- friend Mat3 operator + (const Mat3& c1, const Mat3& c2);
- friend Mat3 operator - (const Mat3& c1, const Mat3& c2);
- friend Mat3 operator * (double r, const Mat3& c1);
- friend Mat3 operator * (const Mat3& c1, double r);
- friend Mat3 operator * (const Mat3& c1, const Mat3& c2);
- friend Mat3 operator / (const Mat3& c1, double r);
- friend Vec3 mxvct (const Mat3& m1, Vec3& v1); // multiply vector with matrix
- friend void gpyr (const Mat3& m1, double& p, double& y, double& r); // get p/y/r
- friend void mxevc (const Mat3& m, double& a, Vec3& v); // eigenvector
- friend std::ostream& operator << (std::ostream& os, const Mat3& c);
-};
-
-// ****************************************************************************
-// defining the following functions here seems to make more compilers happy
- Mat3 mxcon (double r); // constant matrix
- Mat3 mxidn (); // identity matrix
- // friend Mat3 mxtrn (const Mat3& m1); // transposed matrix
- Mat3 xrot (double a); // rotation around x-axis
- Mat3 yrot (double a); // rotation around y-axis
- Mat3 zrot (double a); // rotation around z-axis
-
- Mat3 csmx (double p, double y, double r); // pitch/yaw/roll matrix
- void vcpy (Vec3& v, double& p, double& y); // get pitch and yaw from vector
- void vcrp (Vec3& v, double& p, double& r); // get pitch and roll from vector
- Mat3 mxrox (double& a, Vec3& v); // get matrix from eigenvector and angle
-
-#endif // __attlib_h sentry.
-
diff --git a/src/plugins/render/satellites/mex/planetarySats.cpp \
b/src/plugins/render/satellites/mex/planetarySats.cpp deleted file mode 100644
index ec31a20..0000000
--- a/src/plugins/render/satellites/mex/planetarySats.cpp
+++ /dev/null
@@ -1,727 +0,0 @@
-//
-// This file is part of the Marble Virtual Globe.
-//
-// This program is free software licensed under the GNU LGPL. You can
-// find a copy of this license in LICENSE.txt in the top directory of
-// the source code.
-//
-// Copyright 2012 Gerhard HOLTKAMP
-//
-
-/***************************************************************************
-* Calculate Spacecraft around other planets *
-* *
-* *
-* Open Source Code. License: GNU LGPL Version 2+ *
-* *
-* Author: Gerhard HOLTKAMP, 26-AUG-2012 *
-***************************************************************************/
-
-/*------------ include files and definitions -----------------------------*/
-#include <fstream>
-#include <iostream>
-#include <cstdlib>
-#include <cstring>
-#include <cmath>
-#include <ctime>
-using namespace std;
-
-#include "planetarySats.h"
-#include "attlib.h"
-#include "astrolib.h"
-
-// ################ Planetary Sats Class ####################
-
-PlanetarySats::PlanetarySats()
-{
- plsatinit();
-}
-
-PlanetarySats::~PlanetarySats()
-{
-
-}
-
-double PlanetarySats::atan23 (double y, double x)
- {
- // redefine atan2 so that it does'nt crash when both x and y are 0
- double result;
-
- if ((x == 0) && (y == 0)) result = 0;
- else result = atan2 (y, x);
-
- return result;
- }
-
-void PlanetarySats::plsatinit()
-{
- // initialize planetary sat data
- pls_moonflg = false;
- pls_day = 1;
- pls_month = 1;
- pls_year = 2012;
- pls_hour = 0;
- pls_minute = 0;
- pls_second = 0;
- pls_del_auto = 1;
- pls_step = 60.0;
- pls_delta_rt = 0.0;
- getTime();
- getMars();
-
- strcpy (pls_satelmfl, "./planetarysats.txt");
-
-}
-
-void PlanetarySats::getTime () // Get System Time and Date
-{
- time_t tt;
- int hh, mm, ss;
- double jd, hr;
-
- tt = time(NULL);
- jd = 40587.0 + tt/86400.0; // seconds since 1-JAN-1970
-
- jd = jd + pls_delta_rt / 24.0;
- caldat(jd, hh, mm, ss, hr);
- pls_year = ss;
- pls_month = mm;
- pls_day = hh;
-
- dms(hr, hh, mm, jd);
- pls_hour = hh;
- pls_minute = mm;
- pls_second = int(jd);
- if (pls_del_auto) pls_del_tdut = DefTdUt(pls_year);
- setMJD(pls_year, pls_month, pls_day, hh, mm, jd);
-};
-
-void PlanetarySats::setStepWidth(double s)
-{
- // set the step width (in seconds) used for calculations
- if (s < 0.01) pls_step = 0.01;
- else pls_step = s;
-}
-
-void PlanetarySats::setDeltaRT(double drt)
-{
- pls_delta_rt = drt;
-}
-
-void PlanetarySats::setDeltaTAI_UTC(double d)
-{
- // c is the difference between TAI and UTC according to the IERS
- // we have to add 32.184 sec to get to the difference TT - UT
- // which is used in the calculations here
-
- pls_del_tdut = d + 32.184;
- pls_del_auto = 0;
-}
-
-void PlanetarySats::setAutoTAI_UTC()
-{
- // set the difference between TAI and UTC according to the IERS
- pls_del_auto = true;
- pls_del_tdut = DefTdUt(pls_year);
-}
-
-void PlanetarySats::setMJD(int year, int month, int day, int hour, int min, double \
sec)
-{
- // set the (MJD-) time currently used for calculations to year, month, day, \
hour, min, sec
- double jd;
-
- pls_year = year;
- pls_month = month;
- pls_day = day;
- pls_hour = hour;
- pls_minute = min;
- pls_second = sec;
-
- jd = ddd(hour, min, sec);
- jd = mjd(day, month, year, jd);
-
- pls_time = jd;
-
- if (pls_del_auto) pls_del_tdut = DefTdUt(pls_year);
-
-}
-
-void PlanetarySats::getDatefromMJD(double mjd, int &year, int &month, int &day, int \
&hour, int &min, double &sec)
-{
- // convert times given in Modified Julian Date (MJD) into conventional date and \
time
-
- double magn;
-
- caldat((mjd), day, month, year, magn);
- dms (magn,hour,min,sec);
- if (sec>30.0) min++;;
- if (min>59)
- {
- hour++;
- min=0;
- };
-}
-
-void PlanetarySats::setSatFile(char* fname)
-{
- strcpy (pls_satelmfl, fname);
-
-}
-
-void PlanetarySats::setStateVector(double mjd, double x, double y, double z, double \
vx, double vy, double vz)
-{
- pls_rep[0] = x;
- pls_rep[1] = y;
- pls_rep[2] = z;
- pls_vep[0] = vx;
- pls_vep[1] = vy;
- pls_vep[2] = vz;
-
- int year = 0;
- int month = 0;
- int day = 0;
- int hour = 0;
- int min = 0;
- double sec = 0;
- getDatefromMJD(mjd, year, month, day, hour, min, sec);
- setMJD(year, month, day, hour, min, sec);
- pls_tepoch = pls_time;
- //pls_tepoch = pls_time + pls_del_tdut / 86400.0; // epoch in TT
-}
-
-int PlanetarySats::getStateVector(int nsat)
-{
- // read the state vector from the planetary sat file
- // nsat = number of eligible sat to select (1 if first or only sat)
- // RETURN number of eligible sat selected, 0 if no suitable sats of file problems
-
- int fsc, j, k, nst, utc;
- int yr, month, day, hour, min;
- double sec;
- bool searching;
- ifstream cfle;
- char satname[40]; // name of satellite
- char plntname[40]; // name of planet
-
- strcpy(satname, "");
- strcpy(plntname, "");
- nst = 0;
-
- searching = true;
- fsc = 0;
-
- cfle.open(pls_satelmfl, ios::in);
- if (!cfle)
- {
- searching = false;
- cfle.close();
- };
- if(searching)
- {
- while (searching)
- {
- fsc = 1;
-
- if (!cfle.getline(satname,40)) fsc = 0;
- else
- {
- k = strlen(satname);
- if ((k>1) && (satname[0] == '#'))
- {
- for (j=1; j<k; ++j)
- {
- pls_satname[j-1] = satname[j];
- if(pls_satname[j-1] == '\n') pls_satname[j-1] = '\0';
- };
- pls_satname[k-1] = '\0';
- }
- else fsc = 0;
- };
-
- if (cfle.eof())
- {
- fsc = 0;
- searching = false;
- };
-
- if (fsc)
- {
- if (!cfle.getline(plntname, 40)) fsc = 0;
- else
- {
- k = strlen(plntname);
- if (k>0)
- {
- if(plntname[k-1] == '\n') plntname[k-1] = '\0';
- if((k>1) && (plntname[k-2] == '\r')) plntname[k-2] = '\0';
- };
- };
- };
-
- if (fsc)
- {
- cfle >> yr >> month >> day >> hour >> min >> sec >> utc;
- if (cfle.bad()) fsc = 0;
- if (cfle.eof())
- {
- fsc = 0;
- searching = false;
- };
- };
-
- if (fsc)
- {
- cfle >> pls_rep[0] >> pls_rep[1] >> pls_rep[2];
- if (cfle.bad()) fsc = 0;
- if (cfle.eof())
- {
- fsc = 0;
- searching = false;
- };
- };
-
- if (fsc)
- {
- cfle >> pls_vep[0] >> pls_vep[1] >> pls_vep[2];
- if (cfle.bad()) fsc = 0;
- };
-
- if (fsc)
- {
- if (strncmp(pls_plntname, plntname, 4) == 0) nst++;
- if (nst == nsat)
- {
- searching = false;
-
- setMJD(yr, month, day, hour, min, sec);
- pls_tepoch = pls_time;
- if (utc) pls_tepoch = pls_tepoch + pls_del_tdut/86400.0; // epoch in TT
- };
- };
- };
- cfle.close();
- };
-
- if (fsc == 0) nst = 0;
-
- return nst;
-}
-
-void PlanetarySats::setPlanet(char* pname)
-{
- pls_moonflg = false;
- strcpy(pls_plntname, pname);
- if (strncmp("Mars", pname, 4) == 0) getMars();
- if (strncmp("Venus", pname, 4) == 0) getVenus();
- if (strncmp("Mercury", pname, 4) == 0) getMercury();
- if (strncmp("Moon", pname, 4) == 0) getMoon();
-}
-
-void PlanetarySats::stateToKepler()
-{
- // convert state vector (mean equatorial J2000.0) into planetary Kepler elements
-
- double t, dt, ag, gm, re, j2;
- double n, c, w, a, ecc, inc;
- Vec3 r1, v1;
- Mat3 mx;
-
- dt = (pls_tepoch - 51544.5) / 36525.0;
- gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
- re = pls_R0;
- j2 = pls_J2;
-
- // convert into planet equatorial reference frame
- if (pls_moonflg)
- {
- mx = mxidn();
- r1 = mxvct (mx, pls_rep);
- v1 = mxvct (mx, pls_vep);
-
- }
- else
- {
- ag = (pls_axl0 + pls_axl1 * dt) * M_PI / 180.0;
- mx = zrot(ag + M_PI / 2.0);
- r1 = mxvct (mx, pls_rep);
- v1 = mxvct (mx, pls_vep);
-
- ag = (pls_axb0 + pls_axb1 * dt) * M_PI / 180.0;
- mx = xrot(M_PI / 2.0 - ag);
- r1 = mxvct (mx, r1);
- v1 = mxvct (mx, v1);
- };
-
- v1 *= 86400.0; // convert into km / day
-
- oscelm(gm, pls_tepoch, r1, v1, t, pls_m0, pls_a, pls_ecc, pls_ra, pls_per, \
pls_inc);
-
- // now the mean motion
- a = pls_a;
- ecc = pls_ecc;
- inc = pls_inc;
-
- // preliminary n
- if (a == 0) a = 1.0; // just in case
- if (a < 0) a = -a; // just in case
- n = sqrt (gm / (a*a*a));
-
- // correct for J2 term
- w = 1.0 - ecc*ecc;
- if (w > 0)
- {
- w = pow (w, -1.5);
- c = sin (inc*M_PI/180.0);
- n = n*(1.0 + 1.5*j2*re*re/(a*a) * w * (1.0 - 1.5*c*c));
- }
- else n = 1.0; // do something to avoid a domain error
-
- n = n / (2.0*M_PI);
- if (n > 1000.0) n = 1000.0; // avoid possible errors
-
- pls_n0 = n;
-
-}
-
-void PlanetarySats::getKeplerElements(double &perc, double &apoc, double &inc, \
double &ecc, double &ra, double &tano, double &m0, double &a, double \
&n0)
-{
- // get Kepler elements of orbit with regard to the planetary equator and prime \
meridian.
-
- double t, gm;
- Vec3 r1, v1;
- Mat3 mx;
-
-
- if (pls_moonflg) // for the Moon we normally work in J2000. Now get it into \
planetary
- {
- gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
-
- mx = getSelenographic(pls_tepoch);
- r1 = mxvct (mx, pls_rep);
- v1 = mxvct (mx, pls_vep);
- v1 *= 86400.0; // convert into km / day
-
- oscelm(gm, pls_tepoch, r1, v1, t, m0, a, ecc, ra, tano, inc);
-
- // now the mean motion
- if (a == 0) a = 1.0; // just in case
- if (a < 0) a = -a; // just in case
- n0 = sqrt (gm / (a*a*a));
- n0 = n0 / (2.0*M_PI);
- }
- else
- {
- a = pls_a;
- n0 = pls_n0;
- m0 = pls_m0;
- tano = pls_per;
- ra = pls_ra;
- ecc = pls_ecc;
- inc = pls_inc;
- };
-
- perc = pls_a * (1.0 - pls_ecc) - pls_R0;
- apoc = pls_a * (1.0 + pls_ecc) - pls_R0;
-
-}
-
-int PlanetarySats::selectSat(char* sname)
-{
- // select specified satellite
- // RETURN 1 if successful, 0 if no suitable satellite found
-
- int nst, res, sl;
- bool searching;
-
- searching = true;
- nst = 1;
- sl = strlen(sname);
-
- while (searching)
- {
- res = getStateVector(nst);
- if (res)
- {
- if (strncmp(pls_satname, sname, sl) == 0) searching = false;
- }
- else searching = false;
- nst++;
- };
-
- return res;
-}
-
-void PlanetarySats::getSatName(char* sname)
-{
- strcpy (sname, pls_satname);
-}
-
-void PlanetarySats::currentPos()
-{
- getSatPos(pls_time);
-}
-
-void PlanetarySats::nextStep()
-{
- pls_time = pls_time + pls_step / 86400.0;
- getSatPos(pls_time);
-}
-
-double PlanetarySats::getLastMJD()
-{
- return pls_time;
-}
-
-void PlanetarySats::getPlanetographic(double &lng, double &lat, double &height)
-{
- // planetographic coordinates from current state vector
-
- lng = pls_lng;
- lat = pls_lat;
- height = pls_height;
-
-}
-
-void PlanetarySats::getFixedFrame(double &x, double &y, double &z, double &vx, \
double &vy, double &vz)
-{
- // last state vector coordinates in planetary fixed frame
-
- x = pls_r[0];
- y = pls_r[1];
- z = pls_r[2];
- vx = pls_v[0];
- vy = pls_v[1];
- vz = pls_v[2];
-}
-
-void PlanetarySats::getSatPos (double tutc)
-{
- // Get Position of Satellite at MJD-time t (UTC)
-
- const double mp2 = 2.0*M_PI;
-
- double t, dt, m0, ran, aper, inc, a, ecc, n0, re;
- double f, e, c, s, k, sh, j2, gm, fac, b1, b2, b3;
- int j;
-
- Vec3 r1, v1, rg1, s2;
- Mat3 m1, m2;
-
- // prepare orbit calculation
-
- t = tutc + pls_del_tdut/86400.0;
- dt = t - pls_tepoch;
-
- ecc = pls_ecc;
- if (ecc >= 1.0) ecc = 0.999; // to avoid crashes
- a = pls_a;
- n0 = mp2 * pls_n0;
- if (a < 1.0) a = 1.0; // avoid possible crashes later on
- re = pls_R0;
- f = pls_flat;
- j2 = pls_J2;
- gm = pls_GM * 7.4649600000; // convert from m^3/s^2 into km^3/d^2
-
- // get current orbit elements ready
- aper = (re / a) / (1.0 - ecc*ecc);
- aper = 1.5 * j2 * aper*aper * n0;
- m0 = pls_inc*M_PI/180.0;
- ran = -aper * cos(m0) * dt;
- m0 = sin (m0);
- aper = aper * (2.0 - 2.5*m0*m0) * dt;
- ran = pls_ra * M_PI/180.0 + ran;
- aper = pls_per * M_PI/180.0 + aper;
- m0 = pls_m0*M_PI/180.0 + n0*dt;
- inc = pls_inc * M_PI/180.0;
-
- // solve Kepler equation
- if (a < 1.0) a = 1.0; // avoid possible crashes later on
- e = eccanom(m0, ecc);
- fac = sqrt(1.0 - ecc*ecc);
- c = cos(e);
- s = sin(e);
- r1.assign (a*(c - ecc), a*fac*s, 0.0);
-
- m0 = 1.0 - ecc*c;
- k = sqrt(gm/a);
- v1.assign (-k*s/m0, k*fac*c/m0, 0.0);
-
- // convert into reference plane
- m1 = zrot (-aper);
- m2 = xrot (-inc);
- m1 *= m2;
- m2 = zrot (-ran);
- m2 = m2 * m1;
- r1 = mxvct (m2, r1);
- v1 = mxvct (m2, v1);
- v1 /= 86400.0;
-
- // save state vector in planet-fixed frame
-
- if (pls_moonflg) m1 = getSelenographic(t);
- else m1 = zrot((pls_W + pls_Wd*(t-51544.5))*(M_PI/180.0));
- pls_r = mxvct(m1,r1);
- pls_v = mxvct(m1,v1);
- pls_r *= 1000.0;
- pls_v *= 1000.0;
-
- // get the groundtrack coordinates
-
- rg1 = mxvct(m1,r1);
-
- s2 = carpol(rg1);
- pls_lat = s2[2]; // just preliminary
- pls_lng = s2[1]; // measured with the motion of rotation!
- if (pls_lng > mp2) pls_lng -= mp2;
- if (pls_lng < -M_PI) pls_lng += mp2;
- if (pls_lng > M_PI) pls_lng -= mp2;
-
- // get height above ellipsoid and geodetic latitude
- if (abs(r1) > 0.1)
- {
- if (f == 0) pls_height = abs(r1) - re;
- else
- {
- c = f*(2.0 - f);
- s = r1[0]*r1[0] + r1[1]*r1[1];
- sh = c*r1[2];
-
- for (j=0; j<4; ++j)
- {
- b1 = r1[2] + sh;
- b3 = sqrt(s+b1*b1);
- if (b3 < 1e-5) b3 = sin(pls_lat); // just in case
- else b3 = b1 / b3;
- b2 = re / sqrt(1.0 - c*b3*b3);
- sh = b2 * c * b3;
- };
-
- sh = r1[2] + sh;
- pls_lat = atan20(sh,sqrt(s));
- sh = sqrt(s+sh*sh) - b2;
- pls_height = sh;
- }
- }
- else pls_height = 0; // this should never happen
-
- pls_lat = pls_lat * 180.0 / M_PI;
- pls_lng = pls_lng * 180.0 / M_PI;
-
- }
-
-
-void PlanetarySats::getMars() // Mars planetary constants
-{
- pls_J2 = 1.964e-3;
- pls_R0 = 3397.2;
- pls_flat = 0.00647630;
- pls_axl0 = 317.681;
- pls_axl1 = -0.108;
- pls_axb0 = 52.886;
- pls_axb1 = -0.061;
- pls_W = 176.868;
- pls_Wd = 350.8919830;
- pls_GM = 4.282828596416e+13; // 4.282837405582e+13
-}
-
-void PlanetarySats::getVenus() // Venus planetary constants
-{
- pls_J2 = 0.027e-3;
- pls_R0 = 6051.9;
- pls_flat = 0.0;
- pls_axl0 = 272.72;
- pls_axl1 = 0.0;
- pls_axb0 = 67.15;
- pls_axb1 = 0.0;
- pls_W = 160.26;
- pls_Wd = -1.4813596;
- pls_GM = 3.24858761e+14;
-}
-
-void PlanetarySats::getMercury() // Mercury planetary constants
-{
- pls_J2 = 0.0;
- pls_R0 = 2439.7;
- pls_flat = 0.0;
- pls_axl0 = 281.01;
- pls_axl1 = -0.033;
- pls_axb0 = 61.45;
- pls_axb1 = -0.005;
- pls_W = 329.71;
- pls_Wd = 6.1385025;
- pls_GM = 2.20320802e+13;
-}
-
-void PlanetarySats::getMoon() // Moon planetary constants
-{
- pls_moonflg = true;
- pls_J2 = 0.2027e-3;
- pls_R0 = 1738.0;
- pls_flat = 0.0;
- pls_axl0 = 0.0;
- pls_axl1 = 0.0;
- pls_axb0 = 90.0;
- pls_axb1 = 0.0;
- pls_W = 0.0;
- pls_Wd = 13.17635898;
- pls_GM = 4.90279412e+12;
-}
-
-Mat3 PlanetarySats::getSelenographic (double jd)
-{
- // Calculate the Matrix to transform from Mean of J2000 into selenographic
- // coordinates at MJD time jd.
-
- double t, gam, gmp, l, omg, mln;
- double a, b, c, ic, gn, gp, omp;
- double const degrad = M_PI / 180.0;
- Vec3 v1;
- Mat3 m1, m2;
-
- t = (jd - 15019.5) / 36525.0;
- gam = 281.2208333 + ((0.33333333e-5*t + 0.45277778e-3)*t + 1.7191750)*t;
- gmp = 334.3295556 + ((-0.125e-4*t - 0.10325e-1)*t + 4069.0340333)*t;
- l = 279.6966778 + (0.3025e-3*t + 36000.768925)*t;
- omg = 259.1832750 + ((0.22222222e-5*t + 0.20777778e-2)*t - 1934.1420083)*t;
- b = 23.45229444 + ((0.50277778e-6*t -0.16388889e-5)*t - 0.130125e-1)*t;
- mln = 270.4341639 + ((0.1888889e-5*t -0.11333e-2)*t + 481267.8831417)*t;
- ic = 1.535*degrad;
- gn = (l - gam)*degrad;
- gp = (mln - gmp)*degrad;
- omp = (gmp - omg)*degrad;
- a = -107.0*cos(gp) + 37.0*cos(gp+2.0*omp) -11.0*cos(2.0*(gp+omp));
- a = a*0.000277778*degrad + ic;
- c = (-109.0*sin(gp) + 37.0*sin(gp+2.0*omp) - 11.0*sin(2.0*(gp+omp)))/sin(ic);
- c = (c*0.000277778 + omg)*degrad;
- gn = -12.0*sin(gp) + 59.0*sin(gn) + 18.0*sin(2.0*omp);
- gn = gn*0.000277778*degrad; // tau
-
- b *= degrad;
- gam = cos(a)*cos(b) + sin(a)*sin(b)*cos(c);
- gmp = gam*gam;
- if(gmp > 1.0) gmp = 0;
- else gmp = sqrt(1.0 - gmp);
- gam = atan23(gmp,gam); // theta
- t = cos(a)*sin(b) - sin(a)*sin(b)*cos(c);
- l = -sin(a)*sin(c);
-
- gmp = atan23(l,t); // phi
- t = sin(a)*cos(b) - cos(a)*sin(b)*cos(c);
- l = -sin(b)*sin(c);
- a = atan23(l,t); // delta
- c = a + mln*degrad + gn - c; // psi
-
- // libration rotation matrix from Mean equator to true selenographic
- m1 = zrot(gmp);
- m2 = xrot(gam);
- m1 = m2 * m1;
- m2 = zrot(c);
- m1 = m2 * m1;
-
- t = julcent(jd);
- m2 = pmatequ(0,t); // convert from mean of J2000 to mean of epoch
- m1 = m1 * m2;
-
- return m1;
-}
-
diff --git a/src/plugins/render/satellites/mex/planetarySats.h \
b/src/plugins/render/satellites/mex/planetarySats.h deleted file mode 100644
index 0ed34c7..0000000
--- a/src/plugins/render/satellites/mex/planetarySats.h
+++ /dev/null
@@ -1,104 +0,0 @@
-//
-// This file is part of the Marble Virtual Globe.
-//
-// This program is free software licensed under the GNU LGPL. You can
-// find a copy of this license in LICENSE.txt in the top directory of
-// the source code.
-//
-// Copyright 2012 Gerhard HOLTKAMP
-//
-
-#if !defined(__planetarysats_h)
-#define __planetarysats_h
-
-#include "attlib.h"
-
-class PlanetarySats // Calculate spacecraft around other planets
-{
- public:
- PlanetarySats();
- ~PlanetarySats();
-
- void setStepWidth(double s); // set the step width (seconds)
- void setDeltaTAI_UTC(double d); // set IERS Parameter TAI - UTC
- void setAutoTAI_UTC(); // IERS Parameter TAI - UTC to auto
- void getTime (); // Get System Time and Date
- void setDeltaRT(double drt);
- void setMJD(int year, int month, int day, int hour, int min, double sec);// set \
time
- void getDatefromMJD(double mjd, int &year, int &month, int &day, int &hour, int \
&min, double &sec);
- void setSatFile(char* fname);
- void setPlanet(char* pname);
- int selectSat(char* sname);
- void getSatName(char* sname);
- void setStateVector(double mjd, double x, double y, double z, double vx, double \
vy, double vz);
- int getStateVector(int nsat);
- void stateToKepler();
- void getKeplerElements(double &perc, double &apoc, double &inc, double &ecc, \
double &ra, double &tano, double &m0, double &a, double &n0);
- void getPlanetographic(double &lng, double &lat, double &height);
- void getFixedFrame(double &x, double &y, double &z, double &vx, double &vy, \
double &vz);
- void currentPos();
- void nextStep();
- double getLastMJD();
-
- private:
- void plsatinit(); // initialize PlanetarySats
- double atan23 (double y, double x); // atan without singularity for x,y=0
- void getMercury();
- void getVenus();
- void getMoon();
- void getMars();
- void getSatPos (double t);
- Mat3 getSelenographic (double jd);
-
- // data fields
-
- char pls_satelmfl[205]; // name of file for satellite state vectors
- char pls_satname[40]; // name of satellite
- char pls_plntname[40]; // name of planet
-
- bool pls_moonflg; // true if Moon, false if other body
-
- int pls_day; // date
- int pls_month;
- int pls_year;
- int pls_hour;
- int pls_minute;
- int pls_second;
- double pls_time; // current time in MJD (UTC)
- double pls_del_tdut; // TDT - UT in sec
- int pls_del_auto; // 1 = automatic del_tdut, 0 = manual
- double pls_step; // stepwidth in sec
- double pls_delta_rt; // delta time to R/T in hours
-
- double pls_tepoch; // MJD epoch of state vector (TT)
- Vec3 pls_rep; // state vector km and km/s
- Vec3 pls_vep;
-
- double pls_ra; // Right Ascension
- double pls_per; // argument of perige center
- double pls_m0; // mean anomaly
- double pls_ecc; // eccentricity
- double pls_inc; // inclination
- double pls_a; // semi-major axis (km)
- double pls_n0; // mean motion (rev/day)
-
- double pls_GM; // graviatational constant (m^3/s^2)
- double pls_J2; // J2 gravitational term
- double pls_R0; // equatorial radius (km)
- double pls_flat; // flattening factor
- double pls_axl0; // l-direction of rotation axis
- double pls_axl1; // delta of axl0
- double pls_axb0; // b-direction of rotation axis
- double pls_axb1; // delta of axb0
- double pls_W; // location of prime meridian
- double pls_Wd; // daily variation of W.
-
- Vec3 pls_r; // current state vector m and m/s
- Vec3 pls_v;
- double pls_lat; // planetary latitude (decimal degrees)
- double pls_lng; // planetary longitude (decimal degrees)
- double pls_height; // height above reference ellipsoid (km)
-};
-
-#endif // __planetarysats_h sentry.
-
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic