source: trunk/FACT++/pal/palRdplan.c

Last change on this file was 18347, checked in by tbretz, 9 years ago
File size: 5.5 KB
Line 
1/*
2*+
3* Name:
4* palRdplan
5
6* Purpose:
7* Approximate topocentric apparent RA,Dec of a planet
8
9* Language:
10* Starlink ANSI C
11
12* Type of Module:
13* Library routine
14
15* Invocation:
16* void palRdplan( double date, int np, double elong, double phi,
17* double * ra, double * dec, double * diam );
18
19* Arguments:
20* date = double (Given)
21* MJD of observation (JD-2400000.5) in TDB. For all practical
22* purposes TT can be used instead of TDB, and for many applications
23* UT will do (except for the Moon).
24* np = int (Given)
25* Planet: 1 = Mercury
26* 2 = Venus
27* 3 = Moon
28* 4 = Mars
29* 5 = Jupiter
30* 6 = Saturn
31* 7 = Uranus
32* 8 = Neptune
33* else = Sun
34* elong = double (Given)
35* Observer's east longitude (radians)
36* phi = double (Given)
37* Observer's geodetic latitude (radians)
38* ra = double * (Returned)
39* RA (topocentric apparent, radians)
40* dec = double * (Returned)
41* Dec (topocentric apparent, radians)
42* diam = double * (Returned)
43* Angular diameter (equatorial, radians)
44
45* Description:
46* Approximate topocentric apparent RA,Dec of a planet, and its
47* angular diameter.
48
49* Authors:
50* PTW: Patrick T. Wallace
51* TIMJ: Tim Jenness (JAC, Hawaii)
52* {enter_new_authors_here}
53
54* Notes:
55* - Unlike with slaRdplan, Pluto is not supported.
56* - The longitude and latitude allow correction for geocentric
57* parallax. This is a major effect for the Moon, but in the
58* context of the limited accuracy of the present routine its
59* effect on planetary positions is small (negligible for the
60* outer planets). Geocentric positions can be generated by
61* appropriate use of the routines palDmoon and eraPlan94.
62
63* History:
64* 2012-03-07 (TIMJ):
65* Initial version, with some documentation from SLA/F.
66* Adapted with permission from the Fortran SLALIB library.
67* {enter_further_changes_here}
68
69* Copyright:
70* Copyright (C) 1997 Rutherford Appleton Laboratory
71* Copyright (C) 2012 Science and Technology Facilities Council.
72* All Rights Reserved.
73
74* Licence:
75* This program is free software; you can redistribute it and/or
76* modify it under the terms of the GNU General Public License as
77* published by the Free Software Foundation; either version 3 of
78* the License, or (at your option) any later version.
79*
80* This program is distributed in the hope that it will be
81* useful, but WITHOUT ANY WARRANTY; without even the implied
82* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
83* PURPOSE. See the GNU General Public License for more details.
84*
85* You should have received a copy of the GNU General Public License
86* along with this program; if not, write to the Free Software
87* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
88* MA 02110-1301, USA.
89
90* Bugs:
91* {note_any_bugs_here}
92*-
93*/
94
95#include <math.h>
96
97#include "pal.h"
98#include "palmac.h"
99#include "pal1sofa.h"
100
101void palRdplan( double date, int np, double elong, double phi,
102 double * ra, double * dec, double * diam ) {
103
104 /* AU in km */
105 const double AUKM = 1.49597870e8;
106
107 /* Equatorial radii (km) */
108 const double EQRAU[] = {
109 696000.0, /* Sun */
110 2439.7,
111 6051.9,
112 1738,
113 3397,
114 71492,
115 60268,
116 25559,
117 24764
118 };
119
120 /* Local variables */
121 int i, j;
122 double stl;
123 double vgm[6];
124 double v[6];
125 double rmat[3][3];
126 double vse[6];
127 double vsg[6];
128 double vsp[6];
129 double vgo[6];
130 double dx,dy,dz,r,tl;
131
132 /* Classify np */
133 if (np < 0 || np > 8 ) np=0; /* Sun */
134
135 /* Approximate local sidereal time */
136 stl = palGmst( date - palDt( palEpj(date)) / 86400.0) + elong;
137
138 /* Geocentre to Moon (mean of date) */
139 palDmoon( date, v );
140
141 /* Nutation to true of date */
142 palNut( date, rmat );
143 eraRxp( rmat, v, vgm );
144 eraRxp( rmat, &(v[3]), &(vgm[3]) );
145
146 /* Moon? */
147 if (np == 3) {
148
149 /* geocentre to Moon (true of date) */
150 for (i=0; i<6; i++) {
151 v[i] = vgm[i];
152 }
153
154 } else {
155
156 /* Not moon: precession/nutation matrix J2000 to date */
157 palPrenut( 2000.0, date, rmat );
158
159 /* Sun to Earth-Moon Barycentre (J2000) */
160 palPlanet( date, 3, v, &j );
161
162 /* Precession and nutation to date */
163 eraRxp( rmat, v, vse );
164 eraRxp( rmat, &(v[3]), &(vse[3]) );
165
166 /* Sun to geocentre (true of date) */
167 for (i=0; i<6; i++) {
168 vsg[i] = vse[i] - 0.012150581 * vgm[i];
169 }
170
171 /* Sun ? */
172 if (np == 0) {
173
174 /* Geocentre to Sun */
175 for (i=0; i<6; i++) {
176 v[i] = -vsg[i];
177 }
178
179 } else {
180
181 /* Sun to Planet (J2000) */
182 palPlanet( date, np, v, &j );
183
184 /* Precession and nutation to date */
185 eraRxp( rmat, v, vsp );
186 eraRxp( rmat, &(v[3]), &(vsp[3]) );
187
188 /* Geocentre to planet */
189 for (i=0; i<6; i++) {
190 v[i] = vsp[i] - vsg[i];
191 }
192
193 }
194
195 }
196
197 /* Refer to origina at the observer */
198 palPvobs( phi, 0.0, stl, vgo );
199 for (i=0; i<6; i++) {
200 v[i] -= vgo[i];
201 }
202
203 /* Geometric distance (AU) */
204 dx = v[0];
205 dy = v[1];
206 dz = v[2];
207 r = sqrt( dx*dx + dy*dy + dz*dz );
208
209 /* Light time */
210 tl = PAL__CR * r;
211
212 /* Correct position for planetary aberration */
213 for (i=0; i<3; i++) {
214 v[i] -= tl * v[i+3];
215 }
216
217 /* To RA,Dec */
218 eraC2s( v, ra, dec );
219 *ra = eraAnp( *ra );
220
221 /* Angular diameter (radians) */
222 *diam = 2.0 * asin( EQRAU[np] / (r * AUKM ) );
223
224}
Note: See TracBrowser for help on using the repository browser.