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

You are viewing proxied material from sdf.org. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.