diff --git a/generate/codegen.c b/generate/codegen.c index 8ca457a1..3e23121d 100644 --- a/generate/codegen.c +++ b/generate/codegen.c @@ -31,6 +31,8 @@ #include "ephfile.h" #define CG_MAX_LINE_LENGTH 200 +#define MAX_DATA_PER_LINE 20 +#define IAU_DATA_PER_ROW 11 static const double MJD_BASIS = 2400000.5; @@ -645,6 +647,279 @@ fail: return error; } +static int ScanRealArray( + cg_context_t *context, + const char *filename, + int lnum, + char *line, + int numExpected, + double *data) +{ + int i, t, len, inspace; + char *token[MAX_DATA_PER_LINE]; + + if (numExpected < 1 || numExpected > MAX_DATA_PER_LINE) + return LogError(context, "Invalid value for numExpected=%d\n", numExpected); + + /* Split the line into space delimited tokens. */ + len = strlen(line); + inspace = 1; + t = 0; + for (i=0; i < len; ++i) + { + if (line[i] == ' ' || line[i] == '\t' || line[i] == '\r' || line[i] == '\n') + { + if (!inspace) + { + line[i] = '\0'; /* terminate the previous token */ + inspace = 1; + } + } + else + { + if (inspace) + { + /* we just found the beginning of a new token. */ + if (t < IAU_DATA_PER_ROW) + { + token[t++] = &line[i]; + inspace = 0; + } + else + return LogError(context, "ScanRealArray(%s %d): too many data on line.", filename, lnum); + } + } + } + + /* Verify there are the correct number of tokens. */ + if (t != numExpected) + return LogError(context, "ScanRealArray(%s %d): found %d data, but expected %d\n", filename, lnum, t, numExpected); + + /* Parse each token as a floating point number. */ + for (t=0; t < numExpected; ++t) + if (1 != sscanf(token[t], "%lf", &data[t])) + return LogError(context, "ScanRealArray(%s %d): invalid floating point token '%s'\n", filename, lnum, token[t]); + + return 0; /* successful parse */ +} + +static int OptimizeConst(cg_context_t *context, char *buffer, size_t size, double c, const char *v) +{ + int nprinted; + const char *op; + + if (c == 0.0) + { + buffer[0] = '\0'; + return 0; + } + + if (c < 0.0) + { + op = " - "; + c *= -1.0; + } + else + op = " + "; + + if (c == 1.0) + nprinted = snprintf(buffer, size, "%s%s", op, v); + else + nprinted = snprintf(buffer, size, "%s%0.1lf*%s", op, c, v); + + if (nprinted >= (int)size) + return LogError(context, "OptimizeConst: print buffer overflowed."); + + return 0; +} + +static int OptimizeLinear(cg_context_t *context, char *buffer, size_t size, double a, double b) +{ + int nprinted; + + if (b == 0.0) + nprinted = snprintf(buffer, size, "%0.1lf", a); + else if (a == 0.0) + nprinted = snprintf(buffer, size, "%0.1lf*t", b); + else if (b < 0.0) + nprinted = snprintf(buffer, size, "%0.1lf - %0.1lf*t", a, -b); + else + nprinted = snprintf(buffer, size, "%0.1lf + %0.1lf*t", a, b); + + if (nprinted >= size) + return LogError(context, "OptimizeLinear: print buffer overflowed."); + + return 0; +} + +static int OptIauPython(cg_context_t *context, const double *data) +{ + static const char * const nv[] = {"el", "elp", "f", "d", "om"}; /* variable names */ + int first; + double n; + int i; + const char *op; + char dotprod[200]; + char term[40]; + char linear[100]; + char cpart[100]; + int nprinted; + int lonevar; + + fprintf(context->outfile, "\n"); /* must start on new line to maintain correct indentation in Python source. */ + + /* Optimize dot product of data[0]..data[4] with nv[]. */ + first = 1; + dotprod[0] = '\0'; + for (i=0; i<5; ++i) + { + n = data[i]; + if (n != 0.0) + { + if (n < 0.0) + { + n *= -1.0; + op = first ? "-" : " - "; + } + else + op = first ? "" : " + "; + + if (n == 1.0) + nprinted = snprintf(term, sizeof(term), "%s%s", op, nv[i]); + else + nprinted = snprintf(term, sizeof(term), "%s%0.1lf*%s", op, n, nv[i]); + + if (nprinted >= sizeof(term)) + return LogError(context, "Truncated iau2000b term."); + + if (nprinted + strlen(dotprod) >= sizeof(dotprod)) + return LogError(context, "Dot product overflow in iau2000b formula."); + + strcat(dotprod, term); + first = 0; + } + } + + /* Did we print a lone variable, e.g. "elp"? */ + lonevar = 0; + for (i=0; i<5 && !lonevar; ++i) + if (!strcmp(dotprod, nv[i])) + lonevar = 1; + + if (lonevar) + { + fprintf(context->outfile, " sarg = math.sin(%s)\n", dotprod); + fprintf(context->outfile, " carg = math.cos(%s)\n", dotprod); + } + else + { + fprintf(context->outfile, " arg = %s\n", dotprod); + fprintf(context->outfile, " sarg = math.sin(arg)\n"); + fprintf(context->outfile, " carg = math.cos(arg)\n"); + } + + if (OptimizeLinear(context, linear, sizeof(linear), data[5], data[6])) return 1; + if (OptimizeConst(context, cpart, sizeof(cpart), data[7], "carg")) return 1; + fprintf(context->outfile, " dp += (%s)*sarg%s\n", linear, cpart); + + if (OptimizeLinear(context, linear, sizeof(linear), data[8], data[9])) return 1; + if (OptimizeConst(context, cpart, sizeof(cpart), data[10], "sarg")) return 1; + fprintf(context->outfile, " de += (%s)*carg%s\n", linear, cpart); + + fprintf(context->outfile, "\n"); + return 0; +} + +static int OptIauData(cg_context_t *context) +{ + int error = 1; + int lnum; + FILE *infile; + const char *filename; + char line[100]; + double data[IAU_DATA_PER_ROW]; + + filename = "model_data/iau2000b.txt"; + infile = fopen(filename, "rt"); + if (infile == NULL) goto fail; + + lnum = 0; + while (fgets(line, sizeof(line), infile)) + { + ++lnum; + CHECK(ScanRealArray(context, filename, lnum, line, IAU_DATA_PER_ROW, data)); + switch (context->language) + { + case CODEGEN_LANGUAGE_PYTHON: + CHECK(OptIauPython(context, data)); + break; + + default: + error = LogError(context, "OptIauData: Unsupported language %d", context->language); + goto fail; + } + } + + error = 0; +fail: + if (infile == NULL) + error = LogError(context, "Cannot open input file: %s", filename); + else + fclose(infile); + + return error; +} + +static int OptAddSol(cg_context_t *context) +{ + int nscanned; + double cl, cs, cg, cp, p, q, r, s; + const char *op; + + nscanned = sscanf(context->args, "%lf , %lf , %lf , %lf , %lf , %lf , %lf , %lf", &cl, &cs, &cg, &cp, &p, &q, &r, &s); + if (nscanned != 8) + return LogError(context, "OptAddSol: invalid arguments: '%s'", context->args); + + fprintf(context->outfile, "\n # AddSol(%s)\n", context->args); + + op = ""; + fprintf(context->outfile, " z = "); + if (p != 0.0) + { + fprintf(context->outfile, "ex[%0.0lf][1]", p); + op = " * "; + } + if (q != 0.0) + { + fprintf(context->outfile, "%sex[%0.0lf][2]", op, q); + op = " * "; + } + if (r != 0.0) + { + fprintf(context->outfile, "%sex[%0.0lf][3]", op, r); + op = " * "; + } + if (s != 0.0) + { + fprintf(context->outfile, "%sex[%0.0lf][4]", op, s); + } + fprintf(context->outfile, "\n"); + + if (cl != 0.0) + fprintf(context->outfile, " DLAM += %0.3lf * z.imag\n", cl); + + if (cs != 0.0) + fprintf(context->outfile, " DS += %0.2lf * z.imag\n", cs); + + if (cg != 0.0) + fprintf(context->outfile, " GAM1C += %0.3lf * z.real\n", cg); + + if (cp != 0.0) + fprintf(context->outfile, " SINPI += %0.4lf * z.real\n", cp); + + return 0; +} + static int LogError(const cg_context_t *context, const char *format, ...) { va_list v; @@ -663,6 +938,8 @@ static const cg_directive_entry DirectiveTable[] = { "LIST_CHEBYSHEV", ListChebyshev }, { "C_CHEBYSHEV", CChebyshev }, { "DELTA_T", GenDeltaT }, + { "IAU_DATA", OptIauData }, + { "ADDSOL", OptAddSol }, { NULL, NULL } }; diff --git a/generate/model_data/iau2000b.txt b/generate/model_data/iau2000b.txt new file mode 100644 index 00000000..0c712f1b --- /dev/null +++ b/generate/model_data/iau2000b.txt @@ -0,0 +1,77 @@ + 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 diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index e9980935..d70870fe 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -228,464 +228,7 @@ class _iau2000b: om = math.fmod((450160.398036 - t*6962890.5431), _ASEC360) * _ASEC2RAD dp = 0 de = 0 - sarg = math.sin(om) - carg = math.cos(om) - dp += (-172064161.0 - 174666.0*t)*sarg + 33386.0*carg - de += (92052331.0 + 9086.0*t)*carg + 15377.0*sarg - - arg = 2.0*f - 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-13170906.0 - 1675.0*t)*sarg - 13696.0*carg - de += (5730336.0 - 3015.0*t)*carg - 4587.0*sarg - - arg = 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-2276413.0 - 234.0*t)*sarg + 2796.0*carg - de += (978459.0 - 485.0*t)*carg + 1374.0*sarg - - arg = 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (2074554.0 + 207.0*t)*sarg - 698.0*carg - de += (-897492.0 + 470.0*t)*carg - 291.0*sarg - - sarg = math.sin(elp) - carg = math.cos(elp) - dp += (1475877.0 - 3633.0*t)*sarg + 11817.0*carg - de += (73871.0 - 184.0*t)*carg - 1924.0*sarg - - arg = elp + 2.0*f - 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-516821.0 + 1226.0*t)*sarg - 524.0*carg - de += (224386.0 - 677.0*t)*carg - 174.0*sarg - - sarg = math.sin(el) - carg = math.cos(el) - dp += (711159.0 + 73.0*t)*sarg - 872.0*carg - de += (-6750.0)*carg + 358.0*sarg - - arg = 2.0*f + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-387298.0 - 367.0*t)*sarg + 380.0*carg - de += (200728.0 + 18.0*t)*carg + 318.0*sarg - - arg = el + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-301461.0 - 36.0*t)*sarg + 816.0*carg - de += (129025.0 - 63.0*t)*carg + 367.0*sarg - - arg = -elp + 2.0*f - 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (215829.0 - 494.0*t)*sarg + 111.0*carg - de += (-95929.0 + 299.0*t)*carg + 132.0*sarg - - arg = 2.0*f - 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (128227.0 + 137.0*t)*sarg + 181.0*carg - de += (-68982.0 - 9.0*t)*carg + 39.0*sarg - - arg = -el + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (123457.0 + 11.0*t)*sarg + 19.0*carg - de += (-53311.0 + 32.0*t)*carg - 4.0*sarg - - arg = -el + 2.0*d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (156994.0 + 10.0*t)*sarg - 168.0*carg - de += (-1235.0)*carg + 82.0*sarg - - arg = el + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (63110.0 + 63.0*t)*sarg + 27.0*carg - de += (-33228.0)*carg - 9.0*sarg - - arg = -el + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-57976.0 - 63.0*t)*sarg - 189.0*carg - de += (31429.0)*carg - 75.0*sarg - - arg = -el + 2.0*f + 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-59641.0 - 11.0*t)*sarg + 149.0*carg - de += (25543.0 - 11.0*t)*carg + 66.0*sarg - - arg = el + 2.0*f + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-51613.0 - 42.0*t)*sarg + 129.0*carg - de += (26366.0)*carg + 78.0*sarg - - arg = -2.0*el + 2.0*f + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (45893.0 + 50.0*t)*sarg + 31.0*carg - de += (-24236.0 - 10.0*t)*carg + 20.0*sarg - - arg = 2.0*d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (63384.0 + 11.0*t)*sarg - 150.0*carg - de += (-1220.0)*carg + 29.0*sarg - - arg = 2.0*f + 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-38571.0 - 1.0*t)*sarg + 158.0*carg - de += (16452.0 - 11.0*t)*carg + 68.0*sarg - - arg = -2.0*elp + 2.0*f - 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (32481.0)*sarg - de += (-13870.0)*carg - - arg = -2.0*el + 2.0*d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-47722.0)*sarg - 18.0*carg - de += (477.0)*carg - 25.0*sarg - - arg = 2.0*el + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-31046.0 - 1.0*t)*sarg + 131.0*carg - de += (13238.0 - 11.0*t)*carg + 59.0*sarg - - arg = el + 2.0*f - 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (28593.0)*sarg - carg - de += (-12338.0 + 10.0*t)*carg - 3.0*sarg - - arg = -el + 2.0*f + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (20441.0 + 21.0*t)*sarg + 10.0*carg - de += (-10758.0)*carg - 3.0*sarg - - arg = 2.0*el - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (29243.0)*sarg - 74.0*carg - de += (-609.0)*carg + 13.0*sarg - - arg = 2.0*f - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (25887.0)*sarg - 66.0*carg - de += (-550.0)*carg + 11.0*sarg - - arg = elp + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-14053.0 - 25.0*t)*sarg + 79.0*carg - de += (8551.0 - 2.0*t)*carg - 45.0*sarg - - arg = -el + 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (15164.0 + 10.0*t)*sarg + 11.0*carg - de += (-8001.0)*carg - sarg - - arg = 2.0*elp + 2.0*f - 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-15794.0 + 72.0*t)*sarg - 16.0*carg - de += (6850.0 - 42.0*t)*carg - 5.0*sarg - - arg = -2.0*f + 2.0*d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (21783.0)*sarg + 13.0*carg - de += (-167.0)*carg + 13.0*sarg - - arg = el - 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-12873.0 - 10.0*t)*sarg - 37.0*carg - de += (6953.0)*carg - 14.0*sarg - - arg = -elp + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-12654.0 + 11.0*t)*sarg + 63.0*carg - de += (6415.0)*carg + 26.0*sarg - - arg = -el + 2.0*f + 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-10204.0)*sarg + 25.0*carg - de += (5222.0)*carg + 15.0*sarg - - arg = 2.0*elp - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (16707.0 - 85.0*t)*sarg - 10.0*carg - de += (168.0 - 1.0*t)*carg + 10.0*sarg - - arg = el + 2.0*f + 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-7691.0)*sarg + 44.0*carg - de += (3268.0)*carg + 19.0*sarg - - arg = -2.0*el + 2.0*f - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-11024.0)*sarg - 14.0*carg - de += (104.0)*carg + 2.0*sarg - - arg = elp + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (7566.0 - 21.0*t)*sarg - 11.0*carg - de += (-3250.0)*carg - 5.0*sarg - - arg = 2.0*f + 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-6637.0 - 11.0*t)*sarg + 25.0*carg - de += (3353.0)*carg + 14.0*sarg - - arg = -elp + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-7141.0 + 21.0*t)*sarg + 8.0*carg - de += (3070.0)*carg + 4.0*sarg - - arg = 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-6302.0 - 11.0*t)*sarg + 2.0*carg - de += (3272.0)*carg + 4.0*sarg - - arg = el + 2.0*f - 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (5800.0 + 10.0*t)*sarg + 2.0*carg - de += (-3045.0)*carg - sarg - - arg = 2.0*el + 2.0*f - 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (6443.0)*sarg - 7.0*carg - de += (-2768.0)*carg - 4.0*sarg - - arg = -2.0*el + 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-5774.0 - 11.0*t)*sarg - 15.0*carg - de += (3041.0)*carg - 5.0*sarg - - arg = 2.0*el + 2.0*f + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-5350.0)*sarg + 21.0*carg - de += (2695.0)*carg + 12.0*sarg - - arg = -elp + 2.0*f - 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-4752.0 - 11.0*t)*sarg - 3.0*carg - de += (2719.0)*carg - 3.0*sarg - - arg = -2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-4940.0 - 11.0*t)*sarg - 21.0*carg - de += (2720.0)*carg - 9.0*sarg - - arg = -el - elp + 2.0*d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (7350.0)*sarg - 8.0*carg - de += (-51.0)*carg + 4.0*sarg - - arg = 2.0*el - 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (4065.0)*sarg + 6.0*carg - de += (-2206.0)*carg + sarg - - arg = el + 2.0*d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (6579.0)*sarg - 24.0*carg - de += (-199.0)*carg + 2.0*sarg - - arg = elp + 2.0*f - 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (3579.0)*sarg + 5.0*carg - de += (-1900.0)*carg + sarg - - arg = el - elp - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (4725.0)*sarg - 6.0*carg - de += (-41.0)*carg + 3.0*sarg - - arg = -2.0*el + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-3075.0)*sarg - 2.0*carg - de += (1313.0)*carg - sarg - - arg = 3.0*el + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-2904.0)*sarg + 15.0*carg - de += (1233.0)*carg + 7.0*sarg - - arg = -elp + 2.0*d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (4348.0)*sarg - 10.0*carg - de += (-81.0)*carg + 2.0*sarg - - arg = el - elp + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-2878.0)*sarg + 8.0*carg - de += (1232.0)*carg + 4.0*sarg - - sarg = math.sin(d) - carg = math.cos(d) - dp += (-4230.0)*sarg + 5.0*carg - de += (-20.0)*carg - 2.0*sarg - - arg = -el - elp + 2.0*f + 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-2819.0)*sarg + 7.0*carg - de += (1207.0)*carg + 3.0*sarg - - arg = -el + 2.0*f - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-4056.0)*sarg + 5.0*carg - de += (40.0)*carg - 2.0*sarg - - arg = -elp + 2.0*f + 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-2647.0)*sarg + 11.0*carg - de += (1129.0)*carg + 5.0*sarg - - arg = -2.0*el + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-2294.0)*sarg - 10.0*carg - de += (1266.0)*carg - 4.0*sarg - - arg = el + elp + 2.0*f + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (2481.0)*sarg - 7.0*carg - de += (-1062.0)*carg - 3.0*sarg - - arg = 2.0*el + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (2179.0)*sarg - 2.0*carg - de += (-1129.0)*carg - 2.0*sarg - - arg = -el + elp + d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (3276.0)*sarg + carg - de += (-9.0)*carg - - arg = el + elp - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-3389.0)*sarg + 5.0*carg - de += (35.0)*carg - 2.0*sarg - - arg = el + 2.0*f - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (3339.0)*sarg - 13.0*carg - de += (-107.0)*carg + sarg - - arg = -el + 2.0*f - 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-1987.0)*sarg - 6.0*carg - de += (1073.0)*carg - 2.0*sarg - - arg = el + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-1981.0)*sarg - de += (854.0)*carg - - arg = -el + d - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (4026.0)*sarg - 353.0*carg - de += (-553.0)*carg - 139.0*sarg - - arg = 2.0*f + d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (1660.0)*sarg - 5.0*carg - de += (-710.0)*carg - 2.0*sarg - - arg = -el + 2.0*f + 4.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-1521.0)*sarg + 9.0*carg - de += (647.0)*carg + 4.0*sarg - - arg = -el + elp + d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (1314.0)*sarg - de += (-700.0)*carg - - arg = -2.0*elp + 2.0*f - 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-1283.0)*sarg - de += (672.0)*carg - - arg = el + 2.0*f + 2.0*d + om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (-1331.0)*sarg + 8.0*carg - de += (663.0)*carg + 4.0*sarg - - arg = -2.0*el + 2.0*f + 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (1383.0)*sarg - 2.0*carg - de += (-594.0)*carg - 2.0*sarg - - arg = -el + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (1405.0)*sarg + 4.0*carg - de += (-610.0)*carg + 2.0*sarg - - arg = el + elp + 2.0*f - 2.0*d + 2.0*om - sarg = math.sin(arg) - carg = math.cos(arg) - dp += (1290.0)*sarg - de += (-556.0)*carg - + $ASTRO_IAU_DATA() self.dpsi = -0.000135 + (dp * 1.0e-7) self.deps = +0.000388 + (de * 1.0e-7) @@ -988,694 +531,112 @@ def _CalcMoon(time): I += 1 - - # AddSol(13.902, 14.06, -0.001, 0.2607, 0, 0, 0, 4) - z = ex[4][4] - DLAM += 13.902 * z.imag - DS += 14.06 * z.imag - GAM1C += -0.001 * z.real - SINPI += 0.2607 * z.real - - # AddSol(0.403, -4.01, 0.394, 0.0023, 0, 0, 0, 3) - z = ex[3][4] - DLAM += 0.403 * z.imag - DS += -4.01 * z.imag - GAM1C += 0.394 * z.real - SINPI += 0.0023 * z.real - - # AddSol(2369.912, 2373.36, 0.601, 28.2333, 0, 0, 0, 2) - z = ex[2][4] - DLAM += 2369.912 * z.imag - DS += 2373.36 * z.imag - GAM1C += 0.601 * z.real - SINPI += 28.2333 * z.real - - # AddSol(-125.154, -112.79, -0.725, -0.9781, 0, 0, 0, 1) - z = ex[1][4] - DLAM += -125.154 * z.imag - DS += -112.79 * z.imag - GAM1C += -0.725 * z.real - SINPI += -0.9781 * z.real - - # AddSol(1.979, 6.98, -0.445, 0.0433, 1, 0, 0, 4) - z = ex[1][1] * ex[4][4] - DLAM += 1.979 * z.imag - DS += 6.98 * z.imag - GAM1C += -0.445 * z.real - SINPI += 0.0433 * z.real - - # AddSol(191.953, 192.72, 0.029, 3.0861, 1, 0, 0, 2) - z = ex[1][1] * ex[2][4] - DLAM += 191.953 * z.imag - DS += 192.72 * z.imag - GAM1C += 0.029 * z.real - SINPI += 3.0861 * z.real - - # AddSol(-8.466, -13.51, 0.455, -0.1093, 1, 0, 0, 1) - z = ex[1][1] * ex[1][4] - DLAM += -8.466 * z.imag - DS += -13.51 * z.imag - GAM1C += 0.455 * z.real - SINPI += -0.1093 * z.real - - # AddSol(22639.5, 22609.07, 0.079, 186.5398, 1, 0, 0, 0) - z = ex[1][1] - DLAM += 22639.5 * z.imag - DS += 22609.07 * z.imag - GAM1C += 0.079 * z.real - SINPI += 186.5398 * z.real - - # AddSol(18.609, 3.59, -0.094, 0.0118, 1, 0, 0, -1) - z = ex[1][1] * ex[-1][4] - DLAM += 18.609 * z.imag - DS += 3.59 * z.imag - GAM1C += -0.094 * z.real - SINPI += 0.0118 * z.real - - # AddSol(-4586.465, -4578.13, -0.077, 34.3117, 1, 0, 0, -2) - z = ex[1][1] * ex[-2][4] - DLAM += -4586.465 * z.imag - DS += -4578.13 * z.imag - GAM1C += -0.077 * z.real - SINPI += 34.3117 * z.real - - # AddSol(3.215, 5.44, 0.192, -0.0386, 1, 0, 0, -3) - z = ex[1][1] * ex[-3][4] - DLAM += 3.215 * z.imag - DS += 5.44 * z.imag - GAM1C += 0.192 * z.real - SINPI += -0.0386 * z.real - - # AddSol(-38.428, -38.64, 0.001, 0.6008, 1, 0, 0, -4) - z = ex[1][1] * ex[-4][4] - DLAM += -38.428 * z.imag - DS += -38.64 * z.imag - GAM1C += 0.001 * z.real - SINPI += 0.6008 * z.real - - # AddSol(-0.393, -1.43, -0.092, 0.0086, 1, 0, 0, -6) - z = ex[1][1] * ex[-6][4] - DLAM += -0.393 * z.imag - DS += -1.43 * z.imag - GAM1C += -0.092 * z.real - SINPI += 0.0086 * z.real - - # AddSol(-0.289, -1.59, 0.123, -0.0053, 0, 1, 0, 4) - z = ex[1][2] * ex[4][4] - DLAM += -0.289 * z.imag - DS += -1.59 * z.imag - GAM1C += 0.123 * z.real - SINPI += -0.0053 * z.real - - # AddSol(-24.42, -25.1, 0.04, -0.3, 0, 1, 0, 2) - z = ex[1][2] * ex[2][4] - DLAM += -24.42 * z.imag - DS += -25.1 * z.imag - GAM1C += 0.04 * z.real - SINPI += -0.3 * z.real - - # AddSol(18.023, 17.93, 0.007, 0.1494, 0, 1, 0, 1) - z = ex[1][2] * ex[1][4] - DLAM += 18.023 * z.imag - DS += 17.93 * z.imag - GAM1C += 0.007 * z.real - SINPI += 0.1494 * z.real - - # AddSol(-668.146, -126.98, -1.302, -0.3997, 0, 1, 0, 0) - z = ex[1][2] - DLAM += -668.146 * z.imag - DS += -126.98 * z.imag - GAM1C += -1.302 * z.real - SINPI += -0.3997 * z.real - - # AddSol(0.56, 0.32, -0.001, -0.0037, 0, 1, 0, -1) - z = ex[1][2] * ex[-1][4] - DLAM += 0.56 * z.imag - DS += 0.32 * z.imag - GAM1C += -0.001 * z.real - SINPI += -0.0037 * z.real - - # AddSol(-165.145, -165.06, 0.054, 1.9178, 0, 1, 0, -2) - z = ex[1][2] * ex[-2][4] - DLAM += -165.145 * z.imag - DS += -165.06 * z.imag - GAM1C += 0.054 * z.real - SINPI += 1.9178 * z.real - - # AddSol(-1.877, -6.46, -0.416, 0.0339, 0, 1, 0, -4) - z = ex[1][2] * ex[-4][4] - DLAM += -1.877 * z.imag - DS += -6.46 * z.imag - GAM1C += -0.416 * z.real - SINPI += 0.0339 * z.real - - # AddSol(0.213, 1.02, -0.074, 0.0054, 2, 0, 0, 4) - z = ex[2][1] * ex[4][4] - DLAM += 0.213 * z.imag - DS += 1.02 * z.imag - GAM1C += -0.074 * z.real - SINPI += 0.0054 * z.real - - # AddSol(14.387, 14.78, -0.017, 0.2833, 2, 0, 0, 2) - z = ex[2][1] * ex[2][4] - DLAM += 14.387 * z.imag - DS += 14.78 * z.imag - GAM1C += -0.017 * z.real - SINPI += 0.2833 * z.real - - # AddSol(-0.586, -1.2, 0.054, -0.01, 2, 0, 0, 1) - z = ex[2][1] * ex[1][4] - DLAM += -0.586 * z.imag - DS += -1.2 * z.imag - GAM1C += 0.054 * z.real - SINPI += -0.01 * z.real - - # AddSol(769.016, 767.96, 0.107, 10.1657, 2, 0, 0, 0) - z = ex[2][1] - DLAM += 769.016 * z.imag - DS += 767.96 * z.imag - GAM1C += 0.107 * z.real - SINPI += 10.1657 * z.real - - # AddSol(1.75, 2.01, -0.018, 0.0155, 2, 0, 0, -1) - z = ex[2][1] * ex[-1][4] - DLAM += 1.75 * z.imag - DS += 2.01 * z.imag - GAM1C += -0.018 * z.real - SINPI += 0.0155 * z.real - - # AddSol(-211.656, -152.53, 5.679, -0.3039, 2, 0, 0, -2) - z = ex[2][1] * ex[-2][4] - DLAM += -211.656 * z.imag - DS += -152.53 * z.imag - GAM1C += 5.679 * z.real - SINPI += -0.3039 * z.real - - # AddSol(1.225, 0.91, -0.03, -0.0088, 2, 0, 0, -3) - z = ex[2][1] * ex[-3][4] - DLAM += 1.225 * z.imag - DS += 0.91 * z.imag - GAM1C += -0.03 * z.real - SINPI += -0.0088 * z.real - - # AddSol(-30.773, -34.07, -0.308, 0.3722, 2, 0, 0, -4) - z = ex[2][1] * ex[-4][4] - DLAM += -30.773 * z.imag - DS += -34.07 * z.imag - GAM1C += -0.308 * z.real - SINPI += 0.3722 * z.real - - # AddSol(-0.57, -1.4, -0.074, 0.0109, 2, 0, 0, -6) - z = ex[2][1] * ex[-6][4] - DLAM += -0.57 * z.imag - DS += -1.4 * z.imag - GAM1C += -0.074 * z.real - SINPI += 0.0109 * z.real - - # AddSol(-2.921, -11.75, 0.787, -0.0484, 1, 1, 0, 2) - z = ex[1][1] * ex[1][2] * ex[2][4] - DLAM += -2.921 * z.imag - DS += -11.75 * z.imag - GAM1C += 0.787 * z.real - SINPI += -0.0484 * z.real - - # AddSol(1.267, 1.52, -0.022, 0.0164, 1, 1, 0, 1) - z = ex[1][1] * ex[1][2] * ex[1][4] - DLAM += 1.267 * z.imag - DS += 1.52 * z.imag - GAM1C += -0.022 * z.real - SINPI += 0.0164 * z.real - - # AddSol(-109.673, -115.18, 0.461, -0.949, 1, 1, 0, 0) - z = ex[1][1] * ex[1][2] - DLAM += -109.673 * z.imag - DS += -115.18 * z.imag - GAM1C += 0.461 * z.real - SINPI += -0.949 * z.real - - # AddSol(-205.962, -182.36, 2.056, 1.4437, 1, 1, 0, -2) - z = ex[1][1] * ex[1][2] * ex[-2][4] - DLAM += -205.962 * z.imag - DS += -182.36 * z.imag - GAM1C += 2.056 * z.real - SINPI += 1.4437 * z.real - - # AddSol(0.233, 0.36, 0.012, -0.0025, 1, 1, 0, -3) - z = ex[1][1] * ex[1][2] * ex[-3][4] - DLAM += 0.233 * z.imag - DS += 0.36 * z.imag - GAM1C += 0.012 * z.real - SINPI += -0.0025 * z.real - - # AddSol(-4.391, -9.66, -0.471, 0.0673, 1, 1, 0, -4) - z = ex[1][1] * ex[1][2] * ex[-4][4] - DLAM += -4.391 * z.imag - DS += -9.66 * z.imag - GAM1C += -0.471 * z.real - SINPI += 0.0673 * z.real - - - # AddSol(0.283, 1.53, -0.111, 0.006, 1, -1, 0, 4) - z = ex[1][1] * ex[-1][2] * ex[4][4] - DLAM += 0.283 * z.imag - DS += 1.53 * z.imag - GAM1C += -0.111 * z.real - SINPI += 0.006 * z.real - - # AddSol(14.577, 31.7, -1.54, 0.2302, 1, -1, 0, 2) - z = ex[1][1] * ex[-1][2] * ex[2][4] - DLAM += 14.577 * z.imag - DS += 31.7 * z.imag - GAM1C += -1.54 * z.real - SINPI += 0.2302 * z.real - - # AddSol(147.687, 138.76, 0.679, 1.1528, 1, -1, 0, 0) - z = ex[1][1] * ex[-1][2] - DLAM += 147.687 * z.imag - DS += 138.76 * z.imag - GAM1C += 0.679 * z.real - SINPI += 1.1528 * z.real - - # AddSol(-1.089, 0.55, 0.021, 0.0, 1, -1, 0, -1) - z = ex[1][1] * ex[-1][2] * ex[-1][4] - DLAM += -1.089 * z.imag - DS += 0.55 * z.imag - GAM1C += 0.021 * z.real - - # AddSol(28.475, 23.59, -0.443, -0.2257, 1, -1, 0, -2) - z = ex[1][1] * ex[-1][2] * ex[-2][4] - DLAM += 28.475 * z.imag - DS += 23.59 * z.imag - GAM1C += -0.443 * z.real - SINPI += -0.2257 * z.real - - # AddSol(-0.276, -0.38, -0.006, -0.0036, 1, -1, 0, -3) - z = ex[1][1] * ex[-1][2] * ex[-3][4] - DLAM += -0.276 * z.imag - DS += -0.38 * z.imag - GAM1C += -0.006 * z.real - SINPI += -0.0036 * z.real - - # AddSol(0.636, 2.27, 0.146, -0.0102, 1, -1, 0, -4) - z = ex[1][1] * ex[-1][2] * ex[-4][4] - DLAM += 0.636 * z.imag - DS += 2.27 * z.imag - GAM1C += 0.146 * z.real - SINPI += -0.0102 * z.real - - # AddSol(-0.189, -1.68, 0.131, -0.0028, 0, 2, 0, 2) - z = ex[2][2] * ex[2][4] - DLAM += -0.189 * z.imag - DS += -1.68 * z.imag - GAM1C += 0.131 * z.real - SINPI += -0.0028 * z.real - - # AddSol(-7.486, -0.66, -0.037, -0.0086, 0, 2, 0, 0) - z = ex[2][2] - DLAM += -7.486 * z.imag - DS += -0.66 * z.imag - GAM1C += -0.037 * z.real - SINPI += -0.0086 * z.real - - # AddSol(-8.096, -16.35, -0.74, 0.0918, 0, 2, 0, -2) - z = ex[2][2] * ex[-2][4] - DLAM += -8.096 * z.imag - DS += -16.35 * z.imag - GAM1C += -0.74 * z.real - SINPI += 0.0918 * z.real - - # AddSol(-5.741, -0.04, 0.0, -0.0009, 0, 0, 2, 2) - z = ex[2][3] * ex[2][4] - DLAM += -5.741 * z.imag - DS += -0.04 * z.imag - SINPI += -0.0009 * z.real - - # AddSol(0.255, 0.0, 0.0, 0.0, 0, 0, 2, 1) - z = ex[2][3] * ex[1][4] - DLAM += 0.255 * z.imag - - # AddSol(-411.608, -0.2, 0.0, -0.0124, 0, 0, 2, 0) - z = ex[2][3] - DLAM += -411.608 * z.imag - DS += -0.2 * z.imag - SINPI += -0.0124 * z.real - - # AddSol(0.584, 0.84, 0.0, 0.0071, 0, 0, 2, -1) - z = ex[2][3] * ex[-1][4] - DLAM += 0.584 * z.imag - DS += 0.84 * z.imag - SINPI += 0.0071 * z.real - - # AddSol(-55.173, -52.14, 0.0, -0.1052, 0, 0, 2, -2) - z = ex[2][3] * ex[-2][4] - DLAM += -55.173 * z.imag - DS += -52.14 * z.imag - SINPI += -0.1052 * z.real - - # AddSol(0.254, 0.25, 0.0, -0.0017, 0, 0, 2, -3) - z = ex[2][3] * ex[-3][4] - DLAM += 0.254 * z.imag - DS += 0.25 * z.imag - SINPI += -0.0017 * z.real - - # AddSol(0.025, -1.67, 0.0, 0.0031, 0, 0, 2, -4) - z = ex[2][3] * ex[-4][4] - DLAM += 0.025 * z.imag - DS += -1.67 * z.imag - SINPI += 0.0031 * z.real - - # AddSol(1.06, 2.96, -0.166, 0.0243, 3, 0, 0, 2) - z = ex[3][1] * ex[2][4] - DLAM += 1.06 * z.imag - DS += 2.96 * z.imag - GAM1C += -0.166 * z.real - SINPI += 0.0243 * z.real - - # AddSol(36.124, 50.64, -1.3, 0.6215, 3, 0, 0, 0) - z = ex[3][1] - DLAM += 36.124 * z.imag - DS += 50.64 * z.imag - GAM1C += -1.3 * z.real - SINPI += 0.6215 * z.real - - # AddSol(-13.193, -16.4, 0.258, -0.1187, 3, 0, 0, -2) - z = ex[3][1] * ex[-2][4] - DLAM += -13.193 * z.imag - DS += -16.4 * z.imag - GAM1C += 0.258 * z.real - SINPI += -0.1187 * z.real - - # AddSol(-1.187, -0.74, 0.042, 0.0074, 3, 0, 0, -4) - z = ex[3][1] * ex[-4][4] - DLAM += -1.187 * z.imag - DS += -0.74 * z.imag - GAM1C += 0.042 * z.real - SINPI += 0.0074 * z.real - - # AddSol(-0.293, -0.31, -0.002, 0.0046, 3, 0, 0, -6) - z = ex[3][1] * ex[-6][4] - DLAM += -0.293 * z.imag - DS += -0.31 * z.imag - GAM1C += -0.002 * z.real - SINPI += 0.0046 * z.real - - # AddSol(-0.29, -1.45, 0.116, -0.0051, 2, 1, 0, 2) - z = ex[2][1] * ex[1][2] * ex[2][4] - DLAM += -0.29 * z.imag - DS += -1.45 * z.imag - GAM1C += 0.116 * z.real - SINPI += -0.0051 * z.real - - # AddSol(-7.649, -10.56, 0.259, -0.1038, 2, 1, 0, 0) - z = ex[2][1] * ex[1][2] - DLAM += -7.649 * z.imag - DS += -10.56 * z.imag - GAM1C += 0.259 * z.real - SINPI += -0.1038 * z.real - - # AddSol(-8.627, -7.59, 0.078, -0.0192, 2, 1, 0, -2) - z = ex[2][1] * ex[1][2] * ex[-2][4] - DLAM += -8.627 * z.imag - DS += -7.59 * z.imag - GAM1C += 0.078 * z.real - SINPI += -0.0192 * z.real - - # AddSol(-2.74, -2.54, 0.022, 0.0324, 2, 1, 0, -4) - z = ex[2][1] * ex[1][2] * ex[-4][4] - DLAM += -2.74 * z.imag - DS += -2.54 * z.imag - GAM1C += 0.022 * z.real - SINPI += 0.0324 * z.real - - # AddSol(1.181, 3.32, -0.212, 0.0213, 2, -1, 0, 2) - z = ex[2][1] * ex[-1][2] * ex[2][4] - DLAM += 1.181 * z.imag - DS += 3.32 * z.imag - GAM1C += -0.212 * z.real - SINPI += 0.0213 * z.real - - # AddSol(9.703, 11.67, -0.151, 0.1268, 2, -1, 0, 0) - z = ex[2][1] * ex[-1][2] - DLAM += 9.703 * z.imag - DS += 11.67 * z.imag - GAM1C += -0.151 * z.real - SINPI += 0.1268 * z.real - - # AddSol(-0.352, -0.37, 0.001, -0.0028, 2, -1, 0, -1) - z = ex[2][1] * ex[-1][2] * ex[-1][4] - DLAM += -0.352 * z.imag - DS += -0.37 * z.imag - GAM1C += 0.001 * z.real - SINPI += -0.0028 * z.real - - # AddSol(-2.494, -1.17, -0.003, -0.0017, 2, -1, 0, -2) - z = ex[2][1] * ex[-1][2] * ex[-2][4] - DLAM += -2.494 * z.imag - DS += -1.17 * z.imag - GAM1C += -0.003 * z.real - SINPI += -0.0017 * z.real - - # AddSol(0.36, 0.2, -0.012, -0.0043, 2, -1, 0, -4) - z = ex[2][1] * ex[-1][2] * ex[-4][4] - DLAM += 0.36 * z.imag - DS += 0.2 * z.imag - GAM1C += -0.012 * z.real - SINPI += -0.0043 * z.real - - # AddSol(-1.167, -1.25, 0.008, -0.0106, 1, 2, 0, 0) - z = ex[1][1] * ex[2][2] - DLAM += -1.167 * z.imag - DS += -1.25 * z.imag - GAM1C += 0.008 * z.real - SINPI += -0.0106 * z.real - - # AddSol(-7.412, -6.12, 0.117, 0.0484, 1, 2, 0, -2) - z = ex[1][1] * ex[2][2] * ex[-2][4] - DLAM += -7.412 * z.imag - DS += -6.12 * z.imag - GAM1C += 0.117 * z.real - SINPI += 0.0484 * z.real - - # AddSol(-0.311, -0.65, -0.032, 0.0044, 1, 2, 0, -4) - z = ex[1][1] * ex[2][2] * ex[-4][4] - DLAM += -0.311 * z.imag - DS += -0.65 * z.imag - GAM1C += -0.032 * z.real - SINPI += 0.0044 * z.real - - # AddSol(0.757, 1.82, -0.105, 0.0112, 1, -2, 0, 2) - z = ex[1][1] * ex[-2][2] * ex[2][4] - DLAM += 0.757 * z.imag - DS += 1.82 * z.imag - GAM1C += -0.105 * z.real - SINPI += 0.0112 * z.real - - # AddSol(2.58, 2.32, 0.027, 0.0196, 1, -2, 0, 0) - z = ex[1][1] * ex[-2][2] - DLAM += 2.58 * z.imag - DS += 2.32 * z.imag - GAM1C += 0.027 * z.real - SINPI += 0.0196 * z.real - - # AddSol(2.533, 2.4, -0.014, -0.0212, 1, -2, 0, -2) - z = ex[1][1] * ex[-2][2] * ex[-2][4] - DLAM += 2.533 * z.imag - DS += 2.4 * z.imag - GAM1C += -0.014 * z.real - SINPI += -0.0212 * z.real - - # AddSol(-0.344, -0.57, -0.025, 0.0036, 0, 3, 0, -2) - z = ex[3][2] * ex[-2][4] - DLAM += -0.344 * z.imag - DS += -0.57 * z.imag - GAM1C += -0.025 * z.real - SINPI += 0.0036 * z.real - - # AddSol(-0.992, -0.02, 0.0, 0.0, 1, 0, 2, 2) - z = ex[1][1] * ex[2][3] * ex[2][4] - DLAM += -0.992 * z.imag - DS += -0.02 * z.imag - - # AddSol(-45.099, -0.02, 0.0, -0.001, 1, 0, 2, 0) - z = ex[1][1] * ex[2][3] - DLAM += -45.099 * z.imag - DS += -0.02 * z.imag - SINPI += -0.001 * z.real - - # AddSol(-0.179, -9.52, 0.0, -0.0833, 1, 0, 2, -2) - z = ex[1][1] * ex[2][3] * ex[-2][4] - DLAM += -0.179 * z.imag - DS += -9.52 * z.imag - SINPI += -0.0833 * z.real - - # AddSol(-0.301, -0.33, 0.0, 0.0014, 1, 0, 2, -4) - z = ex[1][1] * ex[2][3] * ex[-4][4] - DLAM += -0.301 * z.imag - DS += -0.33 * z.imag - SINPI += 0.0014 * z.real - - # AddSol(-6.382, -3.37, 0.0, -0.0481, 1, 0, -2, 2) - z = ex[1][1] * ex[-2][3] * ex[2][4] - DLAM += -6.382 * z.imag - DS += -3.37 * z.imag - SINPI += -0.0481 * z.real - - # AddSol(39.528, 85.13, 0.0, -0.7136, 1, 0, -2, 0) - z = ex[1][1] * ex[-2][3] - DLAM += 39.528 * z.imag - DS += 85.13 * z.imag - SINPI += -0.7136 * z.real - - # AddSol(9.366, 0.71, 0.0, -0.0112, 1, 0, -2, -2) - z = ex[1][1] * ex[-2][3] * ex[-2][4] - DLAM += 9.366 * z.imag - DS += 0.71 * z.imag - SINPI += -0.0112 * z.real - - # AddSol(0.202, 0.02, 0.0, 0.0, 1, 0, -2, -4) - z = ex[1][1] * ex[-2][3] * ex[-4][4] - DLAM += 0.202 * z.imag - DS += 0.02 * z.imag - - - # AddSol(0.415, 0.1, 0.0, 0.0013, 0, 1, 2, 0) - z = ex[1][2] * ex[2][3] - DLAM += 0.415 * z.imag - DS += 0.1 * z.imag - SINPI += 0.0013 * z.real - - # AddSol(-2.152, -2.26, 0.0, -0.0066, 0, 1, 2, -2) - z = ex[1][2] * ex[2][3] * ex[-2][4] - DLAM += -2.152 * z.imag - DS += -2.26 * z.imag - SINPI += -0.0066 * z.real - - # AddSol(-1.44, -1.3, 0.0, 0.0014, 0, 1, -2, 2) - z = ex[1][2] * ex[-2][3] * ex[2][4] - DLAM += -1.44 * z.imag - DS += -1.3 * z.imag - SINPI += 0.0014 * z.real - - # AddSol(0.384, -0.04, 0.0, 0.0, 0, 1, -2, -2) - z = ex[1][2] * ex[-2][3] * ex[-2][4] - DLAM += 0.384 * z.imag - DS += -0.04 * z.imag - - # AddSol(1.938, 3.6, -0.145, 0.0401, 4, 0, 0, 0) - z = ex[4][1] - DLAM += 1.938 * z.imag - DS += 3.6 * z.imag - GAM1C += -0.145 * z.real - SINPI += 0.0401 * z.real - - # AddSol(-0.952, -1.58, 0.052, -0.013, 4, 0, 0, -2) - z = ex[4][1] * ex[-2][4] - DLAM += -0.952 * z.imag - DS += -1.58 * z.imag - GAM1C += 0.052 * z.real - SINPI += -0.013 * z.real - - # AddSol(-0.551, -0.94, 0.032, -0.0097, 3, 1, 0, 0) - z = ex[3][1] * ex[1][2] - DLAM += -0.551 * z.imag - DS += -0.94 * z.imag - GAM1C += 0.032 * z.real - SINPI += -0.0097 * z.real - - # AddSol(-0.482, -0.57, 0.005, -0.0045, 3, 1, 0, -2) - z = ex[3][1] * ex[1][2] * ex[-2][4] - DLAM += -0.482 * z.imag - DS += -0.57 * z.imag - GAM1C += 0.005 * z.real - SINPI += -0.0045 * z.real - - # AddSol(0.681, 0.96, -0.026, 0.0115, 3, -1, 0, 0) - z = ex[3][1] * ex[-1][2] - DLAM += 0.681 * z.imag - DS += 0.96 * z.imag - GAM1C += -0.026 * z.real - SINPI += 0.0115 * z.real - - # AddSol(-0.297, -0.27, 0.002, -0.0009, 2, 2, 0, -2) - z = ex[2][1] * ex[2][2] * ex[-2][4] - DLAM += -0.297 * z.imag - DS += -0.27 * z.imag - GAM1C += 0.002 * z.real - SINPI += -0.0009 * z.real - - # AddSol(0.254, 0.21, -0.003, 0.0, 2, -2, 0, -2) - z = ex[2][1] * ex[-2][2] * ex[-2][4] - DLAM += 0.254 * z.imag - DS += 0.21 * z.imag - GAM1C += -0.003 * z.real - - # AddSol(-0.25, -0.22, 0.004, 0.0014, 1, 3, 0, -2) - z = ex[1][1] * ex[3][2] * ex[-2][4] - DLAM += -0.25 * z.imag - DS += -0.22 * z.imag - GAM1C += 0.004 * z.real - SINPI += 0.0014 * z.real - - # AddSol(-3.996, 0.0, 0.0, 0.0004, 2, 0, 2, 0) - z = ex[2][1] * ex[2][3] - DLAM += -3.996 * z.imag - SINPI += 0.0004 * z.real - - # AddSol(0.557, -0.75, 0.0, -0.009, 2, 0, 2, -2) - z = ex[2][1] * ex[2][3] * ex[-2][4] - DLAM += 0.557 * z.imag - DS += -0.75 * z.imag - SINPI += -0.009 * z.real - - # AddSol(-0.459, -0.38, 0.0, -0.0053, 2, 0, -2, 2) - z = ex[2][1] * ex[-2][3] * ex[2][4] - DLAM += -0.459 * z.imag - DS += -0.38 * z.imag - SINPI += -0.0053 * z.real - - # AddSol(-1.298, 0.74, 0.0, 0.0004, 2, 0, -2, 0) - z = ex[2][1] * ex[-2][3] - DLAM += -1.298 * z.imag - DS += 0.74 * z.imag - SINPI += 0.0004 * z.real - - # AddSol(0.538, 1.14, 0.0, -0.0141, 2, 0, -2, -2) - z = ex[2][1] * ex[-2][3] * ex[-2][4] - DLAM += 0.538 * z.imag - DS += 1.14 * z.imag - SINPI += -0.0141 * z.real - - # AddSol(0.263, 0.02, 0.0, 0.0, 1, 1, 2, 0) - z = ex[1][1] * ex[1][2] * ex[2][3] - DLAM += 0.263 * z.imag - DS += 0.02 * z.imag - - # AddSol(0.426, 0.07, 0.0, -0.0006, 1, 1, -2, -2) - z = ex[1][1] * ex[1][2] * ex[-2][3] * ex[-2][4] - DLAM += 0.426 * z.imag - DS += 0.07 * z.imag - SINPI += -0.0006 * z.real - - # AddSol(-0.304, 0.03, 0.0, 0.0003, 1, -1, 2, 0) - z = ex[1][1] * ex[-1][2] * ex[2][3] - DLAM += -0.304 * z.imag - DS += 0.03 * z.imag - SINPI += 0.0003 * z.real - - # AddSol(-0.372, -0.19, 0.0, -0.0027, 1, -1, -2, 2) - z = ex[1][1] * ex[-1][2] * ex[-2][3] * ex[2][4] - DLAM += -0.372 * z.imag - DS += -0.19 * z.imag - SINPI += -0.0027 * z.real - - # AddSol(0.418, 0.0, 0.0, 0.0, 0, 0, 4, 0) - z = ex[4][3] - DLAM += 0.418 * z.imag - - # AddSol(-0.33, -0.04, 0.0, 0.0, 3, 0, 2, 0) - z = ex[3][1] * ex[2][3] - DLAM += -0.33 * z.imag - DS += -0.04 * z.imag + $ASTRO_ADDSOL( 13.902, 14.06,-0.001, 0.2607,0, 0, 0, 4) + $ASTRO_ADDSOL( 0.403, -4.01,+0.394, 0.0023,0, 0, 0, 3) + $ASTRO_ADDSOL( 2369.912, 2373.36,+0.601, 28.2333,0, 0, 0, 2) + $ASTRO_ADDSOL( -125.154, -112.79,-0.725, -0.9781,0, 0, 0, 1) + $ASTRO_ADDSOL( 1.979, 6.98,-0.445, 0.0433,1, 0, 0, 4) + $ASTRO_ADDSOL( 191.953, 192.72,+0.029, 3.0861,1, 0, 0, 2) + $ASTRO_ADDSOL( -8.466, -13.51,+0.455, -0.1093,1, 0, 0, 1) + $ASTRO_ADDSOL( 22639.500,22609.07,+0.079, 186.5398,1, 0, 0, 0) + $ASTRO_ADDSOL( 18.609, 3.59,-0.094, 0.0118,1, 0, 0,-1) + $ASTRO_ADDSOL( -4586.465,-4578.13,-0.077, 34.3117,1, 0, 0,-2) + $ASTRO_ADDSOL( +3.215, 5.44,+0.192, -0.0386,1, 0, 0,-3) + $ASTRO_ADDSOL( -38.428, -38.64,+0.001, 0.6008,1, 0, 0,-4) + $ASTRO_ADDSOL( -0.393, -1.43,-0.092, 0.0086,1, 0, 0,-6) + $ASTRO_ADDSOL( -0.289, -1.59,+0.123, -0.0053,0, 1, 0, 4) + $ASTRO_ADDSOL( -24.420, -25.10,+0.040, -0.3000,0, 1, 0, 2) + $ASTRO_ADDSOL( 18.023, 17.93,+0.007, 0.1494,0, 1, 0, 1) + $ASTRO_ADDSOL( -668.146, -126.98,-1.302, -0.3997,0, 1, 0, 0) + $ASTRO_ADDSOL( 0.560, 0.32,-0.001, -0.0037,0, 1, 0,-1) + $ASTRO_ADDSOL( -165.145, -165.06,+0.054, 1.9178,0, 1, 0,-2) + $ASTRO_ADDSOL( -1.877, -6.46,-0.416, 0.0339,0, 1, 0,-4) + $ASTRO_ADDSOL( 0.213, 1.02,-0.074, 0.0054,2, 0, 0, 4) + $ASTRO_ADDSOL( 14.387, 14.78,-0.017, 0.2833,2, 0, 0, 2) + $ASTRO_ADDSOL( -0.586, -1.20,+0.054, -0.0100,2, 0, 0, 1) + $ASTRO_ADDSOL( 769.016, 767.96,+0.107, 10.1657,2, 0, 0, 0) + $ASTRO_ADDSOL( +1.750, 2.01,-0.018, 0.0155,2, 0, 0,-1) + $ASTRO_ADDSOL( -211.656, -152.53,+5.679, -0.3039,2, 0, 0,-2) + $ASTRO_ADDSOL( +1.225, 0.91,-0.030, -0.0088,2, 0, 0,-3) + $ASTRO_ADDSOL( -30.773, -34.07,-0.308, 0.3722,2, 0, 0,-4) + $ASTRO_ADDSOL( -0.570, -1.40,-0.074, 0.0109,2, 0, 0,-6) + $ASTRO_ADDSOL( -2.921, -11.75,+0.787, -0.0484,1, 1, 0, 2) + $ASTRO_ADDSOL( +1.267, 1.52,-0.022, 0.0164,1, 1, 0, 1) + $ASTRO_ADDSOL( -109.673, -115.18,+0.461, -0.9490,1, 1, 0, 0) + $ASTRO_ADDSOL( -205.962, -182.36,+2.056, +1.4437,1, 1, 0,-2) + $ASTRO_ADDSOL( 0.233, 0.36, 0.012, -0.0025,1, 1, 0,-3) + $ASTRO_ADDSOL( -4.391, -9.66,-0.471, 0.0673,1, 1, 0,-4) + + $ASTRO_ADDSOL( 0.283, 1.53,-0.111, +0.0060,1,-1, 0,+4) + $ASTRO_ADDSOL( 14.577, 31.70,-1.540, +0.2302,1,-1, 0, 2) + $ASTRO_ADDSOL( 147.687, 138.76,+0.679, +1.1528,1,-1, 0, 0) + $ASTRO_ADDSOL( -1.089, 0.55,+0.021, 0.0 ,1,-1, 0,-1) + $ASTRO_ADDSOL( 28.475, 23.59,-0.443, -0.2257,1,-1, 0,-2) + $ASTRO_ADDSOL( -0.276, -0.38,-0.006, -0.0036,1,-1, 0,-3) + $ASTRO_ADDSOL( 0.636, 2.27,+0.146, -0.0102,1,-1, 0,-4) + $ASTRO_ADDSOL( -0.189, -1.68,+0.131, -0.0028,0, 2, 0, 2) + $ASTRO_ADDSOL( -7.486, -0.66,-0.037, -0.0086,0, 2, 0, 0) + $ASTRO_ADDSOL( -8.096, -16.35,-0.740, 0.0918,0, 2, 0,-2) + $ASTRO_ADDSOL( -5.741, -0.04, 0.0 , -0.0009,0, 0, 2, 2) + $ASTRO_ADDSOL( 0.255, 0.0 , 0.0 , 0.0 ,0, 0, 2, 1) + $ASTRO_ADDSOL( -411.608, -0.20, 0.0 , -0.0124,0, 0, 2, 0) + $ASTRO_ADDSOL( 0.584, 0.84, 0.0 , +0.0071,0, 0, 2,-1) + $ASTRO_ADDSOL( -55.173, -52.14, 0.0 , -0.1052,0, 0, 2,-2) + $ASTRO_ADDSOL( 0.254, 0.25, 0.0 , -0.0017,0, 0, 2,-3) + $ASTRO_ADDSOL( +0.025, -1.67, 0.0 , +0.0031,0, 0, 2,-4) + $ASTRO_ADDSOL( 1.060, 2.96,-0.166, 0.0243,3, 0, 0,+2) + $ASTRO_ADDSOL( 36.124, 50.64,-1.300, 0.6215,3, 0, 0, 0) + $ASTRO_ADDSOL( -13.193, -16.40,+0.258, -0.1187,3, 0, 0,-2) + $ASTRO_ADDSOL( -1.187, -0.74,+0.042, 0.0074,3, 0, 0,-4) + $ASTRO_ADDSOL( -0.293, -0.31,-0.002, 0.0046,3, 0, 0,-6) + $ASTRO_ADDSOL( -0.290, -1.45,+0.116, -0.0051,2, 1, 0, 2) + $ASTRO_ADDSOL( -7.649, -10.56,+0.259, -0.1038,2, 1, 0, 0) + $ASTRO_ADDSOL( -8.627, -7.59,+0.078, -0.0192,2, 1, 0,-2) + $ASTRO_ADDSOL( -2.740, -2.54,+0.022, 0.0324,2, 1, 0,-4) + $ASTRO_ADDSOL( 1.181, 3.32,-0.212, 0.0213,2,-1, 0,+2) + $ASTRO_ADDSOL( 9.703, 11.67,-0.151, 0.1268,2,-1, 0, 0) + $ASTRO_ADDSOL( -0.352, -0.37,+0.001, -0.0028,2,-1, 0,-1) + $ASTRO_ADDSOL( -2.494, -1.17,-0.003, -0.0017,2,-1, 0,-2) + $ASTRO_ADDSOL( 0.360, 0.20,-0.012, -0.0043,2,-1, 0,-4) + $ASTRO_ADDSOL( -1.167, -1.25,+0.008, -0.0106,1, 2, 0, 0) + $ASTRO_ADDSOL( -7.412, -6.12,+0.117, 0.0484,1, 2, 0,-2) + $ASTRO_ADDSOL( -0.311, -0.65,-0.032, 0.0044,1, 2, 0,-4) + $ASTRO_ADDSOL( +0.757, 1.82,-0.105, 0.0112,1,-2, 0, 2) + $ASTRO_ADDSOL( +2.580, 2.32,+0.027, 0.0196,1,-2, 0, 0) + $ASTRO_ADDSOL( +2.533, 2.40,-0.014, -0.0212,1,-2, 0,-2) + $ASTRO_ADDSOL( -0.344, -0.57,-0.025, +0.0036,0, 3, 0,-2) + $ASTRO_ADDSOL( -0.992, -0.02, 0.0 , 0.0 ,1, 0, 2, 2) + $ASTRO_ADDSOL( -45.099, -0.02, 0.0 , -0.0010,1, 0, 2, 0) + $ASTRO_ADDSOL( -0.179, -9.52, 0.0 , -0.0833,1, 0, 2,-2) + $ASTRO_ADDSOL( -0.301, -0.33, 0.0 , 0.0014,1, 0, 2,-4) + $ASTRO_ADDSOL( -6.382, -3.37, 0.0 , -0.0481,1, 0,-2, 2) + $ASTRO_ADDSOL( 39.528, 85.13, 0.0 , -0.7136,1, 0,-2, 0) + $ASTRO_ADDSOL( 9.366, 0.71, 0.0 , -0.0112,1, 0,-2,-2) + $ASTRO_ADDSOL( 0.202, 0.02, 0.0 , 0.0 ,1, 0,-2,-4) + + $ASTRO_ADDSOL( 0.415, 0.10, 0.0 , 0.0013,0, 1, 2, 0) + $ASTRO_ADDSOL( -2.152, -2.26, 0.0 , -0.0066,0, 1, 2,-2) + $ASTRO_ADDSOL( -1.440, -1.30, 0.0 , +0.0014,0, 1,-2, 2) + $ASTRO_ADDSOL( 0.384, -0.04, 0.0 , 0.0 ,0, 1,-2,-2) + $ASTRO_ADDSOL( +1.938, +3.60,-0.145, +0.0401,4, 0, 0, 0) + $ASTRO_ADDSOL( -0.952, -1.58,+0.052, -0.0130,4, 0, 0,-2) + $ASTRO_ADDSOL( -0.551, -0.94,+0.032, -0.0097,3, 1, 0, 0) + $ASTRO_ADDSOL( -0.482, -0.57,+0.005, -0.0045,3, 1, 0,-2) + $ASTRO_ADDSOL( 0.681, 0.96,-0.026, 0.0115,3,-1, 0, 0) + $ASTRO_ADDSOL( -0.297, -0.27, 0.002, -0.0009,2, 2, 0,-2) + $ASTRO_ADDSOL( 0.254, +0.21,-0.003, 0.0 ,2,-2, 0,-2) + $ASTRO_ADDSOL( -0.250, -0.22, 0.004, 0.0014,1, 3, 0,-2) + $ASTRO_ADDSOL( -3.996, 0.0 , 0.0 , +0.0004,2, 0, 2, 0) + $ASTRO_ADDSOL( 0.557, -0.75, 0.0 , -0.0090,2, 0, 2,-2) + $ASTRO_ADDSOL( -0.459, -0.38, 0.0 , -0.0053,2, 0,-2, 2) + $ASTRO_ADDSOL( -1.298, 0.74, 0.0 , +0.0004,2, 0,-2, 0) + $ASTRO_ADDSOL( 0.538, 1.14, 0.0 , -0.0141,2, 0,-2,-2) + $ASTRO_ADDSOL( 0.263, 0.02, 0.0 , 0.0 ,1, 1, 2, 0) + $ASTRO_ADDSOL( 0.426, +0.07, 0.0 , -0.0006,1, 1,-2,-2) + $ASTRO_ADDSOL( -0.304, +0.03, 0.0 , +0.0003,1,-1, 2, 0) + $ASTRO_ADDSOL( -0.372, -0.19, 0.0 , -0.0027,1,-1,-2, 2) + $ASTRO_ADDSOL( +0.418, 0.0 , 0.0 , 0.0 ,0, 0, 4, 0) + $ASTRO_ADDSOL( -0.330, -0.04, 0.0 , 0.0 ,3, 0, 2, 0) def ADDN(coeffn, p, q, r, s): return coeffn * (ex[p][1] * ex[q][2] * ex[r][3] * ex[s][4]).imag diff --git a/generate/template/old.py b/generate/template/old.py deleted file mode 100644 index 595e9d9a..00000000 --- a/generate/template/old.py +++ /dev/null @@ -1,1766 +0,0 @@ -#!/usr/bin/env python3 - -import math -import datetime - -_PI2 = 2.0 * math.pi -_EPOCH = datetime.datetime(2000, 1, 1, 12) -_T0 = 2451545.0 -_MJD_BASIS = 2400000.5 -_Y2000_IN_MJD = _T0 - _MJD_BASIS -_DEG2RAD = 0.017453292519943296 -_RAD2DEG = 57.295779513082321 -_ASEC360 = 1296000.0 -_ASEC2RAD = 4.848136811095359935899141e-6 -_ARC = 3600.0 * 180.0 / math.pi # arcseconds per radian -_C_AUDAY = 173.1446326846693 # speed of light in AU/day -_ERAD = 6378136.6 # mean earth radius in meters -_AU = 1.4959787069098932e+11 # astronomical unit in meters -_KM_PER_AU = 1.4959787069098932e+8 -_ANGVEL = 7.2921150e-5 -_SECONDS_PER_DAY = 24.0 * 3600.0 -_SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592 -_MEAN_SYNODIC_MONTH = 29.530588 -_EARTH_ORBITAL_PERIOD = 365.256 -_REFRACTION_NEAR_HORIZON = 34.0 / 60.0 -_SUN_RADIUS_AU = 4.6505e-3 -_MOON_RADIUS_AU = 1.15717e-5 -_ASEC180 = 180.0 * 60.0 * 60.0 -_AU_PER_PARSEC = _ASEC180 / math.pi - -def _LongitudeOffset(diff): - offset = diff - while offset <= -180.0: - offset += 360.0 - while offset > 180.0: - offset -= 360.0 - return offset - -def _NormalizeLongitude(lon): - while lon < 0.0: - lon += 360.0 - while lon >= 360.0: - lon -= 360.0 - return lon - -class Vector: - def __init__(self, x, y, z, t): - self.x = x - self.y = y - self.z = z - self.t = t - - def Length(self): - return math.sqrt(self.x**2 + self.y**2 + self.z**2) - -BODY_INVALID = -1 -BODY_MERCURY = 0 -BODY_VENUS = 1 -BODY_EARTH = 2 -BODY_MARS = 3 -BODY_JUPITER = 4 -BODY_SATURN = 5 -BODY_URANUS = 6 -BODY_NEPTUNE = 7 -BODY_PLUTO = 8 -BODY_SUN = 9 -BODY_MOON = 10 - -BodyName = [ - 'Mercury', - 'Venus', - 'Earth', - 'Mars', - 'Jupiter', - 'Saturn', - 'Uranus', - 'Neptune', - 'Pluto', - 'Sun', - 'Moon', -] - -def BodyCode(name): - return BodyName.index(name) - -def _IsSuperiorPlanet(body): - return body in [BODY_MARS, BODY_JUPITER, BODY_SATURN, BODY_URANUS, BODY_NEPTUNE, BODY_PLUTO] - -_PlanetOrbitalPeriod = [ - 87.969, - 224.701, - _EARTH_ORBITAL_PERIOD, - 686.980, - 4332.589, - 10759.22, - 30685.4, - 60189.0, - 90560.0 -] - -class Error(Exception): - def __init__(self, message): - Exception.__init__(self, message) - -class EarthNotAllowedError(Error): - def __init__(self): - Error.__init__(self, 'The Earth is not allowed as the body.') - -class InvalidBodyError(Error): - def __init__(self): - Error.__init__(self, 'Invalid astronomical body.') - -class BadVectorError(Error): - def __init__(self): - Error.__init__(self, 'Vector is too small to have a direction.') - -class InternalError(Error): - def __init__(self): - Error.__init__(self, 'Internal error - please report issue at https://github.com/cosinekitty/astronomy/issues') - -class NoConvergeError(Error): - def __init__(self): - Error.__init__(self, 'Numeric solver did not converge - please report issue at https://github.com/cosinekitty/astronomy/issues') - -def _SynodicPeriod(body): - if body == BODY_EARTH: - raise EarthNotAllowedError() - if body < 0 or body >= len(_PlanetOrbitalPeriod): - raise InvalidBodyError() - if body == BODY_MOON: - return _MEAN_SYNODIC_MONTH - return abs(_EARTH_ORBITAL_PERIOD / (_EARTH_ORBITAL_PERIOD/_PlanetOrbitalPeriod[body] - 1.0)) - -def _AngleBetween(a, b): - r = a.Length() * b.Length() - if r < 1.0e-8: - return BadVectorError() - dot = (a.x*b.x + a.y*b.y + a.z*b.z) / r - if dot <= -1.0: - return 180.0 - if dot >= +1.0: - return 0.0 - return _RAD2DEG * math.acos(dot) - -class _delta_t_entry_t: - def __init__(self, mjd, dt): - self.mjd = mjd - self.dt = dt - -_DT = $ASTRO_DELTA_T() - -def _DeltaT(mjd): - if mjd <= _DT[0].mjd: - return _DT[0].dt - if mjd >= _DT[-1].mjd: - return _DT[-1].dt - # Do a binary search to find the pair of indexes this mjd lies between. - lo = 0 - hi = len(_DT) - 2 # Make sure there is always an array element after the one we are looking at. - while True: - if lo > hi: - # This should never happen unless there is a bug in the binary search. - raise Error('Could not find delta-t value.') - c = (lo + hi) // 2 - if mjd < _DT[c].mjd: - hi = c-1 - elif mjd > _DT[c+1].mjd: - lo = c+1 - else: - frac = (mjd - _DT[c].mjd) / (_DT[c+1].mjd - _DT[c].mjd) - return _DT[c].dt + frac*(_DT[c+1].dt - _DT[c].dt) - -def _TerrestrialTime(ut): - return ut + _DeltaT(ut + _Y2000_IN_MJD) / 86400.0 - -class Time: - def __init__(self, ut): - self.ut = ut - self.tt = _TerrestrialTime(ut) - self.etilt = None - - @staticmethod - def Make(year, month, day, hour, minute, second): - micro = round(math.fmod(second, 1.0) * 1000000) - second = math.floor(second - micro/1000000) - d = datetime.datetime(year, month, day, hour, minute, second, micro) - ut = (d - _EPOCH).total_seconds() / 86400 - return Time(ut) - - @staticmethod - def Now(): - ut = (datetime.datetime.utcnow() - _EPOCH).total_seconds() / 86400.0 - return Time(ut) - - def AddDays(self, days): - return Time(self.ut + days) - - def __str__(self): - millis = round(self.ut * 86400000.0) - n = _EPOCH + datetime.timedelta(milliseconds=millis) - return '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}.{:03d}Z'.format(n.year, n.month, n.day, n.hour, n.minute, n.second, math.floor(n.microsecond / 1000)) - - def Utc(self): - return _EPOCH + datetime.timedelta(days=self.ut) - - def _etilt(self): - # Calculates precession and nutation of the Earth's axis. - # The calculations are very expensive, so lazy-evaluate and cache - # the result inside this Time object. - if self.etilt is None: - self.etilt = _e_tilt(self) - return self.etilt - - -class Observer: - def __init__(self, latitude, longitude, height=0): - self.latitude = latitude - self.longitude = longitude - self.height = height - -class _iau2000b: - def __init__(self, time): - t = time.tt / 36525.0 - el = math.fmod((485868.249036 + t*1717915923.2178), _ASEC360) * _ASEC2RAD - elp = math.fmod((1287104.79305 + t*129596581.0481), _ASEC360) * _ASEC2RAD - f = math.fmod((335779.526232 + t*1739527262.8478), _ASEC360) * _ASEC2RAD - d = math.fmod((1072260.70369 + t*1602961601.2090), _ASEC360) * _ASEC2RAD - om = math.fmod((450160.398036 - t*6962890.5431), _ASEC360) * _ASEC2RAD - dp = 0 - de = 0 - $ASTRO_IAULOOP() - self.dpsi = -0.000135 + (dp * 1.0e-7) - self.deps = +0.000388 + (de * 1.0e-7) - -def _mean_obliq(tt): - t = tt / 36525 - asec = ( - (((( - 0.0000000434 * t - - 0.000000576 ) * t - + 0.00200340 ) * t - - 0.0001831 ) * t - - 46.836769 ) * t + 84381.406 - ) - return asec / 3600.0 - -class _e_tilt: - def __init__(self, time): - e = _iau2000b(time) - self.dpsi = e.dpsi - self.deps = e.deps - self.mobl = _mean_obliq(time.tt) - self.tobl = self.mobl + (e.deps / 3600.0) - self.tt = time.tt - self.ee = e.dpsi * math.cos(self.mobl * _DEG2RAD) / 15.0 - -def _ecl2equ_vec(time, ecl): - obl = _mean_obliq(time.tt) * _DEG2RAD - cos_obl = math.cos(obl) - sin_obl = math.sin(obl) - return [ - ecl[0], - ecl[1]*cos_obl - ecl[2]*sin_obl, - ecl[1]*sin_obl + ecl[2]*cos_obl - ] - -def _precession(tt1, pos1, tt2): - eps0 = 84381.406 - if tt1 != 0 and tt2 != 0: - raise Error('One of (tt1, tt2) must be zero.') - t = (tt2 - tt1) / 36525 - if tt2 == 0: - t = -t - - psia = (((((- 0.0000000951 * t - + 0.000132851 ) * t - - 0.00114045 ) * t - - 1.0790069 ) * t - + 5038.481507 ) * t) - - omegaa = (((((+ 0.0000003337 * t - - 0.000000467 ) * t - - 0.00772503 ) * t - + 0.0512623 ) * t - - 0.025754 ) * t + eps0) - - chia = (((((- 0.0000000560 * t - + 0.000170663 ) * t - - 0.00121197 ) * t - - 2.3814292 ) * t - + 10.556403 ) * t) - - eps0 *= _ASEC2RAD - psia *= _ASEC2RAD - omegaa *= _ASEC2RAD - chia *= _ASEC2RAD - - sa = math.sin(eps0) - ca = math.cos(eps0) - sb = math.sin(-psia) - cb = math.cos(-psia) - sc = math.sin(-omegaa) - cc = math.cos(-omegaa) - sd = math.sin(chia) - cd = math.cos(chia) - - xx = cd * cb - sb * sd * cc - yx = cd * sb * ca + sd * cc * cb * ca - sa * sd * sc - zx = cd * sb * sa + sd * cc * cb * sa + ca * sd * sc - xy = -sd * cb - sb * cd * cc - yy = -sd * sb * ca + cd * cc * cb * ca - sa * cd * sc - zy = -sd * sb * sa + cd * cc * cb * sa + ca * cd * sc - xz = sb * sc - yz = -sc * cb * ca - sa * cc - zz = -sc * cb * sa + cc * ca - if tt2 == 0.0: - # Perform rotation from other epoch to J2000.0. - return [ - xx * pos1[0] + xy * pos1[1] + xz * pos1[2], - yx * pos1[0] + yy * pos1[1] + yz * pos1[2], - zx * pos1[0] + zy * pos1[1] + zz * pos1[2] - ] - - # Perform rotation from J2000.0 to other epoch. - return [ - xx * pos1[0] + yx * pos1[1] + zx * pos1[2], - xy * pos1[0] + yy * pos1[1] + zy * pos1[2], - xz * pos1[0] + yz * pos1[1] + zz * pos1[2] - ] - -class Equatorial: - def __init__(self, ra, dec, dist): - self.ra = ra - self.dec = dec - self.dist = dist - -def _vector2radec(pos): - xyproj = pos[0]*pos[0] + pos[1]*pos[1] - dist = math.sqrt(xyproj + pos[2]*pos[2]) - if xyproj == 0.0: - if pos[2] == 0.0: - # Indeterminate coordinates: pos vector has zero length. - raise Error('Cannot convert vector to polar coordinates') - ra = 0.0 - if pos[2] < 0.0: - dec = -90.0 - else: - dec = +90.0 - else: - ra = math.atan2(pos[1], pos[0]) / (_DEG2RAD * 15) - if ra < 0: - ra += 24 - dec = _RAD2DEG * math.atan2(pos[2], math.sqrt(xyproj)) - return Equatorial(ra, dec, dist) - - -def _nutation(time, direction, inpos): - tilt = time._etilt() - oblm = tilt.mobl * _DEG2RAD - oblt = tilt.tobl * _DEG2RAD - psi = tilt.dpsi * _ASEC2RAD - cobm = math.cos(oblm) - sobm = math.sin(oblm) - cobt = math.cos(oblt) - sobt = math.sin(oblt) - cpsi = math.cos(psi) - spsi = math.sin(psi) - - xx = cpsi - yx = -spsi * cobm - zx = -spsi * sobm - xy = spsi * cobt - yy = cpsi * cobm * cobt + sobm * sobt - zy = cpsi * sobm * cobt - cobm * sobt - xz = spsi * sobt - yz = cpsi * cobm * sobt - sobm * cobt - zz = cpsi * sobm * sobt + cobm * cobt - - if direction == 0: - # forward rotation - return [ - xx * inpos[0] + yx * inpos[1] + zx * inpos[2], - xy * inpos[0] + yy * inpos[1] + zy * inpos[2], - xz * inpos[0] + yz * inpos[1] + zz * inpos[2] - ] - - # inverse rotation - return [ - xx * inpos[0] + xy * inpos[1] + xz * inpos[2], - yx * inpos[0] + yy * inpos[1] + yz * inpos[2], - zx * inpos[0] + zy * inpos[1] + zz * inpos[2] - ] - -def _era(time): # Earth Rotation Angle - thet1 = 0.7790572732640 + 0.00273781191135448 * time.ut - thet3 = math.fmod(time.ut, 1.0) - theta = 360.0 * math.fmod((thet1 + thet3), 1.0) - if theta < 0.0: - theta += 360.0 - return theta - - -def _sidereal_time(time): - t = time.tt / 36525.0 - eqeq = 15.0 * time._etilt().ee # Replace with eqeq=0 to get GMST instead of GAST (if we ever need it) - theta = _era(time) - st = (eqeq + 0.014506 + - (((( - 0.0000000368 * t - - 0.000029956 ) * t - - 0.00000044 ) * t - + 1.3915817 ) * t - + 4612.156534 ) * t) - gst = math.fmod((st/3600.0 + theta), 360.0) / 15.0 - if gst < 0.0: - gst += 24.0 - return gst - - -def _terra(observer, st): - erad_km = _ERAD / 1000.0 - df = 1.0 - 0.003352819697896 # flattening of the Earth - df2 = df * df - phi = observer.latitude * _DEG2RAD - sinphi = math.sin(phi) - cosphi = math.cos(phi) - c = 1.0 / math.sqrt(cosphi*cosphi + df2*sinphi*sinphi) - s = df2 * c - ht_km = observer.height / 1000.0 - ach = erad_km*c + ht_km - ash = erad_km*s + ht_km - stlocl = (15.0*st + observer.longitude) * _DEG2RAD - sinst = math.sin(stlocl) - cosst = math.cos(stlocl) - return [ - ach * cosphi * cosst / _KM_PER_AU, - ach * cosphi * sinst / _KM_PER_AU, - ash * sinphi / _KM_PER_AU - ] - -def _geo_pos(time, observer): - gast = _sidereal_time(time) - pos1 = _terra(observer, gast) - pos2 = _nutation(time, -1, pos1) - outpos = _precession(time.tt, pos2, 0.0) - return outpos - -def _spin(angle, pos1): - angr = angle * _DEG2RAD - cosang = math.cos(angr) - sinang = math.sin(angr) - return [ - +cosang*pos1[0] + sinang*pos1[1], - -sinang*pos1[0] + cosang*pos1[1], - pos1[2] - ] - -#---------------------------------------------------------------------------- -# BEGIN CalcMoon - -def _Array1(xmin, xmax): - return dict((key, 0j) for key in range(xmin, 1+xmax)) - -def _Array2(xmin, xmax, ymin, ymax): - return dict((key, _Array1(ymin, ymax)) for key in range(xmin, 1+xmax)) - -class _moonpos: - def __init__(self, lon, lat, dist): - self.geo_eclip_lon = lon - self.geo_eclip_lat = lat - self.distance_au = dist - -def _CalcMoon(time): - T = time.tt / 36525 - ex = _Array2(-6, 6, 1, 4) - - def Sine(phi): - return math.sin(_PI2 * phi) - - def Frac(x): - return x - math.floor(x) - - T2 = T*T - DLAM = 0 - DS = 0 - GAM1C = 0 - SINPI = 3422.7000 - S1 = Sine(0.19833+0.05611*T) - S2 = Sine(0.27869+0.04508*T) - S3 = Sine(0.16827-0.36903*T) - S4 = Sine(0.34734-5.37261*T) - S5 = Sine(0.10498-5.37899*T) - S6 = Sine(0.42681-0.41855*T) - S7 = Sine(0.14943-5.37511*T) - DL0 = 0.84*S1+0.31*S2+14.27*S3+ 7.26*S4+ 0.28*S5+0.24*S6 - DL = 2.94*S1+0.31*S2+14.27*S3+ 9.34*S4+ 1.12*S5+0.83*S6 - DLS =-6.40*S1 -1.89*S6 - DF = 0.21*S1+0.31*S2+14.27*S3-88.70*S4-15.30*S5+0.24*S6-1.86*S7 - DD = DL0-DLS - DGAM = ((-3332E-9 * Sine(0.59734-5.37261*T) - -539E-9 * Sine(0.35498-5.37899*T) - -64E-9 * Sine(0.39943-5.37511*T))) - - L0 = _PI2*Frac(0.60643382+1336.85522467*T-0.00000313*T2) + DL0/_ARC - L = _PI2*Frac(0.37489701+1325.55240982*T+0.00002565*T2) + DL /_ARC - LS = _PI2*Frac(0.99312619+ 99.99735956*T-0.00000044*T2) + DLS/_ARC - F = _PI2*Frac(0.25909118+1342.22782980*T-0.00000892*T2) + DF /_ARC - D = _PI2*Frac(0.82736186+1236.85308708*T-0.00000397*T2) + DD /_ARC - - I = 1 - while I <= 4: - if I == 1: - ARG=L; MAX=4; FAC=1.000002208 - elif I == 2: - ARG=LS; MAX=3; FAC=0.997504612-0.002495388*T - elif I == 3: - ARG=F; MAX=4; FAC=1.000002708+139.978*DGAM - else: - ARG=D; MAX=6; FAC=1.0 - - ex[0][I] = complex(1, 0) - ex[1][I] = complex(FAC * math.cos(ARG), FAC * math.sin(ARG)) - - J = 2 - while J <= MAX: - ex[J][I] = ex[J-1][I] * ex[1][I] - J += 1 - - J = 1 - while J <= MAX: - ex[-J][I] = ex[J][I].conjugate() - J += 1 - - I += 1 - - $ASTRO_ADDSOL( 13.902, 14.06,-0.001, 0.2607,0, 0, 0, 4) - $ASTRO_ADDSOL( 0.403, -4.01,+0.394, 0.0023,0, 0, 0, 3) - $ASTRO_ADDSOL( 2369.912, 2373.36,+0.601, 28.2333,0, 0, 0, 2) - $ASTRO_ADDSOL( -125.154, -112.79,-0.725, -0.9781,0, 0, 0, 1) - $ASTRO_ADDSOL( 1.979, 6.98,-0.445, 0.0433,1, 0, 0, 4) - $ASTRO_ADDSOL( 191.953, 192.72,+0.029, 3.0861,1, 0, 0, 2) - $ASTRO_ADDSOL( -8.466, -13.51,+0.455, -0.1093,1, 0, 0, 1) - $ASTRO_ADDSOL( 22639.500,22609.07,+0.079, 186.5398,1, 0, 0, 0) - $ASTRO_ADDSOL( 18.609, 3.59,-0.094, 0.0118,1, 0, 0,-1) - $ASTRO_ADDSOL( -4586.465,-4578.13,-0.077, 34.3117,1, 0, 0,-2) - $ASTRO_ADDSOL( +3.215, 5.44,+0.192, -0.0386,1, 0, 0,-3) - $ASTRO_ADDSOL( -38.428, -38.64,+0.001, 0.6008,1, 0, 0,-4) - $ASTRO_ADDSOL( -0.393, -1.43,-0.092, 0.0086,1, 0, 0,-6) - $ASTRO_ADDSOL( -0.289, -1.59,+0.123, -0.0053,0, 1, 0, 4) - $ASTRO_ADDSOL( -24.420, -25.10,+0.040, -0.3000,0, 1, 0, 2) - $ASTRO_ADDSOL( 18.023, 17.93,+0.007, 0.1494,0, 1, 0, 1) - $ASTRO_ADDSOL( -668.146, -126.98,-1.302, -0.3997,0, 1, 0, 0) - $ASTRO_ADDSOL( 0.560, 0.32,-0.001, -0.0037,0, 1, 0,-1) - $ASTRO_ADDSOL( -165.145, -165.06,+0.054, 1.9178,0, 1, 0,-2) - $ASTRO_ADDSOL( -1.877, -6.46,-0.416, 0.0339,0, 1, 0,-4) - $ASTRO_ADDSOL( 0.213, 1.02,-0.074, 0.0054,2, 0, 0, 4) - $ASTRO_ADDSOL( 14.387, 14.78,-0.017, 0.2833,2, 0, 0, 2) - $ASTRO_ADDSOL( -0.586, -1.20,+0.054, -0.0100,2, 0, 0, 1) - $ASTRO_ADDSOL( 769.016, 767.96,+0.107, 10.1657,2, 0, 0, 0) - $ASTRO_ADDSOL( +1.750, 2.01,-0.018, 0.0155,2, 0, 0,-1) - $ASTRO_ADDSOL( -211.656, -152.53,+5.679, -0.3039,2, 0, 0,-2) - $ASTRO_ADDSOL( +1.225, 0.91,-0.030, -0.0088,2, 0, 0,-3) - $ASTRO_ADDSOL( -30.773, -34.07,-0.308, 0.3722,2, 0, 0,-4) - $ASTRO_ADDSOL( -0.570, -1.40,-0.074, 0.0109,2, 0, 0,-6) - $ASTRO_ADDSOL( -2.921, -11.75,+0.787, -0.0484,1, 1, 0, 2) - $ASTRO_ADDSOL( +1.267, 1.52,-0.022, 0.0164,1, 1, 0, 1) - $ASTRO_ADDSOL( -109.673, -115.18,+0.461, -0.9490,1, 1, 0, 0) - $ASTRO_ADDSOL( -205.962, -182.36,+2.056, +1.4437,1, 1, 0,-2) - $ASTRO_ADDSOL( 0.233, 0.36, 0.012, -0.0025,1, 1, 0,-3) - $ASTRO_ADDSOL( -4.391, -9.66,-0.471, 0.0673,1, 1, 0,-4) - - $ASTRO_ADDSOL( 0.283, 1.53,-0.111, +0.0060,1,-1, 0,+4) - $ASTRO_ADDSOL( 14.577, 31.70,-1.540, +0.2302,1,-1, 0, 2) - $ASTRO_ADDSOL( 147.687, 138.76,+0.679, +1.1528,1,-1, 0, 0) - $ASTRO_ADDSOL( -1.089, 0.55,+0.021, 0.0 ,1,-1, 0,-1) - $ASTRO_ADDSOL( 28.475, 23.59,-0.443, -0.2257,1,-1, 0,-2) - $ASTRO_ADDSOL( -0.276, -0.38,-0.006, -0.0036,1,-1, 0,-3) - $ASTRO_ADDSOL( 0.636, 2.27,+0.146, -0.0102,1,-1, 0,-4) - $ASTRO_ADDSOL( -0.189, -1.68,+0.131, -0.0028,0, 2, 0, 2) - $ASTRO_ADDSOL( -7.486, -0.66,-0.037, -0.0086,0, 2, 0, 0) - $ASTRO_ADDSOL( -8.096, -16.35,-0.740, 0.0918,0, 2, 0,-2) - $ASTRO_ADDSOL( -5.741, -0.04, 0.0 , -0.0009,0, 0, 2, 2) - $ASTRO_ADDSOL( 0.255, 0.0 , 0.0 , 0.0 ,0, 0, 2, 1) - $ASTRO_ADDSOL( -411.608, -0.20, 0.0 , -0.0124,0, 0, 2, 0) - $ASTRO_ADDSOL( 0.584, 0.84, 0.0 , +0.0071,0, 0, 2,-1) - $ASTRO_ADDSOL( -55.173, -52.14, 0.0 , -0.1052,0, 0, 2,-2) - $ASTRO_ADDSOL( 0.254, 0.25, 0.0 , -0.0017,0, 0, 2,-3) - $ASTRO_ADDSOL( +0.025, -1.67, 0.0 , +0.0031,0, 0, 2,-4) - $ASTRO_ADDSOL( 1.060, 2.96,-0.166, 0.0243,3, 0, 0,+2) - $ASTRO_ADDSOL( 36.124, 50.64,-1.300, 0.6215,3, 0, 0, 0) - $ASTRO_ADDSOL( -13.193, -16.40,+0.258, -0.1187,3, 0, 0,-2) - $ASTRO_ADDSOL( -1.187, -0.74,+0.042, 0.0074,3, 0, 0,-4) - $ASTRO_ADDSOL( -0.293, -0.31,-0.002, 0.0046,3, 0, 0,-6) - $ASTRO_ADDSOL( -0.290, -1.45,+0.116, -0.0051,2, 1, 0, 2) - $ASTRO_ADDSOL( -7.649, -10.56,+0.259, -0.1038,2, 1, 0, 0) - $ASTRO_ADDSOL( -8.627, -7.59,+0.078, -0.0192,2, 1, 0,-2) - $ASTRO_ADDSOL( -2.740, -2.54,+0.022, 0.0324,2, 1, 0,-4) - $ASTRO_ADDSOL( 1.181, 3.32,-0.212, 0.0213,2,-1, 0,+2) - $ASTRO_ADDSOL( 9.703, 11.67,-0.151, 0.1268,2,-1, 0, 0) - $ASTRO_ADDSOL( -0.352, -0.37,+0.001, -0.0028,2,-1, 0,-1) - $ASTRO_ADDSOL( -2.494, -1.17,-0.003, -0.0017,2,-1, 0,-2) - $ASTRO_ADDSOL( 0.360, 0.20,-0.012, -0.0043,2,-1, 0,-4) - $ASTRO_ADDSOL( -1.167, -1.25,+0.008, -0.0106,1, 2, 0, 0) - $ASTRO_ADDSOL( -7.412, -6.12,+0.117, 0.0484,1, 2, 0,-2) - $ASTRO_ADDSOL( -0.311, -0.65,-0.032, 0.0044,1, 2, 0,-4) - $ASTRO_ADDSOL( +0.757, 1.82,-0.105, 0.0112,1,-2, 0, 2) - $ASTRO_ADDSOL( +2.580, 2.32,+0.027, 0.0196,1,-2, 0, 0) - $ASTRO_ADDSOL( +2.533, 2.40,-0.014, -0.0212,1,-2, 0,-2) - $ASTRO_ADDSOL( -0.344, -0.57,-0.025, +0.0036,0, 3, 0,-2) - $ASTRO_ADDSOL( -0.992, -0.02, 0.0 , 0.0 ,1, 0, 2, 2) - $ASTRO_ADDSOL( -45.099, -0.02, 0.0 , -0.0010,1, 0, 2, 0) - $ASTRO_ADDSOL( -0.179, -9.52, 0.0 , -0.0833,1, 0, 2,-2) - $ASTRO_ADDSOL( -0.301, -0.33, 0.0 , 0.0014,1, 0, 2,-4) - $ASTRO_ADDSOL( -6.382, -3.37, 0.0 , -0.0481,1, 0,-2, 2) - $ASTRO_ADDSOL( 39.528, 85.13, 0.0 , -0.7136,1, 0,-2, 0) - $ASTRO_ADDSOL( 9.366, 0.71, 0.0 , -0.0112,1, 0,-2,-2) - $ASTRO_ADDSOL( 0.202, 0.02, 0.0 , 0.0 ,1, 0,-2,-4) - - $ASTRO_ADDSOL( 0.415, 0.10, 0.0 , 0.0013,0, 1, 2, 0) - $ASTRO_ADDSOL( -2.152, -2.26, 0.0 , -0.0066,0, 1, 2,-2) - $ASTRO_ADDSOL( -1.440, -1.30, 0.0 , +0.0014,0, 1,-2, 2) - $ASTRO_ADDSOL( 0.384, -0.04, 0.0 , 0.0 ,0, 1,-2,-2) - $ASTRO_ADDSOL( +1.938, +3.60,-0.145, +0.0401,4, 0, 0, 0) - $ASTRO_ADDSOL( -0.952, -1.58,+0.052, -0.0130,4, 0, 0,-2) - $ASTRO_ADDSOL( -0.551, -0.94,+0.032, -0.0097,3, 1, 0, 0) - $ASTRO_ADDSOL( -0.482, -0.57,+0.005, -0.0045,3, 1, 0,-2) - $ASTRO_ADDSOL( 0.681, 0.96,-0.026, 0.0115,3,-1, 0, 0) - $ASTRO_ADDSOL( -0.297, -0.27, 0.002, -0.0009,2, 2, 0,-2) - $ASTRO_ADDSOL( 0.254, +0.21,-0.003, 0.0 ,2,-2, 0,-2) - $ASTRO_ADDSOL( -0.250, -0.22, 0.004, 0.0014,1, 3, 0,-2) - $ASTRO_ADDSOL( -3.996, 0.0 , 0.0 , +0.0004,2, 0, 2, 0) - $ASTRO_ADDSOL( 0.557, -0.75, 0.0 , -0.0090,2, 0, 2,-2) - $ASTRO_ADDSOL( -0.459, -0.38, 0.0 , -0.0053,2, 0,-2, 2) - $ASTRO_ADDSOL( -1.298, 0.74, 0.0 , +0.0004,2, 0,-2, 0) - $ASTRO_ADDSOL( 0.538, 1.14, 0.0 , -0.0141,2, 0,-2,-2) - $ASTRO_ADDSOL( 0.263, 0.02, 0.0 , 0.0 ,1, 1, 2, 0) - $ASTRO_ADDSOL( 0.426, +0.07, 0.0 , -0.0006,1, 1,-2,-2) - $ASTRO_ADDSOL( -0.304, +0.03, 0.0 , +0.0003,1,-1, 2, 0) - $ASTRO_ADDSOL( -0.372, -0.19, 0.0 , -0.0027,1,-1,-2, 2) - $ASTRO_ADDSOL( +0.418, 0.0 , 0.0 , 0.0 ,0, 0, 4, 0) - $ASTRO_ADDSOL( -0.330, -0.04, 0.0 , 0.0 ,3, 0, 2, 0) - - def ADDN(coeffn, p, q, r, s): - return coeffn * (ex[p][1] * ex[q][2] * ex[r][3] * ex[s][4]).imag - - N = 0 - N += ADDN(-526.069, 0, 0,1,-2) - N += ADDN( -3.352, 0, 0,1,-4) - N += ADDN( +44.297,+1, 0,1,-2) - N += ADDN( -6.000,+1, 0,1,-4) - N += ADDN( +20.599,-1, 0,1, 0) - N += ADDN( -30.598,-1, 0,1,-2) - N += ADDN( -24.649,-2, 0,1, 0) - N += ADDN( -2.000,-2, 0,1,-2) - N += ADDN( -22.571, 0,+1,1,-2) - N += ADDN( +10.985, 0,-1,1,-2) - - DLAM += ( - +0.82*Sine(0.7736 -62.5512*T)+0.31*Sine(0.0466 -125.1025*T) - +0.35*Sine(0.5785 -25.1042*T)+0.66*Sine(0.4591+1335.8075*T) - +0.64*Sine(0.3130 -91.5680*T)+1.14*Sine(0.1480+1331.2898*T) - +0.21*Sine(0.5918+1056.5859*T)+0.44*Sine(0.5784+1322.8595*T) - +0.24*Sine(0.2275 -5.7374*T)+0.28*Sine(0.2965 +2.6929*T) - +0.33*Sine(0.3132 +6.3368*T) - ) - S = F + DS/_ARC - lat_seconds = (1.000002708 + 139.978*DGAM)*(18518.511+1.189+GAM1C)*math.sin(S) - 6.24*math.sin(3*S) + N - return _moonpos( - _PI2 * Frac((L0+DLAM/_ARC) / _PI2), - (math.pi / (180 * 3600)) * lat_seconds, - (_ARC * (_ERAD / _AU)) / (0.999953253 * SINPI) - ) - -def GeoMoon(time): - m = _CalcMoon(time) - - # Convert geocentric ecliptic spherical coordinates to Cartesian coordinates. - dist_cos_lat = m.distance_au * math.cos(m.geo_eclip_lat) - gepos = [ - dist_cos_lat * math.cos(m.geo_eclip_lon), - dist_cos_lat * math.sin(m.geo_eclip_lon), - m.distance_au * math.sin(m.geo_eclip_lat) - ] - - # Convert ecliptic coordinates to equatorial coordinates, both in mean equinox of date. - mpos1 = _ecl2equ_vec(time, gepos) - - # Convert from mean equinox of date to J2000. - mpos2 = _precession(time.tt, mpos1, 0) - return Vector(mpos2[0], mpos2[1], mpos2[2], time) - -# END CalcMoon -#---------------------------------------------------------------------------- -# BEGIN VSOP - -_vsop = [ - # Mercury - $ASTRO_LIST_VSOP(Mercury), - - # Venus - $ASTRO_LIST_VSOP(Venus), - - # Earth - $ASTRO_LIST_VSOP(Earth), - - # Mars - $ASTRO_LIST_VSOP(Mars), - - # Jupiter - $ASTRO_LIST_VSOP(Jupiter), - - # Saturn - $ASTRO_LIST_VSOP(Saturn), - - # Uranus - $ASTRO_LIST_VSOP(Uranus), - - # Neptune - $ASTRO_LIST_VSOP(Neptune), -] - -def _CalcVsop(model, time): - spher = [] - t = time.tt / 365250.0 - for formula in model: - tpower = 1.0 - coord = 0.0 - for series in formula: - coord += tpower * sum(A * math.cos(B + C*t) for (A, B, C) in series) - tpower *= t - spher.append(coord) - - # Convert spherical coordinates to ecliptic cartesian coordinates. - r_coslat = spher[2] * math.cos(spher[1]) - ex = r_coslat * math.cos(spher[0]) - ey = r_coslat * math.sin(spher[0]) - ez = spher[2] * math.sin(spher[1]) - - # Convert ecliptic cartesian coordinates to equatorial cartesian coordinates. - vx = ex + 0.000000440360*ey - 0.000000190919*ez - vy = -0.000000479966*ex + 0.917482137087*ey - 0.397776982902*ez - vz = 0.397776982902*ey + 0.917482137087*ez - return Vector(vx, vy, vz, time) - -def _CalcEarth(time): - return _CalcVsop(_vsop[BODY_EARTH], time) - -# END VSOP -#---------------------------------------------------------------------------- -# BEGIN CHEBYSHEV - -_pluto = $ASTRO_LIST_CHEBYSHEV(8) - -def _ChebScale(t_min, t_max, t): - return (2*t - (t_max + t_min)) / (t_max - t_min) - -def _CalcChebyshev(model, time): - # Search for a record that overlaps the given time value. - for record in model: - x = _ChebScale(record['tt'], record['tt'] + record['ndays'], time.tt) - if -1 <= x <= +1: - coeff = record['coeff'] - pos = [] - for d in range(3): - p0 = 1 - sum = coeff[0][d] - p1 = x - sum += coeff[1][d] * p1 - for k in range(2, len(coeff)): - p2 = (2 * x * p1) - p0 - sum += coeff[k][d] * p2 - p0 = p1 - p1 = p2 - pos.append(sum - coeff[0][d]/2) - return Vector(pos[0], pos[1], pos[2], time) - raise Error('Cannot extrapolate Chebyshev model for given Terrestrial Time: {}'.format(time.tt)) - -# END CHEBYSHEV -#---------------------------------------------------------------------------- -# BEGIN Search - -def _QuadInterp(tm, dt, fa, fm, fb): - Q = (fb + fa)/2 - fm - R = (fb - fa)/2 - S = fm - - if Q == 0: - # This is a line, not a parabola. - if R == 0: - # This is a HORIZONTAL line... can't make progress! - return None - x = -S / R - if not (-1 <= x <= +1): - return None # out of bounds - else: - # It really is a parabola. Find roots x1, x2. - u = R*R - 4*Q*S - if u <= 0: - return None - ru = math.sqrt(u) - x1 = (-R + ru) / (2 * Q) - x2 = (-R - ru) / (2 * Q) - - if -1 <= x1 <= +1: - if -1 <= x2 <= +1: - # Two solutions... so parabola intersects twice. - return None - x = x1 - elif -1 <= x2 <= +1: - x = x2 - else: - return None - - t = tm + x*dt - df_dt = (2*Q*x + R) / dt - return (x, t, df_dt) - -def Search(func, context, t1, t2, dt_tolerance_seconds): - dt_days = abs(dt_tolerance_seconds / _SECONDS_PER_DAY) - f1 = func(context, t1) - f2 = func(context, t2) - iter = 0 - iter_limit = 20 - calc_fmid = True - while True: - iter += 1 - if iter > iter_limit: - raise Error('Excessive iteration in Search') - - dt = (t2.tt - t1.tt) / 2.0 - tmid = t1.AddDays(dt) - if abs(dt) < dt_days: - # We are close enough to the event to stop the search. - return tmid - - if calc_fmid: - fmid = func(context, tmid) - else: - # We already have the correct value of fmid from the previous loop. - calc_fmid = True - - # Quadratic interpolation: - # Try to find a parabola that passes through the 3 points we have sampled: - # (t1,f1), (tmid,fmid), (t2,f2). - q = _QuadInterp(tmid.ut, t2.ut - tmid.ut, f1, fmid, f2) - if q: - (q_x, q_ut, q_df_dt) = q - tq = Time(q_ut) - fq = func(context, tq) - if q_df_dt != 0.0: - dt_guess = abs(fq / q_df_dt) - if dt_guess < dt_days: - # The estimated time error is small enough that we can quit now. - return tq - - # Try guessing a tighter boundary with the interpolated root at the center. - dt_guess *= 1.2 - if dt_guess < dt / 10.0: - tleft = tq.AddDays(-dt_guess) - tright = tq.AddDays(+dt_guess) - if (tleft.ut - t1.ut)*(tleft.ut - t2.ut) < 0.0: - if (tright.ut - t1.ut)*(tright.ut - t2.ut) < 0.0: - fleft = func(context, tleft) - fright = func(context, tright) - if fleft < 0.0 and fright >= 0.0: - f1 = fleft - f2 = fright - t1 = tleft - t2 = tright - fmid = fq - calc_fmid = False - continue - - # Quadratic interpolation attempt did not work out. - # Just divide the region in two parts and pick whichever one appears to contain a root. - if f1 < 0.0 and fmid >= 0.0: - t2 = tmid - f2 = fmid - continue - - if fmid < 0.0 and f2 >= 0.0: - t1 = tmid - f1 = fmid - continue - - # Either there is no ascending zero-crossing in this range - # or the search window is too wide (more than one zero-crossing). - return None - -# END Search -#---------------------------------------------------------------------------- - -def HelioVector(body, time): - if body == BODY_PLUTO: - return _CalcChebyshev(_pluto, time) - - if 0 <= body <= len(_vsop): - return _CalcVsop(_vsop[body], time) - - if body == BODY_SUN: - return Vector(0.0, 0.0, 0.0, time) - - if body == BODY_MOON: - e = _CalcEarth(time) - m = GeoMoon(time) - return Vector(e.x+m.x, e.y+m.y, e.z+m.z, time) - - raise InvalidBodyError() - - -def GeoVector(body, time, aberration): - if body == BODY_MOON: - return GeoMoon(time) - - if body == BODY_EARTH: - return Vector(0.0, 0.0, 0.0, time) - - if not aberration: - # No aberration, so calculate Earth's position once, at the time of observation. - earth = _CalcEarth(time) - - # Correct for light-travel time, to get position of body as seen from Earth's center. - ltime = time - for iter in range(10): - h = HelioVector(body, ltime) - if aberration: - # Include aberration, so make a good first-order approximation - # by backdating the Earth's position also. - # This is confusing, but it works for objects within the Solar System - # because the distance the Earth moves in that small amount of light - # travel time (a few minutes to a few hours) is well approximated - # by a line segment that substends the angle seen from the remote - # body viewing Earth. That angle is pretty close to the aberration - # angle of the moving Earth viewing the remote body. - # In other words, both of the following approximate the aberration angle: - # (transverse distance Earth moves) / (distance to body) - # (transverse speed of Earth) / (speed of light). - earth = _CalcEarth(ltime) - - geo = Vector(h.x-earth.x, h.y-earth.y, h.z-earth.z, time) - if body == BODY_SUN: - # The Sun's heliocentric coordinates are always (0,0,0). No need to correct. - return geo - - ltime2 = time.AddDays(-geo.Length() / _C_AUDAY) - dt = abs(ltime2.tt - ltime.tt) - if dt < 1.0e-9: - return geo - - ltime = ltime2 - - raise Error('Light-travel time solver did not converge: dt={}'.format(dt)) - - -def Equator(body, time, observer, ofdate, aberration): - gc_observer = _geo_pos(time, observer) - gc = GeoVector(body, time, aberration) - j2000 = [ - gc.x - gc_observer[0], - gc.y - gc_observer[1], - gc.z - gc_observer[2] - ] - if not ofdate: - return _vector2radec(j2000) - temp = _precession(0, j2000, time.tt) - datevect = _nutation(time, 0, temp) - return _vector2radec(datevect) - -REFRACTION_NONE = 0 -REFRACTION_NORMAL = 1 -REFRACTION_JPLHOR = 2 - -class HorizontalCoordinates: - def __init__(self, azimuth, altitude, ra, dec): - self.azimuth = azimuth - self.altitude = altitude - self.ra = ra - self.dec = dec - -def Horizon(time, observer, ra, dec, refraction): - if not (REFRACTION_NONE <= refraction <= REFRACTION_JPLHOR): - raise Error('Invalid refraction type: ' + str(refraction)) - - sinlat = math.sin(observer.latitude * _DEG2RAD) - coslat = math.cos(observer.latitude * _DEG2RAD) - sinlon = math.sin(observer.longitude * _DEG2RAD) - coslon = math.cos(observer.longitude * _DEG2RAD) - sindc = math.sin(dec * _DEG2RAD) - cosdc = math.cos(dec * _DEG2RAD) - sinra = math.sin(ra * 15 * _DEG2RAD) - cosra = math.cos(ra * 15 * _DEG2RAD) - - uze = [coslat*coslon, coslat*sinlon, sinlat] - une = [-sinlat*coslon, -sinlat*sinlon, coslat] - uwe = [sinlon, -coslon, 0.0] - - angle = -15.0 * _sidereal_time(time) - uz = _spin(angle, uze) - un = _spin(angle, une) - uw = _spin(angle, uwe) - - p = [cosdc*cosra, cosdc*sinra, sindc] - - pz = p[0]*uz[0] + p[1]*uz[1] + p[2]*uz[2] - pn = p[0]*un[0] + p[1]*un[1] + p[2]*un[2] - pw = p[0]*uw[0] + p[1]*uw[1] + p[2]*uw[2] - - proj = math.sqrt(pn*pn + pw*pw) - az = 0.0 - if proj > 0.0: - az = -math.atan2(pw, pn) * _RAD2DEG - if az < 0: - az += 360 - if az >= 360: - az -= 360 - zd = math.atan2(proj, pz) * _RAD2DEG - hor_ra = ra - hor_dec = dec - - if refraction != REFRACTION_NONE: - zd0 = zd - - # http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf - # JPL Horizons says it uses refraction algorithm from - # Meeus "Astronomical Algorithms", 1991, p. 101-102. - # I found the following Go implementation: - # https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go - # This is a translation from the function "Saemundsson" there. - # I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon. - # This is important because the 'refr' formula below goes crazy near hd = -5.11. - hd = 90.0 - zd - if hd < -1.0: - hd = -1.0 - - refr = (1.02 / math.tan((hd+10.3/(hd+5.11))*_DEG2RAD)) / 60.0 - - if refraction == REFRACTION_NORMAL and zd > 91.0: - # In "normal" mode we gradually reduce refraction toward the nadir - # so that we never get an altitude angle less than -90 degrees. - # When horizon angle is -1 degrees, zd = 91, and the factor is exactly 1. - # As zd approaches 180 (the nadir), the fraction approaches 0 linearly. - refr *= (180.0 - zd) / 89.0 - - zd -= refr - - if refr > 0.0 and zd > 3.0e-4: - sinzd = math.sin(zd * _DEG2RAD) - coszd = math.cos(zd * _DEG2RAD) - sinzd0 = math.sin(zd0 * _DEG2RAD) - coszd0 = math.cos(zd0 * _DEG2RAD) - - pr = [(((p[j] - coszd0 * uz[j]) / sinzd0)*sinzd + uz[j]*coszd) for j in range(3)] - proj = math.sqrt(pr[0]*pr[0] + pr[1]*pr[1]) - if proj > 0: - hor_ra = math.atan2(pr[1], pr[0]) * _RAD2DEG / 15 - if hor_ra < 0: - hor_ra += 24 - if hor_ra >= 24: - hor_ra -= 24 - else: - hor_ra = 0 - hor_dec = math.atan2(pr[2], proj) * _RAD2DEG - - return HorizontalCoordinates(az, 90.0 - zd, hor_ra, hor_dec) - -class EclipticCoordinates: - def __init__(self, ex, ey, ez, elat, elon): - self.ex = ex - self.ey = ey - self.ez = ez - self.elat = elat - self.elon = elon - -def _RotateEquatorialToEcliptic(pos, obliq_radians): - cos_ob = math.cos(obliq_radians) - sin_ob = math.sin(obliq_radians) - ex = +pos[0] - ey = +pos[1]*cos_ob + pos[2]*sin_ob - ez = -pos[1]*sin_ob + pos[2]*cos_ob - xyproj = math.sqrt(ex*ex + ey*ey) - if xyproj > 0.0: - elon = _RAD2DEG * math.atan2(ey, ex) - if elon < 0.0: - elon += 360.0 - else: - elon = 0.0 - elat = _RAD2DEG * math.atan2(ez, xyproj) - return EclipticCoordinates(ex, ey, ez, elat, elon) - -def SunPosition(time): - # Correct for light travel time from the Sun. - # Otherwise season calculations (equinox, solstice) will all be early by about 8 minutes! - adjusted_time = time.AddDays(-1.0 / _C_AUDAY) - earth2000 = _CalcEarth(adjusted_time) - sun2000 = [-earth2000.x, -earth2000.y, -earth2000.z] - - # Convert to equatorial Cartesian coordinates of date. - stemp = _precession(0.0, sun2000, adjusted_time.tt) - sun_ofdate = _nutation(adjusted_time, 0, stemp) - - # Convert equatorial coordinates to ecliptic coordinates. - true_obliq = _DEG2RAD * adjusted_time._etilt().tobl - return _RotateEquatorialToEcliptic(sun_ofdate, true_obliq) - -def Ecliptic(equ): - # Based on NOVAS functions equ2ecl() and equ2ecl_vec(). - ob2000 = 0.40909260059599012 # mean obliquity of the J2000 ecliptic in radians - return _RotateEquatorialToEcliptic([equ.x, equ.y, equ.z], ob2000) - -def EclipticLongitude(body, time): - if body == BODY_SUN: - raise InvalidBodyError() - hv = HelioVector(body, time) - eclip = Ecliptic(hv) - return eclip.elon - -def AngleFromSun(body, time): - if body == BODY_EARTH: - raise EarthNotAllowedError() - sv = GeoVector(BODY_SUN, time, True) - bv = GeoVector(body, time, True) - return _AngleBetween(sv, bv) - -def LongitudeFromSun(body, time): - if body == BODY_EARTH: - raise EarthNotAllowedError() - sv = GeoVector(BODY_SUN, time, True) - se = Ecliptic(sv) - bv = GeoVector(body, time, True) - be = Ecliptic(bv) - return _NormalizeLongitude(be.elon - se.elon) - -class ElongationEvent: - def __init__(self, time, visibility, elongation, ecliptic_separation): - self.time = time - self.visibility = visibility - self.elongation = elongation - self.ecliptic_separation = ecliptic_separation - -def Elongation(body, time): - angle = LongitudeFromSun(body, time) - if angle > 180.0: - visibility = 'morning' - esep = 360.0 - angle - else: - visibility = 'evening' - esep = angle - angle = AngleFromSun(body, time) - return ElongationEvent(time, visibility, angle, esep) - -def _rlon_offset(body, time, direction, targetRelLon): - plon = EclipticLongitude(body, time) - elon = EclipticLongitude(BODY_EARTH, time) - diff = direction * (elon - plon) - return _LongitudeOffset(diff - targetRelLon) - -def SearchRelativeLongitude(body, targetRelLon, startTime): - if body == BODY_EARTH: - raise EarthNotAllowedError() - if body == BODY_MOON or body == BODY_SUN: - raise InvalidBodyError() - syn = _SynodicPeriod(body) - direction = +1 if _IsSuperiorPlanet(body) else -1 - # Iterate until we converge on the desired event. - # Calculate the error angle, which will be a negative number of degrees, - # meaning we are "behind" the target relative longitude. - error_angle = _rlon_offset(body, startTime, direction, targetRelLon) - if error_angle > 0.0: - error_angle -= 360.0 # force searching forward in time - time = startTime - iter = 0 - while iter < 100: - # Estimate how many days in the future (positive) or past (negative) - # we have to go to get closer to the target relative longitude. - day_adjust = (-error_angle/360.0) * syn - time = time.AddDays(day_adjust) - if abs(day_adjust) * _SECONDS_PER_DAY < 1.0: - return time - prev_angle = error_angle - error_angle = _rlon_offset(body, time, direction, targetRelLon) - if abs(prev_angle) < 30.0 and prev_angle != error_angle: - # Improve convergence for Mercury/Mars (eccentric orbits) - # by adjusting the synodic period to more closely match the - # variable speed of both planets in this part of their respective orbits. - ratio = prev_angle / (prev_angle - error_angle) - if 0.5 < ratio < 2.0: - syn *= ratio - iter += 1 - raise NoConvergeError() - -def _neg_elong_slope(body, time): - dt = 0.1 - t1 = time.AddDays(-dt/2.0) - t2 = time.AddDays(+dt/2.0) - e1 = AngleFromSun(body, t1) - e2 = AngleFromSun(body, t2) - return (e1 - e2)/dt - -def SearchMaxElongation(body, startTime): - if body == BODY_MERCURY: - s1 = 50.0 - s2 = 85.0 - elif body == BODY_VENUS: - s1 = 40.0 - s2 = 50.0 - else: - raise InvalidBodyError() - syn = _SynodicPeriod(body) - iter = 1 - while iter <= 2: - plon = EclipticLongitude(body, startTime) - elon = EclipticLongitude(BODY_EARTH, startTime) - rlon = _LongitudeOffset(plon - elon) # clamp to (-180, +180] - - # The slope function is not well-behaved when rlon is near 0 degrees or 180 degrees - # because there is a cusp there that causes a discontinuity in the derivative. - # So we need to guard against searching near such times. - if rlon >= -s1 and rlon < +s1: - # Seek to the window [+s1, +s2]. - adjust_days = 0.0 - # Search forward for the time t1 when rel lon = +s1. - rlon_lo = +s1 - # Search forward for the time t2 when rel lon = +s2. - rlon_hi = +s2 - elif rlon > +s2 or rlon < -s2: - # Seek to the next search window at [-s2, -s1]. - adjust_days = 0.0 - # Search forward for the time t1 when rel lon = -s2. - rlon_lo = -s2 - # Search forward for the time t2 when rel lon = -s1. - rlon_hi = -s1 - elif rlon >= 0.0: - # rlon must be in the middle of the window [+s1, +s2]. - # Search BACKWARD for the time t1 when rel lon = +s1. - adjust_days = -syn / 4.0 - rlon_lo = +s1 - rlon_hi = +s2 - # Search forward from t1 to find t2 such that rel lon = +s2. - else: - # rlon must be in the middle of the window [-s2, -s1]. - # Search BACKWARD for the time t1 when rel lon = -s2. - adjust_days = -syn / 4.0 - rlon_lo = -s2 - # Search forward from t1 to find t2 such that rel lon = -s1. - rlon_hi = -s1 - - t_start = startTime.AddDays(adjust_days) - t1 = SearchRelativeLongitude(body, rlon_lo, t_start) - if t1 is None: - return None - - t2 = SearchRelativeLongitude(body, rlon_hi, t1) - if t2 is None: - return None - - # Now we have a time range [t1,t2] that brackets a maximum elongation event. - # Confirm the bracketing. - m1 = _neg_elong_slope(body, t1) - if m1 >= 0.0: - raise InternalError() # there is a bug in the bracketing algorithm! - - m2 = _neg_elong_slope(body, t2) - if m2 <= 0.0: - raise InternalError() # there is a bug in the bracketing algorithm! - - # Use the generic search algorithm to home in on where the slope crosses from negative to positive. - tx = Search(_neg_elong_slope, body, t1, t2, 10.0) - if tx is None: - return None - - if tx.tt >= startTime.tt: - return Elongation(body, tx) - - # This event is in the past (earlier than startTime). - # We need to search forward from t2 to find the next possible window. - # We never need to search more than twice. - startTime = t2.AddDays(1.0) - iter += 1 - - -def _sun_offset(targetLon, time): - ecl = SunPosition(time) - return _LongitudeOffset(ecl.elon - targetLon) - -def SearchSunLongitude(targetLon, startTime, limitDays): - t2 = startTime.AddDays(limitDays) - return Search(_sun_offset, targetLon, startTime, t2, 1.0) - -def MoonPhase(time): - return LongitudeFromSun(BODY_MOON, time) - -def _moon_offset(targetLon, time): - angle = MoonPhase(time) - return _LongitudeOffset(angle - targetLon) - -def SearchMoonPhase(targetLon, startTime, limitDays): - # To avoid discontinuities in the _moon_offset function causing problems, - # we need to approximate when that function will next return 0. - # We probe it with the start time and take advantage of the fact - # that every lunar phase repeats roughly every 29.5 days. - # There is a surprising uncertainty in the quarter timing, - # due to the eccentricity of the moon's orbit. - # I have seen up to 0.826 days away from the simple prediction. - # To be safe, we take the predicted time of the event and search - # +/-0.9 days around it (a 1.8-day wide window). - # But we must return None if the final result goes beyond limitDays after startTime. - uncertainty = 0.9 - ya = _moon_offset(targetLon, startTime) - if ya > 0.0: - ya -= 360.0 # force searching forward in time, not backward - est_dt = -(_MEAN_SYNODIC_MONTH * ya) / 360.0 - dt1 = est_dt - uncertainty - if dt1 > limitDays: - return None # not possible for moon phase to occur within the specified window - dt2 = min(limitDays, est_dt + uncertainty) - t1 = startTime.AddDays(dt1) - t2 = startTime.AddDays(dt2) - return Search(_moon_offset, targetLon, t1, t2, 1.0) - -class MoonQuarter: - def __init__(self, quarter, time): - self.quarter = quarter - self.time = time - -def SearchMoonQuarter(startTime): - angle = MoonPhase(startTime) - quarter = int(1 + math.floor(angle / 90.0)) % 4 - time = SearchMoonPhase(90.0 * quarter, startTime, 10.0) - if time is None: - # The search should never fail. We should always find another lunar quarter. - raise InternalError() - return MoonQuarter(quarter, time) - -def NextMoonQuarter(mq): - # Skip 6 days past the previous found moon quarter to find the next one. - # This is less than the minimum possible increment. - # So far I have seen the interval well contained by the range (6.5, 8.3) days. - time = mq.time.AddDays(6.0) - next_mq = SearchMoonQuarter(time) - # Verify that we found the expected moon quarter. - if next_mq.quarter != (1 + mq.quarter) % 4: - raise InternalError() - return next_mq - - -class IlluminationInfo: - def __init__(self, time, mag, phase, helio_dist, geo_dist, gc, hc, ring_tilt): - self.time = time - self.mag = mag - self.phase_angle = phase - self.phase_fraction = (1.0 + math.cos(_DEG2RAD * phase)) / 2.0 - self.helio_dist = helio_dist - self.geo_dist = geo_dist - self.gc = gc - self.hc = hc - self.ring_tilt = ring_tilt - -def _MoonMagnitude(phase, helio_dist, geo_dist): - # https://astronomy.stackexchange.com/questions/10246/is-there-a-simple-analytical-formula-for-the-lunar-phase-brightness-curve - rad = phase * _DEG2RAD - mag = -12.717 + 1.49*abs(rad) + 0.0431*(rad**4) - moon_mean_distance_au = 385000.6 / _KM_PER_AU - geo_au = geo_dist / moon_mean_distance_au - mag += 5.0 * math.log10(helio_dist * geo_au) - return mag - -def _SaturnMagnitude(phase, helio_dist, geo_dist, gc, time): - # Based on formulas by Paul Schlyter found here: - # http://www.stjarnhimlen.se/comp/ppcomp.html#15 - - # We must handle Saturn's rings as a major component of its visual magnitude. - # Find geocentric ecliptic coordinates of Saturn. - eclip = Ecliptic(gc) - - ir = _DEG2RAD * 28.06 # tilt of Saturn's rings to the ecliptic, in radians - Nr = _DEG2RAD * (169.51 + (3.82e-5 * time.tt)) # ascending node of Saturn's rings, in radians - - # Find tilt of Saturn's rings, as seen from Earth. - lat = _DEG2RAD * eclip.elat - lon = _DEG2RAD * eclip.elon - tilt = math.asin(math.sin(lat)*math.cos(ir) - math.cos(lat)*math.sin(ir)*math.sin(lon-Nr)) - sin_tilt = math.sin(abs(tilt)) - - mag = -9.0 + 0.044*phase - mag += sin_tilt*(-2.6 + 1.2*sin_tilt) - mag += 5.0 * math.log10(helio_dist * geo_dist) - ring_tilt = _RAD2DEG * tilt - return (mag, ring_tilt) - -def _VisualMagnitude(body, phase, helio_dist, geo_dist): - # For Mercury and Venus, see: https://iopscience.iop.org/article/10.1086/430212 - c0 = c1 = c2 = c3 = 0 - if body == BODY_MERCURY: - c0 = -0.60; c1 = +4.98; c2 = -4.88; c3 = +3.02 - elif body == BODY_VENUS: - if phase < 163.6: - c0 = -4.47; c1 = +1.03; c2 = +0.57; c3 = +0.13 - else: - c0 = +0.98; c1 = -1.02 - elif body == BODY_MARS: - c0 = -1.52; c1 = +1.60 - elif body == BODY_JUPITER: - c0 = -9.40; c1 = +0.50 - elif body == BODY_URANUS: - c0 = -7.19; c1 = +0.25 - elif body == BODY_NEPTUNE: - c0 = -6.87 - elif body == BODY_PLUTO: - c0 = -1.00; c1 = +4.00 - else: - raise InvalidBodyError() - - x = phase / 100.0 - mag = c0 + x*(c1 + x*(c2 + x*c3)) - mag += 5.0 * math.log10(helio_dist * geo_dist) - return mag - -def Illumination(body, time): - if body == BODY_EARTH: - raise EarthNotAllowedError() - earth = _CalcEarth(time) - if body == BODY_SUN: - gc = Vector(-earth.x, -earth.y, -earth.z, time) - hc = Vector(0.0, 0.0, 0.0, time) - phase = 0.0 # placeholder value; the Sun does not have a phase angle. - else: - if body == BODY_MOON: - # For extra numeric precision, use geocentric moon formula directly. - gc = GeoMoon(time) - hc = Vector(earth.x + gc.x, earth.y + gc.y, earth.z + gc.z, time) - else: - # For planets, heliocentric vector is most direct to calculate. - hc = HelioVector(body, time) - gc = Vector(hc.x - earth.x, hc.y - earth.y, hc.z - earth.z, time) - phase = _AngleBetween(gc, hc) - - geo_dist = gc.Length() # distance from body to center of Earth - helio_dist = hc.Length() # distance from body to center of Sun - ring_tilt = None # only reported for Saturn - if body == BODY_SUN: - mag = -0.17 + 5.0*math.log10(geo_dist / _AU_PER_PARSEC) - elif body == BODY_MOON: - mag = _MoonMagnitude(phase, helio_dist, geo_dist) - elif body == BODY_SATURN: - mag, ring_tilt = _SaturnMagnitude(phase, helio_dist, geo_dist, gc, time) - else: - mag = _VisualMagnitude(body, phase, helio_dist, geo_dist) - return IlluminationInfo(time, mag, phase, helio_dist, geo_dist, gc, hc, ring_tilt) - -def _mag_slope(body, time): - # The Search() function finds a transition from negative to positive values. - # The derivative of magnitude y with respect to time t (dy/dt) - # is negative as an object gets brighter, because the magnitude numbers - # get smaller. At peak magnitude dy/dt = 0, then as the object gets dimmer, - # dy/dt > 0. - dt = 0.01 - t1 = time.AddDays(-dt/2) - t2 = time.AddDays(+dt/2) - y1 = Illumination(body, t1) - y2 = Illumination(body, t2) - return (y2.mag - y1.mag) / dt - -def SearchPeakMagnitude(body, startTime): - # s1 and s2 are relative longitudes within which peak magnitude of Venus can occur. - s1 = 10.0 - s2 = 30.0 - if body != BODY_VENUS: - raise InvalidBodyError() - - iter = 1 - while iter <= 2: - # Find current heliocentric relative longitude between the - # inferior planet and the Earth. - plon = EclipticLongitude(body, startTime) - elon = EclipticLongitude(BODY_EARTH, startTime) - rlon = _LongitudeOffset(plon - elon) - # The slope function is not well-behaved when rlon is near 0 degrees or 180 degrees - # because there is a cusp there that causes a discontinuity in the derivative. - # So we need to guard against searching near such times. - if -s1 <= rlon < +s1: - # Seek to the window [+s1, +s2]. - adjust_days = 0.0 - # Search forward for the time t1 when rel lon = +s1. - rlon_lo = +s1 - # Search forward for the time t2 when rel lon = +s2. - rlon_hi = +s2 - elif rlon >= +s2 or rlon < -s2: - # Seek to the next search window at [-s2, -s1]. - adjust_days = 0.0 - # Search forward for the time t1 when rel lon = -s2. - rlon_lo = -s2 - # Search forward for the time t2 when rel lon = -s1. - rlon_hi = -s1 - elif rlon >= 0: - # rlon must be in the middle of the window [+s1, +s2]. - # Search BACKWARD for the time t1 when rel lon = +s1. - syn = _SynodicPeriod(body) - adjust_days = -syn / 4 - rlon_lo = +s1 - # Search forward from t1 to find t2 such that rel lon = +s2. - rlon_hi = +s2 - else: - # rlon must be in the middle of the window [-s2, -s1]. - # Search BACKWARD for the time t1 when rel lon = -s2. - syn = _SynodicPeriod(body) - adjust_days = -syn / 4 - rlon_lo = -s2 - # Search forward from t1 to find t2 such that rel lon = -s1. - rlon_hi = -s1 - - t_start = startTime.AddDays(adjust_days) - t1 = SearchRelativeLongitude(body, rlon_lo, t_start) - t2 = SearchRelativeLongitude(body, rlon_hi, t1) - - # Now we have a time range [t1,t2] that brackets a maximum magnitude event. - # Confirm the bracketing. - m1 = _mag_slope(body, t1) - if m1 >= 0.0: - raise InternalError() - - m2 = _mag_slope(body, t2) - if m2 <= 0.0: - raise InternalError() - - # Use the generic search algorithm to home in on where the slope crosses from negative to positive. - tx = Search(_mag_slope, body, t1, t2, 10.0) - if tx is None: - # The search should have found the ascending root in the interval [t1, t2]. - raise InternalError() - - if tx.tt >= startTime.tt: - return Illumination(body, tx) - - # This event is in the past (earlier than startTime). - # We need to search forward from t2 to find the next possible window. - # We never need to search more than twice. - startTime = t2.AddDays(1.0) - iter += 1 - - # We should have found the peak magnitude in at most 2 iterations. - raise InternalError() - -class HourAngleEvent: - def __init__(self, time, hor): - self.time = time - self.hor = hor - -def SearchHourAngle(body, observer, hourAngle, startTime): - if body == BODY_EARTH: - raise EarthNotAllowedError() - - if hourAngle < 0.0 or hourAngle >= 24.0: - raise Error('Invalid hour angle.') - - iter = 0 - time = startTime - while True: - iter += 1 - # Calculate Greenwich Apparent Sidereal Time (GAST) at the given time. - gast = _sidereal_time(time) - ofdate = Equator(body, time, observer, True, True) - - # Calculate the adjustment needed in sidereal time to bring - # the hour angle to the desired value. - delta_sidereal_hours = math.fmod(((hourAngle + ofdate.ra - observer.longitude/15) - gast), 24.0) - if iter == 1: - # On the first iteration, always search forward in time. - if delta_sidereal_hours < 0.0: - delta_sidereal_hours += 24.0 - else: - # On subsequent iterations, we make the smallest possible adjustment, - # either forward or backward in time. - if delta_sidereal_hours < -12.0: - delta_sidereal_hours += 24.0 - elif delta_sidereal_hours > +12.0: - delta_sidereal_hours -= 24.0 - - # If the error is tolerable (less than 0.1 seconds), stop searching. - if abs(delta_sidereal_hours) * 3600.0 < 0.1: - hor = Horizon(time, observer, ofdate.ra, ofdate.dec, REFRACTION_NORMAL) - return HourAngleEvent(time, hor) - - # We need to loop another time to get more accuracy. - # Update the terrestrial time (in solar days) adjusting by sidereal time. - delta_days = (delta_sidereal_hours / 24.0) * _SOLAR_DAYS_PER_SIDEREAL_DAY - time = time.AddDays(delta_days) - -DIRECTION_RISE = +1 -DIRECTION_SET = -1 - -class _peak_altitude_context: - def __init__(self, body, direction, observer, body_radius_au): - self.body = body - self.direction = direction - self.observer = observer - self.body_radius_au = body_radius_au - -def _peak_altitude(context, time): - # Return the angular altitude above or below the horizon - # of the highest part (the peak) of the given object. - # This is defined as the apparent altitude of the center of the body plus - # the body's angular radius. - # The 'direction' parameter controls whether the angle is measured - # positive above the horizon or positive below the horizon, - # depending on whether the caller wants rise times or set times, respectively. - - ofdate = Equator(context.body, time, context.observer, True, True) - - # We calculate altitude without refraction, then add fixed refraction near the horizon. - # This gives us the time of rise/set without the extra work. - hor = Horizon(time, context.observer, ofdate.ra, ofdate.dec, REFRACTION_NONE) - alt = hor.altitude + _RAD2DEG*(context.body_radius_au / ofdate.dist) - return context.direction * (alt + _REFRACTION_NEAR_HORIZON) - -def SearchRiseSet(body, observer, direction, startTime, limitDays): - if body == BODY_EARTH: - raise EarthNotAllowedError() - elif body == BODY_SUN: - body_radius = _SUN_RADIUS_AU - elif body == BODY_MOON: - body_radius = _MOON_RADIUS_AU - else: - body_radius = 0.0 - - if direction == DIRECTION_RISE: - ha_before = 12.0 # minimum altitude (bottom) happens BEFORE the body rises. - ha_after = 0.0 # maximum altitude (culmination) happens AFTER the body rises. - elif direction == DIRECTION_SET: - ha_before = 0.0 # culmination happens BEFORE the body sets. - ha_after = 12.0 # bottom happens AFTER the body sets. - else: - raise Error('Invalid value for direction parameter') - - context = _peak_altitude_context(body, direction, observer, body_radius) - - # See if the body is currently above/below the horizon. - # If we are looking for next rise time and the body is below the horizon, - # we use the current time as the lower time bound and the next culmination - # as the upper bound. - # If the body is above the horizon, we search for the next bottom and use it - # as the lower bound and the next culmination after that bottom as the upper bound. - # The same logic applies for finding set times, only we swap the hour angles. - time_start = startTime - alt_before = _peak_altitude(context, time_start) - if alt_before > 0.0: - # We are past the sought event, so we have to wait for the next "before" event (culm/bottom). - evt_before = SearchHourAngle(body, observer, ha_before, time_start) - time_before = evt_before.time - alt_before = _peak_altitude(context, time_before) - else: - # We are before or at the sought ebvent, so we find the next "after" event (bottom/culm), - # and use the current time as the "before" event. - time_before = time_start - - evt_after = SearchHourAngle(body, observer, ha_after, time_before) - alt_after = _peak_altitude(context, evt_after.time) - - while True: - if alt_before <= 0.0 and alt_after > 0.0: - # Search between the "before time" and the "after time" for the desired event. - event_time = Search(_peak_altitude, context, time_before, evt_after.time, 1.0) - if event_time is not None: - return event_time - # We didn't find the desired event, so use the "after" time to find the next "before" event. - evt_before = SearchHourAngle(body, observer, ha_before, evt_after.time) - evt_after = SearchHourAngle(body, observer, ha_after, evt_before.time) - if evt_before.time.ut >= time_start.ut + limitDays: - return None - time_before = evt_before.time - alt_before = _peak_altitude(context, evt_before.time) - alt_after = _peak_altitude(context, evt_after.time) - -class SeasonInfo: - def __init__(self, mar_equinox, jun_solstice, sep_equinox, dec_solstice): - self.mar_equinox = mar_equinox - self.jun_solstice = jun_solstice - self.sep_equinox = sep_equinox - self.dec_solstice = dec_solstice - -def _FindSeasonChange(targetLon, year, month, day): - startTime = Time.Make(year, month, day, 0, 0, 0) - time = SearchSunLongitude(targetLon, startTime, 4.0) - if time is None: - # We should always be able to find a season change. - raise InternalError() - return time - -def Seasons(year): - mar_equinox = _FindSeasonChange(0, year, 3, 19) - jun_solstice = _FindSeasonChange(90, year, 6, 19) - sep_equinox = _FindSeasonChange(180, year, 9, 21) - dec_solstice = _FindSeasonChange(270, year, 12, 20) - return SeasonInfo(mar_equinox, jun_solstice, sep_equinox, dec_solstice) - -def _MoonDistance(time): - return _CalcMoon(time).distance_au - -def _distance_slope(direction, time): - dt = 0.001 - t1 = time.AddDays(-dt/2.0) - t2 = time.AddDays(+dt/2.0) - dist1 = _MoonDistance(t1) - dist2 = _MoonDistance(t2) - return direction * (dist2 - dist1) / dt - -APSIS_PERICENTER = 0 -APSIS_APOCENTER = 1 -APSIS_INVALID = 2 - -class Apsis: - def __init__(self, time, kind, dist_au): - self.time = time - self.kind = kind - self.dist_au = dist_au - self.dist_km = dist_au * _KM_PER_AU - -def SearchLunarApsis(startTime): - increment = 5.0 # number of days to skip on each iteration - t1 = startTime - m1 = _distance_slope(+1, t1) - iter = 0 - while iter * increment < 2.0 * _MEAN_SYNODIC_MONTH: - t2 = t1.AddDays(increment) - m2 = _distance_slope(+1, t2) - if m1 * m2 <= 0.0: - # There is a change of slope polarity within the time range [t1, t2]. - # Therefore this time range contains an apsis. - # Figure out whether it is perigee or apogee. - if m1 < 0.0 or m2 > 0.0: - # We found a minimum-distance event: perigee. - # Search the time range for the time when the slope goes from negative to positive. - apsis_time = Search(_distance_slope, +1, t1, t2, 1.0) - kind = APSIS_PERICENTER - elif m1 > 0.0 or m2 < 0.0: - # We found a maximum-distance event: apogee. - # Search the time range for the time when the slope goes from positive to negative. - apsis_time = Search(_distance_slope, -1, t1, t2, 1.0) - kind = APSIS_APOCENTER - else: - # This should never happen. It should not be possible for both slopes to be zero. - raise InternalError() - - if apsis_time is None: - # We should always be able to find a lunar apsis! - raise InternalError() - - dist = _MoonDistance(apsis_time) - return Apsis(apsis_time, kind, dist) - - # We have not yet found a slope polarity change. Keep searching. - t1 = t2 - m1 = m2 - iter += 1 - - # It should not be possible to fail to find an apsis within 2 synodic months. - raise InternalError() - - -def NextLunarApsis(apsis): - skip = 11.0 # number of days to skip to start looking for next apsis event - time = apsis.time.AddDays(skip) - next = SearchLunarApsis(time) - # Verify that we found the opposite apsis from the previous one. - if apsis.kind not in [APSIS_APOCENTER, APSIS_PERICENTER]: - raise Error('Parameter "apsis" contains an invalid "kind" value.') - if next.kind + apsis.kind != 1: - raise InternalError() - return next diff --git a/source/python/astronomy.py b/source/python/astronomy.py index f7bada4d..c0b9c712 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -319,464 +319,542 @@ class _iau2000b: om = math.fmod((450160.398036 - t*6962890.5431), _ASEC360) * _ASEC2RAD dp = 0 de = 0 + sarg = math.sin(om) carg = math.cos(om) dp += (-172064161.0 - 174666.0*t)*sarg + 33386.0*carg de += (92052331.0 + 9086.0*t)*carg + 15377.0*sarg + arg = 2.0*f - 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-13170906.0 - 1675.0*t)*sarg - 13696.0*carg de += (5730336.0 - 3015.0*t)*carg - 4587.0*sarg + arg = 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-2276413.0 - 234.0*t)*sarg + 2796.0*carg de += (978459.0 - 485.0*t)*carg + 1374.0*sarg + arg = 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (2074554.0 + 207.0*t)*sarg - 698.0*carg de += (-897492.0 + 470.0*t)*carg - 291.0*sarg + sarg = math.sin(elp) carg = math.cos(elp) dp += (1475877.0 - 3633.0*t)*sarg + 11817.0*carg de += (73871.0 - 184.0*t)*carg - 1924.0*sarg + arg = elp + 2.0*f - 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-516821.0 + 1226.0*t)*sarg - 524.0*carg de += (224386.0 - 677.0*t)*carg - 174.0*sarg + sarg = math.sin(el) carg = math.cos(el) dp += (711159.0 + 73.0*t)*sarg - 872.0*carg de += (-6750.0)*carg + 358.0*sarg + arg = 2.0*f + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-387298.0 - 367.0*t)*sarg + 380.0*carg de += (200728.0 + 18.0*t)*carg + 318.0*sarg + arg = el + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-301461.0 - 36.0*t)*sarg + 816.0*carg de += (129025.0 - 63.0*t)*carg + 367.0*sarg + arg = -elp + 2.0*f - 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (215829.0 - 494.0*t)*sarg + 111.0*carg de += (-95929.0 + 299.0*t)*carg + 132.0*sarg + arg = 2.0*f - 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (128227.0 + 137.0*t)*sarg + 181.0*carg de += (-68982.0 - 9.0*t)*carg + 39.0*sarg + arg = -el + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (123457.0 + 11.0*t)*sarg + 19.0*carg de += (-53311.0 + 32.0*t)*carg - 4.0*sarg + arg = -el + 2.0*d sarg = math.sin(arg) carg = math.cos(arg) dp += (156994.0 + 10.0*t)*sarg - 168.0*carg de += (-1235.0)*carg + 82.0*sarg + arg = el + om sarg = math.sin(arg) carg = math.cos(arg) dp += (63110.0 + 63.0*t)*sarg + 27.0*carg de += (-33228.0)*carg - 9.0*sarg + arg = -el + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-57976.0 - 63.0*t)*sarg - 189.0*carg de += (31429.0)*carg - 75.0*sarg + arg = -el + 2.0*f + 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-59641.0 - 11.0*t)*sarg + 149.0*carg de += (25543.0 - 11.0*t)*carg + 66.0*sarg + arg = el + 2.0*f + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-51613.0 - 42.0*t)*sarg + 129.0*carg de += (26366.0)*carg + 78.0*sarg + arg = -2.0*el + 2.0*f + om sarg = math.sin(arg) carg = math.cos(arg) dp += (45893.0 + 50.0*t)*sarg + 31.0*carg de += (-24236.0 - 10.0*t)*carg + 20.0*sarg + arg = 2.0*d sarg = math.sin(arg) carg = math.cos(arg) dp += (63384.0 + 11.0*t)*sarg - 150.0*carg de += (-1220.0)*carg + 29.0*sarg + arg = 2.0*f + 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-38571.0 - 1.0*t)*sarg + 158.0*carg de += (16452.0 - 11.0*t)*carg + 68.0*sarg + arg = -2.0*elp + 2.0*f - 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (32481.0)*sarg de += (-13870.0)*carg + arg = -2.0*el + 2.0*d sarg = math.sin(arg) carg = math.cos(arg) dp += (-47722.0)*sarg - 18.0*carg de += (477.0)*carg - 25.0*sarg + arg = 2.0*el + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-31046.0 - 1.0*t)*sarg + 131.0*carg de += (13238.0 - 11.0*t)*carg + 59.0*sarg + arg = el + 2.0*f - 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (28593.0)*sarg - carg de += (-12338.0 + 10.0*t)*carg - 3.0*sarg + arg = -el + 2.0*f + om sarg = math.sin(arg) carg = math.cos(arg) dp += (20441.0 + 21.0*t)*sarg + 10.0*carg de += (-10758.0)*carg - 3.0*sarg + arg = 2.0*el sarg = math.sin(arg) carg = math.cos(arg) dp += (29243.0)*sarg - 74.0*carg de += (-609.0)*carg + 13.0*sarg + arg = 2.0*f sarg = math.sin(arg) carg = math.cos(arg) dp += (25887.0)*sarg - 66.0*carg de += (-550.0)*carg + 11.0*sarg + arg = elp + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-14053.0 - 25.0*t)*sarg + 79.0*carg de += (8551.0 - 2.0*t)*carg - 45.0*sarg + arg = -el + 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (15164.0 + 10.0*t)*sarg + 11.0*carg de += (-8001.0)*carg - sarg + arg = 2.0*elp + 2.0*f - 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-15794.0 + 72.0*t)*sarg - 16.0*carg de += (6850.0 - 42.0*t)*carg - 5.0*sarg + arg = -2.0*f + 2.0*d sarg = math.sin(arg) carg = math.cos(arg) dp += (21783.0)*sarg + 13.0*carg de += (-167.0)*carg + 13.0*sarg + arg = el - 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-12873.0 - 10.0*t)*sarg - 37.0*carg de += (6953.0)*carg - 14.0*sarg + arg = -elp + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-12654.0 + 11.0*t)*sarg + 63.0*carg de += (6415.0)*carg + 26.0*sarg + arg = -el + 2.0*f + 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-10204.0)*sarg + 25.0*carg de += (5222.0)*carg + 15.0*sarg + arg = 2.0*elp sarg = math.sin(arg) carg = math.cos(arg) dp += (16707.0 - 85.0*t)*sarg - 10.0*carg de += (168.0 - 1.0*t)*carg + 10.0*sarg + arg = el + 2.0*f + 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-7691.0)*sarg + 44.0*carg de += (3268.0)*carg + 19.0*sarg + arg = -2.0*el + 2.0*f sarg = math.sin(arg) carg = math.cos(arg) dp += (-11024.0)*sarg - 14.0*carg de += (104.0)*carg + 2.0*sarg + arg = elp + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (7566.0 - 21.0*t)*sarg - 11.0*carg de += (-3250.0)*carg - 5.0*sarg + arg = 2.0*f + 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-6637.0 - 11.0*t)*sarg + 25.0*carg de += (3353.0)*carg + 14.0*sarg + arg = -elp + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-7141.0 + 21.0*t)*sarg + 8.0*carg de += (3070.0)*carg + 4.0*sarg + arg = 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-6302.0 - 11.0*t)*sarg + 2.0*carg de += (3272.0)*carg + 4.0*sarg + arg = el + 2.0*f - 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (5800.0 + 10.0*t)*sarg + 2.0*carg de += (-3045.0)*carg - sarg + arg = 2.0*el + 2.0*f - 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (6443.0)*sarg - 7.0*carg de += (-2768.0)*carg - 4.0*sarg + arg = -2.0*el + 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-5774.0 - 11.0*t)*sarg - 15.0*carg de += (3041.0)*carg - 5.0*sarg + arg = 2.0*el + 2.0*f + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-5350.0)*sarg + 21.0*carg de += (2695.0)*carg + 12.0*sarg + arg = -elp + 2.0*f - 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-4752.0 - 11.0*t)*sarg - 3.0*carg de += (2719.0)*carg - 3.0*sarg + arg = -2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-4940.0 - 11.0*t)*sarg - 21.0*carg de += (2720.0)*carg - 9.0*sarg + arg = -el - elp + 2.0*d sarg = math.sin(arg) carg = math.cos(arg) dp += (7350.0)*sarg - 8.0*carg de += (-51.0)*carg + 4.0*sarg + arg = 2.0*el - 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (4065.0)*sarg + 6.0*carg de += (-2206.0)*carg + sarg + arg = el + 2.0*d sarg = math.sin(arg) carg = math.cos(arg) dp += (6579.0)*sarg - 24.0*carg de += (-199.0)*carg + 2.0*sarg + arg = elp + 2.0*f - 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (3579.0)*sarg + 5.0*carg de += (-1900.0)*carg + sarg + arg = el - elp sarg = math.sin(arg) carg = math.cos(arg) dp += (4725.0)*sarg - 6.0*carg de += (-41.0)*carg + 3.0*sarg + arg = -2.0*el + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-3075.0)*sarg - 2.0*carg de += (1313.0)*carg - sarg + arg = 3.0*el + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-2904.0)*sarg + 15.0*carg de += (1233.0)*carg + 7.0*sarg + arg = -elp + 2.0*d sarg = math.sin(arg) carg = math.cos(arg) dp += (4348.0)*sarg - 10.0*carg de += (-81.0)*carg + 2.0*sarg + arg = el - elp + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-2878.0)*sarg + 8.0*carg de += (1232.0)*carg + 4.0*sarg + sarg = math.sin(d) carg = math.cos(d) dp += (-4230.0)*sarg + 5.0*carg de += (-20.0)*carg - 2.0*sarg + arg = -el - elp + 2.0*f + 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-2819.0)*sarg + 7.0*carg de += (1207.0)*carg + 3.0*sarg + arg = -el + 2.0*f sarg = math.sin(arg) carg = math.cos(arg) dp += (-4056.0)*sarg + 5.0*carg de += (40.0)*carg - 2.0*sarg + arg = -elp + 2.0*f + 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-2647.0)*sarg + 11.0*carg de += (1129.0)*carg + 5.0*sarg + arg = -2.0*el + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-2294.0)*sarg - 10.0*carg de += (1266.0)*carg - 4.0*sarg + arg = el + elp + 2.0*f + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (2481.0)*sarg - 7.0*carg de += (-1062.0)*carg - 3.0*sarg + arg = 2.0*el + om sarg = math.sin(arg) carg = math.cos(arg) dp += (2179.0)*sarg - 2.0*carg de += (-1129.0)*carg - 2.0*sarg + arg = -el + elp + d sarg = math.sin(arg) carg = math.cos(arg) dp += (3276.0)*sarg + carg de += (-9.0)*carg + arg = el + elp sarg = math.sin(arg) carg = math.cos(arg) dp += (-3389.0)*sarg + 5.0*carg de += (35.0)*carg - 2.0*sarg + arg = el + 2.0*f sarg = math.sin(arg) carg = math.cos(arg) dp += (3339.0)*sarg - 13.0*carg de += (-107.0)*carg + sarg + arg = -el + 2.0*f - 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-1987.0)*sarg - 6.0*carg de += (1073.0)*carg - 2.0*sarg + arg = el + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-1981.0)*sarg de += (854.0)*carg + arg = -el + d sarg = math.sin(arg) carg = math.cos(arg) dp += (4026.0)*sarg - 353.0*carg de += (-553.0)*carg - 139.0*sarg + arg = 2.0*f + d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (1660.0)*sarg - 5.0*carg de += (-710.0)*carg - 2.0*sarg + arg = -el + 2.0*f + 4.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (-1521.0)*sarg + 9.0*carg de += (647.0)*carg + 4.0*sarg + arg = -el + elp + d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (1314.0)*sarg de += (-700.0)*carg + arg = -2.0*elp + 2.0*f - 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-1283.0)*sarg de += (672.0)*carg + arg = el + 2.0*f + 2.0*d + om sarg = math.sin(arg) carg = math.cos(arg) dp += (-1331.0)*sarg + 8.0*carg de += (663.0)*carg + 4.0*sarg + arg = -2.0*el + 2.0*f + 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (1383.0)*sarg - 2.0*carg de += (-594.0)*carg - 2.0*sarg + arg = -el + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (1405.0)*sarg + 4.0*carg de += (-610.0)*carg + 2.0*sarg + arg = el + elp + 2.0*f - 2.0*d + 2.0*om sarg = math.sin(arg) carg = math.cos(arg) dp += (1290.0)*sarg de += (-556.0)*carg + self.dpsi = -0.000135 + (dp * 1.0e-7) self.deps = +0.000388 + (de * 1.0e-7) @@ -1079,246 +1157,280 @@ def _CalcMoon(time): I += 1 - - # AddSol(13.902, 14.06, -0.001, 0.2607, 0, 0, 0, 4) + + # AddSol( 13.902, 14.06,-0.001, 0.2607,0, 0, 0, 4) z = ex[4][4] DLAM += 13.902 * z.imag DS += 14.06 * z.imag GAM1C += -0.001 * z.real SINPI += 0.2607 * z.real - # AddSol(0.403, -4.01, 0.394, 0.0023, 0, 0, 0, 3) + + # AddSol( 0.403, -4.01,+0.394, 0.0023,0, 0, 0, 3) z = ex[3][4] DLAM += 0.403 * z.imag DS += -4.01 * z.imag GAM1C += 0.394 * z.real SINPI += 0.0023 * z.real - # AddSol(2369.912, 2373.36, 0.601, 28.2333, 0, 0, 0, 2) + + # AddSol( 2369.912, 2373.36,+0.601, 28.2333,0, 0, 0, 2) z = ex[2][4] DLAM += 2369.912 * z.imag DS += 2373.36 * z.imag GAM1C += 0.601 * z.real SINPI += 28.2333 * z.real - # AddSol(-125.154, -112.79, -0.725, -0.9781, 0, 0, 0, 1) + + # AddSol( -125.154, -112.79,-0.725, -0.9781,0, 0, 0, 1) z = ex[1][4] DLAM += -125.154 * z.imag DS += -112.79 * z.imag GAM1C += -0.725 * z.real SINPI += -0.9781 * z.real - # AddSol(1.979, 6.98, -0.445, 0.0433, 1, 0, 0, 4) + + # AddSol( 1.979, 6.98,-0.445, 0.0433,1, 0, 0, 4) z = ex[1][1] * ex[4][4] DLAM += 1.979 * z.imag DS += 6.98 * z.imag GAM1C += -0.445 * z.real SINPI += 0.0433 * z.real - # AddSol(191.953, 192.72, 0.029, 3.0861, 1, 0, 0, 2) + + # AddSol( 191.953, 192.72,+0.029, 3.0861,1, 0, 0, 2) z = ex[1][1] * ex[2][4] DLAM += 191.953 * z.imag DS += 192.72 * z.imag GAM1C += 0.029 * z.real SINPI += 3.0861 * z.real - # AddSol(-8.466, -13.51, 0.455, -0.1093, 1, 0, 0, 1) + + # AddSol( -8.466, -13.51,+0.455, -0.1093,1, 0, 0, 1) z = ex[1][1] * ex[1][4] DLAM += -8.466 * z.imag DS += -13.51 * z.imag GAM1C += 0.455 * z.real SINPI += -0.1093 * z.real - # AddSol(22639.5, 22609.07, 0.079, 186.5398, 1, 0, 0, 0) + + # AddSol( 22639.500,22609.07,+0.079, 186.5398,1, 0, 0, 0) z = ex[1][1] - DLAM += 22639.5 * z.imag + DLAM += 22639.500 * z.imag DS += 22609.07 * z.imag GAM1C += 0.079 * z.real SINPI += 186.5398 * z.real - # AddSol(18.609, 3.59, -0.094, 0.0118, 1, 0, 0, -1) + + # AddSol( 18.609, 3.59,-0.094, 0.0118,1, 0, 0,-1) z = ex[1][1] * ex[-1][4] DLAM += 18.609 * z.imag DS += 3.59 * z.imag GAM1C += -0.094 * z.real SINPI += 0.0118 * z.real - # AddSol(-4586.465, -4578.13, -0.077, 34.3117, 1, 0, 0, -2) + + # AddSol( -4586.465,-4578.13,-0.077, 34.3117,1, 0, 0,-2) z = ex[1][1] * ex[-2][4] DLAM += -4586.465 * z.imag DS += -4578.13 * z.imag GAM1C += -0.077 * z.real SINPI += 34.3117 * z.real - # AddSol(3.215, 5.44, 0.192, -0.0386, 1, 0, 0, -3) + + # AddSol( +3.215, 5.44,+0.192, -0.0386,1, 0, 0,-3) z = ex[1][1] * ex[-3][4] DLAM += 3.215 * z.imag DS += 5.44 * z.imag GAM1C += 0.192 * z.real SINPI += -0.0386 * z.real - # AddSol(-38.428, -38.64, 0.001, 0.6008, 1, 0, 0, -4) + + # AddSol( -38.428, -38.64,+0.001, 0.6008,1, 0, 0,-4) z = ex[1][1] * ex[-4][4] DLAM += -38.428 * z.imag DS += -38.64 * z.imag GAM1C += 0.001 * z.real SINPI += 0.6008 * z.real - # AddSol(-0.393, -1.43, -0.092, 0.0086, 1, 0, 0, -6) + + # AddSol( -0.393, -1.43,-0.092, 0.0086,1, 0, 0,-6) z = ex[1][1] * ex[-6][4] DLAM += -0.393 * z.imag DS += -1.43 * z.imag GAM1C += -0.092 * z.real SINPI += 0.0086 * z.real - # AddSol(-0.289, -1.59, 0.123, -0.0053, 0, 1, 0, 4) + + # AddSol( -0.289, -1.59,+0.123, -0.0053,0, 1, 0, 4) z = ex[1][2] * ex[4][4] DLAM += -0.289 * z.imag DS += -1.59 * z.imag GAM1C += 0.123 * z.real SINPI += -0.0053 * z.real - # AddSol(-24.42, -25.1, 0.04, -0.3, 0, 1, 0, 2) + + # AddSol( -24.420, -25.10,+0.040, -0.3000,0, 1, 0, 2) z = ex[1][2] * ex[2][4] - DLAM += -24.42 * z.imag - DS += -25.1 * z.imag - GAM1C += 0.04 * z.real - SINPI += -0.3 * z.real + DLAM += -24.420 * z.imag + DS += -25.10 * z.imag + GAM1C += 0.040 * z.real + SINPI += -0.3000 * z.real - # AddSol(18.023, 17.93, 0.007, 0.1494, 0, 1, 0, 1) + + # AddSol( 18.023, 17.93,+0.007, 0.1494,0, 1, 0, 1) z = ex[1][2] * ex[1][4] DLAM += 18.023 * z.imag DS += 17.93 * z.imag GAM1C += 0.007 * z.real SINPI += 0.1494 * z.real - # AddSol(-668.146, -126.98, -1.302, -0.3997, 0, 1, 0, 0) + + # AddSol( -668.146, -126.98,-1.302, -0.3997,0, 1, 0, 0) z = ex[1][2] DLAM += -668.146 * z.imag DS += -126.98 * z.imag GAM1C += -1.302 * z.real SINPI += -0.3997 * z.real - # AddSol(0.56, 0.32, -0.001, -0.0037, 0, 1, 0, -1) + + # AddSol( 0.560, 0.32,-0.001, -0.0037,0, 1, 0,-1) z = ex[1][2] * ex[-1][4] - DLAM += 0.56 * z.imag + DLAM += 0.560 * z.imag DS += 0.32 * z.imag GAM1C += -0.001 * z.real SINPI += -0.0037 * z.real - # AddSol(-165.145, -165.06, 0.054, 1.9178, 0, 1, 0, -2) + + # AddSol( -165.145, -165.06,+0.054, 1.9178,0, 1, 0,-2) z = ex[1][2] * ex[-2][4] DLAM += -165.145 * z.imag DS += -165.06 * z.imag GAM1C += 0.054 * z.real SINPI += 1.9178 * z.real - # AddSol(-1.877, -6.46, -0.416, 0.0339, 0, 1, 0, -4) + + # AddSol( -1.877, -6.46,-0.416, 0.0339,0, 1, 0,-4) z = ex[1][2] * ex[-4][4] DLAM += -1.877 * z.imag DS += -6.46 * z.imag GAM1C += -0.416 * z.real SINPI += 0.0339 * z.real - # AddSol(0.213, 1.02, -0.074, 0.0054, 2, 0, 0, 4) + + # AddSol( 0.213, 1.02,-0.074, 0.0054,2, 0, 0, 4) z = ex[2][1] * ex[4][4] DLAM += 0.213 * z.imag DS += 1.02 * z.imag GAM1C += -0.074 * z.real SINPI += 0.0054 * z.real - # AddSol(14.387, 14.78, -0.017, 0.2833, 2, 0, 0, 2) + + # AddSol( 14.387, 14.78,-0.017, 0.2833,2, 0, 0, 2) z = ex[2][1] * ex[2][4] DLAM += 14.387 * z.imag DS += 14.78 * z.imag GAM1C += -0.017 * z.real SINPI += 0.2833 * z.real - # AddSol(-0.586, -1.2, 0.054, -0.01, 2, 0, 0, 1) + + # AddSol( -0.586, -1.20,+0.054, -0.0100,2, 0, 0, 1) z = ex[2][1] * ex[1][4] DLAM += -0.586 * z.imag - DS += -1.2 * z.imag + DS += -1.20 * z.imag GAM1C += 0.054 * z.real - SINPI += -0.01 * z.real + SINPI += -0.0100 * z.real - # AddSol(769.016, 767.96, 0.107, 10.1657, 2, 0, 0, 0) + + # AddSol( 769.016, 767.96,+0.107, 10.1657,2, 0, 0, 0) z = ex[2][1] DLAM += 769.016 * z.imag DS += 767.96 * z.imag GAM1C += 0.107 * z.real SINPI += 10.1657 * z.real - # AddSol(1.75, 2.01, -0.018, 0.0155, 2, 0, 0, -1) + + # AddSol( +1.750, 2.01,-0.018, 0.0155,2, 0, 0,-1) z = ex[2][1] * ex[-1][4] - DLAM += 1.75 * z.imag + DLAM += 1.750 * z.imag DS += 2.01 * z.imag GAM1C += -0.018 * z.real SINPI += 0.0155 * z.real - # AddSol(-211.656, -152.53, 5.679, -0.3039, 2, 0, 0, -2) + + # AddSol( -211.656, -152.53,+5.679, -0.3039,2, 0, 0,-2) z = ex[2][1] * ex[-2][4] DLAM += -211.656 * z.imag DS += -152.53 * z.imag GAM1C += 5.679 * z.real SINPI += -0.3039 * z.real - # AddSol(1.225, 0.91, -0.03, -0.0088, 2, 0, 0, -3) + + # AddSol( +1.225, 0.91,-0.030, -0.0088,2, 0, 0,-3) z = ex[2][1] * ex[-3][4] DLAM += 1.225 * z.imag DS += 0.91 * z.imag - GAM1C += -0.03 * z.real + GAM1C += -0.030 * z.real SINPI += -0.0088 * z.real - # AddSol(-30.773, -34.07, -0.308, 0.3722, 2, 0, 0, -4) + + # AddSol( -30.773, -34.07,-0.308, 0.3722,2, 0, 0,-4) z = ex[2][1] * ex[-4][4] DLAM += -30.773 * z.imag DS += -34.07 * z.imag GAM1C += -0.308 * z.real SINPI += 0.3722 * z.real - # AddSol(-0.57, -1.4, -0.074, 0.0109, 2, 0, 0, -6) + + # AddSol( -0.570, -1.40,-0.074, 0.0109,2, 0, 0,-6) z = ex[2][1] * ex[-6][4] - DLAM += -0.57 * z.imag - DS += -1.4 * z.imag + DLAM += -0.570 * z.imag + DS += -1.40 * z.imag GAM1C += -0.074 * z.real SINPI += 0.0109 * z.real - # AddSol(-2.921, -11.75, 0.787, -0.0484, 1, 1, 0, 2) + + # AddSol( -2.921, -11.75,+0.787, -0.0484,1, 1, 0, 2) z = ex[1][1] * ex[1][2] * ex[2][4] DLAM += -2.921 * z.imag DS += -11.75 * z.imag GAM1C += 0.787 * z.real SINPI += -0.0484 * z.real - # AddSol(1.267, 1.52, -0.022, 0.0164, 1, 1, 0, 1) + + # AddSol( +1.267, 1.52,-0.022, 0.0164,1, 1, 0, 1) z = ex[1][1] * ex[1][2] * ex[1][4] DLAM += 1.267 * z.imag DS += 1.52 * z.imag GAM1C += -0.022 * z.real SINPI += 0.0164 * z.real - # AddSol(-109.673, -115.18, 0.461, -0.949, 1, 1, 0, 0) + + # AddSol( -109.673, -115.18,+0.461, -0.9490,1, 1, 0, 0) z = ex[1][1] * ex[1][2] DLAM += -109.673 * z.imag DS += -115.18 * z.imag GAM1C += 0.461 * z.real - SINPI += -0.949 * z.real + SINPI += -0.9490 * z.real - # AddSol(-205.962, -182.36, 2.056, 1.4437, 1, 1, 0, -2) + + # AddSol( -205.962, -182.36,+2.056, +1.4437,1, 1, 0,-2) z = ex[1][1] * ex[1][2] * ex[-2][4] DLAM += -205.962 * z.imag DS += -182.36 * z.imag GAM1C += 2.056 * z.real SINPI += 1.4437 * z.real - # AddSol(0.233, 0.36, 0.012, -0.0025, 1, 1, 0, -3) + + # AddSol( 0.233, 0.36, 0.012, -0.0025,1, 1, 0,-3) z = ex[1][1] * ex[1][2] * ex[-3][4] DLAM += 0.233 * z.imag DS += 0.36 * z.imag GAM1C += 0.012 * z.real SINPI += -0.0025 * z.real - # AddSol(-4.391, -9.66, -0.471, 0.0673, 1, 1, 0, -4) + + # AddSol( -4.391, -9.66,-0.471, 0.0673,1, 1, 0,-4) z = ex[1][1] * ex[1][2] * ex[-4][4] DLAM += -4.391 * z.imag DS += -9.66 * z.imag @@ -1326,448 +1438,518 @@ def _CalcMoon(time): SINPI += 0.0673 * z.real - # AddSol(0.283, 1.53, -0.111, 0.006, 1, -1, 0, 4) + + # AddSol( 0.283, 1.53,-0.111, +0.0060,1,-1, 0,+4) z = ex[1][1] * ex[-1][2] * ex[4][4] DLAM += 0.283 * z.imag DS += 1.53 * z.imag GAM1C += -0.111 * z.real - SINPI += 0.006 * z.real + SINPI += 0.0060 * z.real - # AddSol(14.577, 31.7, -1.54, 0.2302, 1, -1, 0, 2) + + # AddSol( 14.577, 31.70,-1.540, +0.2302,1,-1, 0, 2) z = ex[1][1] * ex[-1][2] * ex[2][4] DLAM += 14.577 * z.imag - DS += 31.7 * z.imag - GAM1C += -1.54 * z.real + DS += 31.70 * z.imag + GAM1C += -1.540 * z.real SINPI += 0.2302 * z.real - # AddSol(147.687, 138.76, 0.679, 1.1528, 1, -1, 0, 0) + + # AddSol( 147.687, 138.76,+0.679, +1.1528,1,-1, 0, 0) z = ex[1][1] * ex[-1][2] DLAM += 147.687 * z.imag DS += 138.76 * z.imag GAM1C += 0.679 * z.real SINPI += 1.1528 * z.real - # AddSol(-1.089, 0.55, 0.021, 0.0, 1, -1, 0, -1) + + # AddSol( -1.089, 0.55,+0.021, 0.0 ,1,-1, 0,-1) z = ex[1][1] * ex[-1][2] * ex[-1][4] DLAM += -1.089 * z.imag DS += 0.55 * z.imag GAM1C += 0.021 * z.real - # AddSol(28.475, 23.59, -0.443, -0.2257, 1, -1, 0, -2) + + # AddSol( 28.475, 23.59,-0.443, -0.2257,1,-1, 0,-2) z = ex[1][1] * ex[-1][2] * ex[-2][4] DLAM += 28.475 * z.imag DS += 23.59 * z.imag GAM1C += -0.443 * z.real SINPI += -0.2257 * z.real - # AddSol(-0.276, -0.38, -0.006, -0.0036, 1, -1, 0, -3) + + # AddSol( -0.276, -0.38,-0.006, -0.0036,1,-1, 0,-3) z = ex[1][1] * ex[-1][2] * ex[-3][4] DLAM += -0.276 * z.imag DS += -0.38 * z.imag GAM1C += -0.006 * z.real SINPI += -0.0036 * z.real - # AddSol(0.636, 2.27, 0.146, -0.0102, 1, -1, 0, -4) + + # AddSol( 0.636, 2.27,+0.146, -0.0102,1,-1, 0,-4) z = ex[1][1] * ex[-1][2] * ex[-4][4] DLAM += 0.636 * z.imag DS += 2.27 * z.imag GAM1C += 0.146 * z.real SINPI += -0.0102 * z.real - # AddSol(-0.189, -1.68, 0.131, -0.0028, 0, 2, 0, 2) + + # AddSol( -0.189, -1.68,+0.131, -0.0028,0, 2, 0, 2) z = ex[2][2] * ex[2][4] DLAM += -0.189 * z.imag DS += -1.68 * z.imag GAM1C += 0.131 * z.real SINPI += -0.0028 * z.real - # AddSol(-7.486, -0.66, -0.037, -0.0086, 0, 2, 0, 0) + + # AddSol( -7.486, -0.66,-0.037, -0.0086,0, 2, 0, 0) z = ex[2][2] DLAM += -7.486 * z.imag DS += -0.66 * z.imag GAM1C += -0.037 * z.real SINPI += -0.0086 * z.real - # AddSol(-8.096, -16.35, -0.74, 0.0918, 0, 2, 0, -2) + + # AddSol( -8.096, -16.35,-0.740, 0.0918,0, 2, 0,-2) z = ex[2][2] * ex[-2][4] DLAM += -8.096 * z.imag DS += -16.35 * z.imag - GAM1C += -0.74 * z.real + GAM1C += -0.740 * z.real SINPI += 0.0918 * z.real - # AddSol(-5.741, -0.04, 0.0, -0.0009, 0, 0, 2, 2) + + # AddSol( -5.741, -0.04, 0.0 , -0.0009,0, 0, 2, 2) z = ex[2][3] * ex[2][4] DLAM += -5.741 * z.imag DS += -0.04 * z.imag SINPI += -0.0009 * z.real - # AddSol(0.255, 0.0, 0.0, 0.0, 0, 0, 2, 1) + + # AddSol( 0.255, 0.0 , 0.0 , 0.0 ,0, 0, 2, 1) z = ex[2][3] * ex[1][4] DLAM += 0.255 * z.imag - # AddSol(-411.608, -0.2, 0.0, -0.0124, 0, 0, 2, 0) + + # AddSol( -411.608, -0.20, 0.0 , -0.0124,0, 0, 2, 0) z = ex[2][3] DLAM += -411.608 * z.imag - DS += -0.2 * z.imag + DS += -0.20 * z.imag SINPI += -0.0124 * z.real - # AddSol(0.584, 0.84, 0.0, 0.0071, 0, 0, 2, -1) + + # AddSol( 0.584, 0.84, 0.0 , +0.0071,0, 0, 2,-1) z = ex[2][3] * ex[-1][4] DLAM += 0.584 * z.imag DS += 0.84 * z.imag SINPI += 0.0071 * z.real - # AddSol(-55.173, -52.14, 0.0, -0.1052, 0, 0, 2, -2) + + # AddSol( -55.173, -52.14, 0.0 , -0.1052,0, 0, 2,-2) z = ex[2][3] * ex[-2][4] DLAM += -55.173 * z.imag DS += -52.14 * z.imag SINPI += -0.1052 * z.real - # AddSol(0.254, 0.25, 0.0, -0.0017, 0, 0, 2, -3) + + # AddSol( 0.254, 0.25, 0.0 , -0.0017,0, 0, 2,-3) z = ex[2][3] * ex[-3][4] DLAM += 0.254 * z.imag DS += 0.25 * z.imag SINPI += -0.0017 * z.real - # AddSol(0.025, -1.67, 0.0, 0.0031, 0, 0, 2, -4) + + # AddSol( +0.025, -1.67, 0.0 , +0.0031,0, 0, 2,-4) z = ex[2][3] * ex[-4][4] DLAM += 0.025 * z.imag DS += -1.67 * z.imag SINPI += 0.0031 * z.real - # AddSol(1.06, 2.96, -0.166, 0.0243, 3, 0, 0, 2) + + # AddSol( 1.060, 2.96,-0.166, 0.0243,3, 0, 0,+2) z = ex[3][1] * ex[2][4] - DLAM += 1.06 * z.imag + DLAM += 1.060 * z.imag DS += 2.96 * z.imag GAM1C += -0.166 * z.real SINPI += 0.0243 * z.real - # AddSol(36.124, 50.64, -1.3, 0.6215, 3, 0, 0, 0) + + # AddSol( 36.124, 50.64,-1.300, 0.6215,3, 0, 0, 0) z = ex[3][1] DLAM += 36.124 * z.imag DS += 50.64 * z.imag - GAM1C += -1.3 * z.real + GAM1C += -1.300 * z.real SINPI += 0.6215 * z.real - # AddSol(-13.193, -16.4, 0.258, -0.1187, 3, 0, 0, -2) + + # AddSol( -13.193, -16.40,+0.258, -0.1187,3, 0, 0,-2) z = ex[3][1] * ex[-2][4] DLAM += -13.193 * z.imag - DS += -16.4 * z.imag + DS += -16.40 * z.imag GAM1C += 0.258 * z.real SINPI += -0.1187 * z.real - # AddSol(-1.187, -0.74, 0.042, 0.0074, 3, 0, 0, -4) + + # AddSol( -1.187, -0.74,+0.042, 0.0074,3, 0, 0,-4) z = ex[3][1] * ex[-4][4] DLAM += -1.187 * z.imag DS += -0.74 * z.imag GAM1C += 0.042 * z.real SINPI += 0.0074 * z.real - # AddSol(-0.293, -0.31, -0.002, 0.0046, 3, 0, 0, -6) + + # AddSol( -0.293, -0.31,-0.002, 0.0046,3, 0, 0,-6) z = ex[3][1] * ex[-6][4] DLAM += -0.293 * z.imag DS += -0.31 * z.imag GAM1C += -0.002 * z.real SINPI += 0.0046 * z.real - # AddSol(-0.29, -1.45, 0.116, -0.0051, 2, 1, 0, 2) + + # AddSol( -0.290, -1.45,+0.116, -0.0051,2, 1, 0, 2) z = ex[2][1] * ex[1][2] * ex[2][4] - DLAM += -0.29 * z.imag + DLAM += -0.290 * z.imag DS += -1.45 * z.imag GAM1C += 0.116 * z.real SINPI += -0.0051 * z.real - # AddSol(-7.649, -10.56, 0.259, -0.1038, 2, 1, 0, 0) + + # AddSol( -7.649, -10.56,+0.259, -0.1038,2, 1, 0, 0) z = ex[2][1] * ex[1][2] DLAM += -7.649 * z.imag DS += -10.56 * z.imag GAM1C += 0.259 * z.real SINPI += -0.1038 * z.real - # AddSol(-8.627, -7.59, 0.078, -0.0192, 2, 1, 0, -2) + + # AddSol( -8.627, -7.59,+0.078, -0.0192,2, 1, 0,-2) z = ex[2][1] * ex[1][2] * ex[-2][4] DLAM += -8.627 * z.imag DS += -7.59 * z.imag GAM1C += 0.078 * z.real SINPI += -0.0192 * z.real - # AddSol(-2.74, -2.54, 0.022, 0.0324, 2, 1, 0, -4) + + # AddSol( -2.740, -2.54,+0.022, 0.0324,2, 1, 0,-4) z = ex[2][1] * ex[1][2] * ex[-4][4] - DLAM += -2.74 * z.imag + DLAM += -2.740 * z.imag DS += -2.54 * z.imag GAM1C += 0.022 * z.real SINPI += 0.0324 * z.real - # AddSol(1.181, 3.32, -0.212, 0.0213, 2, -1, 0, 2) + + # AddSol( 1.181, 3.32,-0.212, 0.0213,2,-1, 0,+2) z = ex[2][1] * ex[-1][2] * ex[2][4] DLAM += 1.181 * z.imag DS += 3.32 * z.imag GAM1C += -0.212 * z.real SINPI += 0.0213 * z.real - # AddSol(9.703, 11.67, -0.151, 0.1268, 2, -1, 0, 0) + + # AddSol( 9.703, 11.67,-0.151, 0.1268,2,-1, 0, 0) z = ex[2][1] * ex[-1][2] DLAM += 9.703 * z.imag DS += 11.67 * z.imag GAM1C += -0.151 * z.real SINPI += 0.1268 * z.real - # AddSol(-0.352, -0.37, 0.001, -0.0028, 2, -1, 0, -1) + + # AddSol( -0.352, -0.37,+0.001, -0.0028,2,-1, 0,-1) z = ex[2][1] * ex[-1][2] * ex[-1][4] DLAM += -0.352 * z.imag DS += -0.37 * z.imag GAM1C += 0.001 * z.real SINPI += -0.0028 * z.real - # AddSol(-2.494, -1.17, -0.003, -0.0017, 2, -1, 0, -2) + + # AddSol( -2.494, -1.17,-0.003, -0.0017,2,-1, 0,-2) z = ex[2][1] * ex[-1][2] * ex[-2][4] DLAM += -2.494 * z.imag DS += -1.17 * z.imag GAM1C += -0.003 * z.real SINPI += -0.0017 * z.real - # AddSol(0.36, 0.2, -0.012, -0.0043, 2, -1, 0, -4) + + # AddSol( 0.360, 0.20,-0.012, -0.0043,2,-1, 0,-4) z = ex[2][1] * ex[-1][2] * ex[-4][4] - DLAM += 0.36 * z.imag - DS += 0.2 * z.imag + DLAM += 0.360 * z.imag + DS += 0.20 * z.imag GAM1C += -0.012 * z.real SINPI += -0.0043 * z.real - # AddSol(-1.167, -1.25, 0.008, -0.0106, 1, 2, 0, 0) + + # AddSol( -1.167, -1.25,+0.008, -0.0106,1, 2, 0, 0) z = ex[1][1] * ex[2][2] DLAM += -1.167 * z.imag DS += -1.25 * z.imag GAM1C += 0.008 * z.real SINPI += -0.0106 * z.real - # AddSol(-7.412, -6.12, 0.117, 0.0484, 1, 2, 0, -2) + + # AddSol( -7.412, -6.12,+0.117, 0.0484,1, 2, 0,-2) z = ex[1][1] * ex[2][2] * ex[-2][4] DLAM += -7.412 * z.imag DS += -6.12 * z.imag GAM1C += 0.117 * z.real SINPI += 0.0484 * z.real - # AddSol(-0.311, -0.65, -0.032, 0.0044, 1, 2, 0, -4) + + # AddSol( -0.311, -0.65,-0.032, 0.0044,1, 2, 0,-4) z = ex[1][1] * ex[2][2] * ex[-4][4] DLAM += -0.311 * z.imag DS += -0.65 * z.imag GAM1C += -0.032 * z.real SINPI += 0.0044 * z.real - # AddSol(0.757, 1.82, -0.105, 0.0112, 1, -2, 0, 2) + + # AddSol( +0.757, 1.82,-0.105, 0.0112,1,-2, 0, 2) z = ex[1][1] * ex[-2][2] * ex[2][4] DLAM += 0.757 * z.imag DS += 1.82 * z.imag GAM1C += -0.105 * z.real SINPI += 0.0112 * z.real - # AddSol(2.58, 2.32, 0.027, 0.0196, 1, -2, 0, 0) + + # AddSol( +2.580, 2.32,+0.027, 0.0196,1,-2, 0, 0) z = ex[1][1] * ex[-2][2] - DLAM += 2.58 * z.imag + DLAM += 2.580 * z.imag DS += 2.32 * z.imag GAM1C += 0.027 * z.real SINPI += 0.0196 * z.real - # AddSol(2.533, 2.4, -0.014, -0.0212, 1, -2, 0, -2) + + # AddSol( +2.533, 2.40,-0.014, -0.0212,1,-2, 0,-2) z = ex[1][1] * ex[-2][2] * ex[-2][4] DLAM += 2.533 * z.imag - DS += 2.4 * z.imag + DS += 2.40 * z.imag GAM1C += -0.014 * z.real SINPI += -0.0212 * z.real - # AddSol(-0.344, -0.57, -0.025, 0.0036, 0, 3, 0, -2) + + # AddSol( -0.344, -0.57,-0.025, +0.0036,0, 3, 0,-2) z = ex[3][2] * ex[-2][4] DLAM += -0.344 * z.imag DS += -0.57 * z.imag GAM1C += -0.025 * z.real SINPI += 0.0036 * z.real - # AddSol(-0.992, -0.02, 0.0, 0.0, 1, 0, 2, 2) + + # AddSol( -0.992, -0.02, 0.0 , 0.0 ,1, 0, 2, 2) z = ex[1][1] * ex[2][3] * ex[2][4] DLAM += -0.992 * z.imag DS += -0.02 * z.imag - # AddSol(-45.099, -0.02, 0.0, -0.001, 1, 0, 2, 0) + + # AddSol( -45.099, -0.02, 0.0 , -0.0010,1, 0, 2, 0) z = ex[1][1] * ex[2][3] DLAM += -45.099 * z.imag DS += -0.02 * z.imag - SINPI += -0.001 * z.real + SINPI += -0.0010 * z.real - # AddSol(-0.179, -9.52, 0.0, -0.0833, 1, 0, 2, -2) + + # AddSol( -0.179, -9.52, 0.0 , -0.0833,1, 0, 2,-2) z = ex[1][1] * ex[2][3] * ex[-2][4] DLAM += -0.179 * z.imag DS += -9.52 * z.imag SINPI += -0.0833 * z.real - # AddSol(-0.301, -0.33, 0.0, 0.0014, 1, 0, 2, -4) + + # AddSol( -0.301, -0.33, 0.0 , 0.0014,1, 0, 2,-4) z = ex[1][1] * ex[2][3] * ex[-4][4] DLAM += -0.301 * z.imag DS += -0.33 * z.imag SINPI += 0.0014 * z.real - # AddSol(-6.382, -3.37, 0.0, -0.0481, 1, 0, -2, 2) + + # AddSol( -6.382, -3.37, 0.0 , -0.0481,1, 0,-2, 2) z = ex[1][1] * ex[-2][3] * ex[2][4] DLAM += -6.382 * z.imag DS += -3.37 * z.imag SINPI += -0.0481 * z.real - # AddSol(39.528, 85.13, 0.0, -0.7136, 1, 0, -2, 0) + + # AddSol( 39.528, 85.13, 0.0 , -0.7136,1, 0,-2, 0) z = ex[1][1] * ex[-2][3] DLAM += 39.528 * z.imag DS += 85.13 * z.imag SINPI += -0.7136 * z.real - # AddSol(9.366, 0.71, 0.0, -0.0112, 1, 0, -2, -2) + + # AddSol( 9.366, 0.71, 0.0 , -0.0112,1, 0,-2,-2) z = ex[1][1] * ex[-2][3] * ex[-2][4] DLAM += 9.366 * z.imag DS += 0.71 * z.imag SINPI += -0.0112 * z.real - # AddSol(0.202, 0.02, 0.0, 0.0, 1, 0, -2, -4) + + # AddSol( 0.202, 0.02, 0.0 , 0.0 ,1, 0,-2,-4) z = ex[1][1] * ex[-2][3] * ex[-4][4] DLAM += 0.202 * z.imag DS += 0.02 * z.imag - # AddSol(0.415, 0.1, 0.0, 0.0013, 0, 1, 2, 0) + + # AddSol( 0.415, 0.10, 0.0 , 0.0013,0, 1, 2, 0) z = ex[1][2] * ex[2][3] DLAM += 0.415 * z.imag - DS += 0.1 * z.imag + DS += 0.10 * z.imag SINPI += 0.0013 * z.real - # AddSol(-2.152, -2.26, 0.0, -0.0066, 0, 1, 2, -2) + + # AddSol( -2.152, -2.26, 0.0 , -0.0066,0, 1, 2,-2) z = ex[1][2] * ex[2][3] * ex[-2][4] DLAM += -2.152 * z.imag DS += -2.26 * z.imag SINPI += -0.0066 * z.real - # AddSol(-1.44, -1.3, 0.0, 0.0014, 0, 1, -2, 2) + + # AddSol( -1.440, -1.30, 0.0 , +0.0014,0, 1,-2, 2) z = ex[1][2] * ex[-2][3] * ex[2][4] - DLAM += -1.44 * z.imag - DS += -1.3 * z.imag + DLAM += -1.440 * z.imag + DS += -1.30 * z.imag SINPI += 0.0014 * z.real - # AddSol(0.384, -0.04, 0.0, 0.0, 0, 1, -2, -2) + + # AddSol( 0.384, -0.04, 0.0 , 0.0 ,0, 1,-2,-2) z = ex[1][2] * ex[-2][3] * ex[-2][4] DLAM += 0.384 * z.imag DS += -0.04 * z.imag - # AddSol(1.938, 3.6, -0.145, 0.0401, 4, 0, 0, 0) + + # AddSol( +1.938, +3.60,-0.145, +0.0401,4, 0, 0, 0) z = ex[4][1] DLAM += 1.938 * z.imag - DS += 3.6 * z.imag + DS += 3.60 * z.imag GAM1C += -0.145 * z.real SINPI += 0.0401 * z.real - # AddSol(-0.952, -1.58, 0.052, -0.013, 4, 0, 0, -2) + + # AddSol( -0.952, -1.58,+0.052, -0.0130,4, 0, 0,-2) z = ex[4][1] * ex[-2][4] DLAM += -0.952 * z.imag DS += -1.58 * z.imag GAM1C += 0.052 * z.real - SINPI += -0.013 * z.real + SINPI += -0.0130 * z.real - # AddSol(-0.551, -0.94, 0.032, -0.0097, 3, 1, 0, 0) + + # AddSol( -0.551, -0.94,+0.032, -0.0097,3, 1, 0, 0) z = ex[3][1] * ex[1][2] DLAM += -0.551 * z.imag DS += -0.94 * z.imag GAM1C += 0.032 * z.real SINPI += -0.0097 * z.real - # AddSol(-0.482, -0.57, 0.005, -0.0045, 3, 1, 0, -2) + + # AddSol( -0.482, -0.57,+0.005, -0.0045,3, 1, 0,-2) z = ex[3][1] * ex[1][2] * ex[-2][4] DLAM += -0.482 * z.imag DS += -0.57 * z.imag GAM1C += 0.005 * z.real SINPI += -0.0045 * z.real - # AddSol(0.681, 0.96, -0.026, 0.0115, 3, -1, 0, 0) + + # AddSol( 0.681, 0.96,-0.026, 0.0115,3,-1, 0, 0) z = ex[3][1] * ex[-1][2] DLAM += 0.681 * z.imag DS += 0.96 * z.imag GAM1C += -0.026 * z.real SINPI += 0.0115 * z.real - # AddSol(-0.297, -0.27, 0.002, -0.0009, 2, 2, 0, -2) + + # AddSol( -0.297, -0.27, 0.002, -0.0009,2, 2, 0,-2) z = ex[2][1] * ex[2][2] * ex[-2][4] DLAM += -0.297 * z.imag DS += -0.27 * z.imag GAM1C += 0.002 * z.real SINPI += -0.0009 * z.real - # AddSol(0.254, 0.21, -0.003, 0.0, 2, -2, 0, -2) + + # AddSol( 0.254, +0.21,-0.003, 0.0 ,2,-2, 0,-2) z = ex[2][1] * ex[-2][2] * ex[-2][4] DLAM += 0.254 * z.imag DS += 0.21 * z.imag GAM1C += -0.003 * z.real - # AddSol(-0.25, -0.22, 0.004, 0.0014, 1, 3, 0, -2) + + # AddSol( -0.250, -0.22, 0.004, 0.0014,1, 3, 0,-2) z = ex[1][1] * ex[3][2] * ex[-2][4] - DLAM += -0.25 * z.imag + DLAM += -0.250 * z.imag DS += -0.22 * z.imag GAM1C += 0.004 * z.real SINPI += 0.0014 * z.real - # AddSol(-3.996, 0.0, 0.0, 0.0004, 2, 0, 2, 0) + + # AddSol( -3.996, 0.0 , 0.0 , +0.0004,2, 0, 2, 0) z = ex[2][1] * ex[2][3] DLAM += -3.996 * z.imag SINPI += 0.0004 * z.real - # AddSol(0.557, -0.75, 0.0, -0.009, 2, 0, 2, -2) + + # AddSol( 0.557, -0.75, 0.0 , -0.0090,2, 0, 2,-2) z = ex[2][1] * ex[2][3] * ex[-2][4] DLAM += 0.557 * z.imag DS += -0.75 * z.imag - SINPI += -0.009 * z.real + SINPI += -0.0090 * z.real - # AddSol(-0.459, -0.38, 0.0, -0.0053, 2, 0, -2, 2) + + # AddSol( -0.459, -0.38, 0.0 , -0.0053,2, 0,-2, 2) z = ex[2][1] * ex[-2][3] * ex[2][4] DLAM += -0.459 * z.imag DS += -0.38 * z.imag SINPI += -0.0053 * z.real - # AddSol(-1.298, 0.74, 0.0, 0.0004, 2, 0, -2, 0) + + # AddSol( -1.298, 0.74, 0.0 , +0.0004,2, 0,-2, 0) z = ex[2][1] * ex[-2][3] DLAM += -1.298 * z.imag DS += 0.74 * z.imag SINPI += 0.0004 * z.real - # AddSol(0.538, 1.14, 0.0, -0.0141, 2, 0, -2, -2) + + # AddSol( 0.538, 1.14, 0.0 , -0.0141,2, 0,-2,-2) z = ex[2][1] * ex[-2][3] * ex[-2][4] DLAM += 0.538 * z.imag DS += 1.14 * z.imag SINPI += -0.0141 * z.real - # AddSol(0.263, 0.02, 0.0, 0.0, 1, 1, 2, 0) + + # AddSol( 0.263, 0.02, 0.0 , 0.0 ,1, 1, 2, 0) z = ex[1][1] * ex[1][2] * ex[2][3] DLAM += 0.263 * z.imag DS += 0.02 * z.imag - # AddSol(0.426, 0.07, 0.0, -0.0006, 1, 1, -2, -2) + + # AddSol( 0.426, +0.07, 0.0 , -0.0006,1, 1,-2,-2) z = ex[1][1] * ex[1][2] * ex[-2][3] * ex[-2][4] DLAM += 0.426 * z.imag DS += 0.07 * z.imag SINPI += -0.0006 * z.real - # AddSol(-0.304, 0.03, 0.0, 0.0003, 1, -1, 2, 0) + + # AddSol( -0.304, +0.03, 0.0 , +0.0003,1,-1, 2, 0) z = ex[1][1] * ex[-1][2] * ex[2][3] DLAM += -0.304 * z.imag DS += 0.03 * z.imag SINPI += 0.0003 * z.real - # AddSol(-0.372, -0.19, 0.0, -0.0027, 1, -1, -2, 2) + + # AddSol( -0.372, -0.19, 0.0 , -0.0027,1,-1,-2, 2) z = ex[1][1] * ex[-1][2] * ex[-2][3] * ex[2][4] DLAM += -0.372 * z.imag DS += -0.19 * z.imag SINPI += -0.0027 * z.real - # AddSol(0.418, 0.0, 0.0, 0.0, 0, 0, 4, 0) + + # AddSol( +0.418, 0.0 , 0.0 , 0.0 ,0, 0, 4, 0) z = ex[4][3] DLAM += 0.418 * z.imag - # AddSol(-0.33, -0.04, 0.0, 0.0, 3, 0, 2, 0) + + # AddSol( -0.330, -0.04, 0.0 , 0.0 ,3, 0, 2, 0) z = ex[3][1] * ex[2][3] - DLAM += -0.33 * z.imag + DLAM += -0.330 * z.imag DS += -0.04 * z.imag + def ADDN(coeffn, p, q, r, s): return coeffn * (ex[p][1] * ex[q][2] * ex[r][3] * ex[s][4]).imag