sys/src/cmd/astro/satel.c

   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 }