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 }