123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201 |
- #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;
- }
|