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

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

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