You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
405 lines
13 KiB
405 lines
13 KiB
/***************************************************************************
|
|
ksnumbers.cpp - description
|
|
-------------------
|
|
begin : Sun Jan 13 2002
|
|
copyright : (C) 2002-2005 by Jason Harris
|
|
email : kstars@30doradus.org
|
|
copyright : (C) 2004-2005 by Pablo de Vicente
|
|
email : p.devicente@wanadoo.es
|
|
***************************************************************************/
|
|
|
|
/***************************************************************************
|
|
* *
|
|
* 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 "ksnumbers.h"
|
|
|
|
// 63 elements
|
|
const int KSNumbers::arguments[NUTTERMS][5] = {
|
|
{ 0, 0, 0, 0, 1},
|
|
{-2, 0, 0, 2, 2},
|
|
{ 0, 0, 0, 2, 2},
|
|
{ 0, 0, 0, 0, 2},
|
|
{ 0, 1, 0, 0, 0},
|
|
{ 0, 0, 1, 0, 0},
|
|
{-2, 1, 0, 2, 2},
|
|
{ 0, 0, 0, 2, 1},
|
|
{ 0, 0, 1, 2, 2},
|
|
{-2,-1, 0, 2, 2},
|
|
{-2, 0, 1, 0, 0},
|
|
{-2, 0, 0, 2, 1},
|
|
{ 0, 0,-1, 2, 2},
|
|
{ 2, 0, 0, 0, 0},
|
|
{ 0, 0, 1, 0, 1},
|
|
{ 2, 0,-1, 2, 2},
|
|
{ 0, 0,-1, 0, 1},
|
|
{ 0, 0, 1, 2, 1},
|
|
{-2, 0, 2, 0, 0},
|
|
{ 0, 0,-2, 2, 1},
|
|
{ 2, 0, 0, 2, 2},
|
|
{ 0, 0, 2, 2, 2},
|
|
{ 0, 0, 2, 0, 0},
|
|
{-2, 0, 1, 2, 2},
|
|
{ 0, 0, 0, 2, 0},
|
|
{-2, 0, 0, 2, 0},
|
|
{ 0, 0,-1, 2, 1},
|
|
{ 0, 2, 0, 0, 0},
|
|
{ 2, 0,-1, 0, 1},
|
|
{-2, 2, 0, 2, 2},
|
|
{ 0, 1, 0, 0, 1},
|
|
{-2, 0, 1, 0, 1},
|
|
{ 0,-1, 0, 0, 1},
|
|
{ 0, 0, 2,-2, 0},
|
|
{ 2, 0,-1, 2, 1},
|
|
{ 2, 0, 1, 2, 2},
|
|
{ 0, 1, 0, 2, 2},
|
|
{-2, 1, 1, 0, 0},
|
|
{ 0,-1, 0, 2, 2},
|
|
{ 2, 0, 0, 2, 1},
|
|
{ 2, 0, 1, 0, 0},
|
|
{-2, 0, 2, 2, 2},
|
|
{-2, 0, 1, 2, 1},
|
|
{ 2, 0,-2, 0, 1},
|
|
{ 2, 0, 0, 0, 1},
|
|
{ 0,-1, 1, 0, 0},
|
|
{-2,-1, 0, 2, 1},
|
|
{-2, 0, 0, 0, 1},
|
|
{ 0, 0, 2, 2, 1},
|
|
{-2, 0, 2, 0, 1},
|
|
{-2, 1, 0, 2, 1},
|
|
{ 0, 0, 1,-2, 0},
|
|
{-1, 0, 1, 0, 0},
|
|
{-2, 1, 0, 0, 0},
|
|
{ 1, 0, 0, 0, 0},
|
|
{ 0, 0, 1, 2, 0},
|
|
{ 0, 0,-2, 2, 2},
|
|
{-1,-1, 1, 0, 0},
|
|
{ 0, 1, 1, 0, 0},
|
|
{ 0,-1, 1, 2, 2},
|
|
{ 2,-1,-1, 2, 2},
|
|
{ 0, 0, 3, 2, 2},
|
|
{ 2,-1, 0, 2, 2}
|
|
};
|
|
|
|
const int KSNumbers::amp[NUTTERMS][4] = {
|
|
{-171996,-1742, 92025, 89},
|
|
{ -13187, -16, 5736,-31},
|
|
{ -2274, -2, 977, -5},
|
|
{ 2062, 2, -895, 5},
|
|
{ 1426, -34, 54, -1},
|
|
{ 712, 1, -7, 0},
|
|
{ -517, 12, 224, -6},
|
|
{ -386, -4, 200, 0},
|
|
{ -301, 0, 129, -1},
|
|
{ 217, -5, -95, 3},
|
|
{ -158, 0, 0, 0},
|
|
{ 129, 1, -70, 0},
|
|
{ 123, 0, -53, 0},
|
|
{ 63, 0, 0, 0},
|
|
{ 63, 1, -33, 0},
|
|
{ -59, 0, 26, 0},
|
|
{ -58, -1, 32, 0},
|
|
{ -51, 0, 27, 0},
|
|
{ 48, 0, 0, 0},
|
|
{ 46, 0, -24, 0},
|
|
{ -38, 0, 16, 0},
|
|
{ -31, 0, 13, 0},
|
|
{ 29, 0, 0, 0},
|
|
{ 29, 0, -12, 0},
|
|
{ 26, 0, 0, 0},
|
|
{ -22, 0, 0, 0},
|
|
{ 21, 0, -10, 0},
|
|
{ 17, -1, 0, 0},
|
|
{ 16, 0, -8, 0},
|
|
{ -16, 1, 7, 0},
|
|
{ -15, 0, 9, 0},
|
|
{ -13, 0, 7, 0},
|
|
{ -12, 0, 6, 0},
|
|
{ 11, 0, 0, 0},
|
|
{ -10, 0, 5, 0},
|
|
{ -8, 0, 3, 0},
|
|
{ 7, 0, -3, 0},
|
|
{ -7, 0, 0, 0},
|
|
{ -7, 0, 3, 0},
|
|
{ -7, 0, 3, 0},
|
|
{ 6, 0, 0, 0},
|
|
{ 6, 0, -3, 0},
|
|
{ 6, 0, -3, 0},
|
|
{ -6, 0, 3, 0},
|
|
{ -6, 0, 3, 0},
|
|
{ 5, 0, 0, 0},
|
|
{ -5, 0, 3, 0},
|
|
{ -5, 0, 3, 0},
|
|
{ -5, 0, 3, 0},
|
|
{ 4, 0, 0, 0},
|
|
{ 4, 0, 0, 0},
|
|
{ 4, 0, 0, 0},
|
|
{ -4, 0, 0, 0},
|
|
{ -4, 0, 0, 0},
|
|
{ -4, 0, 0, 0},
|
|
{ 3, 0, 0, 0},
|
|
{ -3, 0, 0, 0},
|
|
{ -3, 0, 0, 0},
|
|
{ -3, 0, 0, 0},
|
|
{ -3, 0, 0, 0},
|
|
{ -3, 0, 0, 0},
|
|
{ -3, 0, 0, 0},
|
|
{ -3, 0, 0, 0}
|
|
};
|
|
|
|
|
|
KSNumbers::KSNumbers( long double jd ){
|
|
K.setD( 20.49552 / 3600. ); //set the constant of aberration
|
|
updateValues( jd );
|
|
}
|
|
|
|
KSNumbers::~KSNumbers(){
|
|
}
|
|
|
|
void KSNumbers::updateValues( long double jd ) {
|
|
dms arg;
|
|
double args, argc;
|
|
|
|
days = jd;
|
|
|
|
//Julian Centuries since J2000.0
|
|
T = ( jd - J2000 ) / 36525.;
|
|
|
|
// Julian Millenia since J2000.0
|
|
jm = T / 10.0;
|
|
|
|
double T2 = T*T;
|
|
double T3 = T2*T;
|
|
|
|
//Sun's Mean Longitude
|
|
L.setD( 280.46645 + 36000.76983*T + 0.0003032*T2 );
|
|
|
|
//Mean elongation of the Moon from the Sun
|
|
D.setD( 297.85036 + 445267.111480*T - 0.0019142*T2 + T3/189474.);
|
|
|
|
//Sun's Mean Anomaly
|
|
M.setD( 357.52910 + 35999.05030*T - 0.0001559*T2 - 0.00000048*T3);
|
|
|
|
//Moon's Mean Anomaly
|
|
MM.setD( 134.96298 + 477198.867398*T + 0.0086972*T2 + T3/56250.0 );
|
|
|
|
//Moon's Mean Longitude
|
|
LM.setD( 218.3164591 + 481267.88134236*T - 0.0013268*T2 + T3/538841. - T*T*T*T/6519400.);
|
|
|
|
//Moon's argument of latitude
|
|
F.setD( 93.27191 + 483202.017538*T - 0.0036825*T2 + T3/327270.);
|
|
|
|
//Longitude of Moon's Ascending Node
|
|
O.setD( 125.04452 - 1934.136261*T + 0.0020708*T2 + T3/450000.0 );
|
|
|
|
//Earth's orbital eccentricity
|
|
e = 0.016708617 - 0.000042037*T - 0.0000001236*T2;
|
|
|
|
double C = ( 1.914600 - 0.004817*T - 0.000014*T2 ) * sin( M.radians() )
|
|
+ ( 0.019993 - 0.000101*T ) * sin( 2.0* M.radians() )
|
|
+ 0.000290 * sin( 3.0* M.radians() );
|
|
|
|
//Sun's True Longitude
|
|
L0.setD( L.Degrees() + C );
|
|
|
|
//Sun's True Anomaly
|
|
M0.setD( M.Degrees() + C );
|
|
|
|
//Obliquity of the Ecliptic
|
|
double U = T/100.0;
|
|
double dObliq = -4680.93*U - 1.55*U*U + 1999.25*U*U*U
|
|
- 51.38*U*U*U*U - 249.67*U*U*U*U*U
|
|
- 39.05*U*U*U*U*U*U + 7.12*U*U*U*U*U*U*U
|
|
+ 27.87*U*U*U*U*U*U*U*U + 5.79*U*U*U*U*U*U*U*U*U
|
|
+ 2.45*U*U*U*U*U*U*U*U*U*U;
|
|
Obliquity.setD( 23.43929111 + dObliq/3600.0);
|
|
|
|
//Nutation parameters
|
|
dms L2, M2, O2;
|
|
double sin2L, cos2L, sin2M, cos2M;
|
|
double sinO, cosO, sin2O, cos2O;
|
|
|
|
O2.setD( 2.0*O.Degrees() );
|
|
L2.setD( 2.0*L.Degrees() ); //twice mean ecl. long. of Sun
|
|
M2.setD( 2.0*LM.Degrees() ); //twice mean ecl. long. of Moon
|
|
|
|
O.SinCos( sinO, cosO );
|
|
O2.SinCos( sin2O, cos2O );
|
|
L2.SinCos( sin2L, cos2L );
|
|
M2.SinCos( sin2M, cos2M );
|
|
|
|
// deltaEcLong = ( -17.2*sinO - 1.32*sin2L - 0.23*sin2M + 0.21*sin2O)/3600.0; //Ecl. long. correction
|
|
// deltaObliquity = ( 9.2*cosO + 0.57*cos2L + 0.10*cos2M - 0.09*cos2O)/3600.0; //Obliq. correction
|
|
|
|
deltaEcLong = 0.;
|
|
deltaObliquity = 0.;
|
|
|
|
for (unsigned int i=0; i < NUTTERMS; i++) {
|
|
|
|
arg.setD ( arguments[i][0]*D.Degrees() + arguments[i][1]*M.Degrees() +
|
|
arguments[i][2]*MM.Degrees() + arguments[i][3]*F.Degrees() + arguments[i][4]*O.Degrees() );
|
|
arg.SinCos( args, argc );
|
|
|
|
deltaEcLong += (amp[i][0] + amp[i][1]/10. * T ) * args * 1e-4 ;
|
|
deltaObliquity += (amp[i][2] + amp[i][3]/10. * T ) * argc * 1e-4 ;
|
|
}
|
|
|
|
deltaEcLong/= 3600.0;
|
|
deltaObliquity /= 3600.0;
|
|
|
|
//Compute Precession Matrices:
|
|
XP.setD( 0.6406161*T + 0.0000839*T2 + 0.0000050*T3 );
|
|
YP.setD( 0.5567530*T - 0.0001185*T2 - 0.0000116*T3 );
|
|
ZP.setD( 0.6406161*T + 0.0003041*T2 + 0.0000051*T3 );
|
|
|
|
XP.SinCos( SX, CX );
|
|
YP.SinCos( SY, CY );
|
|
ZP.SinCos( SZ, CZ );
|
|
|
|
//P1 is used to precess from any epoch to J2000
|
|
P1[0][0] = CX*CY*CZ - SX*SZ;
|
|
P1[1][0] = CX*CY*SZ + SX*CZ;
|
|
P1[2][0] = CX*SY;
|
|
P1[0][1] = -1.0*SX*CY*CZ - CX*SZ;
|
|
P1[1][1] = -1.0*SX*CY*SZ + CX*CZ;
|
|
P1[2][1] = -1.0*SX*SY;
|
|
P1[0][2] = -1.0*SY*CZ;
|
|
P1[1][2] = -1.0*SY*SZ;
|
|
P1[2][2] = CY;
|
|
|
|
//P2 is used to precess from J2000 to any other epoch (it is the transpose of P1)
|
|
P2[0][0] = CX*CY*CZ - SX*SZ;
|
|
P2[1][0] = -1.0*SX*CY*CZ - CX*SZ;
|
|
P2[2][0] = -1.0*SY*CZ;
|
|
P2[0][1] = CX*CY*SZ + SX*CZ;
|
|
P2[1][1] = -1.0*SX*CY*SZ + CX*CZ;
|
|
P2[2][1] = -1.0*SY*SZ;
|
|
P2[0][2] = CX*SY;
|
|
P2[1][2] = -1.0*SX*SY;
|
|
P2[2][2] = CY;
|
|
|
|
//Compute Precession Matrices from B1950 to 1984 using Newcomb formulae
|
|
|
|
XB.setD( 0.217697 );
|
|
YB.setD( 0.189274 );
|
|
ZB.setD( 0.217722 );
|
|
|
|
XB.SinCos( SXB, CXB );
|
|
YB.SinCos( SYB, CYB );
|
|
ZB.SinCos( SZB, CZB );
|
|
|
|
//P1B is used to precess from 1984 to B1950:
|
|
|
|
P1B[0][0] = CXB*CYB*CZB - SXB*SZB;
|
|
P1B[1][0] = CXB*CYB*SZB + SXB*CZB;
|
|
P1B[2][0] = CXB*SYB;
|
|
P1B[0][1] = -1.0*SXB*CYB*CZB - CXB*SZB;
|
|
P1B[1][1] = -1.0*SXB*CYB*SZB + CXB*CZB;
|
|
P1B[2][1] = -1.0*SXB*SYB;
|
|
P1B[0][2] = -1.0*SYB*CZB;
|
|
P1B[1][2] = -1.0*SYB*SZB;
|
|
P1B[2][2] = CYB;
|
|
|
|
//P2 is used to precess from B1950 to 1984 (it is the transpose of P1)
|
|
P2B[0][0] = CXB*CYB*CZB - SXB*SZB;
|
|
P2B[1][0] = -1.0*SXB*CYB*CZB - CXB*SZB;
|
|
P2B[2][0] = -1.0*SYB*CZB;
|
|
P2B[0][1] = CXB*CYB*SZB + SXB*CZB;
|
|
P2B[1][1] = -1.0*SXB*CYB*SZB + CXB*CZB;
|
|
P2B[2][1] = -1.0*SYB*SZB;
|
|
P2B[0][2] = CXB*SYB;
|
|
P2B[1][2] = -1.0*SXB*SYB;
|
|
P2B[2][2] = CYB;
|
|
|
|
|
|
// Mean longitudes for the planets. radians
|
|
//
|
|
|
|
// TODO Pasar a grados
|
|
double LVenus = 3.1761467+1021.3285546*T; // Venus
|
|
double LMars = 1.7534703+ 628.3075849*T; // Mars
|
|
double LEarth = 6.2034809+ 334.0612431*T; // Earth
|
|
double LJupiter = 0.5995465+ 52.9690965*T; // Jupiter
|
|
double LSaturn = 0.8740168+ 21.3299095*T; // Saturn
|
|
double LNeptune = 5.3118863+ 3.8133036*T; // Neptune
|
|
double LUranus = 5.4812939+ 7.4781599*T; // Uranus
|
|
|
|
double LMRad = 3.8103444+8399.6847337*T; // Moon
|
|
double DRad = 5.1984667+7771.3771486*T;
|
|
double MMRad = 2.3555559+8328.6914289*T; // Moon
|
|
double FRad = 1.6279052+8433.4661601*T;
|
|
|
|
/** Contibutions to the velocity of the Earth referred to the barycenter of the solar system
|
|
in the J2000 equatorial system
|
|
Velocities 10^{-8} AU/day
|
|
Ron & Vondrak method
|
|
**/
|
|
|
|
double vondrak[36][7] = {
|
|
{LMars, -1719914-2*T, -25, 25-13*T,1578089+156*T, 10+32*T,684185-358*T},
|
|
{2*LMars, 6434+141*T,28007-107*T,25697-95*T, -5904-130*T,11141-48*T, -2559-55*T},
|
|
{LJupiter, 715, 0, 6, -657, -15, -282},
|
|
{LMRad, 715, 0, 0, -656, 0, -285},
|
|
{3*LMars, 486-5*T, -236-4*T, -216-4*T, -446+5*T, -94, -193},
|
|
{LSaturn, 159, 0, 2, -147, -6, -61},
|
|
{FRad, 0, 0, 0, 26, 0, -59},
|
|
{LMRad+MMRad, 39, 0, 0, -36, 0, -16},
|
|
{2*LJupiter, 33, -10, -9, -30, -5, -13},
|
|
{2*LMars-LJupiter, 31, 1, 1, -28, 0, -12},
|
|
{3*LMars-8*LEarth+3*LJupiter, 8, -28, 25, 8, 11, 3},
|
|
{5*LMars-8*LEarth+3*LJupiter, 8, -28, -25, -8, -11, -3},
|
|
{2*LVenus-LMars, 21, 0, 0, -19, 0, -8},
|
|
{LVenus, -19, 0, 0, 17, 0, 8},
|
|
{LNeptune, 17, 0, 0, -16, 0, -7},
|
|
{LMars-2*LJupiter, 16, 0, 0, 15, 1, 7},
|
|
{LUranus, 16, 0, 1, -15, -3, -6},
|
|
{LMars+LJupiter, 11, -1, -1, -10, -1, -5},
|
|
{2*LVenus-2*LMars, 0, -11, -10, 0, -4, 0},
|
|
{LMars-LJupiter, -11, -2, -2, 9, -1, 4},
|
|
{4*LMars, -7, -8, -8, 6, -3, 3},
|
|
{3*LMars-2*LJupiter, -10, 0, 0, 9, 0, 4},
|
|
{LVenus-2*LMars, -9, 0, 0, -9, 0, -4},
|
|
{2*LVenus-3*LMars, -9, 0, 0, -8, 0, -4},
|
|
{2*LSaturn, 0, -9, -8, 0, -3, 0},
|
|
{2*LVenus-4*LMars, 0, -9, 8, 0, 3, 0},
|
|
{3*LMars-2*LEarth, 8, 0, 0, -8, 0, -3},
|
|
{LMRad+2*DRad-MMRad, 8, 0, 0, -7, 0, -3},
|
|
{8*LVenus-12*LMars, -4, -7, -6, 4, -3, 2},
|
|
{8*LVenus-14*LMars, -4, -7, 6, -4, 3, -2},
|
|
{2*LEarth, -6, -5, -4, 5, -2, 2},
|
|
{3*LVenus-4*LMars, -1, -1, -2, -7, 1, -4},
|
|
{2*LMars-2*LJupiter, 4, -6, -5, -4, -2, -2},
|
|
{3*LVenus-3*LMars, 0, -7, -6, 0, -3, 0},
|
|
{2*LMars-2*LEarth, 5, -5, -4, -5, -2, -2},
|
|
{LMRad-2*DRad, 5, 0, 0, -5, 0, -2}
|
|
};
|
|
|
|
dms anglev;
|
|
double sa, ca;
|
|
// Vearth X component
|
|
vearth[0] = 0.;
|
|
// Vearth Y component
|
|
vearth[1] = 0.;
|
|
// Vearth Z component
|
|
vearth[2] = 0.;
|
|
|
|
for (unsigned int i=0; i<36; i++) {
|
|
anglev.setRadians(vondrak[i][0]);
|
|
anglev.SinCos(sa,ca);
|
|
for (unsigned int j=0; j<3; j++) {
|
|
vearth[j] += vondrak[i][2*j+1]*sa +vondrak[i][2*j+2]*ca;
|
|
}
|
|
}
|
|
|
|
const double UA2km = 1.49597870/86400.; // 10^{-8}*UA/dia -> km/s
|
|
|
|
for (unsigned int j=0; j<3; j++) {
|
|
vearth[j] = vearth[j] * UA2km;
|
|
}
|
|
}
|