File Coverage

palsrc/palMapqk.c
Criterion Covered Total %
statement 30 30 100.0
branch 9 10 90.0
condition n/a
subroutine n/a
pod n/a
total 39 40 97.5


line stmt bran cond sub pod time code
1             /*
2             *+
3             * Name:
4             * palMapqk
5              
6             * Purpose:
7             * Quick mean to apparent place
8              
9             * Language:
10             * Starlink ANSI C
11              
12             * Type of Module:
13             * Library routine
14              
15             * Invocation:
16             * void palMapqk ( double rm, double dm, double pr, double pd,
17             * double px, double rv, double amprms[21],
18             * double *ra, double *da );
19              
20             * Arguments:
21             * rm = double (Given)
22             * Mean RA (radians)
23             * dm = double (Given)
24             * Mean declination (radians)
25             * pr = double (Given)
26             * RA proper motion, changes per Julian year (radians)
27             * pd = double (Given)
28             * Dec proper motion, changes per Julian year (radians)
29             * px = double (Given)
30             * Parallax (arcsec)
31             * rv = double (Given)
32             * Radial velocity (km/s, +ve if receding)
33             * amprms = double [21] (Given)
34             * Star-independent mean-to-apparent parameters (see palMappa).
35             * ra = double * (Returned)
36             * Apparent RA (radians)
37             * dec = double * (Returned)
38             * Apparent dec (radians)
39              
40             * Description:
41             * Quick mean to apparent place: transform a star RA,Dec from
42             * mean place to geocentric apparent place, given the
43             * star-independent parameters.
44             *
45             * Use of this routine is appropriate when efficiency is important
46             * and where many star positions, all referred to the same equator
47             * and equinox, are to be transformed for one epoch. The
48             * star-independent parameters can be obtained by calling the
49             * palMappa routine.
50             *
51             * If the parallax and proper motions are zero the palMapqkz
52             * routine can be used instead.
53              
54             * Notes:
55             * - The reference frames and timescales used are post IAU 2006.
56             * - The mean place rm, dm and the vectors amprms[1-3] and amprms[4-6]
57             * are referred to the mean equinox and equator of the epoch
58             * specified when generating the precession/nutation matrix
59             * amprms[12-20]. In the call to palMappa (q.v.) normally used
60             * to populate amprms, this epoch is the first argument (eq).
61             * - Strictly speaking, the routine is not valid for solar-system
62             * sources, though the error will usually be extremely small.
63             * However, to prevent gross errors in the case where the
64             * position of the Sun is specified, the gravitational
65             * deflection term is restrained within about 920 arcsec of the
66             * centre of the Sun's disc. The term has a maximum value of
67             * about 1.85 arcsec at this radius, and decreases to zero as
68             * the centre of the disc is approached.
69              
70             * Authors:
71             * PTW: Patrick T. Wallace
72             * TIMJ: Tim Jenness (JAC, Hawaii)
73             * {enter_new_authors_here}
74              
75             * History:
76             * 2012-03-01 (TIMJ):
77             * Initial version with documentation from SLA/F
78             * Adapted with permission from the Fortran SLALIB library.
79             * {enter_further_changes_here}
80              
81             * Copyright:
82             * Copyright (C) 2000 Rutherford Appleton Laboratory
83             * Copyright (C) 2012 Science and Technology Facilities Council.
84             * All Rights Reserved.
85              
86             * Licence:
87             * This program is free software; you can redistribute it and/or
88             * modify it under the terms of the GNU General Public License as
89             * published by the Free Software Foundation; either version 3 of
90             * the License, or (at your option) any later version.
91             *
92             * This program is distributed in the hope that it will be
93             * useful, but WITHOUT ANY WARRANTY; without even the implied
94             * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
95             * PURPOSE. See the GNU General Public License for more details.
96             *
97             * You should have received a copy of the GNU General Public License
98             * along with this program; if not, write to the Free Software
99             * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
100             * MA 02110-1301, USA.
101              
102             * Bugs:
103             * {note_any_bugs_here}
104             *-
105             */
106              
107             #include "pal.h"
108             #include "palmac.h"
109             #include "pal1sofa.h"
110              
111 1           void palMapqk ( double rm, double dm, double pr, double pd,
112             double px, double rv, double amprms[21],
113             double *ra, double *da ) {
114              
115             /* local constants */
116             const double VF = 0.210945028; /* Km/s to AU/year */
117              
118             /* Local Variables: */
119             int i;
120             double ab1, abv[3], p[3], w, p1dv, p2[3], p3[3];
121             double pmt, gr2e, eb[3], q[3], pxr, em[3];
122             double pde, pdep1, p1[3], ehn[3], pn[3];
123              
124             /* Unpack scalar and vector parameters. */
125 1           pmt = amprms[0];
126 1           gr2e = amprms[7];
127 1           ab1 = amprms[11];
128 4 100         for( i = 0; i < 3; i++ ) {
129 3           eb[i] = amprms[i+1];
130 3           ehn[i] = amprms[i+4];
131 3           abv[i] = amprms[i+8];
132             }
133              
134             /* Spherical to x,y,z. */
135 1           eraS2c( rm, dm, q);
136              
137             /* Space motion (radians per year) */
138 1           pxr = px * PAL__DAS2R;
139 1           w = VF * rv * pxr;
140 1           em[0] = -pr * q[1] - pd * cos(rm) * sin(dm) + w * q[0];
141 1           em[1] = pr * q[0] - pd * sin(rm) * sin(dm) + w * q[1];
142 1           em[2] = pd * cos(dm) + w * q[2];
143              
144             /* Geocentric direction of star (normalised) */
145 4 100         for( i = 0; i < 3; i++ ) {
146 3           p[i] = q[i] + pmt * em[i] - pxr * eb[i];
147             }
148 1           eraPn( p, &w, pn );
149              
150             /* Light deflection (restrained within the Sun's disc) */
151 1           pde = eraPdp( pn, ehn );
152 1           pdep1 = pde + 1.0;
153 1 50         w = gr2e / ( pdep1 > 1.0e-5 ? pdep1 : 1.0e-5 );
154 4 100         for( i = 0; i < 3; i++) {
155 3           p1[i] = pn[i] + w * ( ehn[i] - pde * pn[i] );
156             }
157              
158             /* Aberration (normalisation omitted). */
159 1           p1dv = eraPdp( p, abv );
160 1           w = 1.0 + p1dv / ( ab1 + 1.0 );
161 4 100         for( i = 0; i < 3; i++ ) {
162 3           p2[i] = ( ab1 * p1[i] ) + ( w * abv[i] );
163             }
164              
165             /* Precession and nutation. */
166 1           eraRxp( (double(*)[3]) &rms[12], p2, p3 );
167              
168             /* Geocentric apparent RA,dec. */
169 1           eraC2s( p3, ra, da );
170 1           *ra = eraAnp( *ra );
171              
172 1           }