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

Last change on this file since 2411 was 1719, checked in by cmv, 24 years ago

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