00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "ksnumbers.h"
00021
00022
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. );
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
00171 T = ( jd - J2000 ) / 36525.;
00172
00173
00174 jm = T / 10.0;
00175
00176 double T2 = T*T;
00177 double T3 = T2*T;
00178
00179
00180 L.setD( 280.46645 + 36000.76983*T + 0.0003032*T2 );
00181
00182
00183 D.setD( 297.85036 + 445267.111480*T - 0.0019142*T2 + T3/189474.);
00184
00185
00186 M.setD( 357.52910 + 35999.05030*T - 0.0001559*T2 - 0.00000048*T3);
00187
00188
00189 MM.setD( 134.96298 + 477198.867398*T + 0.0086972*T2 + T3/56250.0 );
00190
00191
00192 LM.setD( 218.3164591 + 481267.88134236*T - 0.0013268*T2 + T3/538841. - T*T*T*T/6519400.);
00193
00194
00195 F.setD( 93.27191 + 483202.017538*T - 0.0036825*T2 + T3/327270.);
00196
00197
00198 O.setD( 125.04452 - 1934.136261*T + 0.0020708*T2 + T3/450000.0 );
00199
00200
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
00208 L0.setD( L.Degrees() + C );
00209
00210
00211 M0.setD( M.Degrees() + C );
00212
00213
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
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() );
00229 M2.setD( 2.0*LM.Degrees() );
00230
00231 O.SinCos( sinO, cosO );
00232 O2.SinCos( sin2O, cos2O );
00233 L2.SinCos( sin2L, cos2L );
00234 M2.SinCos( sin2M, cos2M );
00235
00236
00237
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
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
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
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
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
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
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
00321
00322
00323
00324 double LVenus = 3.1761467+1021.3285546*T;
00325 double LMars = 1.7534703+ 628.3075849*T;
00326 double LEarth = 6.2034809+ 334.0612431*T;
00327 double LJupiter = 0.5995465+ 52.9690965*T;
00328 double LSaturn = 0.8740168+ 21.3299095*T;
00329 double LNeptune = 5.3118863+ 3.8133036*T;
00330 double LUranus = 5.4812939+ 7.4781599*T;
00331
00332 double LMRad = 3.8103444+8399.6847337*T;
00333 double DRad = 5.1984667+7771.3771486*T;
00334 double MMRad = 2.3555559+8328.6914289*T;
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
00385 vearth[0] = 0.;
00386
00387 vearth[1] = 0.;
00388
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.;
00400
00401 for (unsigned int j=0; j<3; j++) {
00402 vearth[j] = vearth[j] * UA2km;
00403 }
00404 }