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