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