[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 | */
|
---|
| 13 | static ObjPl plobj[NOBJ];
|
---|
| 14 | static Now plnow[NOBJ];
|
---|
| 15 |
|
---|
| 16 | /* public builtin storage
|
---|
| 17 | */
|
---|
| 18 | static Obj builtin[NBUILTIN];
|
---|
| 19 |
|
---|
| 20 | static char *moondir;
|
---|
| 21 |
|
---|
| 22 | static void setMoon (Now *np, Obj *moonop, Obj *planop, MoonData *mdp);
|
---|
| 23 | static void init1BI (int idx, int pl, int moon, char *name);
|
---|
| 24 | static void initPlobj(void);
|
---|
| 25 | static 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 | */
|
---|
| 30 | void
|
---|
| 31 | setMoonDir (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 | */
|
---|
| 39 | int
|
---|
| 40 | getBuiltInObjs (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 |
|
---|
| 89 | static void
|
---|
| 90 | init1BI (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 | */
|
---|
| 101 | int
|
---|
| 102 | plmoon_cir (Now *np, Obj *moonop)
|
---|
| 103 | {
|
---|
| 104 | Obj *sunop = (Obj*)&plobj[SUN];
|
---|
| 105 | MoonData md[X_MAXNMOONS];
|
---|
| 106 | double sz, t1, t2;
|
---|
[2643] | 107 | double pra, pdec;
|
---|
[2551] | 108 | MoonData *mdp;
|
---|
| 109 | Obj *planop;
|
---|
| 110 |
|
---|
| 111 | /* init plobj[] */
|
---|
[2818] | 112 | if (!((Obj *)&plobj[0])->o_type)
|
---|
[2551] | 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 */
|
---|
[2643] | 136 | marsm_data (mjd, moondir, sunop, planop, &sz, &pra, &pdec, md);
|
---|
[2551] | 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 */
|
---|
[2643] | 154 | jupiter_data (mjd,moondir,sunop,planop,&sz,&t1,&t2,&pra,&pdec,md);
|
---|
[2551] | 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 */
|
---|
[2643] | 178 | saturn_data (mjd,moondir,sunop,planop,&sz,&t1,&t2,&pra,&pdec,md);
|
---|
[2551] | 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 */
|
---|
[2643] | 199 | uranus_data (mjd, moondir, sunop, planop, &sz, &pra, &pdec, md);
|
---|
[2551] | 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 |
|
---|
| 216 | static void
|
---|
| 217 | initPlobj()
|
---|
| 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 | */
|
---|
| 230 | static void
|
---|
| 231 | setMoon (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 */
|
---|
| 275 | static void
|
---|
| 276 | rotate (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 | }
|
---|