source: Sophya/trunk/SophyaExt/XephemAstroLib/satmoon.c@ 3654

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

mise a niveau Xephem 3.7.4, cmv 16/07/2009
File size: 17.9 KB
RevLine 
[2551]1/* saturn moon info */
2
3#include <stdio.h>
4#include <stdlib.h>
5#include <string.h>
[3111]6#include <errno.h>
[2551]7#include <math.h>
8
9#include "astro.h"
10#include "bdl.h"
11
12static int use_bdl (double JD, char *dir, MoonData md[S_NMOONS]);
13static void bruton_saturn (Obj *sop, double JD, MoonData md[S_NMOONS]);
14static void moonradec (double satsize, MoonData md[S_NMOONS]);
15static void moonSVis (Obj *eop, Obj *sop, MoonData md[S_NMOONS]);
16static void moonEVis (MoonData md[S_NMOONS]);
[2643]17static void moonPShad (Obj *eop, Obj *sop, MoonData md[S_NMOONS]);
18static void moonTrans (MoonData md[S_NMOONS]);
[2551]19
20/* moon table and a few other goodies and when it was last computed */
21static double mdmjd = -123456;
22static MoonData smd[S_NMOONS] = {
23 {"Saturn", NULL},
24 {"Mimas", "I"},
25 {"Enceladus","II"},
26 {"Tethys", "III"},
27 {"Dione", "IV"},
28 {"Rhea", "V"},
29 {"Titan", "VI"},
30 {"Hyperion","VII"},
31 {"Iapetus", "VIII"},
32};
33static double sizemjd;
34static double etiltmjd;
35static double stiltmjd;
36
[2643]37/* These values are from the Explanatory Supplement.
38 * Precession degrades them gradually over time.
39 */
40#define POLE_RA degrad(40.58) /* RA of Saturn's north pole */
41#define POLE_DEC degrad(83.54) /* Dec of Saturn's north pole */
42
43
[2551]44/* get saturn info in md[0], moon info in md[1..S_NMOONS-1].
45 * if !dir always use bruton model.
46 * if !sop caller just wants md[] for names
47 * N.B. we assume eop and sop are updated.
48 */
49void
50saturn_data (
51double Mjd, /* mjd */
52char dir[], /* dir in which to look for helper files */
53Obj *eop, /* earth == Sun */
54Obj *sop, /* saturn */
55double *sizep, /* saturn's angular diam, rads */
56double *etiltp, double *stiltp, /* earth and sun tilts -- +S */
[2643]57double *polera, double *poledec,/* pole location */
[2551]58MoonData md[S_NMOONS]) /* return info */
59{
60 double JD;
61
62 /* always copy back at least for name */
63 memcpy (md, smd, sizeof(smd));
64
[2643]65 /* pole */
66 if (polera) *polera = POLE_RA;
67 if (poledec) *poledec = POLE_DEC;
68
[2551]69 /* nothing else if repeat call or just want names */
70 if (Mjd == mdmjd || !sop) {
71 if (sop) {
72 *sizep = sizemjd;
73 *etiltp = etiltmjd;
74 *stiltp = stiltmjd;
75 }
76 return;
77 }
78 JD = Mjd + MJD0;
79
80 /* planet in [0] */
81 md[0].ra = sop->s_ra;
82 md[0].dec = sop->s_dec;
83 md[0].mag = get_mag(sop);
84 md[0].x = 0;
85 md[0].y = 0;
86 md[0].z = 0;
87 md[0].evis = 1;
88 md[0].svis = 1;
89
90 /* size is straight from sop */
91 *sizep = degrad(sop->s_size/3600.0);
92
93 /* Visual Magnitude of the Satellites */
94
95 md[1].mag = 13; md[2].mag = 11.8; md[3].mag = 10.3; md[4].mag = 10.2;
96 md[5].mag = 9.8; md[6].mag = 8.4; md[7].mag = 14.3; md[8].mag = 11.2;
97
98 /* get tilts from sky and tel code */
99 satrings (sop->s_hlat, sop->s_hlong, sop->s_sdist, eop->s_hlong,
100 eop->s_edist, JD, etiltp, stiltp);
101
102 /* get moon x,y,z from BDL if possible, else Bruton's model */
[3477]103 if (!dir || use_bdl (JD, dir, md) < 0)
[2551]104 bruton_saturn (sop, JD, md);
105
106 /* set visibilities */
107 moonSVis (eop, sop, md);
[2643]108 moonPShad (eop, sop, md);
[2551]109 moonEVis (md);
[2643]110 moonTrans (md);
[2551]111
112 /* fill in moon ra and dec */
113 moonradec (*sizep, md);
114
115 /* save */
116 mdmjd = Mjd;
117 etiltmjd = *etiltp;
118 stiltmjd = *stiltp;
119 sizemjd = *sizep;
120 memcpy (smd, md, sizeof(smd));
121}
122
123/* hunt for BDL file in dir[] and use if possible.
124 * return 0 if ok, else -1
125 */
126static int
127use_bdl (
128double JD, /* julian date */
129char dir[], /* directory */
130MoonData md[S_NMOONS]) /* fill md[1..NM-1].x/y/z for each moon */
131{
132#define SATRAU .0004014253 /* saturn radius, AU */
133 double x[S_NMOONS], y[S_NMOONS], z[S_NMOONS];
134 char buf[1024];
135 FILE *fp;
[3111]136 char *fn;
[2551]137 int i;
138
[3111]139 /* check ranges and appropriate data file */
140 if (JD < 2451179.50000) /* Jan 1 1999 UTC */
[2551]141 return (-1);
[3111]142 if (JD < 2455562.5) /* Jan 1 2011 UTC */
143 fn = "saturne.9910";
144 else if (JD < 2459215.5) /* Jan 1 2021 UTC */
145 fn = "saturne.1020";
146 else
147 return (-1);
[2551]148
149 /* open */
[3111]150 (void) sprintf (buf, "%s/%s", dir, fn);
[2551]151 fp = fopen (buf, "r");
[3111]152 if (!fp) {
153 fprintf (stderr, "%s: %s\n", fn, strerror(errno));
[2551]154 return (-1);
[3111]155 }
[2551]156
157 /* use it */
158 if ((i = read_bdl (fp, JD, x, y, z, buf)) < 0) {
[3111]159 fprintf (stderr, "%s: %s\n", fn, buf);
[2551]160 fclose (fp);
161 return (-1);
162 }
163 if (i != S_NMOONS-1) {
[3111]164 fprintf (stderr, "%s: BDL says %d moons, code expects %d", fn,
[2551]165 i, S_NMOONS-1);
166 fclose (fp);
167 return (-1);
168 }
169
170 /* copy into md[1..NM-1] with our scale and sign conventions */
171 for (i = 1; i < S_NMOONS; i++) {
172 md[i].x = x[i-1]/SATRAU; /* we want sat radii +E */
173 md[i].y = -y[i-1]/SATRAU; /* we want sat radii +S */
174 md[i].z = -z[i-1]/SATRAU; /* we want sat radii +front */
175 }
176
177 /* ok */
178 fclose (fp);
179 return (0);
180}
181
182/* */
183/* SS2TXT.BAS Dan Bruton, astro@tamu.edu */
184/* */
185/* This is a text version of SATSAT2.BAS. It is smaller, */
186/* making it easier to convert other languages (250 lines */
187/* compared to 850 lines). */
188/* */
189/* This BASIC program computes and displays the locations */
190/* of Saturn's Satellites for a given date and time. See */
191/* "Practical Astronomy with your Calculator" by Peter */
192/* Duffett-Smith and the Astronomical Almanac for explanations */
193/* of some of the calculations here. The code is included so */
194/* that users can make changes or convert to other languages. */
195/* This code was made using QBASIC (comes with DOS 5.0). */
196/* */
197/* ECD: merged with Sky and Tel, below, for better earth and sun ring tilt */
198/* */
199
200/* ECD: BASICeze */
201#define FOR for
202#define IF if
203#define ELSE else
204#define COS cos
205#define SIN sin
206#define TAN tan
207#define ATN atan
208#define ABS fabs
209#define SQR sqrt
210
211/* find saturn moon data from Bruton's model */
212/* this originally computed +X:East +Y:North +Z:behind in [1..8] indeces.
213 * and +tilt:front south, rads
214 * then we adjust things in md[].x/y/z/mag to fit into our MoonData format.
215 */
216static void
217bruton_saturn (
218Obj *sop, /* saturn */
219double JD, /* julian date */
220MoonData md[S_NMOONS]) /* fill md[1..NM-1].x/y/z for each moon */
221{
222 /* ECD: code does not use [0].
223 * ECD and why 11 here? seems like 9 would do
224 */
225 double SMA[11], U[11], U0[11], PD[11];
226 double X[S_NMOONS], Y[S_NMOONS], Z[S_NMOONS];
227
228 double P,TP,TE,EP,EE,RE0,RP0,RS;
229 double JDE,LPE,LPP,LEE,LEP;
230 double NN,ME,MP,VE,VP;
231 double LE,LP,RE,RP,DT,II,F,F1;
232 double RA,DECL;
233 double TVA,PVA,TVC,PVC,DOT1,INC,TVB,PVB,DOT2,INCI;
234 double TRIP,GAM,TEMPX,TEMPY,TEMPZ;
235 int I;
236
237 /* saturn */
238 RA = sop->s_ra;
239 DECL = sop->s_dec;
240
241 /* ******************************************************************** */
242 /* * * */
243 /* * Constants * */
244 /* * * */
245 /* ******************************************************************** */
246 P = PI / 180;
247 /* Orbital Rate of Saturn in Radians per Days */
248 TP = 2 * PI / (29.45771 * 365.2422);
249 /* Orbital Rate of Earth in Radians per Day */
250 TE = 2 * PI / (1.00004 * 365.2422);
251 /* Eccentricity of Saturn's Orbit */
252 EP = .0556155;
253 /* Eccentricity of Earth's Orbit */
254 EE = .016718;
255 /* Semimajor axis of Earth's and Saturn's orbit in Astronomical Units */
256 RE0 = 1; RP0 = 9.554747;
257 /* Semimajor Axis of the Satellites' Orbit in Kilometers */
258 SMA[1] = 185600; SMA[2] = 238100; SMA[3] = 294700; SMA[4] = 377500;
259 SMA[5] = 527200; SMA[6] = 1221600; SMA[7] = 1483000; SMA[8] = 3560100;
260 /* Eccentricity of Satellites' Orbit [Program uses 0] */
261 /* Synodic Orbital Period of Moons in Days */
262 PD[1] = .9425049;
263 PD[2] = 1.3703731;
264 PD[3] = 1.8880926;
265 PD[4] = 2.7375218;
266 PD[5] = 4.5191631;
267 PD[6] = 15.9669028;
268 PD[7] = 21.3174647;
269 PD[8] = 79.9190206; /* personal mail 1/14/95 */
270 RS = 60330; /* Radius of planet in kilometers */
271
272 /* ******************************************************************** */
273 /* * * */
274 /* * Epoch Information * */
275 /* * * */
276 /* ******************************************************************** */
277 JDE = 2444238.5; /* Epoch Jan 0.0 1980 = December 31,1979 0:0:0 UT */
278 LPE = 165.322242 * P; /* Longitude of Saturn at Epoch */
279 LPP = 92.6653974 * P; /* Longitude of Saturn`s Perihelion */
280 LEE = 98.83354 * P; /* Longitude of Earth at Epoch */
281 LEP = 102.596403 * P; /* Longitude of Earth's Perihelion */
282 /* U0[I] = Angle from inferior geocentric conjuction */
283 /* measured westward along the orbit at epoch */
284 U0[1] = 18.2919 * P;
285 U0[2] = 174.2135 * P;
286 U0[3] = 172.8546 * P;
287 U0[4] = 76.8438 * P;
288 U0[5] = 37.2555 * P;
289 U0[6] = 57.7005 * P;
290 U0[7] = 266.6977 * P;
291 U0[8] = 195.3513 * P; /* from personal mail 1/14/1995 */
292
293 /* ******************************************************************** */
294 /* * * */
295 /* * Orbit Calculations * */
296 /* * * */
297 /* ******************************************************************** */
298 /* ****************** FIND MOON ORBITAL ANGLES ************************ */
299 NN = JD - JDE; /* NN = Number of days since epoch */
300 ME = ((TE * NN) + LEE - LEP); /* Mean Anomoly of Earth */
301 MP = ((TP * NN) + LPE - LPP); /* Mean Anomoly of Saturn */
302 VE = ME; VP = MP; /* True Anomolies - Solve Kepler's Equation */
303 FOR (I = 1; I <= 3; I++) {
304 VE = VE - (VE - (EE * SIN(VE)) - ME) / (1 - (EE * COS(VE)));
305 VP = VP - (VP - (EP * SIN(VP)) - MP) / (1 - (EP * COS(VP)));
306 }
307 VE = 2 * ATN(SQR((1 + EE) / (1 - EE)) * TAN(VE / 2));
308 IF (VE < 0) VE = (2 * PI) + VE;
309 VP = 2 * ATN(SQR((1 + EP) / (1 - EP)) * TAN(VP / 2));
310 IF (VP < 0) VP = (2 * PI) + VP;
311 /* Heliocentric Longitudes of Earth and Saturn */
312 LE = VE + LEP; IF (LE > (2 * PI)) LE = LE - (2 * PI);
313 LP = VP + LPP; IF (LP > (2 * PI)) LP = LP - (2 * PI);
314 /* Distances of Earth and Saturn from the Sun in AU's */
315 RE = RE0 * (1 - EE * EE) / (1 + EE * COS(VE));
316 RP = RP0 * (1 - EP * EP) / (1 + EP * COS(VP));
317 /* DT = Distance from Saturn to Earth in AU's - Law of Cosines */
318 DT = SQR((RE * RE) + (RP * RP) - (2 * RE * RP * COS(LE - LP)));
319 /* II = Angle between Earth and Sun as seen from Saturn */
320 II = RE * SIN(LE - LP) / DT;
321 II = ATN(II / SQR(1 - II * II)); /* ArcSIN and Law of Sines */
322 /* F = NN - (Light Time to Earth in days) */
323 F = NN - (DT / 173.83);
324 F1 = II + MP - VP;
325 /* U(I) = Angle from inferior geocentric conjuction measured westward */
326 FOR (I = 1; I < S_NMOONS; I++) {
327 U[I] = U0[I] + (F * 2 * PI / PD[I]) + F1;
328 U[I] = ((U[I] / (2 * PI)) - (int)(U[I] / (2 * PI))) * 2 * PI;
329
330 }
331
332 /* **************** FIND INCLINATION OF RINGS ************************* */
333 /* Use dot product of Earth-Saturn vector and Saturn's rotation axis */
334 TVA = (90 - 83.51) * P; /* Theta coordinate of Saturn's axis */
335 PVA = 40.27 * P; /* Phi coordinate of Saturn's axis */
336 TVC = (PI / 2) - DECL;
337 PVC = RA;
338 DOT1 = SIN(TVA) * COS(PVA) * SIN(TVC) * COS(PVC);
339 DOT1 = DOT1 + SIN(TVA) * SIN(PVA) * SIN(TVC) * SIN(PVC);
340 DOT1 = DOT1 + COS(TVA) * COS(TVC);
341 INC = ATN(SQR(1 - DOT1 * DOT1) / DOT1); /* ArcCOS */
342 IF (INC > 0) INC = (PI / 2) - INC; ELSE INC = -(PI / 2) - INC;
343
344 /* ************* FIND INCLINATION OF IAPETUS' ORBIT ******************* */
345 /* Use dot product of Earth-Saturn vector and Iapetus' orbit axis */
346 /* Vector B */
347 TVB = (90 - 75.6) * P; /* Theta coordinate of Iapetus' orbit axis (estimate) */
348 PVB = 21.34 * 2 * PI / 24; /* Phi coordinate of Iapetus' orbit axis (estimate) */
349 DOT2 = SIN(TVB) * COS(PVB) * SIN(TVC) * COS(PVC);
350 DOT2 = DOT2 + SIN(TVB) * SIN(PVB) * SIN(TVC) * SIN(PVC);
351 DOT2 = DOT2 + COS(TVB) * COS(TVC);
352 INCI = ATN(SQR(1 - DOT2 * DOT2) / DOT2); /* ArcCOS */
353 IF (INCI > 0) INCI = (PI / 2) - INCI; ELSE INCI = -(PI / 2) - INCI;
354
355 /* ************* FIND ROTATION ANGLE OF IAPETUS' ORBIT **************** */
356 /* Use inclination of Iapetus' orbit with respect to ring plane */
357 /* Triple Product */
358 TRIP = SIN(TVC) * COS(PVC) * SIN(TVA) * SIN(PVA) * COS(TVB);
359 TRIP = TRIP - SIN(TVC) * COS(PVC) * SIN(TVB) * SIN(PVB) * COS(TVA);
360 TRIP = TRIP + SIN(TVC) * SIN(PVC) * SIN(TVB) * COS(PVB) * COS(TVA);
361 TRIP = TRIP - SIN(TVC) * SIN(PVC) * SIN(TVA) * COS(PVA) * COS(TVB);
362 TRIP = TRIP + COS(TVC) * SIN(TVA) * COS(PVA) * SIN(TVB) * SIN(PVB);
363 TRIP = TRIP - COS(TVC) * SIN(TVB) * COS(PVB) * SIN(TVA) * SIN(PVA);
364 GAM = -1 * ATN(TRIP / SQR(1 - TRIP * TRIP)); /* ArcSIN */
365
366 /* ******************************************************************** */
367 /* * * */
368 /* * Compute Moon Positions * */
369 /* * * */
370 /* ******************************************************************** */
371 FOR (I = 1; I < S_NMOONS - 1; I++) {
372 X[I] = -1 * SMA[I] * SIN(U[I]) / RS;
373 Z[I] = -1 * SMA[I] * COS(U[I]) / RS; /* ECD */
374 Y[I] = SMA[I] * COS(U[I]) * SIN(INC) / RS;
375 }
376 /* ************************* Iapetus' Orbit *************************** */
377 TEMPX = -1 * SMA[8] * SIN(U[8]) / RS;
378 TEMPZ = -1 * SMA[8] * COS(U[8]) / RS;
379 TEMPY = SMA[8] * COS(U[8]) * SIN(INCI) / RS;
380 X[8] = TEMPX * COS(GAM) + TEMPY * SIN(GAM); /* Rotation */
381 Z[8] = TEMPZ * COS(GAM) + TEMPY * SIN(GAM);
382 Y[8] = -1 * TEMPX * SIN(GAM) + TEMPY * COS(GAM);
383
384#ifdef SHOWALL
385 /* ******************************************************************** */
386 /* * * */
387 /* * Show Results * */
388 /* * * */
389 /* ******************************************************************** */
390 printf (" Julian Date : %g\n", JD);
391 printf (" Right Ascension of Saturn : %g Hours\n", RA * 24 / (2 * PI));
392 printf (" Declination of Saturn : %g\n", DECL / P);
393 printf (" Ring Inclination as seen from Earth : %g\n", -1 * INC / P);
394 printf (" Heliocentric Longitude of Saturn : %g\n", LP / P);
395 printf (" Heliocentric Longitude of Earth : %g\n", LE / P);
396 printf (" Sun-Saturn-Earth Angle : %g\n", II / P);
397 printf (" Distance between Saturn and Earth : %g AU = %g million miles\n", DT, (DT * 93));
398 printf (" Light time from Saturn to Earth : %g minutes\n", DT * 8.28);
399 TEMP = 2 * ATN(TAN(165.6 * P / (2 * 3600)) / DT) * 3600 / P;
400 printf (" Angular Size of Saturn : %g arcsec\n", TEMP);
401 printf (" Major Angular Size of Saturn's Rings : %g arcsec\n", RS4 * TEMP / RS);
402 printf (" Minor Angular Size of Saturn's Rings : %g arcsec\n", ABS(RS4 * TEMP * SIN(INC) / RS));
403#endif
404
405 /* copy into md[1..S_NMOONS-1] with our sign conventions */
406 for (I = 1; I < S_NMOONS; I++) {
407 md[I].x = X[I]; /* we want +E */
408 md[I].y = -Y[I]; /* we want +S */
409 md[I].z = -Z[I]; /* we want +front */
410 }
411}
412
413/* given saturn loc in md[0].ra/dec and size, and location of each moon in
414 * md[1..NM-1].x/y in sat radii, find ra/dec of each moon in md[1..NM-1].ra/dec.
415 */
416static void
417moonradec (
418double satsize, /* sat diameter, rads */
419MoonData md[S_NMOONS]) /* fill in RA and Dec */
420{
421 double satrad = satsize/2;
422 double satra = md[0].ra;
423 double satdec = md[0].dec;
424 int i;
425
426 for (i = 1; i < S_NMOONS; i++) {
427 double dra = satrad * md[i].x;
428 double ddec = satrad * md[i].y;
429 md[i].ra = satra + dra;
430 md[i].dec = satdec - ddec;
431 }
432}
433
434/* set svis according to whether moon is in sun light */
435static void
436moonSVis(
437Obj *eop, /* earth == SUN */
438Obj *sop, /* saturn */
439MoonData md[S_NMOONS])
440{
441 double esd = eop->s_edist;
442 double eod = sop->s_edist;
443 double sod = sop->s_sdist;
444 double soa = degrad(sop->s_elong);
445 double esa = asin(esd*sin(soa)/sod);
446 double h = sod*sop->s_hlat;
447 double nod = h*(1./eod - 1./sod);
448 double sca = cos(esa), ssa = sin(esa);
449 int i;
450
451 for (i = 1; i < S_NMOONS; i++) {
452 MoonData *mdp = &md[i];
453 double xp = sca*mdp->x + ssa*mdp->z;
454 double yp = mdp->y;
455 double zp = -ssa*mdp->x + sca*mdp->z;
456 double ca = cos(nod), sa = sin(nod);
457 double xpp = xp;
458 double ypp = ca*yp - sa*zp;
459 double zpp = sa*yp + ca*zp;
460 int outside = xpp*xpp + ypp*ypp > 1.0;
461 int infront = zpp > 0.0;
462 mdp->svis = outside || infront;
463 }
464}
465
466/* set evis according to whether moon is geometrically visible from earth */
467static void
468moonEVis (MoonData md[S_NMOONS])
469{
470 int i;
471
472 for (i = 1; i < S_NMOONS; i++) {
473 MoonData *mdp = &md[i];
474 int outside = mdp->x*mdp->x + mdp->y*mdp->y > 1.0;
475 int infront = mdp->z > 0.0;
476 mdp->evis = outside || infront;
477 }
478}
479
[2643]480/* set pshad and sx,sy shadow info */
481static void
482moonPShad(
483Obj *eop, /* earth == SUN */
484Obj *sop, /* saturn */
485MoonData md[S_NMOONS])
486{
487 int i;
488
489 for (i = 1; i < S_NMOONS; i++) {
490 MoonData *mdp = &md[i];
491 mdp->pshad = !plshadow (sop, eop, POLE_RA, POLE_DEC, mdp->x,
492 mdp->y, mdp->z, &mdp->sx, &mdp->sy);
493 }
494}
495
496
497/* set whether moons are transiting */
498static void
499moonTrans (MoonData md[S_NMOONS])
500{
501 int i;
502
503 for (i = 1; i < S_NMOONS; i++) {
504 MoonData *mdp = &md[i];
505 mdp->trans = mdp->z > 0 && mdp->x*mdp->x + mdp->y*mdp->y < 1;
506 }
507}
508
[2551]509/* For RCS Only -- Do Not Edit */
[3654]510static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: satmoon.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.