source: Sophya/trunk/SophyaExt/XephemAstroLib/umoon.c@ 3660

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

mise a niveau Xephem 3.7.4, cmv 16/07/2009
File size: 6.5 KB
Line 
1/* uranus moon info */
2
3#include <stdio.h>
4#include <stdlib.h>
5#include <string.h>
6#include <errno.h>
7#include <math.h>
8
9#include "astro.h"
10#include "bdl.h"
11
12static int use_bdl (double jd, char *dir, MoonData md[U_NMOONS]);
13static void moonradec (double usize, MoonData md[U_NMOONS]);
14static void moonSVis (Obj *sop, Obj *uop, MoonData md[U_NMOONS]);
15static void moonEVis (MoonData md[U_NMOONS]);
16static void moonPShad (Obj *sop, Obj *uop, MoonData md[U_NMOONS]);
17static void moonTrans (MoonData md[U_NMOONS]);
18
19/* moon table and a few other goodies and when it was last computed */
20static double mdmjd = -123456;
21static MoonData umd[U_NMOONS] = {
22 {"Uranus", NULL},
23 {"Ariel", "I"},
24 {"Umbriel", "II"},
25 {"Titania", "III"},
26 {"Oberon", "IV"},
27 {"Miranda", "V"},
28};
29static double sizemjd; /* size at last mjd */
30
31/* These values are from the Explanatory Supplement.
32 * Precession degrades them gradually over time.
33 */
34#define POLE_RA degrad(257.43) /* RA of Uranus' north pole */
35#define POLE_DEC degrad(-15.10) /* Dec of Uranus' north pole */
36
37
38/* get uranus info in md[0], moon info in md[1..U_NMOONS-1].
39 * if !uop caller just wants md[] for names
40 * N.B. we assume sop and uop are updated.
41 */
42void
43uranus_data (
44double Mjd, /* mjd */
45char dir[], /* dir in which to look for helper files */
46Obj *sop, /* Sun */
47Obj *uop, /* uranus */
48double *sizep, /* u angular diam, rads */
49double *polera, double *poledec, /* pole location */
50MoonData md[U_NMOONS]) /* return info */
51{
52 double JD;
53
54 /* always copy back at least for name */
55 memcpy (md, umd, sizeof(umd));
56
57 /* pole */
58 if (polera) *polera = POLE_RA;
59 if (poledec) *poledec = POLE_DEC;
60
61 /* nothing else if repeat call or just want names */
62 if (Mjd == mdmjd || !uop) {
63 if (uop) {
64 *sizep = sizemjd;
65 }
66 return;
67 }
68 JD = Mjd + MJD0;
69
70 /* planet in [0] */
71 md[0].ra = uop->s_ra;
72 md[0].dec = uop->s_dec;
73 md[0].mag = get_mag(uop);
74 md[0].x = 0;
75 md[0].y = 0;
76 md[0].z = 0;
77 md[0].evis = 1;
78 md[0].svis = 1;
79
80 /* size is straight from uop */
81 *sizep = degrad(uop->s_size/3600.0);
82
83 /* from Pasachoff/Menzel */
84
85 md[1].mag = 14.2;
86 md[2].mag = 14.8;
87 md[3].mag = 13.7;
88 md[4].mag = 14.0;
89 md[5].mag = 16.3;
90
91 /* get moon x,y,z from BDL if possible */
92 if (!dir || use_bdl (JD, dir, md) < 0) {
93 int i;
94 for (i = 1; i < U_NMOONS; i++)
95 md[i].x = md[i].y = md[i].z = 0.0;
96 fprintf (stderr, "No uranus model available\n");
97 }
98
99 /* set visibilities */
100 moonSVis (sop, uop, md);
101 moonPShad (sop, uop, md);
102 moonEVis (md);
103 moonTrans (md);
104
105 /* fill in moon ra and dec */
106 moonradec (*sizep, md);
107
108 /* save */
109 mdmjd = Mjd;
110 sizemjd = *sizep;
111 memcpy (umd, md, sizeof(umd));
112}
113
114/* hunt for BDL file in dir[] and use if possible
115 * return 0 if ok, else -1
116 */
117static int
118use_bdl (
119double JD, /* julian date */
120char dir[], /* directory */
121MoonData md[U_NMOONS]) /* fill md[1..NM-1].x/y/z for each moon */
122{
123#define URAU .0001597 /* Uranus radius, AU */
124 double x[U_NMOONS], y[U_NMOONS], z[U_NMOONS];
125 char buf[1024];
126 FILE *fp;
127 char *fn;
128 int i;
129
130 /* check ranges and appropriate data file */
131 if (JD < 2451179.50000) /* Jan 1 1999 UTC */
132 return (-1);
133 if (JD < 2455562.5) /* Jan 1 2011 UTC */
134 fn = "uranus.9910";
135 else if (JD < 2459215.5) /* Jan 1 2021 UTC */
136 fn = "uranus.1020";
137 else
138 return (-1);
139
140 /* open */
141 (void) sprintf (buf, "%s/%s", dir, fn);
142 fp = fopen (buf, "r");
143 if (!fp) {
144 fprintf (stderr, "%s: %s\n", fn, strerror(errno));
145 return (-1);
146 }
147
148 /* use it */
149 if ((i = read_bdl (fp, JD, x, y, z, buf)) < 0) {
150 fprintf (stderr, "%s: %s\n", fn, buf);
151 fclose (fp);
152 return (-1);
153 }
154 if (i != U_NMOONS-1) {
155 fprintf (stderr, "%s: BDL says %d moons, code expects %d", fn,
156 i, U_NMOONS-1);
157 fclose (fp);
158 return (-1);
159 }
160
161 /* copy into md[1..NM-1] with our scale and sign conventions */
162 for (i = 1; i < U_NMOONS; i++) {
163 md[i].x = x[i-1]/URAU; /* we want u radii +E */
164 md[i].y = -y[i-1]/URAU; /* we want u radii +S */
165 md[i].z = -z[i-1]/URAU; /* we want u radii +front */
166 }
167
168 /* ok */
169 fclose (fp);
170 return (0);
171}
172
173/* given uranus loc in md[0].ra/dec and size, and location of each moon in
174 * md[1..NM-1].x/y in ura radii, find ra/dec of each moon in md[1..NM-1].ra/dec.
175 */
176static void
177moonradec (
178double usize, /* ura diameter, rads */
179MoonData md[U_NMOONS]) /* fill in RA and Dec */
180{
181 double urad = usize/2;
182 double ura = md[0].ra;
183 double udec = md[0].dec;
184 int i;
185
186 for (i = 1; i < U_NMOONS; i++) {
187 double dra = urad * md[i].x;
188 double ddec = urad * md[i].y;
189 md[i].ra = ura + dra;
190 md[i].dec = udec - ddec;
191 }
192}
193
194/* set svis according to whether moon is in sun light */
195static void
196moonSVis(
197Obj *sop, /* SUN */
198Obj *uop, /* uranus */
199MoonData md[U_NMOONS])
200{
201 double esd = sop->s_edist;
202 double eod = uop->s_edist;
203 double sod = uop->s_sdist;
204 double soa = degrad(uop->s_elong);
205 double esa = asin(esd*sin(soa)/sod);
206 double h = sod*uop->s_hlat;
207 double nod = h*(1./eod - 1./sod);
208 double sca = cos(esa), ssa = sin(esa);
209 int i;
210
211 for (i = 1; i < U_NMOONS; i++) {
212 MoonData *mdp = &md[i];
213 double xp = sca*mdp->x + ssa*mdp->z;
214 double yp = mdp->y;
215 double zp = -ssa*mdp->x + sca*mdp->z;
216 double ca = cos(nod), sa = sin(nod);
217 double xpp = xp;
218 double ypp = ca*yp - sa*zp;
219 double zpp = sa*yp + ca*zp;
220 int outside = xpp*xpp + ypp*ypp > 1.0;
221 int infront = zpp > 0.0;
222 mdp->svis = outside || infront;
223 }
224}
225
226/* set evis according to whether moon is geometrically visible from earth */
227static void
228moonEVis (MoonData md[U_NMOONS])
229{
230 int i;
231
232 for (i = 1; i < U_NMOONS; i++) {
233 MoonData *mdp = &md[i];
234 int outside = mdp->x*mdp->x + mdp->y*mdp->y > 1.0;
235 int infront = mdp->z > 0.0;
236 mdp->evis = outside || infront;
237 }
238}
239
240/* set pshad and sx,sy shadow info */
241static void
242moonPShad(
243Obj *sop, /* SUN */
244Obj *uop, /* uranus */
245MoonData md[U_NMOONS])
246{
247 int i;
248
249 for (i = 1; i < U_NMOONS; i++) {
250 MoonData *mdp = &md[i];
251 mdp->pshad = !plshadow (uop, sop, POLE_RA, POLE_DEC, mdp->x,
252 mdp->y, mdp->z, &mdp->sx, &mdp->sy);
253 }
254}
255
256/* set whether moons are transiting */
257static void
258moonTrans (MoonData md[U_NMOONS])
259{
260 int i;
261
262 for (i = 1; i < U_NMOONS; i++) {
263 MoonData *mdp = &md[i];
264 mdp->trans = mdp->z > 0 && mdp->x*mdp->x + mdp->y*mdp->y < 1;
265 }
266}
267
268
269/* For RCS Only -- Do Not Edit */
270static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: umoon.c,v $ $Date: 2009-07-16 10:34:39 $ $Revision: 1.6 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.