sys/src/cmd/astro/occ.c

   1 #include "astro.h"
   2
   3 Occ      o1, o2;
   4 Obj2     xo1, xo2;
   5
   6 void
   7 occult(Obj2 *p1, Obj2 *p2, double)
   8 {
   9         int i, i1, N;
  10         double d1, d2, d3, d4;
  11         double x, di, dx, x1;
  12
  13         d3 = 0;
  14         d2 = 0;
  15         occ.t1 = -100;
  16         occ.t2 = -100;
  17         occ.t3 = -100;
  18         occ.t4 = -100;
  19         occ.t5 = -100;
  20         for(i=0; i<=NPTS+1; i++) {
  21                 d1 = d2;
  22                 d2 = d3;
  23                 d3 = dist(&p1->point[i], &p2->point[i]);
  24                 if(i >= 2 && d2 <= d1 && d2 <= d3)
  25                         goto found;
  26         }
  27         return;
  28
  29 found:
  30         N = 2880*PER/NPTS; /* 1 min steps */
  31         i -= 2;
  32         set3pt(p1, i, &o1);
  33         set3pt(p2, i, &o2);
  34
  35         di = i;
  36         x = 0;
  37         dx = 2./N;
  38         for(i=0; i<=N; i++) {
  39                 setpt(&o1, x);
  40                 setpt(&o2, x);
  41                 d1 = d2;
  42                 d2 = d3;
  43                 d3 = dist(&o1.act, &o2.act);
  44                 if(i >= 2 && d2 <= d1 && d2 <= d3)
  45                         goto found1;
  46                 x += dx;
  47         }
  48         fprint(2, "bad 1 \n");
  49         return;
  50
  51 found1:
  52         if(d2 > o1.act.semi2+o2.act.semi2+50)
  53                 return;
  54         di += x - 3*dx;
  55         x = 0;
  56         for(i=0; i<3; i++) {
  57                 setime(day + deld*(di + x));
  58                 (*p1->obj)();
  59                 setobj(&xo1.point[i]);
  60                 (*p2->obj)();
  61                 setobj(&xo2.point[i]);
  62                 x += 2*dx;
  63         }
  64         dx /= 60;
  65         x = 0;
  66         set3pt(&xo1, 0, &o1);
  67         set3pt(&xo2, 0, &o2);
  68         for(i=0; i<=240; i++) {
  69                 setpt(&o1, x);
  70                 setpt(&o2, x);
  71                 d1 = d2;
  72                 d2 = d3;
  73                 d3 = dist(&o1.act, &o2.act);
  74                 if(i >= 2 && d2 <= d1 && d2 <= d3)
  75                         goto found2;
  76                 x += 1./120.;
  77         }
  78         fprint(2, "bad 2 \n");
  79         return;
  80
  81 found2:
  82         if(d2 > o1.act.semi2 + o2.act.semi2)
  83                 return;
  84         i1 = i-1;
  85         x1 = x - 1./120.;
  86         occ.t3 = di + i1 * dx;
  87         occ.e3 = o1.act.el;
  88         d3 = o1.act.semi2 - o2.act.semi2;
  89         if(d3 < 0)
  90                 d3 = -d3;
  91         d4 = o1.act.semi2 + o2.act.semi2;
  92         for(i=i1,x=x1;; i++) {
  93                 setpt(&o1, x);
  94                 setpt(&o2, x);
  95                 d1 = d2;
  96                 d2 = dist(&o1.act, &o2.act);
  97                 if(i != i1) {
  98                         if(d1 <= d3 && d2 > d3) {
  99                                 occ.t4 = di + (i-.5) * dx;
 100                                 occ.e4 = o1.act.el;
 101                         }
 102                         if(d2 > d4) {
 103                                 if(d1 <= d4) {
 104                                         occ.t5 = di + (i-.5) * dx;
 105                                         occ.e5 = o1.act.el;
 106                                 }
 107                                 break;
 108                         }
 109                 }
 110                 x += 1./120.;
 111         }
 112         for(i=i1,x=x1;; i--) {
 113                 setpt(&o1, x);
 114                 setpt(&o2, x);
 115                 d1 = d2;
 116                 d2 = dist(&o1.act, &o2.act);
 117                 if(i != i1) {
 118                         if(d1 <= d3 && d2 > d3) {
 119                                 occ.t2 = di + (i+.5) * dx;
 120                                 occ.e2 = o1.act.el;
 121                         }
 122                         if(d2 > d4) {
 123                                 if(d1 <= d4) {
 124                                         occ.t1 = di + (i+.5) * dx;
 125                                         occ.e1 = o1.act.el;
 126                                 }
 127                                 break;
 128                         }
 129                 }
 130                 x -= 1./120.;
 131         }
 132 }
 133
 134 void
 135 setpt(Occ *o, double x)
 136 {
 137         double y;
 138
 139         y = x * (x-1);
 140         o->act.ra = o->del0.ra +
 141                 x*o->del1.ra + y*o->del2.ra;
 142         o->act.decl2 = o->del0.decl2 +
 143                 x*o->del1.decl2 + y*o->del2.decl2;
 144         o->act.semi2 = o->del0.semi2 +
 145                 x*o->del1.semi2 + y*o->del2.semi2;
 146         o->act.el = o->del0.el +
 147                 x*o->del1.el + y*o->del2.el;
 148 }
 149
 150
 151 double
 152 pinorm(double a)
 153 {
 154
 155         while(a < -pi)
 156                 a += pipi;
 157         while(a > pi)
 158                 a -= pipi;
 159         return a;
 160 }
 161
 162 void
 163 set3pt(Obj2 *p, int i, Occ *o)
 164 {
 165         Obj1 *p1, *p2, *p3;
 166         double a;
 167
 168         p1 = p->point+i;
 169         p2 = p1+1;
 170         p3 = p2+1;
 171
 172         o->del0.ra = p1->ra;
 173         o->del0.decl2 = p1->decl2;
 174         o->del0.semi2 = p1->semi2;
 175         o->del0.el = p1->el;
 176
 177         a = p2->ra - p1->ra;
 178         o->del1.ra = pinorm(a);
 179         a = p2->decl2 - p1->decl2;
 180         o->del1.decl2 = pinorm(a);
 181         o->del1.semi2 = p2->semi2 - p1->semi2;
 182         o->del1.el = p2->el - p1->el;
 183
 184         a = p1->ra + p3->ra - 2*p2->ra;
 185         o->del2.ra = pinorm(a)/2;
 186         a = p1->decl2 + p3->decl2 - 2*p2->decl2;
 187         o->del2.decl2 = pinorm(a)/2;
 188         o->del2.semi2 = (p1->semi2 + p3->semi2 - 2*p2->semi2) / 2;
 189         o->del2.el = (p1->el + p3->el - 2*p2->el) / 2;
 190 }