tocc.c - plan9port - [fork] Plan 9 from user space
HTML git clone git://src.adamsgaard.dk/plan9port
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tocc.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 }