| 1 | /*
|
|---|
| 2 | *+
|
|---|
| 3 | * Name:
|
|---|
| 4 | * palPlantu
|
|---|
| 5 |
|
|---|
| 6 | * Purpose:
|
|---|
| 7 | * Topocentric RA,Dec of a Solar-System object from universal elements
|
|---|
| 8 |
|
|---|
| 9 | * Language:
|
|---|
| 10 | * Starlink ANSI C
|
|---|
| 11 |
|
|---|
| 12 | * Type of Module:
|
|---|
| 13 | * Library routine
|
|---|
| 14 |
|
|---|
| 15 | * Invocation:
|
|---|
| 16 | * void palPlantu ( double date, double elong, double phi, const double u[13],
|
|---|
| 17 | * double *ra, double *dec, double *r, int *jstat ) {
|
|---|
| 18 |
|
|---|
| 19 | * Description:
|
|---|
| 20 | * Topocentric apparent RA,Dec of a Solar-System object whose
|
|---|
| 21 | * heliocentric universal elements are known.
|
|---|
| 22 |
|
|---|
| 23 | * Arguments:
|
|---|
| 24 | * date = double (Given)
|
|---|
| 25 | * TT MJD of observation (JD-2400000.5)
|
|---|
| 26 | * elong = double (Given)
|
|---|
| 27 | * Observer's east longitude (radians)
|
|---|
| 28 | * phi = double (Given)
|
|---|
| 29 | * Observer's geodetic latitude (radians)
|
|---|
| 30 | * u = const double [13] (Given)
|
|---|
| 31 | * Universal orbital elements
|
|---|
| 32 | * - (0) combined mass (M+m)
|
|---|
| 33 | * - (1) total energy of the orbit (alpha)
|
|---|
| 34 | * - (2) reference (osculating) epoch (t0)
|
|---|
| 35 | * - (3-5) position at reference epoch (r0)
|
|---|
| 36 | * - (6-8) velocity at reference epoch (v0)
|
|---|
| 37 | * - (9) heliocentric distance at reference epoch
|
|---|
| 38 | * - (10) r0.v0
|
|---|
| 39 | * - (11) date (t)
|
|---|
| 40 | * - (12) universal eccentric anomaly (psi) of date, approx
|
|---|
| 41 | * ra = double * (Returned)
|
|---|
| 42 | * Topocentric apparent RA (radians)
|
|---|
| 43 | * dec = double * (Returned)
|
|---|
| 44 | * Topocentric apparent Dec (radians)
|
|---|
| 45 | * r = double * (Returned)
|
|---|
| 46 | * Distance from observer (AU)
|
|---|
| 47 | * jstat = int * (Returned)
|
|---|
| 48 | * status: 0 = OK
|
|---|
| 49 | * - -1 = radius vector zero
|
|---|
| 50 | * - -2 = failed to converge
|
|---|
| 51 |
|
|---|
| 52 | * Authors:
|
|---|
| 53 | * PTW: Pat Wallace (STFC)
|
|---|
| 54 | * TIMJ: Tim Jenness (JAC, Hawaii)
|
|---|
| 55 | * {enter_new_authors_here}
|
|---|
| 56 |
|
|---|
| 57 | * Notes:
|
|---|
| 58 | * - DATE is the instant for which the prediction is required. It is
|
|---|
| 59 | * in the TT timescale (formerly Ephemeris Time, ET) and is a
|
|---|
| 60 | * Modified Julian Date (JD-2400000.5).
|
|---|
| 61 | * - The longitude and latitude allow correction for geocentric
|
|---|
| 62 | * parallax. This is usually a small effect, but can become
|
|---|
| 63 | * important for near-Earth asteroids. Geocentric positions can be
|
|---|
| 64 | * generated by appropriate use of routines palEpv (or palEvp) and
|
|---|
| 65 | * palUe2pv.
|
|---|
| 66 | * - The "universal" elements are those which define the orbit for the
|
|---|
| 67 | * purposes of the method of universal variables (see reference 2).
|
|---|
| 68 | * They consist of the combined mass of the two bodies, an epoch,
|
|---|
| 69 | * and the position and velocity vectors (arbitrary reference frame)
|
|---|
| 70 | * at that epoch. The parameter set used here includes also various
|
|---|
| 71 | * quantities that can, in fact, be derived from the other
|
|---|
| 72 | * information. This approach is taken to avoiding unnecessary
|
|---|
| 73 | * computation and loss of accuracy. The supplementary quantities
|
|---|
| 74 | * are (i) alpha, which is proportional to the total energy of the
|
|---|
| 75 | * orbit, (ii) the heliocentric distance at epoch, (iii) the
|
|---|
| 76 | * outwards component of the velocity at the given epoch, (iv) an
|
|---|
| 77 | * estimate of psi, the "universal eccentric anomaly" at a given
|
|---|
| 78 | * date and (v) that date.
|
|---|
| 79 | * - The universal elements are with respect to the J2000 equator and
|
|---|
| 80 | * equinox.
|
|---|
| 81 |
|
|---|
| 82 | * See Also:
|
|---|
| 83 | * - Sterne, Theodore E., "An Introduction to Celestial Mechanics",
|
|---|
| 84 | * Interscience Publishers Inc., 1960. Section 6.7, p199.
|
|---|
| 85 | * - Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
|
|---|
| 86 |
|
|---|
| 87 | * History:
|
|---|
| 88 | * 2012-03-12 (TIMJ):
|
|---|
| 89 | * Initial version direct conversion of SLA/F.
|
|---|
| 90 | * Adapted with permission from the Fortran SLALIB library.
|
|---|
| 91 | * {enter_further_changes_here}
|
|---|
| 92 |
|
|---|
| 93 | * Copyright:
|
|---|
| 94 | * Copyright (C) 2005 Patrick T. Wallace
|
|---|
| 95 | * Copyright (C) 2012 Science and Technology Facilities Council.
|
|---|
| 96 | * All Rights Reserved.
|
|---|
| 97 |
|
|---|
| 98 | * Licence:
|
|---|
| 99 | * This program is free software; you can redistribute it and/or
|
|---|
| 100 | * modify it under the terms of the GNU General Public License as
|
|---|
| 101 | * published by the Free Software Foundation; either version 3 of
|
|---|
| 102 | * the License, or (at your option) any later version.
|
|---|
| 103 | *
|
|---|
| 104 | * This program is distributed in the hope that it will be
|
|---|
| 105 | * useful, but WITHOUT ANY WARRANTY; without even the implied
|
|---|
| 106 | * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|---|
| 107 | * PURPOSE. See the GNU General Public License for more details.
|
|---|
| 108 | *
|
|---|
| 109 | * You should have received a copy of the GNU General Public License
|
|---|
| 110 | * along with this program; if not, write to the Free Software
|
|---|
| 111 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|---|
| 112 | * MA 02110-1301, USA.
|
|---|
| 113 |
|
|---|
| 114 | * Bugs:
|
|---|
| 115 | * {note_any_bugs_here}
|
|---|
| 116 | *-
|
|---|
| 117 | */
|
|---|
| 118 |
|
|---|
| 119 | #include <math.h>
|
|---|
| 120 |
|
|---|
| 121 | #include "pal.h"
|
|---|
| 122 | #include "palmac.h"
|
|---|
| 123 |
|
|---|
| 124 | #include "pal1sofa.h"
|
|---|
| 125 |
|
|---|
| 126 | void palPlantu ( double date, double elong, double phi, const double u[13],
|
|---|
| 127 | double *ra, double *dec, double *r, int *jstat ) {
|
|---|
| 128 |
|
|---|
| 129 | int i;
|
|---|
| 130 | double dvb[3], dpb[3], vsg[6], vsp[6], v[6], rmat[3][3],
|
|---|
| 131 | vgp[6], stl, vgo[6], dx, dy, dz, d, tl;
|
|---|
| 132 |
|
|---|
| 133 | double ucp[13];
|
|---|
| 134 |
|
|---|
| 135 | /* To retain the stated const API and conform to the documentation
|
|---|
| 136 | we must copy the contents of the u array as palUe2pv updates
|
|---|
| 137 | the final two elements */
|
|---|
| 138 | for (i=0;i<13;i++) {
|
|---|
| 139 | ucp[i] = u[i];
|
|---|
| 140 | }
|
|---|
| 141 |
|
|---|
| 142 | /* Sun to geocentre (J2000, velocity in AU/s) */
|
|---|
| 143 | palEpv( date, vsg, &(vsg[3]), dpb, dvb );
|
|---|
| 144 | for (i=3; i < 6; i++) {
|
|---|
| 145 | vsg[i] /= PAL__SPD;
|
|---|
| 146 | }
|
|---|
| 147 |
|
|---|
| 148 | /* Sun to planet (J2000) */
|
|---|
| 149 | palUe2pv( date, ucp, vsp, jstat );
|
|---|
| 150 |
|
|---|
| 151 | /* Geocentre to planet (J2000) */
|
|---|
| 152 | for (i=0; i<6; i++) {
|
|---|
| 153 | v[i] = vsp[i] - vsg[i];
|
|---|
| 154 | }
|
|---|
| 155 |
|
|---|
| 156 | /* Precession and nutation to date */
|
|---|
| 157 | palPrenut( 2000.0, date, rmat );
|
|---|
| 158 | eraRxp(rmat, v, vgp);
|
|---|
| 159 | eraRxp( rmat, &(v[3]), &(vgp[3]) );
|
|---|
| 160 |
|
|---|
| 161 | /* Geocentre to observer (date) */
|
|---|
| 162 | stl = palGmst( date - palDt( palEpj(date) ) / PAL__SPD ) + elong;
|
|---|
| 163 | palPvobs( phi, 0.0, stl, vgo );
|
|---|
| 164 |
|
|---|
| 165 | /* Observer to planet (date) */
|
|---|
| 166 | for (i=0; i<6; i++) {
|
|---|
| 167 | v[i] = vgp[i] - vgo[i];
|
|---|
| 168 | }
|
|---|
| 169 |
|
|---|
| 170 | /* Geometric distance (AU) */
|
|---|
| 171 | dx = v[0];
|
|---|
| 172 | dy = v[1];
|
|---|
| 173 | dz = v[2];
|
|---|
| 174 | d = sqrt( dx*dx + dy*dy + dz*dz );
|
|---|
| 175 |
|
|---|
| 176 | /* Light time (sec) */
|
|---|
| 177 | tl = PAL__CR * d;
|
|---|
| 178 |
|
|---|
| 179 | /* Correct position for planetary aberration */
|
|---|
| 180 | for (i=0; i<3; i++) {
|
|---|
| 181 | v[i] -= tl * v[i+3];
|
|---|
| 182 | }
|
|---|
| 183 |
|
|---|
| 184 | /* To RA,Dec */
|
|---|
| 185 | eraC2s( v, ra, dec );
|
|---|
| 186 | *ra = eraAnp( *ra );
|
|---|
| 187 | *r = d;
|
|---|
| 188 | }
|
|---|
| 189 |
|
|---|
