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

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

update pour version 3.6.3 Xephem cmv 17/01/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;
[2643]108 double pra, pdec;
[2551]109 MoonData *mdp;
110 Obj *planop;
111
112 /* init plobj[] */
113 if (!((Obj *)&plobj[0])->o_type) {
114 initPlobj();
115 initDir();
116 }
117
118 /* get sun @ np */
119 if (memcmp (&plnow[SUN], np, sizeof(Now))) {
120 obj_cir (np, (Obj*)&plobj[SUN]);
121 memcpy (&plnow[SUN], np, sizeof(Now));
122 }
123
124 /* get parent planet and moon info @ np */
125 switch (moonop->pl_code) {
126
127 case MARS:
128 case PHOBOS:
129 case DEIMOS:
130
131 planop = (Obj*)&plobj[MARS];
132
133 if (memcmp (&plnow[MARS], np, sizeof(Now))) {
134 obj_cir (np, planop);
135 memcpy (&plnow[MARS], np, sizeof(Now));
136 }
137
138 /* don't worry, this already caches based on same mjd */
[2643]139 marsm_data (mjd, moondir, sunop, planop, &sz, &pra, &pdec, md);
[2551]140 mdp = &md[moonop->pl_moon];
141 break;
142
143 case JUPITER:
144 case IO:
145 case EUROPA:
146 case GANYMEDE:
147 case CALLISTO:
148
149 planop = (Obj*)&plobj[JUPITER];
150
151 if (memcmp (&plnow[JUPITER], np, sizeof(Now))) {
152 obj_cir (np, planop);
153 memcpy (&plnow[JUPITER], np, sizeof(Now));
154 }
155
156 /* don't worry, this already caches based on same mjd */
[2643]157 jupiter_data (mjd,moondir,sunop,planop,&sz,&t1,&t2,&pra,&pdec,md);
[2551]158 mdp = &md[moonop->pl_moon];
159 moonop->pl_aux1 = t1;
160 moonop->pl_aux2 = t2;
161 break;
162
163 case SATURN:
164 case MIMAS:
165 case ENCELADUS:
166 case TETHYS:
167 case DIONE:
168 case RHEA:
169 case TITAN:
170 case HYPERION:
171 case IAPETUS:
172
173 planop = (Obj*)&plobj[SATURN];
174
175 if (memcmp (&plnow[SATURN], np, sizeof(Now))) {
176 obj_cir (np, planop);
177 memcpy (&plnow[SATURN], np, sizeof(Now));
178 }
179
180 /* don't worry, this already caches based on same mjd */
[2643]181 saturn_data (mjd,moondir,sunop,planop,&sz,&t1,&t2,&pra,&pdec,md);
[2551]182 mdp = &md[moonop->pl_moon];
183 moonop->pl_aux1 = t1;
184 moonop->pl_aux2 = t2;
185 break;
186
187 case URANUS:
188 case ARIEL:
189 case UMBRIEL:
190 case TITANIA:
191 case OBERON:
192 case MIRANDA:
193
194 planop = (Obj*)&plobj[URANUS];
195
196 if (memcmp (&plnow[URANUS], np, sizeof(Now))) {
197 obj_cir (np, planop);
198 memcpy (&plnow[URANUS], np, sizeof(Now));
199 }
200
201 /* don't worry, this already caches based on same mjd */
[2643]202 uranus_data (mjd, moondir, sunop, planop, &sz, &pra, &pdec, md);
[2551]203 mdp = &md[moonop->pl_moon];
204 break;
205
206 default:
207
208 printf ("Called plmoon_cir with bad code: %d\n",moonop->pl_code);
209 return (-1);
210
211 }
212
213 /* set moonop */
214 setMoon (np, moonop, planop, mdp);
215
216 return (0);
217}
218
219static void
220initPlobj()
221{
222 int i;
223
224 for (i = 0; i < NOBJ; i++) {
225 ((Obj*)&plobj[i])->o_type = PLANET;
226 ((Obj*)&plobj[i])->pl_code = i;
227 }
228}
229
230/* install a default if no moondir set yet */
231static void
232initDir()
233{
234 if (!moondir)
235 moondir = ".";
236}
237
238/* set moonop->s_* fields.
239 * np is needed to get local parallactic angle.
240 */
241static void
242setMoon (Now *np, Obj *moonop, Obj *planop, MoonData *mdp)
243{
244 double pa, dra, ddec;
245
246 /* just copy most fields from planet for now */
247 moonop->s_gaera = planop->s_gaera; /* TODO */
248 moonop->s_gaedec = planop->s_gaedec; /* TODO */
249 moonop->s_elong = planop->s_elong; /* TODO */
250 moonop->s_size = 0; /* TODO */
251 moonop->s_sdist = planop->s_sdist; /* TODO */
252 moonop->s_edist = planop->s_edist; /* TODO */
253 moonop->s_hlat = planop->s_hlat; /* TODO */
254 moonop->s_hlong = planop->s_hlong; /* TODO */
255 moonop->s_phase = planop->s_phase; /* TODO */
256
257 /* new ra/dec directly from mdp */
258 moonop->s_ra = mdp->ra;
259 moonop->s_dec = mdp->dec;
260
261 /* geoemtry info */
262 moonop->pl_x = mdp->x;
263 moonop->pl_y = mdp->y;
264 moonop->pl_z = mdp->z;
265 moonop->pl_evis = mdp->evis;
266 moonop->pl_svis = mdp->svis;
267
268 /* tweak alt/az by change in ra/dec rotated by pa */
269 pa = parallacticLDA (lat, planop->s_dec, planop->s_alt);
270 if (planop->s_az < PI)
271 pa = -pa; /* rotation radec to altaz */
272 dra = (moonop->s_ra - planop->s_ra)*cos(planop->s_dec);
273 ddec = moonop->s_dec - planop->s_dec;
274 rotate (pa, &dra, &ddec);
275 moonop->s_alt = planop->s_alt + ddec;
276 moonop->s_az = planop->s_az - dra/cos(planop->s_alt);
277
278 /* new mag directly from mdp */
279 set_smag (moonop, mdp->mag);
280
281 /* name */
282 strcpy (moonop->o_name, mdp->full);
283}
284
285/* rotate ccw by a */
286static void
287rotate (double a, double *x, double *y)
288{
289 double sa = sin(a);
290 double ca = cos(a);
291 double xp = (*x)*ca - (*y)*sa;
292 double yp = (*x)*sa + (*y)*ca;
293 *x = xp;
294 *y = yp;
295}
Note: See TracBrowser for help on using the repository browser.