satel.c (3832B)
1 #include "astro.h" 2 3 char* satlst[] = 4 { 5 0 6 }; 7 8 struct 9 { 10 double time; 11 double tilt; 12 double pnni; 13 double psi; 14 double ppi; 15 double d1pp; 16 double peri; 17 double d1per; 18 double e0; 19 double deo; 20 double rdp; 21 double st; 22 double ct; 23 double rot; 24 double semi; 25 } satl; 26 27 void 28 satels(void) 29 { 30 double ifa[10], t, t1, t2, tinc; 31 char **satp; 32 int flag, f, i, n; 33 34 satp = satlst; 35 36 loop: 37 if(*satp == 0) 38 return; 39 f = open(*satp, 0); 40 if(f < 0) { 41 fprint(2, "cannot open %s\n", *satp); 42 satp += 2; 43 goto loop; 44 } 45 satp++; 46 rline(f); 47 tinc = atof(skip(6)); 48 rline(f); 49 rline(f); 50 for(i=0; i<9; i++) 51 ifa[i] = atof(skip(i)); 52 n = ifa[0]; 53 i = ifa[1]; 54 t = dmo[i-1] + ifa[2] - 1.; 55 if(n%4 == 0 && i > 2) 56 t += 1.; 57 for(i=1970; i<n; i++) { 58 t += 365.; 59 if(i%4 == 0) 60 t += 1.; 61 } 62 t = (t * 24. + ifa[3]) * 60. + ifa[4]; 63 satl.time = t * 60.; 64 satl.tilt = ifa[5] * radian; 65 satl.pnni = ifa[6] * radian; 66 satl.psi = ifa[7]; 67 satl.ppi = ifa[8] * radian; 68 rline(f); 69 for(i=0; i<5; i++) 70 ifa[i] = atof(skip(i)); 71 satl.d1pp = ifa[0] * radian; 72 satl.peri = ifa[1]; 73 satl.d1per = ifa[2]; 74 satl.e0 = ifa[3]; 75 satl.deo = 0; 76 satl.rdp = ifa[4]; 77 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); 82 83 n = PER*288.; /* 5 min steps */ 84 t = day; 85 for(i=0; i<n; i++) { 86 if(sunel((t-day)/deld) > 0.) 87 goto out; 88 satel(t); 89 if(el > 0) { 90 t1 = t; 91 flag = 0; 92 do { 93 if(el > 30.) 94 flag++; 95 t -= tinc/(24.*3600.); 96 satel(t); 97 } while(el > 0.); 98 t2 = (t - day)/deld; 99 t = t1; 100 do { 101 t += tinc/(24.*3600.); 102 satel(t); 103 if(el > 30.) 104 flag++; 105 } while(el > 0.); 106 if(flag) 107 if((*satp)[0] == '-') 108 event("%s pass at ", (*satp)+1, "", 109 t2, SIGNIF+PTIME+DARK); else 110 event("%s pass at ", *satp, "", 111 t2, PTIME+DARK); 112 } 113 out: 114 t += 5./(24.*60.); 115 } 116 close(f); 117 satp++; 118 goto loop; 119 } 120 121 void 122 satel(double time) 123 { 124 int i; 125 double amean, an, coc, csl, d, de, enom, eo; 126 double pp, q, rdp, slong, ssl, t, tp; 127 128 i = 500; 129 el = -1; 130 time = (time-25567.5) * 86400; 131 t = (time - satl.time)/60; 132 if(t < 0) 133 return; /* too early for satelites */ 134 an = floor(t/satl.peri); 135 while(an*satl.peri + an*an*satl.d1per/2. <= t) { 136 an += 1; 137 if(--i == 0) 138 return; 139 } 140 while((tp = an*satl.peri + an*an*satl.d1per/2.) > t) { 141 an -= 1; 142 if(--i == 0) 143 return; 144 } 145 amean = (t-tp)/(satl.peri+(an+.5)*satl.d1per); 146 pp = satl.ppi+(an+amean)*satl.d1pp; 147 amean *= pipi; 148 eo = satl.e0+satl.deo*an; 149 rdp = satl.semi/(1+eo); 150 enom = amean+eo*sin(amean); 151 do { 152 de = (amean-enom+eo*sin(enom))/(1.0-eo*cos(enom)); 153 enom += de; 154 if(--i == 0) 155 return; 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); 162 csl = pyth(ssl); 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)); 166 el /= radian; 167 } 168 } 169 170 int 171 vis(double t, double slat, double slong, double q) 172 { 173 double t0, t1, t2, d; 174 175 d = t/86400 - .005375 + 3653; 176 t0 = 6.238030674 + .01720196977*d; 177 t2 = t0 + .0167253303*sin(t0); 178 do { 179 t1 = (t0 - t2 + .0167259152*sin(t2)) / 180 (1 - .0167259152*cos(t2)); 181 t2 = t2 + t1; 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)); 187 if(t0 < -pi/2) 188 t0 = t0 + pipi; 189 d = 4.88097876 + 6.30038809*d - t0; 190 t0 = 0.43366079 * t1; 191 t1 = pyth(t0); 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) 195 return 0; 196 } 197 t2 = t1*cos(glat)*cos(d-wlong) - t0*sin(glat); 198 if(t2 < .1) 199 return 0; 200 return 1; 201 }