source: Sophya/trunk/SophyaExt/XephemAstroLib/umoon.c@ 2641

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