source: Sophya/trunk/SophyaExt/XephemAstroLib/riset_cir.c@ 3451

Last change on this file since 3451 was 3111, checked in by cmv, 19 years ago

mise en conformite xephem 3.7.2 cmv 22/11/2006
File size: 10.0 KB
Line 
1/* find rise and set circumstances, ie, riset_cir() and related functions. */
2
3#include <stdio.h>
4#include <math.h>
5#include <stdlib.h>
6#include <string.h>
7
8#include "astro.h"
9
10#define TMACC (10./3600./24.0) /* convergence accuracy, days */
11
12static void e_riset_cir (Now *np, Obj *op, double dis, RiseSet *rp);
13static int find_0alt (double dt, double dis, Now *np, Obj *op);
14static int find_transit (double dt, Now *np, Obj *op);
15static int find_max (Now *np, Obj *op, double tr, double ts, double *tp,
16 double *alp);
17
18/* find where and when an object, op, will rise and set and
19 * it's transit circumstances. all times are utc mjd, angles rads e of n.
20 * dis is the angle down from an ideal horizon, in rads (see riset()).
21 * N.B. dis should NOT include refraction, we do that here.
22 */
23void
24riset_cir (Now *np, Obj *op, double dis, RiseSet *rp)
25{
26 double mjdn; /* mjd of local noon */
27 double lstn; /* lst at local noon */
28 double lr, ls; /* lst rise/set times */
29 double ar, as; /* az of rise/set */
30 double ran; /* RA at noon */
31 Now n; /* copy to move time around */
32 Obj o; /* copy to get circumstances at n */
33 int rss; /* temp status */
34
35 /* work with local copies so we can move the time around */
36 (void) memcpy ((void *)&n, (void *)np, sizeof(n));
37 (void) memcpy ((void *)&o, (void *)op, sizeof(o));
38
39 /* fast Earth satellites need a different approach.
40 * "fast" here is pretty arbitrary -- just too fast to work with the
41 * iterative approach based on refining the times for a "fixed" object.
42 */
43 if (op->o_type == EARTHSAT && op->es_n > FAST_SAT_RPD) {
44 e_riset_cir (&n, &o, dis, rp);
45 return;
46 }
47
48 /* assume no problems initially */
49 rp->rs_flags = 0;
50
51 /* start the iteration at local noon */
52 mjdn = mjd_day(mjd - tz/24.0) + tz/24.0 + 0.5;
53 n.n_mjd = mjdn;
54 now_lst (&n, &lstn);
55
56 /* first approximation is to find rise/set times of a fixed object
57 * at the current epoch in its position at local noon.
58 * N.B. add typical refraction if dis is above horizon for initial
59 * go/no-go test. if it passes, real code does refraction rigorously.
60 */
61 n.n_mjd = mjdn;
62 if (obj_cir (&n, &o) < 0) {
63 rp->rs_flags = RS_ERROR;
64 return;
65 }
66 ran = o.s_gaera;
67 riset (o.s_gaera, o.s_gaedec, lat, dis+(dis>.01 ? 0 : .01), &lr, &ls,
68 &ar, &as, &rss);
69 switch (rss) {
70 case 0: break;
71 case 1: rp->rs_flags = RS_NEVERUP; return;
72 case -1: rp->rs_flags = RS_CIRCUMPOLAR; goto dotransit;
73 default: rp->rs_flags = RS_ERROR; return;
74 }
75
76 /* iterate to find better rise time */
77 n.n_mjd = mjdn;
78 switch (find_0alt ((lr - lstn)/SIDRATE, dis, &n, &o)) {
79 case 0: /* ok */
80 rp->rs_risetm = n.n_mjd;
81 rp->rs_riseaz = o.s_az;
82 break;
83 case -1: /* obj_cir error */
84 rp->rs_flags |= RS_RISERR;
85 break;
86 case -2: /* converged but not today */ /* FALLTHRU */
87 case -3: /* probably never up */
88 rp->rs_flags |= RS_NORISE;
89 break;
90 }
91
92 /* iterate to find better set time */
93 n.n_mjd = mjdn;
94 switch (find_0alt ((ls - lstn)/SIDRATE, dis, &n, &o)) {
95 case 0: /* ok */
96 rp->rs_settm = n.n_mjd;
97 rp->rs_setaz = o.s_az;
98 break;
99 case -1: /* obj_cir error */
100 rp->rs_flags |= RS_SETERR;
101 break;
102 case -2: /* converged but not today */ /* FALLTHRU */
103 case -3: /* probably circumpolar */
104 rp->rs_flags |= RS_NOSET;
105 break;
106 }
107
108 /* can try transit even if rise or set failed */
109 dotransit:
110 n.n_mjd = mjdn;
111 switch (find_transit ((radhr(ran) - lstn)/SIDRATE, &n, &o)) {
112 case 0: /* ok */
113 rp->rs_trantm = n.n_mjd;
114 rp->rs_tranalt = o.s_alt;
115 break;
116 case -1: /* did not converge */
117 rp->rs_flags |= RS_TRANSERR;
118 break;
119 case -2: /* converged but not today */
120 rp->rs_flags |= RS_NOTRANS;
121 break;
122 }
123}
124
125/* find local times when sun is dis rads below horizon.
126 */
127void
128twilight_cir (Now *np, double dis, double *dawn, double *dusk, int *status)
129{
130 RiseSet rs;
131 Obj o;
132
133 memset (&o, 0, sizeof(o));
134 o.o_type = PLANET;
135 o.pl_code = SUN;
136 (void) strcpy (o.o_name, "Sun");
137 riset_cir (np, &o, dis, &rs);
138 *dawn = rs.rs_risetm;
139 *dusk = rs.rs_settm;
140 *status = rs.rs_flags;
141}
142
143/* find where and when a fast-moving Earth satellite, op, will rise and set and
144 * it's transit circumstances. all times are mjd, angles rads e of n.
145 * dis is the angle down from the local topo horizon, in rads (see riset()).
146 * idea is to walk forward in time looking for alt+dis==0 crossings.
147 * initial time step is a few degrees (based on average daily motion).
148 * we stop as soon as we see both a rise and set.
149 * N.B. we assume *np and *op are working copies we can mess up.
150 */
151static void
152e_riset_cir (Now *np, Obj *op, double dis, RiseSet *rp)
153{
154#define DEGSTEP 5 /* time step is about this many degrees */
155 int steps; /* max number of time steps */
156 double dt; /* time change per step, days */
157 double t0, t1; /* current and next mjd values */
158 double a0, a1; /* altitude at t0 and t1 */
159 int rise, set; /* flags to check when we find these events */
160 int i;
161
162 dt = DEGSTEP * (1.0/360.0/op->es_n);
163 steps = (int)(1.0/dt);
164 rise = set = 0;
165 rp->rs_flags = 0;
166
167 if (obj_cir (np, op) < 0) {
168 rp->rs_flags |= RS_ERROR;
169 return;
170 }
171
172 t0 = mjd;
173 a0 = op->s_alt + dis;
174
175 for (i = 0; i < steps && (!rise || !set); i++) {
176 mjd = t1 = t0 + dt;
177 if (obj_cir (np, op) < 0) {
178 rp->rs_flags |= RS_ERROR;
179 return;
180 }
181 a1 = op->s_alt + dis;
182
183 if (a0 < 0 && a1 > 0 && !rise) {
184 /* found a rise event -- interate to refine */
185 switch (find_0alt (0.0, dis, np, op)) {
186 case 0: /* ok */
187 rp->rs_risetm = np->n_mjd;
188 rp->rs_riseaz = op->s_az;
189 rise = 1;
190 break;
191 case -1: /* obj_cir error */
192 rp->rs_flags |= RS_RISERR;
193 return;
194 case -2: /* converged but not today */ /* FALLTHRU */
195 case -3: /* probably never up */
196 rp->rs_flags |= RS_NORISE;
197 return;
198 }
199 } else if (a0 > 0 && a1 < 0 && !set) {
200 /* found a setting event -- interate to refine */
201 switch (find_0alt (0.0, dis, np, op)) {
202 case 0: /* ok */
203 rp->rs_settm = np->n_mjd;
204 rp->rs_setaz = op->s_az;
205 set = 1;
206 break;
207 case -1: /* obj_cir error */
208 rp->rs_flags |= RS_SETERR;
209 return;
210 case -2: /* converged but not today */ /* FALLTHRU */
211 case -3: /* probably circumpolar */
212 rp->rs_flags |= RS_NOSET;
213 return;
214 }
215 }
216
217 t0 = t1;
218 a0 = a1;
219 }
220
221 /* instead of transit, for satellites we find time of maximum
222 * altitude, if we know both the rise and set times and the former
223 * occurs before the latter.
224 */
225 if (rise && set && rp->rs_risetm < rp->rs_settm) {
226 double tt, al;
227 if (find_max (np, op, rp->rs_risetm, rp->rs_settm, &tt, &al) < 0) {
228 rp->rs_flags |= RS_TRANSERR;
229 return;
230 }
231 rp->rs_trantm = tt;
232 rp->rs_tranalt = al;
233 } else
234 rp->rs_flags |= RS_NOTRANS;
235
236 /* check for some bad conditions */
237 if (!rise) {
238 if (a0 > 0)
239 rp->rs_flags |= RS_CIRCUMPOLAR;
240 else
241 rp->rs_flags |= RS_NORISE;
242 }
243 if (!set) {
244 if (a0 < 0)
245 rp->rs_flags |= RS_NEVERUP;
246 else
247 rp->rs_flags |= RS_NOSET;
248 }
249}
250
251/* given a Now at noon and a dt from noon, in hours, for a first approximation
252 * to a rise or set event, refine the event by searching for when alt+dis = 0.
253 * return 0: if find one within 12 hours of noon with np and op set to the
254 * better time and circumstances;
255 * return -1: if error from obj_cir;
256 * return -2: if converges but not today;
257 * return -3: if does not converge at all (probably circumpolar or never up);
258 */
259static int
260find_0alt (
261double dt, /* hours from noon to first guess at event */
262double dis, /* horizon displacement, rads */
263Now *np, /* working Now -- starts with mjd is noon, returns as answer */
264Obj *op) /* working object -- returns as answer */
265{
266#define MAXPASSES 20 /* max iterations to try */
267#define FIRSTSTEP (1.0/60.0/24.0) /* first time step, days */
268
269 double a0 = 0;
270 double mjdn = mjd;
271 int npasses;
272
273 /* insure initial guess is today -- if not, move by 24 hours */
274 if (dt < -12.0 && !find_0alt (dt+24, dis, np, op))
275 return (0);
276 mjd = mjdn;
277 if (dt > 12.0 && !find_0alt (dt-24, dis, np, op))
278 return (0);
279 mjd = mjdn;
280
281 /* convert dt to days for remainder of algorithm */
282 dt /= 24.0;
283
284 /* use secant method to look for s_alt + dis == 0 */
285 npasses = 0;
286 do {
287 double a1;
288
289 mjd += dt;
290 if (obj_cir (np, op) < 0)
291 return (-1);
292 a1 = op->s_alt;
293
294 dt = (npasses == 0) ? FIRSTSTEP : (dis+a1)*dt/(a0-a1);
295 a0 = a1;
296
297 } while (++npasses < MAXPASSES && fabs(dt) > TMACC);
298
299 /* return codes */
300 if (npasses == MAXPASSES)
301 return (-3);
302 return (fabs(mjdn-mjd) < .5 ? 0 : -2);
303
304#undef MAXPASSES
305#undef FIRSTSTEP
306}
307
308/* find when the given object transits. start the search when LST matches the
309 * object's RA at noon.
310 * if ok, return 0 with np and op set to the transit conditions; if can't
311 * converge return -1; if converges ok but not today return -2.
312 * N.B. we assume np is passed set to local noon.
313 */
314static int
315find_transit (double dt, Now *np, Obj *op)
316{
317#define MAXLOOPS 10
318#define MAXERR (0.25/60.) /* hours */
319 double mjdn = mjd;
320 double lst;
321 int i;
322
323 /* insure initial guess is today -- if not, move by 24 hours */
324 if (dt < -12.0)
325 dt += 24.0;
326 if (dt > 12.0)
327 dt -= 24.0;
328
329 i = 0;
330 do {
331 mjd += dt/24.0;
332 if (obj_cir (np, op) < 0)
333 return (-1);
334 now_lst (np, &lst);
335 dt = (radhr(op->s_gaera) - lst);
336 if (dt < -12.0)
337 dt += 24.0;
338 if (dt > 12.0)
339 dt -= 24.0;
340 } while (++i < MAXLOOPS && fabs(dt) > MAXERR);
341
342 /* return codes */
343 if (i == MAXLOOPS)
344 return (-1);
345 return (fabs(mjd - mjdn) < 0.5 ? 0 : -2);
346
347#undef MAXLOOPS
348#undef MAXERR
349}
350
351/* find the mjd time of max altitude between the given rise and set times.
352 * N.B. we assume *np and *op are working copies we can mess up.
353 * N.B. we just assume max occurs at the center time.
354 * return 0 if ok, else -1.
355 */
356static int
357find_max (
358Now *np,
359Obj *op,
360double tr, double ts, /* times of rise and set */
361double *tp, double *alp) /* time of max altitude, and that altitude */
362{
363 mjd = (ts + tr)/2;
364 if (obj_cir (np, op) < 0)
365 return (-1);
366 *tp = mjd;
367 *alp = op->s_alt;
368 return (0);
369}
370
371/* For RCS Only -- Do Not Edit */
372static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: riset_cir.c,v $ $Date: 2006-11-22 13:53:30 $ $Revision: 1.6 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.