ref: 4821c261c4cea895ddcbd7c70bc71103c9a2763b
dir: /sys/src/cmd/astro/satel.c/
#include "astro.h" char* satlst[] = { 0, }; struct { double time; double tilt; double pnni; double psi; double ppi; double d1pp; double peri; double d1per; double e0; double deo; double rdp; double st; double ct; double rot; double semi; } satl; void satels(void) { double ifa[10], t, t1, t2, tinc; char **satp; int flag, f, i, n; satp = satlst; loop: if(*satp == 0) return; f = open(*satp, 0); if(f < 0) { fprint(2, "cannot open %s\n", *satp); satp += 2; goto loop; } satp++; rline(f); tinc = atof(skip(6)); rline(f); rline(f); for(i=0; i<9; i++) ifa[i] = atof(skip(i)); n = ifa[0]; i = ifa[1]; t = dmo[i-1] + ifa[2] - 1.; if(n%4 == 0 && i > 2) t += 1.; for(i=1970; i<n; i++) { t += 365.; if(i%4 == 0) t += 1.; } t = (t * 24. + ifa[3]) * 60. + ifa[4]; satl.time = t * 60.; satl.tilt = ifa[5] * radian; satl.pnni = ifa[6] * radian; satl.psi = ifa[7]; satl.ppi = ifa[8] * radian; rline(f); for(i=0; i<5; i++) ifa[i] = atof(skip(i)); satl.d1pp = ifa[0] * radian; satl.peri = ifa[1]; satl.d1per = ifa[2]; satl.e0 = ifa[3]; satl.deo = 0; satl.rdp = ifa[4]; satl.st = sin(satl.tilt); satl.ct = cos(satl.tilt); satl.rot = pipi / (1440. + satl.psi); satl.semi = satl.rdp * (1. + satl.e0); n = PER*288.; /* 5 min steps */ t = day; for(i=0; i<n; i++) { if(sunel((t-day)/deld) > 0.) goto out; satel(t); if(el > 0) { t1 = t; flag = 0; do { if(el > 30.) flag++; t -= tinc/(24.*3600.); satel(t); } while(el > 0.); t2 = (t - day)/deld; t = t1; do { t += tinc/(24.*3600.); satel(t); if(el > 30.) flag++; } while(el > 0.); if(flag) if((*satp)[0] == '-') event("%s pass at ", (*satp)+1, "", t2, SIGNIF+PTIME+DARK); else event("%s pass at ", *satp, "", t2, PTIME+DARK); } out: t += 5./(24.*60.); } close(f); satp++; goto loop; } void satel(double time) { int i; double amean, an, coc, csl, d, de, enom, eo; double pp, q, rdp, slong, ssl, t, tp; i = 500; el = -1; time = (time-25567.5) * 86400; t = (time - satl.time)/60; if(t < 0) return; /* too early for satelites */ an = floor(t/satl.peri); while(an*satl.peri + an*an*satl.d1per/2. <= t) { an += 1; if(--i == 0) return; } while((tp = an*satl.peri + an*an*satl.d1per/2.) > t) { an -= 1; if(--i == 0) return; } amean = (t-tp)/(satl.peri+(an+.5)*satl.d1per); pp = satl.ppi+(an+amean)*satl.d1pp; amean *= pipi; eo = satl.e0+satl.deo*an; rdp = satl.semi/(1+eo); enom = amean+eo*sin(amean); do { de = (amean-enom+eo*sin(enom))/(1.0-eo*cos(enom)); enom += de; if(--i == 0) return; } while(fabs(de) >= 1.0e-6); q = 3963.35*erad/(rdp*(1-eo*cos(enom))/(1-eo)); d = pp + 2*atan2(sqrt((1+eo)/(1-eo))*sin(enom/2),cos(enom/2)); slong = satl.pnni + t*satl.rot - atan2(satl.ct*sin(d), cos(d)); ssl = satl.st*sin(d); csl = pyth(ssl); if(vis(time, atan2(ssl,csl), slong, q)) { coc = ssl*sin(glat) + csl*cos(glat)*cos(wlong-slong); el = atan2(coc-q, pyth(coc)); el /= radian; } } int vis(double t, double slat, double slong, double q) { double t0, t1, t2, d; d = t/86400 - .005375 + 3653; t0 = 6.238030674 + .01720196977*d; t2 = t0 + .0167253303*sin(t0); do { t1 = (t0 - t2 + .0167259152*sin(t2)) / (1 - .0167259152*cos(t2)); t2 = t2 + t1; } while(fabs(t1) >= 1.e-4); t0 = 2*atan2(1.01686816*sin(t2/2), cos(t2/2)); t0 = 4.926234925 + 8.214985538e-7*d + t0; t1 = 0.91744599 * sin(t0); t0 = atan2(t1, cos(t0)); if(t0 < -pi/2) t0 = t0 + pipi; d = 4.88097876 + 6.30038809*d - t0; t0 = 0.43366079 * t1; t1 = pyth(t0); t2 = t1*cos(slat)*cos(d-slong) - t0*sin(slat); if(t2 > 0.46949322e-2) { if(0.46949322e-2*t2 + 0.999988979*pyth(t2) < q) return 0; } t2 = t1*cos(glat)*cos(d-wlong) - t0*sin(glat); if(t2 < .1) return 0; return 1; }