Blame

Date:
Message:
Astro with some minor changes to placate Unix.
Actions:
History | Blob | Raw File
  1 cd5bae78 2004-04-21 devnull #include "astro.h"

  2 cd5bae78 2004-04-21 devnull 

  3 cd5bae78 2004-04-21 devnull Occ	 o1, o2;

  4 cd5bae78 2004-04-21 devnull Obj2	 xo1, xo2;

  5 cd5bae78 2004-04-21 devnull 

  6 cd5bae78 2004-04-21 devnull void

  7 cd5bae78 2004-04-21 devnull occult(Obj2 *p1, Obj2 *p2, double d)

  8 cd5bae78 2004-04-21 devnull {

  9 cd5bae78 2004-04-21 devnull 	int i, i1, N;

 10 cd5bae78 2004-04-21 devnull 	double d1, d2, d3, d4;

 11 cd5bae78 2004-04-21 devnull 	double x, di, dx, x1;

 12 cd5bae78 2004-04-21 devnull 

 13 cd5bae78 2004-04-21 devnull 	USED(d);

 14 cd5bae78 2004-04-21 devnull 

 15 cd5bae78 2004-04-21 devnull 	d3 = 0;

 16 cd5bae78 2004-04-21 devnull 	d2 = 0;

 17 cd5bae78 2004-04-21 devnull 	occ.t1 = -100;

 18 cd5bae78 2004-04-21 devnull 	occ.t2 = -100;

 19 cd5bae78 2004-04-21 devnull 	occ.t3 = -100;

 20 cd5bae78 2004-04-21 devnull 	occ.t4 = -100;

 21 cd5bae78 2004-04-21 devnull 	occ.t5 = -100;

 22 cd5bae78 2004-04-21 devnull 	for(i=0; i<=NPTS+1; i++) {

 23 cd5bae78 2004-04-21 devnull 		d1 = d2;

 24 cd5bae78 2004-04-21 devnull 		d2 = d3;

 25 cd5bae78 2004-04-21 devnull 		d3 = dist(&p1->point[i], &p2->point[i]);

 26 cd5bae78 2004-04-21 devnull 		if(i >= 2 && d2 <= d1 && d2 <= d3)

 27 cd5bae78 2004-04-21 devnull 			goto found;

 28 cd5bae78 2004-04-21 devnull 	}

 29 cd5bae78 2004-04-21 devnull 	return;

 30 cd5bae78 2004-04-21 devnull 

 31 cd5bae78 2004-04-21 devnull found:

 32 cd5bae78 2004-04-21 devnull 	N = 2880*PER/NPTS; /* 1 min steps */

 33 cd5bae78 2004-04-21 devnull 	i -= 2;

 34 cd5bae78 2004-04-21 devnull 	set3pt(p1, i, &o1);

 35 cd5bae78 2004-04-21 devnull 	set3pt(p2, i, &o2);

 36 cd5bae78 2004-04-21 devnull 

 37 cd5bae78 2004-04-21 devnull 	di = i;

 38 cd5bae78 2004-04-21 devnull 	x = 0;

 39 cd5bae78 2004-04-21 devnull 	dx = 2./N;

 40 cd5bae78 2004-04-21 devnull 	for(i=0; i<=N; i++) {

 41 cd5bae78 2004-04-21 devnull 		setpt(&o1, x);

 42 cd5bae78 2004-04-21 devnull 		setpt(&o2, x);

 43 cd5bae78 2004-04-21 devnull 		d1 = d2;

 44 cd5bae78 2004-04-21 devnull 		d2 = d3;

 45 cd5bae78 2004-04-21 devnull 		d3 = dist(&o1.act, &o2.act);

 46 cd5bae78 2004-04-21 devnull 		if(i >= 2 && d2 <= d1 && d2 <= d3)

 47 cd5bae78 2004-04-21 devnull 			goto found1;

 48 cd5bae78 2004-04-21 devnull 		x += dx;

 49 cd5bae78 2004-04-21 devnull 	}

 50 cd5bae78 2004-04-21 devnull 	fprint(2, "bad 1 \n");

 51 cd5bae78 2004-04-21 devnull 	return;

 52 cd5bae78 2004-04-21 devnull 

 53 cd5bae78 2004-04-21 devnull found1:

 54 cd5bae78 2004-04-21 devnull 	if(d2 > o1.act.semi2+o2.act.semi2+50)

 55 cd5bae78 2004-04-21 devnull 		return;

 56 cd5bae78 2004-04-21 devnull 	di += x - 3*dx;

 57 cd5bae78 2004-04-21 devnull 	x = 0;

 58 cd5bae78 2004-04-21 devnull 	for(i=0; i<3; i++) {

 59 cd5bae78 2004-04-21 devnull 		setime(day + deld*(di + x));

 60 cd5bae78 2004-04-21 devnull 		(*p1->obj)();

 61 cd5bae78 2004-04-21 devnull 		setobj(&xo1.point[i]);

 62 cd5bae78 2004-04-21 devnull 		(*p2->obj)();

 63 cd5bae78 2004-04-21 devnull 		setobj(&xo2.point[i]);

 64 cd5bae78 2004-04-21 devnull 		x += 2*dx;

 65 cd5bae78 2004-04-21 devnull 	}

 66 cd5bae78 2004-04-21 devnull 	dx /= 60;

 67 cd5bae78 2004-04-21 devnull 	x = 0;

 68 cd5bae78 2004-04-21 devnull 	set3pt(&xo1, 0, &o1);

 69 cd5bae78 2004-04-21 devnull 	set3pt(&xo2, 0, &o2);

 70 cd5bae78 2004-04-21 devnull 	for(i=0; i<=240; i++) {

 71 cd5bae78 2004-04-21 devnull 		setpt(&o1, x);

 72 cd5bae78 2004-04-21 devnull 		setpt(&o2, x);

 73 cd5bae78 2004-04-21 devnull 		d1 = d2;

 74 cd5bae78 2004-04-21 devnull 		d2 = d3;

 75 cd5bae78 2004-04-21 devnull 		d3 = dist(&o1.act, &o2.act);

 76 cd5bae78 2004-04-21 devnull 		if(i >= 2 && d2 <= d1 && d2 <= d3)

 77 cd5bae78 2004-04-21 devnull 			goto found2;

 78 cd5bae78 2004-04-21 devnull 		x += 1./120.;

 79 cd5bae78 2004-04-21 devnull 	}

 80 cd5bae78 2004-04-21 devnull 	fprint(2, "bad 2 \n");

 81 cd5bae78 2004-04-21 devnull 	return;

 82 cd5bae78 2004-04-21 devnull 

 83 cd5bae78 2004-04-21 devnull found2:

 84 cd5bae78 2004-04-21 devnull 	if(d2 > o1.act.semi2 + o2.act.semi2)

 85 cd5bae78 2004-04-21 devnull 		return;

 86 cd5bae78 2004-04-21 devnull 	i1 = i-1;

 87 cd5bae78 2004-04-21 devnull 	x1 = x - 1./120.;

 88 cd5bae78 2004-04-21 devnull 	occ.t3 = di + i1 * dx;

 89 cd5bae78 2004-04-21 devnull 	occ.e3 = o1.act.el;

 90 cd5bae78 2004-04-21 devnull 	d3 = o1.act.semi2 - o2.act.semi2;

 91 cd5bae78 2004-04-21 devnull 	if(d3 < 0)

 92 cd5bae78 2004-04-21 devnull 		d3 = -d3;

 93 cd5bae78 2004-04-21 devnull 	d4 = o1.act.semi2 + o2.act.semi2;

 94 cd5bae78 2004-04-21 devnull 	for(i=i1,x=x1;; i++) {

 95 cd5bae78 2004-04-21 devnull 		setpt(&o1, x);

 96 cd5bae78 2004-04-21 devnull 		setpt(&o2, x);

 97 cd5bae78 2004-04-21 devnull 		d1 = d2;

 98 cd5bae78 2004-04-21 devnull 		d2 = dist(&o1.act, &o2.act);

 99 cd5bae78 2004-04-21 devnull 		if(i != i1) {

100 cd5bae78 2004-04-21 devnull 			if(d1 <= d3 && d2 > d3) {

101 cd5bae78 2004-04-21 devnull 				occ.t4 = di + (i-.5) * dx;

102 cd5bae78 2004-04-21 devnull 				occ.e4 = o1.act.el;

103 cd5bae78 2004-04-21 devnull 			}

104 cd5bae78 2004-04-21 devnull 			if(d2 > d4) {

105 cd5bae78 2004-04-21 devnull 				if(d1 <= d4) {

106 cd5bae78 2004-04-21 devnull 					occ.t5 = di + (i-.5) * dx;

107 cd5bae78 2004-04-21 devnull 					occ.e5 = o1.act.el;

108 cd5bae78 2004-04-21 devnull 				}

109 cd5bae78 2004-04-21 devnull 				break;

110 cd5bae78 2004-04-21 devnull 			}

111 cd5bae78 2004-04-21 devnull 		}

112 cd5bae78 2004-04-21 devnull 		x += 1./120.;

113 cd5bae78 2004-04-21 devnull 	}

114 cd5bae78 2004-04-21 devnull 	for(i=i1,x=x1;; i--) {

115 cd5bae78 2004-04-21 devnull 		setpt(&o1, x);

116 cd5bae78 2004-04-21 devnull 		setpt(&o2, x);

117 cd5bae78 2004-04-21 devnull 		d1 = d2;

118 cd5bae78 2004-04-21 devnull 		d2 = dist(&o1.act, &o2.act);

119 cd5bae78 2004-04-21 devnull 		if(i != i1) {

120 cd5bae78 2004-04-21 devnull 			if(d1 <= d3 && d2 > d3) {

121 cd5bae78 2004-04-21 devnull 				occ.t2 = di + (i+.5) * dx;

122 cd5bae78 2004-04-21 devnull 				occ.e2 = o1.act.el;

123 cd5bae78 2004-04-21 devnull 			}

124 cd5bae78 2004-04-21 devnull 			if(d2 > d4) {

125 cd5bae78 2004-04-21 devnull 				if(d1 <= d4) {

126 cd5bae78 2004-04-21 devnull 					occ.t1 = di + (i+.5) * dx;

127 cd5bae78 2004-04-21 devnull 					occ.e1 = o1.act.el;

128 cd5bae78 2004-04-21 devnull 				}

129 cd5bae78 2004-04-21 devnull 				break;

130 cd5bae78 2004-04-21 devnull 			}

131 cd5bae78 2004-04-21 devnull 		}

132 cd5bae78 2004-04-21 devnull 		x -= 1./120.;

133 cd5bae78 2004-04-21 devnull 	}

134 cd5bae78 2004-04-21 devnull }

135 cd5bae78 2004-04-21 devnull 

136 cd5bae78 2004-04-21 devnull void

137 cd5bae78 2004-04-21 devnull setpt(Occ *o, double x)

138 cd5bae78 2004-04-21 devnull {

139 cd5bae78 2004-04-21 devnull 	double y;

140 cd5bae78 2004-04-21 devnull 

141 cd5bae78 2004-04-21 devnull 	y = x * (x-1);

142 cd5bae78 2004-04-21 devnull 	o->act.ra = o->del0.ra +

143 cd5bae78 2004-04-21 devnull 		x*o->del1.ra + y*o->del2.ra;

144 cd5bae78 2004-04-21 devnull 	o->act.decl2 = o->del0.decl2 +

145 cd5bae78 2004-04-21 devnull 		x*o->del1.decl2 + y*o->del2.decl2;

146 cd5bae78 2004-04-21 devnull 	o->act.semi2 = o->del0.semi2 +

147 cd5bae78 2004-04-21 devnull 		x*o->del1.semi2 + y*o->del2.semi2;

148 cd5bae78 2004-04-21 devnull 	o->act.el = o->del0.el +

149 cd5bae78 2004-04-21 devnull 		x*o->del1.el + y*o->del2.el;

150 cd5bae78 2004-04-21 devnull }

151 cd5bae78 2004-04-21 devnull 

152 cd5bae78 2004-04-21 devnull 

153 cd5bae78 2004-04-21 devnull double

154 cd5bae78 2004-04-21 devnull pinorm(double a)

155 cd5bae78 2004-04-21 devnull {

156 cd5bae78 2004-04-21 devnull 

157 cd5bae78 2004-04-21 devnull 	while(a < -pi)

158 cd5bae78 2004-04-21 devnull 		a += pipi;

159 cd5bae78 2004-04-21 devnull 	while(a > pi)

160 cd5bae78 2004-04-21 devnull 		a -= pipi;

161 cd5bae78 2004-04-21 devnull 	return a;

162 cd5bae78 2004-04-21 devnull }

163 cd5bae78 2004-04-21 devnull 

164 cd5bae78 2004-04-21 devnull void

165 cd5bae78 2004-04-21 devnull set3pt(Obj2 *p, int i, Occ *o)

166 cd5bae78 2004-04-21 devnull {

167 cd5bae78 2004-04-21 devnull 	Obj1 *p1, *p2, *p3;

168 cd5bae78 2004-04-21 devnull 	double a;

169 cd5bae78 2004-04-21 devnull 

170 cd5bae78 2004-04-21 devnull 	p1 = p->point+i;

171 cd5bae78 2004-04-21 devnull 	p2 = p1+1;

172 cd5bae78 2004-04-21 devnull 	p3 = p2+1;

173 cd5bae78 2004-04-21 devnull 

174 cd5bae78 2004-04-21 devnull 	o->del0.ra = p1->ra;

175 cd5bae78 2004-04-21 devnull 	o->del0.decl2 = p1->decl2;

176 cd5bae78 2004-04-21 devnull 	o->del0.semi2 = p1->semi2;

177 cd5bae78 2004-04-21 devnull 	o->del0.el = p1->el;

178 cd5bae78 2004-04-21 devnull 

179 cd5bae78 2004-04-21 devnull 	a = p2->ra - p1->ra;

180 cd5bae78 2004-04-21 devnull 	o->del1.ra = pinorm(a);

181 cd5bae78 2004-04-21 devnull 	a = p2->decl2 - p1->decl2;

182 cd5bae78 2004-04-21 devnull 	o->del1.decl2 = pinorm(a);

183 cd5bae78 2004-04-21 devnull 	o->del1.semi2 = p2->semi2 - p1->semi2;

184 cd5bae78 2004-04-21 devnull 	o->del1.el = p2->el - p1->el;

185 cd5bae78 2004-04-21 devnull 

186 cd5bae78 2004-04-21 devnull 	a = p1->ra + p3->ra - 2*p2->ra;

187 cd5bae78 2004-04-21 devnull 	o->del2.ra = pinorm(a)/2;

188 cd5bae78 2004-04-21 devnull 	a = p1->decl2 + p3->decl2 - 2*p2->decl2;

189 cd5bae78 2004-04-21 devnull 	o->del2.decl2 = pinorm(a)/2;

190 cd5bae78 2004-04-21 devnull 	o->del2.semi2 = (p1->semi2 + p3->semi2 - 2*p2->semi2) / 2;

191 cd5bae78 2004-04-21 devnull 	o->del2.el = (p1->el + p3->el - 2*p2->el) / 2;

192 cd5bae78 2004-04-21 devnull }