FACT++  1.0
palPlantu.c
Go to the documentation of this file.
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 
static const double PAL__CR
Definition: palmac.h:99
void palPrenut(double epoch, double date, double rmatpn[3][3])
Definition: palPrenut.c:70
double palEpj(double date)
Definition: palOne2One.c:1169
int i
Definition: db_dim_client.c:21
double palDt(double epoch)
Definition: palDt.c:97
void palEpv(double date, double ph[3], double vh[3], double pb[3], double vb[3])
Definition: palEpv.c:79
void eraC2s(double p[3], double *theta, double *phi)
Definition: c2s.c:3
void eraRxp(double r[3][3], double p[3], double rp[3])
Definition: rxp.c:3
void palPlantu(double date, double elong, double phi, const double u[13], double *ra, double *dec, double *r, int *jstat)
Definition: palPlantu.c:126
void palPvobs(double p, double h, double stl, double pv[6])
Definition: palPvobs.c:77
double eraAnp(double a)
Definition: anp.c:3
static const double PAL__SPD
Definition: palmac.h:102
double palGmst(double ut1)
Definition: palOne2One.c:1332
void palUe2pv(double date, double u[13], double pv[], int *jstat)