source: Sophya/trunk/SophyaExt/XephemAstroLib/plans.c@ 2900

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

Update de Xephem 3.7 cmv 21/08/2005
File size: 7.2 KB
Line 
1/* rewritten for Bureau des Longitude theories by Bretagnon and Chapront
2 * Michael Sternberg <sternberg@physik.tu-chemnitz.de>
3 */
4#include <stdio.h>
5#include <math.h>
6
7#include "astro.h"
8#include "vsop87.h"
9#include "chap95.h"
10
11static void pluto_ell (double mj, double *ret);
12static void chap_trans (double mj, double *ret);
13static void planpos (double mj, int obj, double prec, double *ret);
14
15/* coordinate transformation
16 * from:
17 * J2000.0 rectangular equatoreal ret[{0,1,2}] = {x,y,z}
18 * to:
19 * mean equinox of date spherical ecliptical ret[{0,1,2}] = {l,b,r}
20 */
21static void
22chap_trans (
23double mj, /* destination epoch */
24double *ret) /* vector to be transformed _IN PLACE_ */
25{
26 double ra, dec, r, eps;
27 double sr, cr, sd, cd, se, ce;
28
29 cartsph(ret[0], ret[1], ret[2], &ra, &dec, &r);
30 precess(J2000, mj, &ra, &dec);
31 obliquity(mj, &eps);
32 sr = sin(ra); cr = cos(ra);
33 sd = sin(dec); cd = cos(dec);
34 se = sin(eps); ce = cos(eps);
35 ret[0] = atan2( sr * ce + sd/cd * se, cr); /* long */
36 ret[1] = asin( sd * ce - cd * se * sr); /* lat */
37 ret[2] = r; /* radius */
38}
39
40/* low precision ecliptic coordinates of Pluto from mean orbit.
41 * Only for sake of completeness outside available perturbation theories.
42 */
43static void
44pluto_ell (
45double mj, /* epoch */
46double *ret) /* ecliptic coordinates {l,b,r} at equinox of date */
47{
48 /* mean orbital elements of Pluto.
49 * The origin of these is somewhat obscure.
50 */
51 double a = 39.543, /* semimajor axis, au */
52 e = 0.2490, /* excentricity */
53 inc0 = 17.140, /* inclination, deg */
54 Om0 = 110.307, /* long asc node, deg */
55 omeg0 = 113.768, /* arg of perihel, deg */
56 mjp = 2448045.539 - MJD0, /* epoch of perihel */
57 mjeq = J2000, /* equinox of elements */
58 n = 144.9600/36525.; /* daily motion, deg */
59
60 double inc, Om, omeg; /* orbital elements at epoch of date */
61 double ma, ea, nu; /* mean, excentric and true anomaly */
62 double lo, slo, clo; /* longitude in orbit from asc node */
63
64 reduce_elements(mjeq, mj, degrad(inc0), degrad(omeg0), degrad(Om0),
65 &inc, &omeg, &Om);
66 ma = degrad((mj - mjp) * n);
67 anomaly(ma, e, &nu, &ea);
68 ret[2] = a * (1.0 - e*cos(ea)); /* r */
69 lo = omeg + nu;
70 slo = sin(lo);
71 clo = cos(lo);
72 ret[1] = asin(slo * sin(inc)); /* b */
73 ret[0] = atan2(slo * cos(inc), clo) + Om; /* l */
74}
75
76/*************************************************************/
77
78/* geometric heliocentric position of planet, mean ecliptic of date
79 * (not corrected for light-time)
80 */
81static void
82planpos (double mj, int obj, double prec, double *ret)
83{
84 if (mj >= CHAP_BEGIN && mj <= CHAP_END) {
85 if (obj >= JUPITER) { /* prefer Chapront */
86 chap95(mj, obj, prec, ret);
87 chap_trans (mj, ret);
88 } else { /* VSOP for inner planets */
89 vsop87(mj, obj, prec, ret);
90 }
91 } else { /* outside Chapront time: */
92 if (obj != PLUTO) { /* VSOP for all but Pluto */
93 vsop87(mj, obj, prec, ret);
94 } else { /* Pluto mean elliptic orbit */
95 pluto_ell(mj, ret);
96 }
97 }
98}
99
100/*************************************************************/
101
102/* visual elements of planets
103 * [planet][0] = angular size at 1 AU
104 * [planet][1] = magnitude at 1 AU from sun and earth and 0 deg phase angle
105 * [planet][2] = A
106 * [planet][3] = B
107 * [planet][4] = C
108 * where mag correction = A*(i/100) + B*(i/100)^2 + C*(i/100)^3
109 * i = angle between sun and earth from planet, degrees
110 * from Explanatory Supplement, 1992
111 */
112static double vis_elements[8][5] = {
113 /* Mercury */ { 6.74, -0.36, 3.8, -2.73, 2.00},
114 /* Venus */ { 16.92, -4.29, 0.09, 2.39, -.65},
115 /* Mars */ { 9.36, -1.52, 1.60, 0., 0.},
116 /* Jupiter */ { 196.74, -9.25, 0.50, 0., 0.},
117 /* Saturn */ { 165.6, -8.88, 4.40, 0., 0.},
118 /* Uranus */ { 65.8, -7.19, 0.28, 0., 0.},
119 /* Neptune */ { 62.2, -6.87, 0., 0., 0.},
120 /* Pluto */ { 8.2, -1.01, 4.1, 0., 0.}
121};
122
123/* given a modified Julian date, mj, and a planet, p, find:
124 * lpd0: heliocentric longitude,
125 * psi0: heliocentric latitude,
126 * rp0: distance from the sun to the planet,
127 * rho0: distance from the Earth to the planet,
128 * none corrected for light time, ie, they are the true values for the
129 * given instant.
130 * lam: geocentric ecliptic longitude,
131 * bet: geocentric ecliptic latitude,
132 * each corrected for light time, ie, they are the apparent values as
133 * seen from the center of the Earth for the given instant.
134 * dia: angular diameter in arcsec at 1 AU,
135 * mag: visual magnitude
136 *
137 * all angles are in radians, all distances in AU.
138 *
139 * corrections for nutation and abberation must be made by the caller. The RA
140 * and DEC calculated from the fully-corrected ecliptic coordinates are then
141 * the apparent geocentric coordinates. Further corrections can be made, if
142 * required, for atmospheric refraction and geocentric parallax.
143 */
144void
145plans (double mj, PLCode p, double *lpd0, double *psi0, double *rp0,
146double *rho0, double *lam, double *bet, double *dia, double *mag)
147{
148 static double lastmj = -10000;
149 static double lsn, bsn, rsn; /* geocentric coords of sun */
150 static double xsn, ysn, zsn; /* cartesian " */
151 double lp, bp, rp; /* heliocentric coords of planet */
152 double xp, yp, zp, rho; /* rect. coords and geocentric dist. */
153 double dt; /* light time */
154 double *vp; /* vis_elements[p] */
155 double ci, i; /* sun/earth angle: cos, degrees */
156 int pass;
157
158 /* get sun cartesian; needed only once at mj */
159 if (mj != lastmj) {
160 sunpos (mj, &lsn, &rsn, &bsn);
161 sphcart (lsn, bsn, rsn, &xsn, &ysn, &zsn);
162 lastmj = mj;
163 }
164
165 /* first find the true position of the planet at mj.
166 * then repeat a second time for a slightly different time based
167 * on the position found in the first pass to account for light-travel
168 * time.
169 */
170 dt = 0.0;
171 for (pass = 0; pass < 2; pass++) {
172 double ret[6];
173
174 /* get spherical coordinates of planet from precision routines,
175 * retarded for light time in second pass;
176 * alternative option: vsop allows calculating rates.
177 */
178 planpos(mj - dt, p, 0.0, ret);
179
180 lp = ret[0];
181 bp = ret[1];
182 rp = ret[2];
183
184 sphcart (lp, bp, rp, &xp, &yp, &zp);
185 cartsph (xp + xsn, yp + ysn, zp + zsn, lam, bet, &rho);
186
187 if (pass == 0) {
188 /* save heliocentric coordinates at first pass since, being
189 * true, they are NOT to be corrected for light-travel time.
190 */
191 *lpd0 = lp;
192 range (lpd0, 2.*PI);
193 *psi0 = bp;
194 *rp0 = rp;
195 *rho0 = rho;
196 }
197
198 /* when we view a planet we see it in the position it occupied
199 * dt days ago, where rho is the distance between it and earth,
200 * in AU. use this as the new time for the next pass.
201 */
202 dt = rho * 5.7755183e-3;
203 }
204
205 vp = vis_elements[p];
206 *dia = vp[0];
207
208 /* solve plane triangle, assume sun/earth dist == 1 */
209 ci = (rp*rp + rho*rho - 1)/(2*rp*rho);
210
211 /* expl supp equation for mag */
212 if (ci < -1) ci = -1;
213 if (ci > 1) ci = 1;
214 i = raddeg(acos(ci))/100.;
215 *mag = vp[1] + 5*log10(rho*rp) + i*(vp[2] + i*(vp[3] + i*vp[4]));
216
217 /* rings contribution if SATURN */
218 if (p == SATURN) {
219 double et, st, set;
220 satrings (bp, lp, rp, lsn+PI, rsn, mj+MJD0, &et, &st);
221 set = sin(fabs(et));
222 *mag += (-2.60 + 1.25*set)*set;
223 }
224}
225
226/* For RCS Only -- Do Not Edit */
227static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: plans.c,v $ $Date: 2005-08-21 10:02:39 $ $Revision: 1.5 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.