source: Sophya/trunk/SophyaExt/XephemAstroLib/plmoon.c@ 2551

Last change on this file since 2551 was 2551, checked in by cmv, 21 years ago

nouvelle version de xephem/libastro (3.6) cmv 15/6/04
File size: 7.2 KB
RevLine 
[2551]1/* compute Obj fields for natural satellites.
2 */
3
4#include <stdio.h>
5#include <string.h>
6#include <math.h>
7
8#include "astro.h"
9
10/* private cache of planet ephemerides and when they were computed
11 * N.B. don't use ones in builtin[] -- they are the user's responsibility.
12 */
13static ObjPl plobj[NOBJ];
14static Now plnow[NOBJ];
15
16/* public builtin storage
17 */
18static Obj builtin[NBUILTIN];
19
20static char *moondir;
21
22static void setMoon (Now *np, Obj *moonop, Obj *planop, MoonData *mdp);
23static void init1BI (int idx, int pl, int moon, char *name);
24static void initPlobj(void);
25static void initDir(void);
26static void rotate (double a, double *x, double *y);
27
28/* directory in which to look for auxil moon data files.
29 * N.B. caller must supply persistent storage.
30 */
31void
32setMoonDir (char *dir)
33{
34 moondir = dir;
35}
36
37/* return set of builtin objects.
38 * caller can use this storage but should never try to free anything.
39 */
40int
41getBuiltInObjs (Obj **opp)
42{
43 if (!builtin[MERCURY].o_name[0]) {
44 /* first time only */
45
46 init1BI (MERCURY, MERCURY, X_PLANET, "Mercury");
47
48 init1BI (VENUS, VENUS, X_PLANET, "Venus");
49
50 init1BI (MARS, MARS, X_PLANET, "Mars");
51 init1BI (PHOBOS, MARS, M_PHOBOS, "Phobos");
52 init1BI (DEIMOS, MARS, M_DEIMOS, "Deimos");
53
54 init1BI (JUPITER, JUPITER, X_PLANET, "Jupiter");
55 init1BI (IO, JUPITER, J_IO, "Io");
56 init1BI (EUROPA, JUPITER, J_EUROPA, "Europa");
57 init1BI (GANYMEDE, JUPITER, J_GANYMEDE, "Ganymede");
58 init1BI (CALLISTO, JUPITER, J_CALLISTO, "Callisto");
59
60 init1BI (SATURN, SATURN, X_PLANET, "Saturn");
61 init1BI (MIMAS, SATURN, S_MIMAS, "Mimas");
62 init1BI (ENCELADUS, SATURN, S_ENCELADUS, "Enceladus");
63 init1BI (TETHYS, SATURN, S_TETHYS, "Tethys");
64 init1BI (DIONE, SATURN, S_DIONE, "Dione");
65 init1BI (RHEA, SATURN, S_RHEA, "Rhea");
66 init1BI (TITAN, SATURN, S_TITAN, "Titan");
67 init1BI (HYPERION, SATURN, S_HYPERION, "Hyperion");
68 init1BI (IAPETUS, SATURN, S_IAPETUS, "Iapetus");
69
70 init1BI (URANUS, URANUS, X_PLANET, "Uranus");
71 init1BI (ARIEL, URANUS, U_ARIEL, "Ariel");
72 init1BI (UMBRIEL, URANUS, U_UMBRIEL, "Umbriel");
73 init1BI (TITANIA, URANUS, U_TITANIA, "Titania");
74 init1BI (OBERON, URANUS, U_OBERON, "Oberon");
75 init1BI (MIRANDA, URANUS, U_MIRANDA, "Miranda");
76
77 init1BI (NEPTUNE, NEPTUNE, X_PLANET, "Neptune");
78
79 init1BI (PLUTO, PLUTO, X_PLANET, "Pluto");
80
81 init1BI (SUN, SUN, X_PLANET, "Sun");
82
83 init1BI (MOON, MOON, X_PLANET, "Moon");
84 }
85
86 *opp = builtin;
87 return (NBUILTIN);
88}
89
90static void
91init1BI (int idx, int pl, int moon, char *name)
92{
93 strcpy (builtin[idx].o_name, name);
94 builtin[idx].o_type = PLANET;
95 builtin[idx].pl_code = pl;
96 builtin[idx].pl_moon = moon;
97}
98
99/* find the circumstances for natural satellite object op at np.
100 * TODO: distances and helio coords just copied from parent planet.
101 */
102int
103plmoon_cir (Now *np, Obj *moonop)
104{
105 Obj *sunop = (Obj*)&plobj[SUN];
106 MoonData md[X_MAXNMOONS];
107 double sz, t1, t2;
108 MoonData *mdp;
109 Obj *planop;
110
111 /* init plobj[] */
112 if (!((Obj *)&plobj[0])->o_type) {
113 initPlobj();
114 initDir();
115 }
116
117 /* get sun @ np */
118 if (memcmp (&plnow[SUN], np, sizeof(Now))) {
119 obj_cir (np, (Obj*)&plobj[SUN]);
120 memcpy (&plnow[SUN], np, sizeof(Now));
121 }
122
123 /* get parent planet and moon info @ np */
124 switch (moonop->pl_code) {
125
126 case MARS:
127 case PHOBOS:
128 case DEIMOS:
129
130 planop = (Obj*)&plobj[MARS];
131
132 if (memcmp (&plnow[MARS], np, sizeof(Now))) {
133 obj_cir (np, planop);
134 memcpy (&plnow[MARS], np, sizeof(Now));
135 }
136
137 /* don't worry, this already caches based on same mjd */
138 marsm_data (mjd, moondir, sunop, planop, &sz, md);
139 mdp = &md[moonop->pl_moon];
140 break;
141
142 case JUPITER:
143 case IO:
144 case EUROPA:
145 case GANYMEDE:
146 case CALLISTO:
147
148 planop = (Obj*)&plobj[JUPITER];
149
150 if (memcmp (&plnow[JUPITER], np, sizeof(Now))) {
151 obj_cir (np, planop);
152 memcpy (&plnow[JUPITER], np, sizeof(Now));
153 }
154
155 /* don't worry, this already caches based on same mjd */
156 jupiter_data (mjd, moondir, sunop, planop, &sz, &t1, &t2, md);
157 mdp = &md[moonop->pl_moon];
158 moonop->pl_aux1 = t1;
159 moonop->pl_aux2 = t2;
160 break;
161
162 case SATURN:
163 case MIMAS:
164 case ENCELADUS:
165 case TETHYS:
166 case DIONE:
167 case RHEA:
168 case TITAN:
169 case HYPERION:
170 case IAPETUS:
171
172 planop = (Obj*)&plobj[SATURN];
173
174 if (memcmp (&plnow[SATURN], np, sizeof(Now))) {
175 obj_cir (np, planop);
176 memcpy (&plnow[SATURN], np, sizeof(Now));
177 }
178
179 /* don't worry, this already caches based on same mjd */
180 saturn_data (mjd, moondir, sunop, planop, &sz, &t1, &t2, md);
181 mdp = &md[moonop->pl_moon];
182 moonop->pl_aux1 = t1;
183 moonop->pl_aux2 = t2;
184 break;
185
186 case URANUS:
187 case ARIEL:
188 case UMBRIEL:
189 case TITANIA:
190 case OBERON:
191 case MIRANDA:
192
193 planop = (Obj*)&plobj[URANUS];
194
195 if (memcmp (&plnow[URANUS], np, sizeof(Now))) {
196 obj_cir (np, planop);
197 memcpy (&plnow[URANUS], np, sizeof(Now));
198 }
199
200 /* don't worry, this already caches based on same mjd */
201 uranus_data (mjd, moondir, sunop, planop, &sz, md);
202 mdp = &md[moonop->pl_moon];
203 break;
204
205 default:
206
207 printf ("Called plmoon_cir with bad code: %d\n",moonop->pl_code);
208 return (-1);
209
210 }
211
212 /* set moonop */
213 setMoon (np, moonop, planop, mdp);
214
215 return (0);
216}
217
218static void
219initPlobj()
220{
221 int i;
222
223 for (i = 0; i < NOBJ; i++) {
224 ((Obj*)&plobj[i])->o_type = PLANET;
225 ((Obj*)&plobj[i])->pl_code = i;
226 }
227}
228
229/* install a default if no moondir set yet */
230static void
231initDir()
232{
233 if (!moondir)
234 moondir = ".";
235}
236
237/* set moonop->s_* fields.
238 * np is needed to get local parallactic angle.
239 */
240static void
241setMoon (Now *np, Obj *moonop, Obj *planop, MoonData *mdp)
242{
243 double pa, dra, ddec;
244
245 /* just copy most fields from planet for now */
246 moonop->s_gaera = planop->s_gaera; /* TODO */
247 moonop->s_gaedec = planop->s_gaedec; /* TODO */
248 moonop->s_elong = planop->s_elong; /* TODO */
249 moonop->s_size = 0; /* TODO */
250 moonop->s_sdist = planop->s_sdist; /* TODO */
251 moonop->s_edist = planop->s_edist; /* TODO */
252 moonop->s_hlat = planop->s_hlat; /* TODO */
253 moonop->s_hlong = planop->s_hlong; /* TODO */
254 moonop->s_phase = planop->s_phase; /* TODO */
255
256 /* new ra/dec directly from mdp */
257 moonop->s_ra = mdp->ra;
258 moonop->s_dec = mdp->dec;
259
260 /* geoemtry info */
261 moonop->pl_x = mdp->x;
262 moonop->pl_y = mdp->y;
263 moonop->pl_z = mdp->z;
264 moonop->pl_evis = mdp->evis;
265 moonop->pl_svis = mdp->svis;
266
267 /* tweak alt/az by change in ra/dec rotated by pa */
268 pa = parallacticLDA (lat, planop->s_dec, planop->s_alt);
269 if (planop->s_az < PI)
270 pa = -pa; /* rotation radec to altaz */
271 dra = (moonop->s_ra - planop->s_ra)*cos(planop->s_dec);
272 ddec = moonop->s_dec - planop->s_dec;
273 rotate (pa, &dra, &ddec);
274 moonop->s_alt = planop->s_alt + ddec;
275 moonop->s_az = planop->s_az - dra/cos(planop->s_alt);
276
277 /* new mag directly from mdp */
278 set_smag (moonop, mdp->mag);
279
280 /* name */
281 strcpy (moonop->o_name, mdp->full);
282}
283
284/* rotate ccw by a */
285static void
286rotate (double a, double *x, double *y)
287{
288 double sa = sin(a);
289 double ca = cos(a);
290 double xp = (*x)*ca - (*y)*sa;
291 double yp = (*x)*sa + (*y)*ca;
292 *x = xp;
293 *y = yp;
294}
Note: See TracBrowser for help on using the repository browser.