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

Last change on this file since 3572 was 2818, checked in by cmv, 20 years ago

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