#include #include #include #include #include #include #include "Aurinkokunta.h" #include "Cities.h" #include "GetCity.h" #ifndef PI #define PI 3.14159265358979323846 #endif #define LE 6 static struct BODY defaultBody = { {0}, {0}, VAL_NA, VAL_NA, VAL_NA, VAL_NA, VAL_NA, VAL_NA, VAL_NA, VAL_NA, VAL_NA, VAL_NA, VAL_NA, VAL_NA, "default" }; static int pv=0, kk=0, vs=0, jp=0, zq=0, p1=0, xx=0, tp=0; static double y1=0, m1=0, a=0, jd=0, t=0, a0=0, m=0, eps=0, e2=0, qq=0, e=0, v=0, r=0, l=0, o=0, u=0, b=0, ls=0, in=0, oa=0, ra=0, lm=0, de=0, be=0, ek=0, eg=0, al=0, dc=0, rr=0, dd=0, ha=0, tt=0, gt=0, gl=0, gn=0, hh=0, ko=0, j1=0, j2=0, j3=0, j4=0, j5=0, j6=0, j7=0, j8=0, j9=0, a1=0, a2=0, a3=0, a4=0, a5, a9=0; static double el[LE][6][4] = { { {0, 1.22173048e-7, 2608.78753, 1.78511194}, {0, 2.24449342e-5, 1021.32835, 3.71062619}, {-5.75958653e-8, -2.61799388e-6, 628.301945, 6.25658352}, {0, 3.15904595e-6, 334.053484, 5.57666084}, {0, -1.26012772e-5, 52.9653675, 3.93272126}, {0, -8.76155284e-6, 21.3200951, 3.06246326} },{ {0, -3e-8, 2.046e-5, .20561421}, {0, 9.1e-8, 4.774e-5, .00682069}, {0, -1.26e-7, -4.18e-5, .01675104}, {0, 7.7e-8, 9.2064e-5, .0933129}, { -1.7e-9, -4.676e-7, 1.6418e-4, .04833475}, {7.4e-10, -7.28e-7, -3.455e-4, .05589232} },{ {0, 0, 0, .3870986}, {0, 0, 0, .7233316}, {0, 0, 0, 1.0000002}, {0, 0, 0, 1.5236883}, {0, 0, 0, 5.202561}, {0, 0, 0, 9.554747} },{ {0, 5.25518638e-6, 2608.81468, 3.10981157}, {0, 5.4052847e-6, 1021.35292, 5.98241364}, {0, 5.27962099e-6, 628.331951, 4.88162797}, {0, 5.42273799e-6, 334.085611, 5.12668362}, {-2.87979327e-8, 5.841617e-6, 52.9934667, 4.15474332}, {-1.01229097e-7, 5.66359342e-6, 21.3542759, 4.65242605} },{ {0, 3.03512757e-6, .0206857872, .822851951}, {0, 7.15584993e-6, .0157053453, 1.32260435}, {8.77900614e-9, -2.86233997e-8, -2.27110969e-4, .409319747}, {-9.30260491e-8, -2.44346095e-8, .0134563437, .851484043}, {-1.48527519e-7, 6.14739869e-6, .0176370757, 1.73561499}, {-9.26769833e-8, -2.65604206e-6, .0152401295, 1.96856409} },{ {0, -3.19395253e-7, 3.24770867e-5, .122223331}, {0, -1.74532925e-8, 1.75545216e-5, .0592300345}, {0, 0, 0, 0}, {0, 2.19911486e-7, -1.17809724e-5, .0322944031}, {0, 6.80678408e-8, -9.94156995e-5, .0228417522}, {6.98131701e-10, -2.70351501e-7, -6.8397708e-5, .0435026632} } }; static char pl[LE][10]= {"Mercury", "Venus", "Sun", "Mars", "Jupiter", "Saturn"}; static const double cc = PI/180; static int fsgn(double x) { return (x ? ((x < 0) ? -1 : 1) : 0); } double FNM(double _x) { return _x-2*PI*((int)(_x/(PI*2.0))); } double FNP(double _t) { return ((el[p1-1][jp-1][1-1]*_t+el[p1-1][jp-1][2-1])*_t+ el[p1-1][jp-1][3-1])*_t+el[p1-1][jp-1][4-1]; } void F5200(void) { a0 = FNP(t); if (p1<2 || p1>3) a0=FNM(a0); } void F5300(void) { e2=m; qq=0; do { e = e2; e2 = e + (m + eps * sin (e) - e) / (1 - eps * cos (e)); } while (qq++ < 100 && fabs (e2 - e) > 1e-10); } void F5400(void) { a9 = sqrt (a1 * a1 + a2 * a2); a9 = 1 / a9; a1 = a9 * a1 ; a2 = a9 * a2; if (fabs(a2) < 1e-9) { a0 = .5*PI; if (a1<0) { a0+=PI; } else { if (a2<0) a0+=PI; } if (a0<0) a0+=2*PI; } else { a0 = atan(a1/a2); if (a2<0) a0+=PI; if (a0<0) a0+=2*PI; } } void F5500(void) { a4 = fsgn(a0) ; a0 = fabs(a0 / cc); a1 = (int)(a0); a0 = 60 * (a0 - a1); a2 = (int)(a0); a1 = a4 * (a1 + a2 / 100); } void F5800(void) { p1=1; F5200(); m=a0; p1=2; F5200(); eps=a0; F5300(); a1 = sqrt ((1 + eps) / (1 - eps)) * sin (.5 * e); a2 = cos(.5 * e); F5400(); v=2*a0; } void F6700(void) { a1 = (int)(gt); gt = 60 * (gt - a1); a2 = (int)(gt); a1 = a1 + a2 / 100; } static char months[2][12] = { {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}, {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31} }; struct BODY *Aurinkokunta(enum BODYNAME bodyName, const struct tm *day, const struct CITY *city) { struct BODY *body; if (!city) city = GetDefaultCity(); vs = day->tm_year+1900; kk = day->tm_mon+1; pv = day->tm_mday; if (vs < 1950 || vs > 2049) { printf("Illegal year %d (range 1950..2049)\n", vs); return NULL; } if (kk < 1 || kk > 12) { printf("Illegal month %d (range 1..12)\n", kk); return NULL; } if (pv < 1 || pv > months[vs&3 ? 0:1][kk-1]) { printf("Illegal day %d (range 1..%d)\n", pv, months[vs&3 ? 0:1][kk-1]); return NULL; } if (!(body = malloc(sizeof(struct BODY)))) { printf("Unable to allocate storage for body\n"); return NULL; } switch(bodyName) { case BN_SUN: jp = 3; break; case BN_MERCURY: jp = 1; break; case BN_VENUS: jp = 2; break; case BN_MARS: jp = 4; break; case BN_JUPITER: jp = 5; break; case BN_SATURN: jp = 6; break; default: printf("Internal error #2\n"); return NULL; break; } *body = defaultBody; body->day = *day; body->city = *city; strcpy(body->name, pl[jp-1]); zq = 0; y1 = vs; m1 = kk + 1; if (!(kk > 2)) { y1 = vs-1; m1=kk+13; } a = (int)(.01*y1); jd = (int)(365.25 * y1) + (int)(30.6001 * m1) + pv - a + (int)(.25 * a) + 1720996.5; t = (jd - 2415020.0) / 36525.0; F5800(); p1 = 3; F5200(); a=a0; r = a * (1 - eps * cos (e)); p1 = 4 ; F5200() ; l = a0; p1 = 5 ; F5200() ; o = a0; u = l + v - m - o; while (u<0) u+= 2*PI; u = FNM(u); p1 = 6 ; F5200(); in = a0; a1 = cos (in) * sin (u); a2 = cos (u); F5400(); ls = FNM(a0 + o); a1 = sin (u) * sin (in); a2 = sqrt (1 - a1 * a1); F5400(); b = a0 ; if (b > PI) b = b - 2 * PI; j1 = jd + 0.5; j2 = (int)(j1); j3 = j1 - (int)(j1); j4 = j2 + 1 + (int)((j2 - 1867216.25) / 36524.25) - (int)(((j2 - 1867216.25) / 36524.25) / 4); j5 = j4 + 1524; j6 = (int)((j5 - 122.1) / 365.25) ; j7 = (int)(365.25 * j6); j8 = (int)((j5 - j7) / 30.6001); j9 = j5 - j7 - (int)(30.6001 * j8) + j3; kk = j8 - ((j8 < 13.5) ? 1 : 13); pv = (int)(j9); a0 = ls; F5500(); a0 = b; F5500(); xx = jp; jp = 3; F5800(); p1 = 4; F5200(); l = a0; oa = l + v - m; if (oa < 0) oa = oa + 2 * PI; oa = FNM(oa); ra = 1.0000002 * (1 - eps * cos (e)); a0 = r * cos (b); a1 = a0 * sin (ls - oa); a3 = a1; a2 = a0 * cos (ls - oa) + ra; a4 = a2; F5400(); lm = FNM(a0 + oa); a5 = r * sin (b); de = sqrt(a3*a3 + a4*a4 + a5*a5); a1 = a5; a2 = de; F5400(); be = a0; if (be > PI) be = be - 2 * PI; p1 = 5; F5200() ; ek = a0; a2 = cos (be) * cos (lm - oa); a1 = sqrt (1 - a2 * a2); F5400(); eg = a0; if (eg > PI) eg = eg - 2 * PI; eg = eg * 180 / PI; eg = (int)(eg); a1 = sin (lm) * cos (ek) - tan (be) * sin (ek); a2 = cos (lm); F5400(); al = a0; a1 = sin (be) * cos (ek) + cos (be) * sin (ek) * sin (lm); a2 = sqrt (1 - a1 * a1); F5400(); dc = a0; if (dc > PI) dc = dc - 2 * PI; a0 = al / 15; F5500(); rr = a1; a0 = dc; F5500(); dd = a1; if (xx == 3) goto c7100; switch(xx) { case 1: ha = 6.68 / de; break; case 2: ha = 16.82 / de; break; case 4: ha = 9.36 / de; break; case 5: ha = 196.94 / de; break; case 6: ha = 166.66 / de; break; default: printf("Internal error #1\n"); exit(10); } ha = (int)(10 * ha) / 10.0; tt = 6.6460656 + 2400.051262 * t + 0.00002581 * t * t; if (tt >= 24) tt = tt - ((int)(tt / 24)) * 24; hh = ( -0.00989 - sin (city->latitude * cc) * sin (dc)) / ( cos (city->latitude * cc) * cos (dc)); if (fabs (hh) > 1) { body->timeOfRise = VAL_LINE; goto ret6400; } hh = -atan (hh / sqrt( - hh * hh + 1)) + PI / 2 ; hh = hh * 180 / PI / 15; gn = 24 - hh + al / 15 / cc; gl = hh + al / 15 / cc; gt = (gn - city->longitude/15.0 - tt) / 1.002737908 - city->timeZone/3600.0; while (gt >= 24) gt -= 24; while (gt < 0) gt += 24; F6700(); body->timeOfRise = a1; gt = (gl - city->longitude/15.0 - tt) / 1.002737908 - city->timeZone/3600.0; while (gt >= 24) gt -= 24; while (gt < 0) gt += 24; F6700(); body->timeOfSet = a1; ret6400: ko = ( (int) (10 * (dd + 90.0-city->latitude))) / 10.0; body->rectascension = rr; body->declination = dd; body->elongation = eg; body->distance = de; body->diameter = ha; body->maxHeight = ko; goto out; c7100: for (tp=1; tp<=6; tp++) { tt = 6.6460656 + 2400.051262 * t + 0.00002581 * t * t; if (tt >= 24) tt = tt - ((int)(tt / 24)) * 24; if (tp == 1 || tp == 6) hh = ( -18 * PI / 180 - sin (city->latitude * cc) /* -18° */ * sin (dc)) / ( cos (city->latitude * cc) * cos (dc)); if (tp == 2 || tp == 5) hh = ( -6 * PI / 180 - sin (city->latitude * cc) /* -6° */ * sin (dc)) / ( cos (city->latitude * cc) * cos (dc)); else if (tp == 3 || tp == 4) hh = ( - 0.01454 - sin (city->latitude * cc) * sin (dc)) / ( cos (city->latitude * cc) * cos (dc)); if (fabs (hh) > 1) goto c7310; hh = -atan (hh / sqrt ( - hh * hh + 1)) + PI / 2 ; hh = hh * 180 / PI / 15; if (tp == 4 || tp == 5 || tp == 6) goto c7180; gn = 24 - hh + al / 15 / cc; gt = (gn - city->longitude/15.0 - tt) / 1.002737908 - city->timeZone/3600.0; goto c7200; c7180: gl = hh + al / 15 / cc; gt = (gl - city->longitude/15.0 - tt) / 1.002737908 - city->timeZone/3600.0; c7200: while (gt >= 24) gt -= 24; while (gt < 0) gt += 24; F6700(); switch(tp) { case 1: body->timeOf18DegRise = a1; break; case 2: body->timeOf6DegRise = a1; break; case 3: body->timeOfRise = a1; break; case 4: body->timeOfSet = a1; break; case 5: body->timeOf6DegSet = a1; break; case 6: body->timeOf18DegSet = a1; break; } goto c7399; c7310: switch(tp) { case 1: body->timeOf18DegRise = VAL_LINE; break; case 2: body->timeOf6DegRise = VAL_LINE; break; case 3: body->timeOfRise = VAL_LINE; break; case 4: body->timeOfSet = VAL_LINE; break; case 5: body->timeOf6DegSet = VAL_LINE; break; case 6: body->timeOf18DegSet = VAL_LINE; break; } c7399: ; /* empty statement */ } ko = ((int)(10 * (dd + 90.0-city->latitude))) / 10.0; body->rectascension = rr; body->declination = dd; body->maxHeight = ko; out: return body; }