/* copyright 2020, Chen Wei version 0.0.3 Implement astronomical algorithms for finding solar terms and moon phases. Truncated IAU2000B from USNO NOVAS c source is used for nutation. The higher accuracy light abberation algorithm from A & A. p156, */ #include #include #include "astro.h" /* 77 terms IAU2000B Luni-Solar nutation table from USNO NOVAS c source Luni-Solar argument multipliers, coefficients, unit 1e-7 arcsec L L' F D Om longitude (sin, t*sin, cos), obliquity (cos, t*cos, sin) */ static double IAU2000BNutationTable[77][11] = { { 0, 0, 0, 0, 1, -172064161, -174666, 33386, 92052331, 9086, 15377 }, { 0, 0, 2, -2, 2, -13170906, -1675, -13696, 5730336, -3015, -4587 }, { 0, 0, 2, 0, 2, -2276413, -234, 2796, 978459, -485, 1374 }, { 0, 0, 0, 0, 2, 2074554, 207, -698, -897492, 470, -291 }, { 0, 1, 0, 0, 0, 1475877, -3633, 11817, 73871, -184, -1924 }, { 0, 1, 2, -2, 2, -516821, 1226, -524, 224386, -677, -174 }, { 1, 0, 0, 0, 0, 711159, 73, -872, -6750, 0, 358 }, { 0, 0, 2, 0, 1, -387298, -367, 380, 200728, 18, 318 }, { 1, 0, 2, 0, 2, -301461, -36, 816, 129025, -63, 367 }, { 0, -1, 2, -2, 2, 215829, -494, 111, -95929, 299, 132 }, { 0, 0, 2, -2, 1, 128227, 137, 181, -68982, -9, 39 }, { -1, 0, 2, 0, 2, 123457, 11, 19, -53311, 32, -4 }, { -1, 0, 0, 2, 0, 156994, 10, -168, -1235, 0, 82 }, { 1, 0, 0, 0, 1, 63110, 63, 27, -33228, 0, -9 }, { -1, 0, 0, 0, 1, -57976, -63, -189, 31429, 0, -75 }, { -1, 0, 2, 2, 2, -59641, -11, 149, 25543, -11, 66 }, { 1, 0, 2, 0, 1, -51613, -42, 129, 26366, 0, 78 }, { -2, 0, 2, 0, 1, 45893, 50, 31, -24236, -10, 20 }, { 0, 0, 0, 2, 0, 63384, 11, -150, -1220, 0, 29 }, { 0, 0, 2, 2, 2, -38571, -1, 158, 16452, -11, 68 }, { 0, -2, 2, -2, 2, 32481, 0, 0, -13870, 0, 0 }, { -2, 0, 0, 2, 0, -47722, 0, -18, 477, 0, -25 }, { 2, 0, 2, 0, 2, -31046, -1, 131, 13238, -11, 59 }, { 1, 0, 2, -2, 2, 28593, 0, -1, -12338, 10, -3 }, { -1, 0, 2, 0, 1, 20441, 21, 10, -10758, 0, -3 }, { 2, 0, 0, 0, 0, 29243, 0, -74, -609, 0, 13 }, { 0, 0, 2, 0, 0, 25887, 0, -66, -550, 0, 11 }, { 0, 1, 0, 0, 1, -14053, -25, 79, 8551, -2, -45 }, { -1, 0, 0, 2, 1, 15164, 10, 11, -8001, 0, -1 }, { 0, 2, 2, -2, 2, -15794, 72, -16, 6850, -42, -5 }, { 0, 0, -2, 2, 0, 21783, 0, 13, -167, 0, 13 }, { 1, 0, 0, -2, 1, -12873, -10, -37, 6953, 0, -14 }, { 0, -1, 0, 0, 1, -12654, 11, 63, 6415, 0, 26 }, { -1, 0, 2, 2, 1, -10204, 0, 25, 5222, 0, 15 }, { 0, 2, 0, 0, 0, 16707, -85, -10, 168, -1, 10 }, { 1, 0, 2, 2, 2, -7691, 0, 44, 3268, 0, 19 }, { -2, 0, 2, 0, 0, -11024, 0, -14, 104, 0, 2 }, { 0, 1, 2, 0, 2, 7566, -21, -11, -3250, 0, -5 }, { 0, 0, 2, 2, 1, -6637, -11, 25, 3353, 0, 14 }, { 0, -1, 2, 0, 2, -7141, 21, 8, 3070, 0, 4 }, { 0, 0, 0, 2, 1, -6302, -11, 2, 3272, 0, 4 }, { 1, 0, 2, -2, 1, 5800, 10, 2, -3045, 0, -1 }, { 2, 0, 2, -2, 2, 6443, 0, -7, -2768, 0, -4 }, { -2, 0, 0, 2, 1, -5774, -11, -15, 3041, 0, -5 }, { 2, 0, 2, 0, 1, -5350, 0, 21, 2695, 0, 12 }, { 0, -1, 2, -2, 1, -4752, -11, -3, 2719, 0, -3 }, { 0, 0, 0, -2, 1, -4940, -11, -21, 2720, 0, -9 }, { -1, -1, 0, 2, 0, 7350, 0, -8, -51, 0, 4 }, { 2, 0, 0, -2, 1, 4065, 0, 6, -2206, 0, 1 }, { 1, 0, 0, 2, 0, 6579, 0, -24, -199, 0, 2 }, { 0, 1, 2, -2, 1, 3579, 0, 5, -1900, 0, 1 }, { 1, -1, 0, 0, 0, 4725, 0, -6, -41, 0, 3 }, { -2, 0, 2, 0, 2, -3075, 0, -2, 1313, 0, -1 }, { 3, 0, 2, 0, 2, -2904, 0, 15, 1233, 0, 7 }, { 0, -1, 0, 2, 0, 4348, 0, -10, -81, 0, 2 }, { 1, -1, 2, 0, 2, -2878, 0, 8, 1232, 0, 4 }, { 0, 0, 0, 1, 0, -4230, 0, 5, -20, 0, -2 }, { -1, -1, 2, 2, 2, -2819, 0, 7, 1207, 0, 3 }, { -1, 0, 2, 0, 0, -4056, 0, 5, 40, 0, -2 }, { 0, -1, 2, 2, 2, -2647, 0, 11, 1129, 0, 5 }, { -2, 0, 0, 0, 1, -2294, 0, -10, 1266, 0, -4 }, { 1, 1, 2, 0, 2, 2481, 0, -7, -1062, 0, -3 }, { 2, 0, 0, 0, 1, 2179, 0, -2, -1129, 0, -2 }, { -1, 1, 0, 1, 0, 3276, 0, 1, -9, 0, 0 }, { 1, 1, 0, 0, 0, -3389, 0, 5, 35, 0, -2 }, { 1, 0, 2, 0, 0, 3339, 0, -13, -107, 0, 1 }, { -1, 0, 2, -2, 1, -1987, 0, -6, 1073, 0, -2 }, { 1, 0, 0, 0, 2, -1981, 0, 0, 854, 0, 0 }, { -1, 0, 0, 1, 0, 4026, 0, -353, -553, 0, -139 }, { 0, 0, 2, 1, 2, 1660, 0, -5, -710, 0, -2 }, { -1, 0, 2, 4, 2, -1521, 0, 9, 647, 0, 4 }, { -1, 1, 0, 1, 1, 1314, 0, 0, -700, 0, 0 }, { 0, -2, 2, -2, 1, -1283, 0, 0, 672, 0, 0 }, { 1, 0, 2, 2, 1, -1331, 0, 8, 663, 0, 4 }, { -2, 0, 2, 2, 2, 1383, 0, -2, -594, 0, -2 }, { -1, 0, 0, 0, 2, 1405, 0, 4, -610, 0, 2 }, { 1, 1, 2, -2, 2, 1290, 0, 0, -556, 0, 0 } }; /* Calculate nutation angles using the IAU 2000B model. The table is from NOVAS(USNO) c source flle. "...the 77-term, IAU-approved truncated nutation series, IAU 2000B, which is accurate to about 0.001 arcsecond in the interval 1995-2050" Arg: jd as jd Return: nutation of longitude in radians */ double nutation(double jd) { double t, L, Lp, F, D, Om; t = (jd - J2000) / 36525.0; /* Mean anomaly of the Moon, in arcsec */ L = 485868.249036 + t * 1717915923.2178; /* Mean anomaly of the Sun. */ Lp = 1287104.79305 + t * 129596581.0481; /* Mean argument of the latitude of the Moon. */ F = 335779.526232 + t * 1739527262.8478; /* Mean elongation of the Moon from the Sun. */ D = 1072260.70369 + t * 1602961601.2090; /* Mean longitude of the ascending node of the Moon. */ Om = 450160.398036 - t * 6962890.5431; double arg, lon, dpplan; lon = 0; int i; for (i = 0; i < 77; i++) { arg = ( IAU2000BNutationTable[i][0] * L + IAU2000BNutationTable[i][1] * Lp + IAU2000BNutationTable[i][2] * F + IAU2000BNutationTable[i][3] * D + IAU2000BNutationTable[i][4] * Om) * ASEC2RAD; lon += ( IAU2000BNutationTable[i][5] + IAU2000BNutationTable[i][6] * t) * sin(arg) + IAU2000BNutationTable[i][7] * cos(arg); } /* unit of longitude is 1.0e-7 arcsec, convert it to arcsec */ lon *= 1.0e-7; /* Constant account for the missing long-period planetary terms in the truncated nutation model, in arcsec */ dpplan = -0.000135; lon += dpplan; lon *= ASEC2RAD; /* lon = fmod(lon, TWOPI); */ return lon; } /* the higher accuracy light abberation algorithm from A & A. p156, * the error will be less than 0.001" */ double lightabbr_high(double jd) { double tm, tm2, tm3, t, t2, t3; t = (jd - J2000) / 36525.0; t2 = t * t; t3 = t * t2; tm = t / 10.0; tm2 = tm * tm; tm3 = tm * tm2; // variation of the Sun's longitude double var_lon; var_lon = 3548.330 + 118.568 * sin(1.527664004990 + 6283.075850600876 * tm) + 2.476 * sin(1.484508993906 + 12566.151699456424 * tm) + 1.376 * sin(0.486077687339 + 77713.771468687730 * tm) + 0.119 * sin(1.276490181677 + 7860.419392621536 * tm) + 0.114 * sin(5.885711004647 + 5753.384884566103 * tm) + 0.086 * sin(3.884055717388 + 11506.769769132206 * tm) + 0.078 * sin(2.841633387025 + 161000.685738008644 * tm) + 0.054 * sin(1.441333038870 + 18849.227550057301 * tm) + 0.052 * sin(2.993569534399 + 3930.209696310768 * tm) + 0.034 * sin(0.529208263814 + 71430.695618086858 * tm) + 0.033 * sin(2.091087703461 + 5884.926847413107 * tm) + 0.023 * sin(4.320419446313 + 5223.693920276658 * tm) + 0.023 * sin(5.674983441420 + 5507.553238331428 * tm) + 0.021 * sin(2.707426294193 + 11790.629088932305 * tm) + 7.311 * tm * sin(5.819826570714 + 6283.075850600876 * tm) + 0.305 * tm * sin(5.776715192860 + 12566.151699456424 * tm) + 0.010 * tm * sin(5.733703298774 + 18849.227550057301 * tm) + 0.309 * tm2 * sin(4.214128894867 + 6283.075850600876 * tm) + 0.021 * tm2 * sin(3.578766215288 + 12566.151699456424 * tm) + 0.004 * tm3 * sin(5.198655163283 + 77713.771468687730 * tm) + 0.010 * tm3 * sin(2.700139544566 + 6283.075850600876 * tm); double M, e, C, v, R, res; // mean anomaly of the Sun M = (357.52910 + 35999.0503 * t - 0.0001559 * t2 - 0.00000048 * t3) * DEG2RAD; // the eccentricity of Earth's orbit e = 0.016708617 - 0.000042037 * t - 0.0000001236 * t2; // Sun's equation of center C = ( (1.9146 - 0.004817 * t - 0.000014 * t2) * sin(M) + (0.019993 - 0.000101 * t) * sin(2 * M) + 0.00029 * sin(3 * M)) * DEG2RAD; // true anomaly v = M + C; // Sun's distance from the Earth, in AU R = (1.000001018 * (1 - e * e)) / (1 + e * cos(v)); res = -0.005775518 * R * var_lon * ASEC2RAD; return res; }