File: palMapqk.c

package info (click to toggle)
starlink-pal 0.9.8-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 1,808 kB
  • sloc: ansic: 6,689; makefile: 128; sh: 81
file content (172 lines) | stat: -rw-r--r-- 5,563 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
/*
*+
*  Name:
*     palMapqk

*  Purpose:
*     Quick mean to apparent place

*  Language:
*     Starlink ANSI C

*  Type of Module:
*     Library routine

*  Invocation:
*     void palMapqk ( double rm, double dm, double pr, double pd,
*                     double px, double rv, double amprms[21],
*                     double *ra, double *da );

*  Arguments:
*     rm = double (Given)
*        Mean RA (radians)
*     dm = double (Given)
*        Mean declination (radians)
*     pr = double (Given)
*        RA proper motion, changes per Julian year (radians)
*     pd = double (Given)
*        Dec proper motion, changes per Julian year (radians)
*     px = double (Given)
*        Parallax (arcsec)
*     rv = double (Given)
*        Radial velocity (km/s, +ve if receding)
*     amprms = double [21] (Given)
*        Star-independent mean-to-apparent parameters (see palMappa).
*     ra = double * (Returned)
*        Apparent RA (radians)
*     dec = double * (Returned)
*        Apparent dec (radians)

*  Description:
*     Quick mean to apparent place:  transform a star RA,Dec from
*     mean place to geocentric apparent place, given the
*     star-independent parameters.
*
*     Use of this routine is appropriate when efficiency is important
*     and where many star positions, all referred to the same equator
*     and equinox, are to be transformed for one epoch.  The
*     star-independent parameters can be obtained by calling the
*     palMappa routine.
*
*     If the parallax and proper motions are zero the palMapqkz
*     routine can be used instead.

*  Notes:
*     - The reference frames and timescales used are post IAU 2006.
*     - The mean place rm, dm and the vectors amprms[1-3] and amprms[4-6]
*       are referred to the mean equinox and equator of the epoch
*       specified when generating the precession/nutation matrix
*       amprms[12-20].  In the call to palMappa (q.v.) normally used
*       to populate amprms, this epoch is the first argument (eq).
*     - Strictly speaking, the routine is not valid for solar-system
*       sources, though the error will usually be extremely small.
*       However, to prevent gross errors in the case where the
*       position of the Sun is specified, the gravitational
*       deflection term is restrained within about 920 arcsec of the
*       centre of the Sun's disc.  The term has a maximum value of
*       about 1.85 arcsec at this radius, and decreases to zero as
*       the centre of the disc is approached.

*  Authors:
*     PTW: Patrick T. Wallace
*     TIMJ: Tim Jenness (JAC, Hawaii)
*     {enter_new_authors_here}

*  History:
*     2012-03-01 (TIMJ):
*        Initial version with documentation from SLA/F
*        Adapted with permission from the Fortran SLALIB library.
*     {enter_further_changes_here}

*  Copyright:
*     Copyright (C) 2000 Rutherford Appleton Laboratory
*     Copyright (C) 2012 Science and Technology Facilities Council.
*     All Rights Reserved.

*  Licence:
*     This program is free software; you can redistribute it and/or
*     modify it under the terms of the GNU General Public License as
*     published by the Free Software Foundation; either version 3 of
*     the License, or (at your option) any later version.
*
*     This program is distributed in the hope that it will be
*     useful, but WITHOUT ANY WARRANTY; without even the implied
*     warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
*     PURPOSE. See the GNU General Public License for more details.
*
*     You should have received a copy of the GNU General Public License
*     along with this program; if not, write to the Free Software
*     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
*     MA 02110-1301, USA.

*  Bugs:
*     {note_any_bugs_here}
*-
*/

#include "pal.h"
#include "palmac.h"
#include "pal1sofa.h"

void palMapqk ( double rm, double dm, double pr, double pd,
                double px, double rv, double amprms[21],
                double *ra, double *da ) {

/* local constants */
   const double VF = 0.210945028; /* Km/s to AU/year */

/* Local Variables: */
   int i;
   double ab1, abv[3], p[3], w, p1dv, p2[3], p3[3];
   double pmt, gr2e, eb[3], q[3], pxr, em[3];
   double pde, pdep1, p1[3], ehn[3], pn[3];

/* Unpack scalar and vector parameters. */
   pmt = amprms[0];
   gr2e = amprms[7];
   ab1 = amprms[11];
   for( i = 0; i < 3; i++ ) {
      eb[i] = amprms[i+1];
      ehn[i] = amprms[i+4];
      abv[i] = amprms[i+8];
   }

/* Spherical to x,y,z. */
   eraS2c( rm, dm, q);

 /* Space motion (radians per year) */
   pxr = px * PAL__DAS2R;
   w = VF * rv * pxr;
   em[0] = -pr * q[1] - pd * cos(rm) * sin(dm) + w * q[0];
   em[1] =  pr * q[0] - pd * sin(rm) * sin(dm) + w * q[1];
   em[2] =              pd * cos(dm)           + w * q[2];

/* Geocentric direction of star (normalised) */
   for( i = 0; i < 3; i++ ) {
      p[i] = q[i] + pmt * em[i] - pxr * eb[i];
   }
   eraPn( p, &w, pn );

/* Light deflection (restrained within the Sun's disc) */
   pde = eraPdp( pn, ehn );
   pdep1 = pde + 1.0;
   w = gr2e / ( pdep1 > 1.0e-5 ? pdep1 : 1.0e-5 );
   for( i = 0; i < 3; i++) {
      p1[i] = pn[i] + w * ( ehn[i] - pde * pn[i] );
   }

/* Aberration (normalisation omitted). */
   p1dv = eraPdp( p, abv );
   w = 1.0 + p1dv / ( ab1 + 1.0 );
   for( i = 0; i < 3; i++ ) {
      p2[i] = ( ab1 * p1[i] ) + ( w * abv[i] );
   }

/* Precession and nutation. */
   eraRxp( (double(*)[3]) &amprms[12], p2, p3 );

/* Geocentric apparent RA,dec. */
   eraC2s( p3, ra, da );
   *ra = eraAnp( *ra );

}