Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 
}