Implemented CustomRotation for Earth using the long period extension to the
P03 precession model.ver1_6_1
Chris Laurel
parent
a045b83ab1
commit
5b7c56b522
|
|
@ -15,6 +15,7 @@
|
|||
#include <celengine/customrotation.h>
|
||||
#include <celengine/rotation.h>
|
||||
#include <celengine/astro.h>
|
||||
#include <celengine/precession.h>
|
||||
|
||||
using namespace std;
|
||||
|
||||
|
|
@ -69,6 +70,67 @@ private:
|
|||
};
|
||||
|
||||
|
||||
|
||||
/******* Earth rotation model *******/
|
||||
|
||||
class EarthRotationModel : public CachingRotationModel
|
||||
{
|
||||
public:
|
||||
EarthRotationModel()
|
||||
{
|
||||
}
|
||||
|
||||
~EarthRotationModel()
|
||||
{
|
||||
}
|
||||
|
||||
Quatd computeSpin(double tjd) const
|
||||
{
|
||||
// TODO: Use a more accurate model for sidereal time
|
||||
double t = tjd - astro::J2000;
|
||||
double theta = 2 * PI * (t * 24.0 / 23.9344694 - 280.5 / 360.0);
|
||||
|
||||
return Quatd::yrotation(-theta);
|
||||
}
|
||||
|
||||
Quatd computeEquatorOrientation(double tjd) const
|
||||
{
|
||||
double T = (tjd - astro::J2000) / 36525.0;
|
||||
astro::PrecessionAngles prec = astro::PrecObliquity_P03LP(T);
|
||||
astro::EclipticPole pole = astro::EclipticPrecession_P03LP(T);
|
||||
|
||||
double obliquity = degToRad(prec.epsA / 3600);
|
||||
double precession = degToRad(prec.pA / 3600);
|
||||
|
||||
// Calculate the angles pi and Pi from the ecliptic pole coordinates
|
||||
// P and Q:
|
||||
// P = sin(pi)*sin(Pi)
|
||||
// Q = sin(pi)*cos(Pi)
|
||||
double P = pole.PA * 2.0 * PI / 1296000;
|
||||
double Q = pole.QA * 2.0 * PI / 1296000;
|
||||
double piA = asin(sqrt(P * P + Q * Q));
|
||||
double PiA = atan2(P, Q);
|
||||
|
||||
// Calculate the rotation from the J2000 ecliptic to the ecliptic
|
||||
// of date.
|
||||
Quatd RPi = Quatd::zrotation(PiA);
|
||||
Quatd rpi = Quatd::xrotation(piA);
|
||||
Quatd eclRotation = ~RPi * rpi * RPi;
|
||||
|
||||
Quatd q = Quatd::xrotation(obliquity) * Quatd::zrotation(-precession) * ~eclRotation;
|
||||
|
||||
// convert to Celestia's coordinate system
|
||||
return Quatd::xrotation(PI / 2.0) * q * Quatd::xrotation(-PI / 2.0);
|
||||
}
|
||||
|
||||
double getPeriod() const
|
||||
{
|
||||
return 24.0 / 23.9344694;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
/******* IAU rotation models for the planets *******/
|
||||
|
||||
|
||||
|
|
@ -477,6 +539,8 @@ GetCustomRotationModel(const std::string& name)
|
|||
{
|
||||
CustomRotationModelsInitialized = true;
|
||||
|
||||
CustomRotationModels["earth-p03lp"] = new EarthRotationModel();
|
||||
|
||||
CustomRotationModels["iau-mercury"] = new IAUPrecessingRotationModel(281.01, -0.033,
|
||||
61.45, -0.005,
|
||||
329.548, 6.1385025);
|
||||
|
|
|
|||
|
|
@ -1,185 +1,378 @@
|
|||
// precession.cpp
|
||||
//
|
||||
// Calculate precession angles for Earth.
|
||||
//
|
||||
// Copyright (C) 2008, the Celestia Development Team
|
||||
// Initial version by Chris Laurel, claurel@gmail.com
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License
|
||||
// as published by the Free Software Foundation; either version 2
|
||||
// of the License, or (at your option) any later version.
|
||||
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <celmath/mathlib.h>
|
||||
#include <celmath/vecmath.h>
|
||||
#include "precession.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
|
||||
// Periodic term for the long-period extension of the P03 precession
|
||||
// model.
|
||||
struct EclipticPrecessionTerm
|
||||
{
|
||||
double Pc;
|
||||
double Qc;
|
||||
double Ps;
|
||||
double Qs;
|
||||
double period;
|
||||
};
|
||||
|
||||
|
||||
static EclipticPrecessionTerm EclipticPrecessionTerms[] =
|
||||
{
|
||||
{ 486.230527, 2559.065245, -2578.462809, 485.116645, 2308.98 },
|
||||
{ -963.825784, 247.582718, -237.405076, -971.375498, 1831.25 },
|
||||
{ -1868.737098, -957.399054, 1007.593090, -1930.464338, 687.52 },
|
||||
{ -1589.172175, 493.021354, -423.035168, -1634.905683, 729.97 },
|
||||
{ 429.442489, -328.301413, 337.266785, 429.594383, 492.21 },
|
||||
{ -2244.742029, -339.969833, 221.240093, -2131.745072, 708.13 },
|
||||
};
|
||||
|
||||
|
||||
// Periodic term for the long-period extension of the P03 precession
|
||||
// model.
|
||||
struct PrecessionTerm
|
||||
{
|
||||
double pc;
|
||||
double epsc;
|
||||
double ps;
|
||||
double epss;
|
||||
double period;
|
||||
};
|
||||
|
||||
|
||||
static PrecessionTerm PrecessionTerms[] =
|
||||
{
|
||||
{ -6180.062400, 807.904635, -2434.845716, -2056.455197, 409.90 },
|
||||
{ -2721.869299, -177.959383, 538.034071, -912.727303, 396.15 },
|
||||
{ 1460.746498, 371.942696, -1245.689351, 447.710000, 536.91 },
|
||||
{ -1838.488899, -176.029134, 529.220775, -611.297411, 402.90 },
|
||||
{ 949.518077, -89.154030, 277.195375, 315.900626, 417.15 },
|
||||
{ 32.701460, -336.048179, 945.979710, 12.390157, 288.92 },
|
||||
{ 598.054819, -17.415730, -955.163661, -15.922155, 4042.97 },
|
||||
{ -293.145284, -28.084479, 93.894079, -102.870153, 304.90 },
|
||||
{ 66.354942, 21.456146, 0.671968, 24.123484, 281.46 },
|
||||
{ 18.894136, 30.917011, -184.663935, 2.512708, 204.38 },
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
/*! Compute the precession of the ecliptic, based on a long-period
|
||||
* extension of the the P03 model, presented in "Long-periodic Precession
|
||||
* Parameters", J. Vondrak (2006)
|
||||
* http://www.astronomy2006.com/files/vondrak.pdf
|
||||
* For an explanation of the angles used in the P03 model, see
|
||||
* "Expressions for IAU2000 precession quantities", N. Capitaine et al,
|
||||
* Astronomy & Astrophysics, v.412, p.567-586 (2003).
|
||||
*
|
||||
* 6 long-periodic terms, plus a cubic polynomial for longer terms.
|
||||
* The terms are fitted to the P03 model withing 1000 years of J2000.
|
||||
*
|
||||
* T is the time in centuries since J2000. The angles returned are
|
||||
* in arcseconds.
|
||||
*/
|
||||
astro::EclipticAngles
|
||||
astro::EclipticPrecession_P03LP(double T)
|
||||
{
|
||||
EclipticAngles angles;
|
||||
|
||||
double T2 = T * T;
|
||||
double T3 = T2 * T;
|
||||
|
||||
angles.PA = (5750.804069
|
||||
+ 0.1948311 * T
|
||||
- 0.00016739 * T2
|
||||
- 4.8e-8 * T3);
|
||||
angles.QA = (-1673.999018
|
||||
+ 0.3474459 * T
|
||||
+ 0.00011243 * T2
|
||||
- 6.4e-8 * T3);
|
||||
|
||||
unsigned int nTerms = sizeof(EclipticPrecessionTerms) / sizeof(EclipticPrecessionTerms[0]);
|
||||
for (unsigned int i = 0; i < nTerms; i++)
|
||||
{
|
||||
const EclipticPrecessionTerm& p = EclipticPrecessionTerms[i];
|
||||
double theta = 2.0 * PI * T / p.period;
|
||||
double s = sin(theta);
|
||||
double c = cos(theta);
|
||||
angles.PA += p.Pc * c + p.Ps * s;
|
||||
angles.QA += p.Qc * c + p.Qs * s;
|
||||
}
|
||||
|
||||
return angles;
|
||||
}
|
||||
|
||||
|
||||
/*! Compute the general precession and obliquity, based on the model
|
||||
* presented in "Long-periodic Precession Parameters", J. Vondrak, 2006
|
||||
* http://www.astronomy2006.com/files/vondrak.pdf
|
||||
*
|
||||
* T is the time in centuries since J2000. The angles returned are
|
||||
* in arcseconds.
|
||||
*/
|
||||
astro::PrecessionAngles
|
||||
astro::PrecObliquity_P03LP(double T)
|
||||
{
|
||||
PrecessionAngles angles;
|
||||
|
||||
double T2 = T * T;
|
||||
double T3 = T2 * T;
|
||||
|
||||
angles.pA = ( 7907.295950
|
||||
+ 5044.374034 * T
|
||||
- 0.00713473 * T2
|
||||
+ 6e-9 * T3);
|
||||
angles.epsA = ( 83973.876448
|
||||
- 0.0425899 * T
|
||||
- 0.00000113 * T2);
|
||||
|
||||
unsigned int nTerms = sizeof(PrecessionTerms) / sizeof(PrecessionTerms[0]);
|
||||
for (unsigned int i = 0; i < nTerms; i++)
|
||||
{
|
||||
const PrecessionTerm& p = PrecessionTerms[i];
|
||||
double theta = 2.0 * PI * T / p.period;
|
||||
double s = sin(theta);
|
||||
double c = cos(theta);
|
||||
angles.pA += p.pc * c + p.ps * s;
|
||||
angles.epsA += p.epsc * c + p.epss * s;
|
||||
}
|
||||
|
||||
return angles;
|
||||
}
|
||||
|
||||
|
||||
|
||||
#define TEST 0
|
||||
|
||||
#if TEST
|
||||
|
||||
int main(int argc, char* argv[])
|
||||
{
|
||||
for (int i = -100; i <= 100; i++)
|
||||
{
|
||||
// Get time in Julian centuries from J2000
|
||||
double T = (double) i * 100.0;
|
||||
|
||||
astro::EclipticAngles ecl = astro::EclipticPrecession_P03LP(T);
|
||||
|
||||
Vec3d v(0.0, 0.0, 0.0);
|
||||
v.x = ecl.PA * 2.0 * PI / 1296000;
|
||||
v.y = -ecl.QA * 2.0 * PI / 1296000;
|
||||
v.z = sqrt(1.0 - v.x * v.x - v.y * v.y);
|
||||
|
||||
// Show the angle between the J2000 ecliptic pole and the ecliptic
|
||||
// pole at T centuries from J2000
|
||||
clog << i << ": " << radToDeg(acos(v.z)) << endl;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
#endif // TEST
|
||||
// precession.cpp
|
||||
//
|
||||
// Calculate precession angles for Earth.
|
||||
//
|
||||
// Copyright (C) 2008, the Celestia Development Team
|
||||
// Initial version by Chris Laurel, claurel@gmail.com
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License
|
||||
// as published by the Free Software Foundation; either version 2
|
||||
// of the License, or (at your option) any later version.
|
||||
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
#include <celmath/mathlib.h>
|
||||
#include <celmath/vecmath.h>
|
||||
#include "precession.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
|
||||
// Periodic term for the long-period extension of the P03 precession
|
||||
// model.
|
||||
struct EclipticPrecessionTerm
|
||||
{
|
||||
double Pc;
|
||||
double Qc;
|
||||
double Ps;
|
||||
double Qs;
|
||||
double period;
|
||||
};
|
||||
|
||||
|
||||
static EclipticPrecessionTerm EclipticPrecessionTerms[] =
|
||||
{
|
||||
{ 486.230527, 2559.065245, -2578.462809, 485.116645, 2308.98 },
|
||||
{ -963.825784, 247.582718, -237.405076, -971.375498, 1831.25 },
|
||||
{ -1868.737098, -957.399054, 1007.593090, -1930.464338, 687.52 },
|
||||
{ -1589.172175, 493.021354, -423.035168, -1634.905683, 729.97 },
|
||||
{ 429.442489, -328.301413, 337.266785, 429.594383, 492.21 },
|
||||
{ -2244.742029, -339.969833, 221.240093, -2131.745072, 708.13 },
|
||||
};
|
||||
|
||||
|
||||
// Periodic term for the long-period extension of the P03 precession
|
||||
// model.
|
||||
struct PrecessionTerm
|
||||
{
|
||||
double pc;
|
||||
double epsc;
|
||||
double ps;
|
||||
double epss;
|
||||
double period;
|
||||
};
|
||||
|
||||
|
||||
static PrecessionTerm PrecessionTerms[] =
|
||||
{
|
||||
{ -6180.062400, 807.904635, -2434.845716, -2056.455197, 409.90 },
|
||||
{ -2721.869299, -177.959383, 538.034071, -912.727303, 396.15 },
|
||||
{ 1460.746498, 371.942696, -1245.689351, 447.710000, 536.91 },
|
||||
{ -1838.488899, -176.029134, 529.220775, -611.297411, 402.90 },
|
||||
{ 949.518077, -89.154030, 277.195375, 315.900626, 417.15 },
|
||||
{ 32.701460, -336.048179, 945.979710, 12.390157, 288.92 },
|
||||
{ 598.054819, -17.415730, -955.163661, -15.922155, 4042.97 },
|
||||
{ -293.145284, -28.084479, 93.894079, -102.870153, 304.90 },
|
||||
{ 66.354942, 21.456146, 0.671968, 24.123484, 281.46 },
|
||||
{ 18.894136, 30.917011, -184.663935, 2.512708, 204.38 },
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
/*! Compute the precession of the ecliptic, based on a long-period
|
||||
* extension of the the P03 model, presented in "Long-periodic Precession
|
||||
* Parameters", J. Vondrak (2006)
|
||||
* http://www.astronomy2006.com/files/vondrak.pdf
|
||||
* For an explanation of the angles used in the P03 model, see
|
||||
* "Expressions for IAU2000 precession quantities", N. Capitaine et al,
|
||||
* Astronomy & Astrophysics, v.412, p.567-586 (2003).
|
||||
*
|
||||
* Also: "Expressions for the Precession Quantities", J. H. Lieske et al,
|
||||
* Astronomy & Astrophysics, v.58, p. 1-16 (1977).
|
||||
*
|
||||
* 6 long-periodic terms, plus a cubic polynomial for longer terms.
|
||||
* The terms are fitted to the P03 model withing 1000 years of J2000.
|
||||
*
|
||||
* T is the time in centuries since J2000. The angles returned are
|
||||
* in arcseconds.
|
||||
*/
|
||||
astro::EclipticPole
|
||||
astro::EclipticPrecession_P03LP(double T)
|
||||
{
|
||||
EclipticPole pole;
|
||||
|
||||
double T2 = T * T;
|
||||
double T3 = T2 * T;
|
||||
|
||||
pole.PA = (5750.804069
|
||||
+ 0.1948311 * T
|
||||
- 0.00016739 * T2
|
||||
- 4.8e-8 * T3);
|
||||
pole.QA = (-1673.999018
|
||||
+ 0.3474459 * T
|
||||
+ 0.00011243 * T2
|
||||
- 6.4e-8 * T3);
|
||||
|
||||
unsigned int nTerms = sizeof(EclipticPrecessionTerms) / sizeof(EclipticPrecessionTerms[0]);
|
||||
for (unsigned int i = 0; i < nTerms; i++)
|
||||
{
|
||||
const EclipticPrecessionTerm& p = EclipticPrecessionTerms[i];
|
||||
double theta = 2.0 * PI * T / p.period;
|
||||
double s = sin(theta);
|
||||
double c = cos(theta);
|
||||
pole.PA += p.Pc * c + p.Ps * s;
|
||||
pole.QA += p.Qc * c + p.Qs * s;
|
||||
}
|
||||
|
||||
return pole;
|
||||
}
|
||||
|
||||
|
||||
/*! Compute the general precession and obliquity, based on the model
|
||||
* presented in "Long-periodic Precession Parameters", J. Vondrak, 2006
|
||||
* http://www.astronomy2006.com/files/vondrak.pdf
|
||||
*
|
||||
* T is the time in centuries since J2000. The angles returned are
|
||||
* in arcseconds.
|
||||
*/
|
||||
astro::PrecessionAngles
|
||||
astro::PrecObliquity_P03LP(double T)
|
||||
{
|
||||
PrecessionAngles angles;
|
||||
|
||||
double T2 = T * T;
|
||||
double T3 = T2 * T;
|
||||
|
||||
angles.pA = ( 7907.295950
|
||||
+ 5044.374034 * T
|
||||
- 0.00713473 * T2
|
||||
+ 6e-9 * T3);
|
||||
angles.epsA = ( 83973.876448
|
||||
- 0.0425899 * T
|
||||
- 0.00000113 * T2);
|
||||
|
||||
unsigned int nTerms = sizeof(PrecessionTerms) / sizeof(PrecessionTerms[0]);
|
||||
for (unsigned int i = 0; i < nTerms; i++)
|
||||
{
|
||||
const PrecessionTerm& p = PrecessionTerms[i];
|
||||
double theta = 2.0 * PI * T / p.period;
|
||||
double s = sin(theta);
|
||||
double c = cos(theta);
|
||||
angles.pA += p.pc * c + p.ps * s;
|
||||
angles.epsA += p.epsc * c + p.epss * s;
|
||||
}
|
||||
|
||||
return angles;
|
||||
}
|
||||
|
||||
|
||||
/*! Compute equatorial precession angles z, zeta, and theta using the P03
|
||||
* precession model.
|
||||
*/
|
||||
astro::EquatorialPrecessionAngles
|
||||
astro::EquatorialPrecessionAngles_P03(double T)
|
||||
{
|
||||
EquatorialPrecessionAngles prec;
|
||||
double T2 = T * T;
|
||||
double T3 = T2 * T;
|
||||
double T4 = T3 * T;
|
||||
double T5 = T4 * T;
|
||||
|
||||
prec.zetaA = ( 2.650545
|
||||
+ 2306.083227 * T
|
||||
+ 0.2988499 * T2
|
||||
+ 0.01801828 * T3
|
||||
- 0.000005971 * T4
|
||||
- 0.0000003173 * T5);
|
||||
prec.zA = ( - 2.650545
|
||||
+ 2306.077181 * T
|
||||
+ 1.0927348 * T2
|
||||
+ 0.01826837 * T3
|
||||
- 0.000028596 * T4
|
||||
- 0.0000002904 * T5);
|
||||
prec.thetaA = ( 2004.191903 * T
|
||||
- 0.4294934 * T2
|
||||
- 0.04182264 * T3
|
||||
- 0.000007089 * T4
|
||||
- 0.0000001274 * T5);
|
||||
|
||||
return prec;
|
||||
}
|
||||
|
||||
|
||||
/*! Compute the ecliptic pole coordinates PA and QA using the P03 precession
|
||||
* model. The quantities PA and QA are coordinates, but they are given in
|
||||
* units of arcseconds in P03. They should be divided by 1296000/2pi.
|
||||
*/
|
||||
astro::EclipticPole
|
||||
astro::EclipticPrecession_P03(double T)
|
||||
{
|
||||
astro::EclipticPole pole;
|
||||
double T2 = T * T;
|
||||
double T3 = T2 * T;
|
||||
double T4 = T3 * T;
|
||||
double T5 = T4 * T;
|
||||
|
||||
pole.PA = ( 4.199094 * T
|
||||
+ 0.1939873 * T2
|
||||
- 0.00022466 * T3
|
||||
- 0.000000912 * T4
|
||||
+ 0.0000000120 * T5);
|
||||
pole.QA = (-46.811015 * T
|
||||
+ 0.0510283 * T2
|
||||
+ 0.00052413 * T3
|
||||
- 0.00000646 * T4
|
||||
- 0.0000000172 * T5);
|
||||
|
||||
return pole;
|
||||
}
|
||||
|
||||
|
||||
/*! Calculate the angles of the ecliptic of date with respect to
|
||||
* the J2000 ecliptic using the P03 precession model.
|
||||
*/
|
||||
astro::EclipticAngles
|
||||
astro::EclipticPrecessionAngles_P03(double T)
|
||||
{
|
||||
astro::EclipticAngles ecl;
|
||||
double T2 = T * T;
|
||||
double T3 = T2 * T;
|
||||
double T4 = T3 * T;
|
||||
double T5 = T4 * T;
|
||||
|
||||
ecl.piA = ( 46.998973 * T
|
||||
- 0.0334926 * T2
|
||||
- 0.00012559 * T3
|
||||
+ 0.000000113 * T4
|
||||
- 0.0000000022 * T5);
|
||||
ecl.PiA = (629546.7936
|
||||
- 867.95758 * T
|
||||
+ 0.157992 * T2
|
||||
- 0.0005371 * T3
|
||||
- 0.00004797 * T4
|
||||
+ 0.000000072 * T5);
|
||||
|
||||
return ecl;
|
||||
}
|
||||
|
||||
|
||||
//static const double eps0 = 23.4392911;
|
||||
static const double eps0 = 23.4392802;
|
||||
|
||||
/*! Compute the general precession and obliquity using the P03
|
||||
* precession model. See PrecObliquity_P03LP for more details.
|
||||
*/
|
||||
astro::PrecessionAngles
|
||||
astro::PrecObliquity_P03(double T)
|
||||
{
|
||||
astro::PrecessionAngles prec;
|
||||
double T2 = T * T;
|
||||
double T3 = T2 * T;
|
||||
double T4 = T3 * T;
|
||||
double T5 = T4 * T;
|
||||
|
||||
prec.epsA = (eps0 * 3600
|
||||
- 46.836769 * T
|
||||
- 0.0001831 * T2
|
||||
+ 0.00200340 * T3
|
||||
- 0.000000576 * T4
|
||||
- 0.0000000434 * T5);
|
||||
prec.pA = ( 5028.796195 * T
|
||||
+ 1.1054348 * T2
|
||||
+ 0.00007964 * T3
|
||||
- 0.000023857 * T4
|
||||
- 0.0000000383 * T5);
|
||||
#if 0
|
||||
prec.chiA = ( 10.556403 * T
|
||||
- 2.3814292 * T2
|
||||
- 0.00121197 * T3
|
||||
+ 0.000170663 * T4
|
||||
- 0.0000000560 * T5);
|
||||
#endif
|
||||
|
||||
return prec;
|
||||
}
|
||||
|
||||
|
||||
#define TEST 0
|
||||
|
||||
#if TEST
|
||||
|
||||
#include <celmath/quaternion.h>
|
||||
|
||||
|
||||
int main(int argc, char* argv[])
|
||||
{
|
||||
double step = 10.0;
|
||||
|
||||
for (int i = -100; i <= 100; i++)
|
||||
{
|
||||
// Get time in Julian centuries from J2000
|
||||
double T = (double) i * step / 100;
|
||||
|
||||
astro::EclipticPole ecl = astro::EclipticPrecession_P03LP(T);
|
||||
astro::EclipticPole eclP03 = astro::EclipticPrecession_P03(T);
|
||||
astro::EclipticAngles eclAnglesP03 = astro::EclipticPrecessionAngles_P03(T);
|
||||
|
||||
//clog << ecl.PA - eclP03.PA << ", " << ecl.QA - eclP03.QA << endl;
|
||||
ecl = eclP03;
|
||||
Vec3d v(0.0, 0.0, 0.0);
|
||||
v.x = ecl.PA * 2.0 * PI / 1296000;
|
||||
v.y = -ecl.QA * 2.0 * PI / 1296000;
|
||||
v.z = sqrt(1.0 - v.x * v.x - v.y * v.y);
|
||||
|
||||
Quatd eclRotation = Quatd::vecToVecRotation(Vec3d(0.0, 0.0, 1.0), v);
|
||||
Quatd eclRotation2 = Quatd::zrotation(-degToRad(eclAnglesP03.PiA / 3600.0)) *
|
||||
Quatd::xrotation(degToRad(eclAnglesP03.piA / 3600.0)) *
|
||||
Quatd::zrotation(degToRad(eclAnglesP03.PiA / 3600.0));
|
||||
Quatd eclRotation3;
|
||||
{
|
||||
// Calculate the angles pi and Pi from the ecliptic pole coordinates
|
||||
// P and Q:
|
||||
// P = sin(pi)*sin(Pi)
|
||||
// Q = sin(pi)*cos(Pi)
|
||||
double P = eclP03.PA * 2.0 * PI / 1296000;
|
||||
double Q = eclP03.QA * 2.0 * PI / 1296000;
|
||||
double piA = asin(sqrt(P * P + Q * Q));
|
||||
double PiA = atan2(P, Q);
|
||||
|
||||
// Calculate the rotation from the J2000 ecliptic to the ecliptic
|
||||
// of date.
|
||||
eclRotation3 = Quatd::zrotation(-PiA) *
|
||||
Quatd::xrotation(piA) *
|
||||
Quatd::zrotation(PiA);
|
||||
|
||||
piA = radToDeg(piA) * 3600;
|
||||
PiA = radToDeg(PiA) * 3600;
|
||||
clog << "Pi: " << PiA << ", " << eclAnglesP03.PiA << endl;
|
||||
clog << "pi: " << piA << ", " << eclAnglesP03.piA << endl;
|
||||
}
|
||||
|
||||
astro::PrecessionAngles prec = astro::PrecObliquity_P03LP(T);
|
||||
Quatd p03lpRot = Quatd::xrotation(degToRad(prec.epsA / 3600.0)) *
|
||||
Quatd::zrotation(-degToRad(prec.pA / 3600.0));
|
||||
p03lpRot = p03lpRot * ~eclRotation3;
|
||||
//p03lpRot = ~eclRotation2 * p03lpRot;
|
||||
|
||||
astro::PrecessionAngles prec2 = astro::PrecObliquity_P03(T);
|
||||
//clog << prec.epsA - prec2.epsA << ", " << prec.pA - prec2.pA << endl;
|
||||
|
||||
astro::EquatorialPrecessionAngles precP03 = astro::EquatorialPrecessionAngles_P03(T);
|
||||
Quatd p03Rot = Quatd::zrotation(-degToRad(precP03.zA / 3600.0)) *
|
||||
Quatd::yrotation( degToRad(precP03.thetaA / 3600.0)) *
|
||||
Quatd::zrotation(-degToRad(precP03.zetaA / 3600.0));
|
||||
p03Rot = p03Rot * Quatd::xrotation(degToRad(eps0));
|
||||
|
||||
Vec3d xaxis(0, 0, 1);
|
||||
Vec3d v0 = xaxis * p03lpRot.toMatrix3();
|
||||
Vec3d v1 = xaxis * p03Rot.toMatrix3();
|
||||
|
||||
// Show the angle between the J2000 ecliptic pole and the ecliptic
|
||||
// pole at T centuries from J2000
|
||||
double theta0 = acos(xaxis * v0);
|
||||
double theta1 = acos(xaxis * v1);
|
||||
|
||||
clog << T * 100 << ": "
|
||||
<< radToDeg(acos(v0 * v1)) * 3600 << ", "
|
||||
<< radToDeg(theta0) << ", "
|
||||
<< radToDeg(theta1) << ", "
|
||||
<< radToDeg(theta1 - theta0) * 3600
|
||||
<< endl;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
#endif // TEST
|
||||
|
|
|
|||
|
|
@ -15,17 +15,32 @@ namespace astro
|
|||
|
||||
// PA and QA are the location of the pole of the ecliptic of date
|
||||
// with respect to the fixed ecliptic of J2000.0
|
||||
struct EclipticAngles
|
||||
struct EclipticPole
|
||||
{
|
||||
double PA;
|
||||
double QA;
|
||||
};
|
||||
|
||||
|
||||
// epsA is the obliquity with respect to the ecliptic of date. pA is
|
||||
// the general precession---the angle between the ascending node of the
|
||||
// equator on the ecliptic of date, and the ascending node of the ecliptic
|
||||
// of date on the J2000.0 ecliptic.
|
||||
// piA and PiA are angles that transform the J2000 ecliptic to the
|
||||
// ecliptic of date. They are related to the ecliptic pole coordinates
|
||||
// PA and QA:
|
||||
// PA = sin(piA)*sin(PiA)
|
||||
// QA = sin(piA)*cos(PiA)
|
||||
//
|
||||
// PiA is the angle along the J2000 ecliptic between the J2000 equinox
|
||||
// and the intersection of the J2000 ecliptic and ecliptic of date.
|
||||
struct EclipticAngles
|
||||
{
|
||||
double piA;
|
||||
double PiA;
|
||||
};
|
||||
|
||||
|
||||
// epsA is the angle between the ecliptic and mean equator of date. pA is the
|
||||
// general precession: the difference between angles L and PiA. L is the angle
|
||||
// along the mean ecliptic of date from the equinox of date to the
|
||||
// intersection of the J2000 ecliptic and ecliptic of date.
|
||||
struct PrecessionAngles
|
||||
{
|
||||
double pA; // precession
|
||||
|
|
@ -33,7 +48,20 @@ struct PrecessionAngles
|
|||
};
|
||||
|
||||
|
||||
extern EclipticAngles EclipticPrecession_P03LP(double T);
|
||||
struct EquatorialPrecessionAngles
|
||||
{
|
||||
double zetaA;
|
||||
double zA;
|
||||
double thetaA;
|
||||
};
|
||||
|
||||
|
||||
extern EclipticPole EclipticPrecession_P03LP(double T);
|
||||
extern PrecessionAngles PrecObliquity_P03LP(double T);
|
||||
|
||||
extern EclipticPole EclipticPrecession_P03(double T);
|
||||
extern EclipticAngles EclipticPrecessionAngles_P03(double T);
|
||||
extern PrecessionAngles PrecObliquity_P03(double T);
|
||||
extern EquatorialPrecessionAngles EquatorialPrecessionAngles_P03(double T);
|
||||
|
||||
};
|
||||
|
|
|
|||
Loading…
Reference in New Issue