1 #include "astro.h"
2
3 Occ o1, o2;
4 Obj2 xo1, xo2;
5
6 void
occult(Obj2 * p1,Obj2 * p2,double)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
setpt(Occ * o,double x)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
pinorm(double a)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
set3pt(Obj2 * p,int i,Occ * o)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 }
191