kstars

ksnumbers.cpp

00001 /***************************************************************************
00002                           ksnumbers.cpp  -  description
00003                              -------------------
00004     begin                : Sun Jan 13 2002
00005     copyright            : (C) 2002-2005 by Jason Harris
00006     email                : kstars@30doradus.org
00007     copyright            : (C) 2004-2005 by Pablo de Vicente
00008     email                : p.devicente@wanadoo.es
00009  ***************************************************************************/
00010 
00011 /***************************************************************************
00012  *                                                                         *
00013  *   This program is free software; you can redistribute it and/or modify  *
00014  *   it under the terms of the GNU General Public License as published by  *
00015  *   the Free Software Foundation; either version 2 of the License, or     *
00016  *   (at your option) any later version.                                   *
00017  *                                                                         *
00018  ***************************************************************************/
00019 
00020 #include "ksnumbers.h"
00021 
00022 // 63 elements
00023 const int KSNumbers::arguments[NUTTERMS][5] = {
00024 { 0, 0, 0, 0, 1},
00025 {-2, 0, 0, 2, 2},
00026 { 0, 0, 0, 2, 2},
00027 { 0, 0, 0, 0, 2},
00028 { 0, 1, 0, 0, 0},
00029 { 0, 0, 1, 0, 0},
00030 {-2, 1, 0, 2, 2},
00031 { 0, 0, 0, 2, 1},
00032 { 0, 0, 1, 2, 2},
00033 {-2,-1, 0, 2, 2},
00034 {-2, 0, 1, 0, 0},
00035 {-2, 0, 0, 2, 1},
00036 { 0, 0,-1, 2, 2},
00037 { 2, 0, 0, 0, 0},
00038 { 0, 0, 1, 0, 1},
00039 { 2, 0,-1, 2, 2},
00040 { 0, 0,-1, 0, 1},
00041 { 0, 0, 1, 2, 1},
00042 {-2, 0, 2, 0, 0},
00043 { 0, 0,-2, 2, 1},
00044 { 2, 0, 0, 2, 2},
00045 { 0, 0, 2, 2, 2},
00046 { 0, 0, 2, 0, 0},
00047 {-2, 0, 1, 2, 2},
00048 { 0, 0, 0, 2, 0},
00049 {-2, 0, 0, 2, 0},
00050 { 0, 0,-1, 2, 1},
00051 { 0, 2, 0, 0, 0},
00052 { 2, 0,-1, 0, 1},
00053 {-2, 2, 0, 2, 2},
00054 { 0, 1, 0, 0, 1},
00055 {-2, 0, 1, 0, 1},
00056 { 0,-1, 0, 0, 1},
00057 { 0, 0, 2,-2, 0},
00058 { 2, 0,-1, 2, 1},
00059 { 2, 0, 1, 2, 2},
00060 { 0, 1, 0, 2, 2},
00061 {-2, 1, 1, 0, 0},
00062 { 0,-1, 0, 2, 2},
00063 { 2, 0, 0, 2, 1},
00064 { 2, 0, 1, 0, 0},
00065 {-2, 0, 2, 2, 2},
00066 {-2, 0, 1, 2, 1},
00067 { 2, 0,-2, 0, 1},
00068 { 2, 0, 0, 0, 1},
00069 { 0,-1, 1, 0, 0},
00070 {-2,-1, 0, 2, 1},
00071 {-2, 0, 0, 0, 1},
00072 { 0, 0, 2, 2, 1},
00073 {-2, 0, 2, 0, 1},
00074 {-2, 1, 0, 2, 1},
00075 { 0, 0, 1,-2, 0},
00076 {-1, 0, 1, 0, 0},
00077 {-2, 1, 0, 0, 0},
00078 { 1, 0, 0, 0, 0},
00079 { 0, 0, 1, 2, 0},
00080 { 0, 0,-2, 2, 2},
00081 {-1,-1, 1, 0, 0},
00082 { 0, 1, 1, 0, 0},
00083 { 0,-1, 1, 2, 2},
00084 { 2,-1,-1, 2, 2},
00085 { 0, 0, 3, 2, 2},
00086 { 2,-1, 0, 2, 2}
00087 };
00088 
00089 const int KSNumbers::amp[NUTTERMS][4] = {
00090 {-171996,-1742, 92025, 89},
00091 { -13187,  -16,  5736,-31},
00092 {  -2274,   -2,   977, -5},
00093 {   2062,    2,  -895,  5},
00094 {   1426,  -34,    54, -1},
00095 {    712,    1,    -7,  0},
00096 {   -517,   12,   224, -6},
00097 {   -386,   -4,   200,  0},
00098 {   -301,    0,   129, -1},
00099 {    217,   -5,   -95,  3},
00100 {   -158,    0,     0,  0},
00101 {    129,    1,   -70,  0},
00102 {    123,    0,   -53,  0},
00103 {     63,    0,     0,  0},
00104 {     63,    1,   -33,  0},
00105 {    -59,    0,    26,  0},
00106 {    -58,   -1,    32,  0},
00107 {    -51,    0,    27,  0},
00108 {     48,    0,     0,  0},
00109 {     46,    0,   -24,  0},
00110 {    -38,    0,    16,  0},
00111 {    -31,    0,    13,  0},
00112 {     29,    0,     0,  0},
00113 {     29,    0,   -12,  0},
00114 {     26,    0,     0,  0},
00115 {    -22,    0,     0,  0},
00116 {     21,    0,   -10,  0},
00117 {     17,   -1,     0,  0},
00118 {     16,    0,    -8,  0},
00119 {    -16,    1,     7,  0},
00120 {    -15,    0,     9,  0},
00121 {    -13,    0,     7,  0},
00122 {    -12,    0,     6,  0},
00123 {     11,    0,     0,  0},
00124 {    -10,    0,     5,  0},
00125 {     -8,    0,     3,  0},
00126 {      7,    0,    -3,  0},
00127 {     -7,    0,     0,  0},
00128 {     -7,    0,     3,  0},
00129 {     -7,    0,     3,  0},
00130 {      6,    0,     0,  0},
00131 {      6,    0,    -3,  0},
00132 {      6,    0,    -3,  0},
00133 {     -6,    0,     3,  0},
00134 {     -6,    0,     3,  0},
00135 {      5,    0,     0,  0},
00136 {     -5,    0,     3,  0},
00137 {     -5,    0,     3,  0},
00138 {     -5,    0,     3,  0},
00139 {      4,    0,     0,  0},
00140 {      4,    0,     0,  0},
00141 {      4,    0,     0,  0},
00142 {     -4,    0,     0,  0},
00143 {     -4,    0,     0,  0},
00144 {     -4,    0,     0,  0},
00145 {      3,    0,     0,  0},
00146 {     -3,    0,     0,  0},
00147 {     -3,    0,     0,  0},
00148 {     -3,    0,     0,  0},
00149 {     -3,    0,     0,  0},
00150 {     -3,    0,     0,  0},
00151 {     -3,    0,     0,  0},
00152 {     -3,    0,     0,  0}
00153 };
00154 
00155 
00156 KSNumbers::KSNumbers( long double jd ){
00157     K.setD( 20.49552 / 3600. );  //set the constant of aberration
00158     updateValues( jd );
00159 }
00160 
00161 KSNumbers::~KSNumbers(){
00162 }
00163 
00164 void KSNumbers::updateValues( long double jd ) {
00165     dms arg;
00166     double args, argc;
00167 
00168     days = jd;
00169 
00170     //Julian Centuries since J2000.0
00171     T = ( jd - J2000 ) / 36525.;
00172 
00173     // Julian Millenia since J2000.0
00174     jm = T / 10.0;
00175 
00176     double T2 = T*T;
00177     double T3 = T2*T;
00178 
00179     //Sun's Mean Longitude
00180     L.setD( 280.46645 + 36000.76983*T + 0.0003032*T2 );
00181 
00182     //Mean elongation of the Moon from the Sun
00183     D.setD( 297.85036 + 445267.111480*T - 0.0019142*T2 + T3/189474.);
00184         
00185     //Sun's Mean Anomaly
00186     M.setD( 357.52910 + 35999.05030*T - 0.0001559*T2 - 0.00000048*T3);
00187     
00188     //Moon's Mean Anomaly
00189     MM.setD( 134.96298 + 477198.867398*T + 0.0086972*T2 + T3/56250.0 );
00190 
00191     //Moon's Mean Longitude
00192     LM.setD( 218.3164591 + 481267.88134236*T - 0.0013268*T2 + T3/538841. - T*T*T*T/6519400.);
00193     
00194     //Moon's argument of latitude
00195     F.setD( 93.27191 + 483202.017538*T - 0.0036825*T2 + T3/327270.);
00196 
00197     //Longitude of Moon's Ascending Node
00198     O.setD( 125.04452 - 1934.136261*T + 0.0020708*T2 + T3/450000.0 );
00199 
00200     //Earth's orbital eccentricity
00201     e = 0.016708617 - 0.000042037*T - 0.0000001236*T2;
00202     
00203     double C = ( 1.914600 - 0.004817*T - 0.000014*T2 ) * sin( M.radians() )
00204                         + ( 0.019993 - 0.000101*T ) * sin( 2.0* M.radians() )
00205                         + 0.000290 * sin( 3.0* M.radians() );
00206 
00207     //Sun's True Longitude
00208     L0.setD( L.Degrees() + C );
00209 
00210     //Sun's True Anomaly
00211     M0.setD( M.Degrees() + C );
00212 
00213     //Obliquity of the Ecliptic
00214     double U = T/100.0;
00215     double dObliq = -4680.93*U - 1.55*U*U + 1999.25*U*U*U
00216                     - 51.38*U*U*U*U - 249.67*U*U*U*U*U
00217                     - 39.05*U*U*U*U*U*U + 7.12*U*U*U*U*U*U*U
00218                     + 27.87*U*U*U*U*U*U*U*U + 5.79*U*U*U*U*U*U*U*U*U
00219                     + 2.45*U*U*U*U*U*U*U*U*U*U;
00220     Obliquity.setD( 23.43929111 + dObliq/3600.0);
00221 
00222     //Nutation parameters
00223     dms  L2, M2, O2;
00224     double sin2L, cos2L, sin2M, cos2M;
00225     double sinO, cosO, sin2O, cos2O;
00226     
00227     O2.setD( 2.0*O.Degrees() );
00228     L2.setD( 2.0*L.Degrees() ); //twice mean ecl. long. of Sun
00229     M2.setD( 2.0*LM.Degrees() );  //twice mean ecl. long. of Moon
00230     
00231     O.SinCos( sinO, cosO );
00232     O2.SinCos( sin2O, cos2O );
00233     L2.SinCos( sin2L, cos2L );
00234     M2.SinCos( sin2M, cos2M );
00235     
00236 //  deltaEcLong = ( -17.2*sinO - 1.32*sin2L - 0.23*sin2M + 0.21*sin2O)/3600.0; //Ecl. long. correction
00237 //  deltaObliquity = ( 9.2*cosO + 0.57*cos2L + 0.10*cos2M - 0.09*cos2O)/3600.0; //Obliq. correction
00238 
00239     deltaEcLong = 0.;
00240     deltaObliquity = 0.;
00241     
00242     for (unsigned int i=0; i < NUTTERMS; i++) {
00243 
00244         arg.setD ( arguments[i][0]*D.Degrees() + arguments[i][1]*M.Degrees() + 
00245             arguments[i][2]*MM.Degrees() + arguments[i][3]*F.Degrees() + arguments[i][4]*O.Degrees() );
00246         arg.SinCos( args, argc );
00247 
00248         deltaEcLong +=  (amp[i][0] + amp[i][1]/10. * T ) * args * 1e-4 ;
00249         deltaObliquity += (amp[i][2] + amp[i][3]/10. * T ) * argc * 1e-4 ;
00250     }
00251 
00252     deltaEcLong/= 3600.0;
00253     deltaObliquity /= 3600.0;
00254 
00255     //Compute Precession Matrices:
00256     XP.setD( 0.6406161*T + 0.0000839*T2 + 0.0000050*T3 );
00257     YP.setD( 0.5567530*T - 0.0001185*T2 - 0.0000116*T3 );
00258     ZP.setD( 0.6406161*T + 0.0003041*T2 + 0.0000051*T3 );
00259 
00260     XP.SinCos( SX, CX );
00261     YP.SinCos( SY, CY );
00262     ZP.SinCos( SZ, CZ );
00263     
00264 //P1 is used to precess from any epoch to J2000
00265     P1[0][0] = CX*CY*CZ - SX*SZ;
00266     P1[1][0] = CX*CY*SZ + SX*CZ;
00267     P1[2][0] = CX*SY;
00268     P1[0][1] = -1.0*SX*CY*CZ - CX*SZ;
00269     P1[1][1] = -1.0*SX*CY*SZ + CX*CZ;
00270     P1[2][1] = -1.0*SX*SY;
00271     P1[0][2] = -1.0*SY*CZ;
00272     P1[1][2] = -1.0*SY*SZ;
00273     P1[2][2] = CY;
00274 
00275 //P2 is used to precess from J2000 to any other epoch (it is the transpose of P1)
00276     P2[0][0] = CX*CY*CZ - SX*SZ;
00277     P2[1][0] = -1.0*SX*CY*CZ - CX*SZ;
00278     P2[2][0] = -1.0*SY*CZ;
00279     P2[0][1] = CX*CY*SZ + SX*CZ;
00280     P2[1][1] = -1.0*SX*CY*SZ + CX*CZ;
00281     P2[2][1] = -1.0*SY*SZ;
00282     P2[0][2] = CX*SY;
00283     P2[1][2] = -1.0*SX*SY;
00284     P2[2][2] = CY;
00285 
00286 //Compute Precession Matrices from B1950 to 1984 using Newcomb formulae
00287 
00288     XB.setD( 0.217697 );
00289     YB.setD( 0.189274 );
00290     ZB.setD( 0.217722 );
00291 
00292     XB.SinCos( SXB, CXB );
00293     YB.SinCos( SYB, CYB );
00294     ZB.SinCos( SZB, CZB );
00295 
00296 //P1B is used to precess from 1984 to B1950:
00297 
00298     P1B[0][0] = CXB*CYB*CZB - SXB*SZB;
00299     P1B[1][0] = CXB*CYB*SZB + SXB*CZB;
00300     P1B[2][0] = CXB*SYB;
00301     P1B[0][1] = -1.0*SXB*CYB*CZB - CXB*SZB;
00302     P1B[1][1] = -1.0*SXB*CYB*SZB + CXB*CZB;
00303     P1B[2][1] = -1.0*SXB*SYB;
00304     P1B[0][2] = -1.0*SYB*CZB;
00305     P1B[1][2] = -1.0*SYB*SZB;
00306     P1B[2][2] = CYB;
00307 
00308 //P2 is used to precess from B1950 to 1984 (it is the transpose of P1)
00309     P2B[0][0] = CXB*CYB*CZB - SXB*SZB;
00310     P2B[1][0] = -1.0*SXB*CYB*CZB - CXB*SZB;
00311     P2B[2][0] = -1.0*SYB*CZB;
00312     P2B[0][1] = CXB*CYB*SZB + SXB*CZB;
00313     P2B[1][1] = -1.0*SXB*CYB*SZB + CXB*CZB;
00314     P2B[2][1] = -1.0*SYB*SZB;
00315     P2B[0][2] = CXB*SYB;
00316     P2B[1][2] = -1.0*SXB*SYB;
00317     P2B[2][2] = CYB;
00318 
00319     
00320     // Mean longitudes for the planets. radians
00321     //
00322 
00323     // TODO Pasar a grados
00324     double LVenus   = 3.1761467+1021.3285546*T; // Venus
00325     double LMars    = 1.7534703+ 628.3075849*T; // Mars
00326     double LEarth   = 6.2034809+ 334.0612431*T; // Earth
00327     double LJupiter = 0.5995465+  52.9690965*T; // Jupiter
00328     double LSaturn  = 0.8740168+  21.3299095*T; // Saturn
00329     double LNeptune = 5.3118863+   3.8133036*T; // Neptune
00330     double LUranus  = 5.4812939+   7.4781599*T; // Uranus
00331     
00332     double LMRad = 3.8103444+8399.6847337*T; // Moon
00333     double DRad  = 5.1984667+7771.3771486*T;
00334     double MMRad = 2.3555559+8328.6914289*T; // Moon
00335     double FRad  = 1.6279052+8433.4661601*T;
00336 
00343     double vondrak[36][7] = {
00344     {LMars,             -1719914-2*T,        -25,   25-13*T,1578089+156*T,   10+32*T,684185-358*T},
00345     {2*LMars,             6434+141*T,28007-107*T,25697-95*T,  -5904-130*T,11141-48*T,  -2559-55*T},
00346     {LJupiter,                   715,           0,        6,         -657,       -15,        -282},
00347     {LMRad,                      715,           0,        0,         -656,         0,        -285},
00348     {3*LMars,                486-5*T,    -236-4*T, -216-4*T,     -446+5*T,       -94,        -193},
00349     {LSaturn,                    159,           0,        2,         -147,        -6,         -61},
00350     {FRad,                         0,           0,        0,           26,         0,         -59},
00351     {LMRad+MMRad,                 39,           0,        0,          -36,         0,         -16},
00352     {2*LJupiter,                  33,         -10,       -9,          -30,        -5,         -13},
00353     {2*LMars-LJupiter,            31,           1,        1,          -28,         0,         -12},
00354     {3*LMars-8*LEarth+3*LJupiter,  8,         -28,       25,            8,        11,           3},
00355     {5*LMars-8*LEarth+3*LJupiter,  8,         -28,      -25,           -8,       -11,          -3},
00356     {2*LVenus-LMars,              21,           0,        0,          -19,         0,          -8},
00357     {LVenus,                     -19,           0,        0,           17,         0,           8},
00358     {LNeptune,                    17,           0,        0,          -16,         0,          -7},
00359     {LMars-2*LJupiter,            16,           0,        0,           15,         1,           7},
00360     {LUranus,                     16,           0,        1,          -15,        -3,          -6},
00361     {LMars+LJupiter,              11,          -1,       -1,          -10,        -1,          -5},
00362     {2*LVenus-2*LMars,             0,         -11,      -10,            0,        -4,           0},
00363     {LMars-LJupiter,             -11,          -2,       -2,            9,        -1,           4},
00364     {4*LMars,                     -7,          -8,       -8,            6,        -3,           3},
00365     {3*LMars-2*LJupiter,         -10,           0,        0,            9,         0,           4},
00366     {LVenus-2*LMars,              -9,           0,        0,           -9,         0,          -4},
00367     {2*LVenus-3*LMars,            -9,           0,        0,           -8,         0,          -4},
00368     {2*LSaturn,                    0,          -9,       -8,            0,        -3,           0},
00369     {2*LVenus-4*LMars,             0,          -9,        8,            0,         3,           0},
00370     {3*LMars-2*LEarth,             8,           0,        0,           -8,         0,          -3},
00371     {LMRad+2*DRad-MMRad,           8,           0,        0,           -7,         0,          -3},
00372     {8*LVenus-12*LMars,           -4,          -7,       -6,            4,        -3,           2},
00373     {8*LVenus-14*LMars,           -4,          -7,        6,           -4,         3,          -2},
00374     {2*LEarth,                    -6,          -5,       -4,            5,        -2,           2},
00375     {3*LVenus-4*LMars,            -1,          -1,       -2,           -7,         1,          -4},
00376     {2*LMars-2*LJupiter,           4,          -6,       -5,           -4,        -2,          -2},
00377     {3*LVenus-3*LMars,             0,          -7,       -6,            0,        -3,           0},
00378     {2*LMars-2*LEarth,             5,          -5,       -4,           -5,        -2,          -2},
00379     {LMRad-2*DRad,                 5,           0,        0,           -5,         0,          -2}
00380     };
00381 
00382     dms anglev;
00383     double sa, ca;
00384     // Vearth X component
00385     vearth[0] = 0.;  
00386     // Vearth Y component
00387     vearth[1] = 0.; 
00388     // Vearth Z component
00389     vearth[2] = 0.;
00390 
00391     for (unsigned int i=0; i<36; i++) {
00392         anglev.setRadians(vondrak[i][0]);
00393         anglev.SinCos(sa,ca);
00394         for (unsigned int j=0; j<3; j++) {
00395             vearth[j] += vondrak[i][2*j+1]*sa +vondrak[i][2*j+2]*ca;
00396         }
00397     }
00398 
00399     const double UA2km  =  1.49597870/86400.;     // 10^{-8}*UA/dia -> km/s
00400 
00401     for (unsigned int j=0; j<3; j++) {
00402         vearth[j] = vearth[j] * UA2km;
00403     }
00404 }
kstars

ksnumbers.cpp

kstars

API Reference
