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
palPlantu(double date,double elong,double phi,const double u[13],double * ra,double * dec,double * r,int * jstat)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
190