File Coverage

palsrc/palRdplan.c
Criterion Covered Total %
statement 39 39 100.0
branch 17 18 94.4
condition n/a
subroutine n/a
pod n/a
total 56 57 98.2


line stmt bran cond sub pod time code
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
96              
97             #include "pal.h"
98             #include "palmac.h"
99             #include "pal1sofa.h"
100              
101 9           void 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 9           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 9 50         if (np < 0 || np > 8 ) np=0; /* Sun */
134              
135             /* Approximate local sidereal time */
136 9           stl = palGmst( date - palDt( palEpj(date)) / 86400.0) + elong;
137              
138             /* Geocentre to Moon (mean of date) */
139 9           palDmoon( date, v );
140              
141             /* Nutation to true of date */
142 9           palNut( date, rmat );
143 9           eraRxp( rmat, v, vgm );
144 9           eraRxp( rmat, &(v[3]), &(vgm[3]) );
145              
146             /* Moon? */
147 9 100         if (np == 3) {
148              
149             /* geocentre to Moon (true of date) */
150 7 100         for (i=0; i<6; i++) {
151 6           v[i] = vgm[i];
152             }
153              
154             } else {
155              
156             /* Not moon: precession/nutation matrix J2000 to date */
157 8           palPrenut( 2000.0, date, rmat );
158              
159             /* Sun to Earth-Moon Barycentre (J2000) */
160 8           palPlanet( date, 3, v, &j );
161              
162             /* Precession and nutation to date */
163 8           eraRxp( rmat, v, vse );
164 8           eraRxp( rmat, &(v[3]), &(vse[3]) );
165              
166             /* Sun to geocentre (true of date) */
167 56 100         for (i=0; i<6; i++) {
168 48           vsg[i] = vse[i] - 0.012150581 * vgm[i];
169             }
170              
171             /* Sun ? */
172 8 100         if (np == 0) {
173              
174             /* Geocentre to Sun */
175 7 100         for (i=0; i<6; i++) {
176 6           v[i] = -vsg[i];
177             }
178              
179             } else {
180              
181             /* Sun to Planet (J2000) */
182 7           palPlanet( date, np, v, &j );
183              
184             /* Precession and nutation to date */
185 7           eraRxp( rmat, v, vsp );
186 7           eraRxp( rmat, &(v[3]), &(vsp[3]) );
187              
188             /* Geocentre to planet */
189 49 100         for (i=0; i<6; i++) {
190 42           v[i] = vsp[i] - vsg[i];
191             }
192              
193             }
194              
195             }
196              
197             /* Refer to origina at the observer */
198 9           palPvobs( phi, 0.0, stl, vgo );
199 63 100         for (i=0; i<6; i++) {
200 54           v[i] -= vgo[i];
201             }
202              
203             /* Geometric distance (AU) */
204 9           dx = v[0];
205 9           dy = v[1];
206 9           dz = v[2];
207 9           r = sqrt( dx*dx + dy*dy + dz*dz );
208              
209             /* Light time */
210 9           tl = PAL__CR * r;
211              
212             /* Correct position for planetary aberration */
213 36 100         for (i=0; i<3; i++) {
214 27           v[i] -= tl * v[i+3];
215             }
216              
217             /* To RA,Dec */
218 9           eraC2s( v, ra, dec );
219 9           *ra = eraAnp( *ra );
220              
221             /* Angular diameter (radians) */
222 9           *diam = 2.0 * asin( EQRAU[np] / (r * AUKM ) );
223              
224 9           }