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

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

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