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