30 double ifa[10], t, t1, t2, tinc;
41 fprint(2, "cannot open %s\n", *satp);
51 ifa[i] = atof(skip(i));
54 t = dmo[i-1] + ifa[2] - 1.;
57 for(i=1970; i<n; i++) {
62 t = (t * 24. + ifa[3]) * 60. + ifa[4];
64 satl.tilt = ifa[5] * radian;
65 satl.pnni = ifa[6] * radian;
67 satl.ppi = ifa[8] * radian;
70 ifa[i] = atof(skip(i));
71 satl.d1pp = ifa[0] * radian;
78 satl.st = sin(satl.tilt);
79 satl.ct = cos(satl.tilt);
80 satl.rot = pipi / (1440. + satl.psi);
81 satl.semi = satl.rdp * (1. + satl.e0);
83 n = PER*288.; /* 5 min steps */
86 if(sunel((t-day)/deld) > 0.)
95 t -= tinc/(24.*3600.);
101 t += tinc/(24.*3600.);
107 if((*satp)[0] == '-')
108 event("%s pass at ", (*satp)+1, "",
109 t2, SIGNIF+PTIME+DARK); else
110 event("%s pass at ", *satp, "",
125 double amean, an, coc, csl, d, de, enom, eo;
126 double pp, q, rdp, slong, ssl, t, tp;
130 time = (time-25567.5) * 86400;
131 t = (time - satl.time)/60;
133 return; /* too early for satelites */
134 an = floor(t/satl.peri);
135 while(an*satl.peri + an*an*satl.d1per/2. <= t) {
140 while((tp = an*satl.peri + an*an*satl.d1per/2.) > t) {
145 amean = (t-tp)/(satl.peri+(an+.5)*satl.d1per);
146 pp = satl.ppi+(an+amean)*satl.d1pp;
148 eo = satl.e0+satl.deo*an;
149 rdp = satl.semi/(1+eo);
150 enom = amean+eo*sin(amean);
152 de = (amean-enom+eo*sin(enom))/(1.0-eo*cos(enom));
156 } while(fabs(de) >= 1.0e-6);
157 q = 3963.35*erad/(rdp*(1-eo*cos(enom))/(1-eo));
158 d = pp + 2*atan2(sqrt((1+eo)/(1-eo))*sin(enom/2),cos(enom/2));
159 slong = satl.pnni + t*satl.rot -
160 atan2(satl.ct*sin(d), cos(d));
161 ssl = satl.st*sin(d);
163 if(vis(time, atan2(ssl,csl), slong, q)) {
164 coc = ssl*sin(glat) + csl*cos(glat)*cos(wlong-slong);
165 el = atan2(coc-q, pyth(coc));
171 vis(double t, double slat, double slong, double q)
173 double t0, t1, t2, d;
175 d = t/86400 - .005375 + 3653;
176 t0 = 6.238030674 + .01720196977*d;
177 t2 = t0 + .0167253303*sin(t0);
179 t1 = (t0 - t2 + .0167259152*sin(t2)) /
180 (1 - .0167259152*cos(t2));
182 } while(fabs(t1) >= 1.e-4);
183 t0 = 2*atan2(1.01686816*sin(t2/2), cos(t2/2));
184 t0 = 4.926234925 + 8.214985538e-7*d + t0;
185 t1 = 0.91744599 * sin(t0);
186 t0 = atan2(t1, cos(t0));
189 d = 4.88097876 + 6.30038809*d - t0;
190 t0 = 0.43366079 * t1;
192 t2 = t1*cos(slat)*cos(d-slong) - t0*sin(slat);
193 if(t2 > 0.46949322e-2) {
194 if(0.46949322e-2*t2 + 0.999988979*pyth(t2) < q)
197 t2 = t1*cos(glat)*cos(d-wlong) - t0*sin(glat);