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

Last change on this file since 3848 was 3654, checked in by cmv, 16 years ago

mise a niveau Xephem 3.7.4, cmv 16/07/2009
File size: 10.3 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, err but give times anyway */
87 rp->rs_risetm = n.n_mjd;
88 rp->rs_riseaz = o.s_az;
89 rp->rs_flags |= RS_NORISE;
90 break;
91 case -3: /* probably never up */
92 rp->rs_flags |= RS_NEVERUP;
93 break;
94 }
95
96 /* iterate to find better set time */
97 n.n_mjd = mjdn;
98 switch (find_0alt ((ls - lstn)/SIDRATE, dis, &n, &o)) {
99 case 0: /* ok */
100 rp->rs_settm = n.n_mjd;
101 rp->rs_setaz = o.s_az;
102 break;
103 case -1: /* obj_cir error */
104 rp->rs_flags |= RS_SETERR;
105 break;
106 case -2: /* converged but not today, err but give times anyway */
107 rp->rs_settm = n.n_mjd;
108 rp->rs_setaz = o.s_az;
109 rp->rs_flags |= RS_NOSET;
110 break;
111 case -3: /* probably circumpolar */
112 rp->rs_flags |= RS_CIRCUMPOLAR;
113 break;
114 }
115
116 /* can try transit even if rise or set failed */
117 dotransit:
118 n.n_mjd = mjdn;
119 switch (find_transit ((radhr(ran) - lstn)/SIDRATE, &n, &o)) {
120 case 0: /* ok */
121 rp->rs_trantm = n.n_mjd;
122 rp->rs_tranalt = o.s_alt;
123 break;
124 case -1: /* did not converge */
125 rp->rs_flags |= RS_TRANSERR;
126 break;
127 case -2: /* converged but not today */
128 rp->rs_flags |= RS_NOTRANS;
129 break;
130 }
131}
132
133/* find local times when sun is dis rads below horizon.
134 */
135void
136twilight_cir (Now *np, double dis, double *dawn, double *dusk, int *status)
137{
138 RiseSet rs;
139 Obj o;
140
141 memset (&o, 0, sizeof(o));
142 o.o_type = PLANET;
143 o.pl_code = SUN;
144 (void) strcpy (o.o_name, "Sun");
145 riset_cir (np, &o, dis, &rs);
146 *dawn = rs.rs_risetm;
147 *dusk = rs.rs_settm;
148 *status = rs.rs_flags;
149}
150
151/* find where and when a fast-moving Earth satellite, op, will rise and set and
152 * it's transit circumstances. all times are mjd, angles rads e of n.
153 * dis is the angle down from the local topo horizon, in rads (see riset()).
154 * idea is to walk forward in time looking for alt+dis==0 crossings.
155 * initial time step is a few degrees (based on average daily motion).
156 * we stop as soon as we see both a rise and set.
157 * N.B. we assume *np and *op are working copies we can mess up.
158 */
159static void
160e_riset_cir (Now *np, Obj *op, double dis, RiseSet *rp)
161{
162#define DEGSTEP 5 /* time step is about this many degrees */
163 int steps; /* max number of time steps */
164 double dt; /* time change per step, days */
165 double t0, t1; /* current and next mjd values */
166 double a0, a1; /* altitude at t0 and t1 */
167 int rise, set; /* flags to check when we find these events */
168 int i;
169
170 dt = DEGSTEP * (1.0/360.0/op->es_n);
171 steps = (int)(1.0/dt);
172 rise = set = 0;
173 rp->rs_flags = 0;
174
175 if (obj_cir (np, op) < 0) {
176 rp->rs_flags |= RS_ERROR;
177 return;
178 }
179
180 t0 = mjd;
181 a0 = op->s_alt + dis;
182
183 for (i = 0; i < steps && (!rise || !set); i++) {
184 mjd = t1 = t0 + dt;
185 if (obj_cir (np, op) < 0) {
186 rp->rs_flags |= RS_ERROR;
187 return;
188 }
189 a1 = op->s_alt + dis;
190
191 if (a0 < 0 && a1 > 0 && !rise) {
192 /* found a rise event -- interate to refine */
193 switch (find_0alt (0.0, dis, np, op)) {
194 case 0: /* ok */
195 rp->rs_risetm = np->n_mjd;
196 rp->rs_riseaz = op->s_az;
197 rise = 1;
198 break;
199 case -1: /* obj_cir error */
200 rp->rs_flags |= RS_RISERR;
201 return;
202 case -2: /* converged but not today */ /* FALLTHRU */
203 case -3: /* probably never up */
204 rp->rs_flags |= RS_NORISE;
205 return;
206 }
207 } else if (a0 > 0 && a1 < 0 && !set) {
208 /* found a setting event -- interate to refine */
209 switch (find_0alt (0.0, dis, np, op)) {
210 case 0: /* ok */
211 rp->rs_settm = np->n_mjd;
212 rp->rs_setaz = op->s_az;
213 set = 1;
214 break;
215 case -1: /* obj_cir error */
216 rp->rs_flags |= RS_SETERR;
217 return;
218 case -2: /* converged but not today */ /* FALLTHRU */
219 case -3: /* probably circumpolar */
220 rp->rs_flags |= RS_NOSET;
221 return;
222 }
223 }
224
225 t0 = t1;
226 a0 = a1;
227 }
228
229 /* instead of transit, for satellites we find time of maximum
230 * altitude, if we know both the rise and set times.
231 */
232 if (rise && set) {
233 double tt, al;
234 if (find_max (np, op, rp->rs_risetm, rp->rs_settm, &tt, &al) < 0) {
235 rp->rs_flags |= RS_TRANSERR;
236 return;
237 }
238 rp->rs_trantm = tt;
239 rp->rs_tranalt = al;
240 } else
241 rp->rs_flags |= RS_NOTRANS;
242
243 /* check for some bad conditions */
244 if (!rise) {
245 if (a0 > 0)
246 rp->rs_flags |= RS_CIRCUMPOLAR;
247 else
248 rp->rs_flags |= RS_NORISE;
249 }
250 if (!set) {
251 if (a0 < 0)
252 rp->rs_flags |= RS_NEVERUP;
253 else
254 rp->rs_flags |= RS_NOSET;
255 }
256}
257
258/* given a Now at noon and a dt from noon, in hours, for a first approximation
259 * to a rise or set event, refine the event by searching for when alt+dis = 0.
260 * return 0: if find one within 12 hours of noon with np and op set to the
261 * better time and circumstances;
262 * return -1: if error from obj_cir;
263 * return -2: if converges but not today;
264 * return -3: if does not converge at all (probably circumpolar or never up);
265 */
266static int
267find_0alt (
268double dt, /* hours from noon to first guess at event */
269double dis, /* horizon displacement, rads */
270Now *np, /* working Now -- starts with mjd is noon, returns as answer */
271Obj *op) /* working object -- returns as answer */
272{
273#define MAXPASSES 20 /* max iterations to try */
274#define FIRSTSTEP (1.0/60.0/24.0) /* first time step, days */
275#define MAXSTEP (12.0/24.0)/* max time step,days (to detect flat)*/
276
277 double a0 = 0;
278 double mjdn = mjd;
279 int npasses;
280
281 /* insure initial guess is today -- if not, move by 24 hours */
282 if (dt < -12.0 && !find_0alt (dt+24, dis, np, op))
283 return (0);
284 mjd = mjdn;
285 if (dt > 12.0 && !find_0alt (dt-24, dis, np, op))
286 return (0);
287 mjd = mjdn;
288
289 /* convert dt to days for remainder of algorithm */
290 dt /= 24.0;
291
292 /* use secant method to look for s_alt + dis == 0 */
293 npasses = 0;
294 do {
295 double a1;
296
297 mjd += dt;
298 if (obj_cir (np, op) < 0)
299 return (-1);
300 a1 = op->s_alt;
301
302 dt = (npasses == 0) ? FIRSTSTEP : (dis+a1)*dt/(a0-a1);
303 a0 = a1;
304
305 if (++npasses > MAXPASSES || fabs(dt) >= MAXSTEP)
306 return (-3);
307
308 } while (fabs(dt)>TMACC);
309
310 /* return codes */
311 return (fabs(mjdn-mjd) < .5 ? 0 : -2);
312
313#undef MAXPASSES
314#undef FIRSTSTEP
315#undef MAXSTEP
316}
317
318/* find when the given object transits. start the search when LST matches the
319 * object's RA at noon.
320 * if ok, return 0 with np and op set to the transit conditions; if can't
321 * converge return -1; if converges ok but not today return -2.
322 * N.B. we assume np is passed set to local noon.
323 */
324static int
325find_transit (double dt, Now *np, Obj *op)
326{
327#define MAXLOOPS 10
328#define MAXERR (0.25/60.) /* hours */
329 double mjdn = mjd;
330 double lst;
331 int i;
332
333 /* insure initial guess is today -- if not, move by 24 hours */
334 if (dt < -12.0)
335 dt += 24.0;
336 if (dt > 12.0)
337 dt -= 24.0;
338
339 i = 0;
340 do {
341 mjd += dt/24.0;
342 if (obj_cir (np, op) < 0)
343 return (-1);
344 now_lst (np, &lst);
345 dt = (radhr(op->s_gaera) - lst);
346 if (dt < -12.0)
347 dt += 24.0;
348 if (dt > 12.0)
349 dt -= 24.0;
350 } while (++i < MAXLOOPS && fabs(dt) > MAXERR);
351
352 /* return codes */
353 if (i == MAXLOOPS)
354 return (-1);
355 return (fabs(mjd - mjdn) < 0.5 ? 0 : -2);
356
357#undef MAXLOOPS
358#undef MAXERR
359}
360
361/* find the mjd time of max altitude between the given rise and set times.
362 * N.B. we assume *np and *op are working copies we can mess up.
363 * N.B. we just assume max occurs at the center time.
364 * return 0 if ok, else -1.
365 */
366static int
367find_max (
368Now *np,
369Obj *op,
370double tr, double ts, /* times of rise and set */
371double *tp, double *alp) /* time of max altitude, and that altitude */
372{
373 /* want rise before set */
374 while (ts < tr)
375 tr -= 1.0/op->es_n;
376 mjd = (ts + tr)/2;
377 if (obj_cir (np, op) < 0)
378 return (-1);
379 *tp = mjd;
380 *alp = op->s_alt;
381 return (0);
382}
383
384/* For RCS Only -- Do Not Edit */
385static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: riset_cir.c,v $ $Date: 2009-07-16 10:34:39 $ $Revision: 1.8 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.