source: Sophya/trunk/SophyaExt/XephemAstroLib/marsmoon.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: 5.2 KB
Line 
1/* mars 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[M_NMOONS]);
12static void moonradec (double msize, MoonData md[M_NMOONS]);
13static void moonSVis (Obj *eop, Obj *mop, MoonData md[M_NMOONS]);
14static void moonEVis (MoonData md[M_NMOONS]);
15
16/* moon table and a few other goodies and when it was last computed */
17static double mdmjd = -123456;
18static MoonData mmd[M_NMOONS] = {
19 {"Mars", NULL},
20 {"Phobos", "I"},
21 {"Deimos", "II"},
22};
23static double sizemjd;
24
25/* file containing BDL coefficients */
26static char mbdlfn[] = "mars.9910";
27
28/* get mars info in md[0], moon info in md[1..M_NMOONS-1].
29 * if !dir always use bruton model.
30 * if !mop caller just wants md[] for names
31 * N.B. we assume eop and mop are updated.
32 */
33void
34marsm_data (
35double Mjd, /* mjd */
36char dir[], /* dir in which to look for helper files */
37Obj *eop, /* earth == Sun */
38Obj *mop, /* mars */
39double *sizep, /* mars's angular diam, rads */
40MoonData md[M_NMOONS]) /* return info */
41{
42 double JD, dmag;
43
44 /* always copy back at least for name */
45 memcpy (md, mmd, sizeof(mmd));
46
47 /* nothing else if repeat call or just want names */
48 if (Mjd == mdmjd || !mop) {
49 if (mop) {
50 *sizep = sizemjd;
51 }
52 return;
53 }
54 JD = Mjd + MJD0;
55
56 /* planet in [0] */
57 md[0].ra = mop->s_ra;
58 md[0].dec = mop->s_dec;
59 md[0].mag = get_mag(mop);
60 md[0].x = 0;
61 md[0].y = 0;
62 md[0].z = 0;
63 md[0].evis = 1;
64 md[0].svis = 1;
65
66 /* size is straight from mop */
67 *sizep = degrad(mop->s_size/3600.0);
68
69 /* from Pasachoff/Menzel: brightest @ .6 AU */
70 dmag = 5.0*log10(mop->s_edist + 0.4);
71 md[1].mag = 11.8 + dmag;
72 md[2].mag = 12.9 + dmag;
73
74 /* get moon x,y,z from BDL if possible */
75 if (dir && use_bdl (JD, dir, md) < 0) {
76 int i;
77 for (i = 1; i < M_NMOONS; i++)
78 md[i].x = md[i].y = md[i].z = 0.0;
79 fprintf (stderr, "No mars model available\n");
80 }
81
82 /* set visibilities */
83 moonSVis (eop, mop, md);
84 moonEVis (md);
85
86 /* fill in moon ra and dec */
87 moonradec (*sizep, md);
88
89 /* save */
90 mdmjd = Mjd;
91 sizemjd = *sizep;
92 memcpy (mmd, md, sizeof(mmd));
93}
94
95/* hunt for BDL file in dir[] and use if possible.
96 * return 0 if ok, else -1
97 */
98static int
99use_bdl (
100double JD, /* julian date */
101char dir[], /* directory */
102MoonData md[M_NMOONS]) /* fill md[1..NM-1].x/y/z for each moon */
103{
104#define MRAU .00002269 /* Mars radius, AU */
105 double x[M_NMOONS], y[M_NMOONS], z[M_NMOONS];
106 char buf[1024];
107 FILE *fp;
108 int i;
109
110 /* only valid 1999 through 2010 */
111 if (JD < 2451179.50000 || JD >= 2455562.5)
112 return (-1);
113
114 /* open */
115 (void) sprintf (buf, "%s/%s", dir, mbdlfn);
116 fp = fopen (buf, "r");
117 if (!fp)
118 return (-1);
119
120 /* use it */
121 if ((i = read_bdl (fp, JD, x, y, z, buf)) < 0) {
122 fprintf (stderr, "%s: %s\n", mbdlfn, buf);
123 fclose (fp);
124 return (-1);
125 }
126 if (i != M_NMOONS-1) {
127 fprintf (stderr, "%s: BDL says %d moons, code expects %d", mbdlfn,
128 i, M_NMOONS-1);
129 fclose (fp);
130 return (-1);
131 }
132
133 /* copy into md[1..NM-1] with our scale and sign conventions */
134 for (i = 1; i < M_NMOONS; i++) {
135 md[i].x = x[i-1]/MRAU; /* we want mars radii +E */
136 md[i].y = -y[i-1]/MRAU; /* we want mars radii +S */
137 md[i].z = -z[i-1]/MRAU; /* we want mars radii +front */
138 }
139
140 /* ok */
141 fclose (fp);
142 return (0);
143}
144
145/* given mars loc in md[0].ra/dec and size, and location of each moon in
146 * md[1..NM-1].x/y in mars radii,find ra/dec of each moon in md[1..NM-1].ra/dec.
147 */
148static void
149moonradec (
150double msize, /* mars diameter, rads */
151MoonData md[M_NMOONS]) /* fill in RA and Dec */
152{
153 double mrad = msize/2;
154 double mra = md[0].ra;
155 double mdec = md[0].dec;
156 int i;
157
158 for (i = 1; i < M_NMOONS; i++) {
159 double dra = mrad * md[i].x;
160 double ddec = mrad * md[i].y;
161 md[i].ra = mra + dra;
162 md[i].dec = mdec - ddec;
163 }
164}
165
166/* set svis according to whether moon is in sun light */
167static void
168moonSVis(
169Obj *eop, /* earth == SUN */
170Obj *mop, /* mars */
171MoonData md[M_NMOONS])
172{
173 double esd = eop->s_edist;
174 double eod = mop->s_edist;
175 double sod = mop->s_sdist;
176 double soa = degrad(mop->s_elong);
177 double esa = asin(esd*sin(soa)/sod);
178 double h = sod*mop->s_hlat;
179 double nod = h*(1./eod - 1./sod);
180 double sca = cos(esa), ssa = sin(esa);
181 int i;
182
183 for (i = 1; i < M_NMOONS; i++) {
184 MoonData *mdp = &md[i];
185 double xp = sca*mdp->x + ssa*mdp->z;
186 double yp = mdp->y;
187 double zp = -ssa*mdp->x + sca*mdp->z;
188 double ca = cos(nod), sa = sin(nod);
189 double xpp = xp;
190 double ypp = ca*yp - sa*zp;
191 double zpp = sa*yp + ca*zp;
192 int outside = xpp*xpp + ypp*ypp > 1.0;
193 int infront = zpp > 0.0;
194 mdp->svis = outside || infront;
195 }
196}
197
198/* set evis according to whether moon is geometrically visible from earth */
199static void
200moonEVis (MoonData md[M_NMOONS])
201{
202 int i;
203
204 for (i = 1; i < M_NMOONS; i++) {
205 MoonData *mdp = &md[i];
206 int outside = mdp->x*mdp->x + mdp->y*mdp->y > 1.0;
207 int infront = mdp->z > 0.0;
208 mdp->evis = outside || infront;
209 }
210}
211
212/* For RCS Only -- Do Not Edit */
213static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: marsmoon.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.