1 /*
2 *+
3 *  Name:
4 *     palDh2e
5 
6 *  Purpose:
7 *     Horizon to equatorial coordinates: Az,El to HA,Dec
8 
9 *  Language:
10 *     Starlink ANSI C
11 
12 *  Type of Module:
13 *     Library routine
14 
15 *  Invocation:
16 *     palDh2e( double az, double el, double phi, double * ha, double * dec );
17 
18 *  Arguments:
19 *     az = double (Given)
20 *        Azimuth (radians)
21 *     el = double (Given)
22 *        Elevation (radians)
23 *     phi = double (Given)
24 *        Observatory latitude (radians)
25 *     ha = double * (Returned)
26 *        Hour angle (radians)
27 *     dec = double * (Returned)
28 *        Declination (radians)
29 
30 *  Description:
31 *     Convert horizon to equatorial coordinates.
32 
33 *  Authors:
34 *     PTW: Pat Wallace (STFC)
35 *     TIMJ: Tim Jenness (JAC, Hawaii)
36 *     {enter_new_authors_here}
37 
38 *  Notes:
39 *     - All the arguments are angles in radians.
40 *     - The sign convention for azimuth is north zero, east +pi/2.
41 *     - HA is returned in the range +/-pi.  Declination is returned
42 *       in the range +/-pi/2.
43 *     - The latitude is (in principle) geodetic.  In critical
44 *       applications, corrections for polar motion should be applied.
45 *     - In some applications it will be important to specify the
46 *       correct type of elevation in order to produce the required
47 *       type of HA,Dec.  In particular, it may be important to
48 *       distinguish between the elevation as affected by refraction,
49 *       which will yield the "observed" HA,Dec, and the elevation
50 *       in vacuo, which will yield the "topocentric" HA,Dec.  If the
51 *       effects of diurnal aberration can be neglected, the
52 *       topocentric HA,Dec may be used as an approximation to the
53 *       "apparent" HA,Dec.
54 *     - No range checking of arguments is done.
55 *     - In applications which involve many such calculations, rather
56 *       than calling the present routine it will be more efficient to
57 *       use inline code, having previously computed fixed terms such
58 *       as sine and cosine of latitude.
59 
60 *  History:
61 *     2012-02-08 (TIMJ):
62 *        Initial version with documentation taken from Fortran SLA
63 *        Adapted with permission from the Fortran SLALIB library.
64 *     {enter_further_changes_here}
65 
66 *  Copyright:
67 *     Copyright (C) 1996 Rutherford Appleton Laboratory
68 *     Copyright (C) 2012 Science and Technology Facilities Council.
69 *     All Rights Reserved.
70 
71 *  Licence:
72 *     This program is free software: you can redistribute it and/or
73 *     modify it under the terms of the GNU Lesser General Public
74 *     License as published by the Free Software Foundation, either
75 *     version 3 of the License, or (at your option) any later
76 *     version.
77 *
78 *     This program is distributed in the hope that it will be useful,
79 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
80 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
81 *     GNU Lesser General Public License for more details.
82 *
83 *     You should have received a copy of the GNU Lesser General
84 *     License along with this program.  If not, see
85 *     <http://www.gnu.org/licenses/>.
86 
87 *  Bugs:
88 *     {note_any_bugs_here}
89 *-
90 */
91 
92 #include "pal.h"
93 #include <math.h>
94 
95 void
palDh2e(double az,double el,double phi,double * ha,double * dec)96 palDh2e ( double az, double el, double phi, double *ha, double *dec) {
97 
98   double sa;
99   double ca;
100   double se;
101   double ce;
102   double sp;
103   double cp;
104 
105   double x;
106   double y;
107   double z;
108   double r;
109 
110   /*  Useful trig functions */
111   sa = sin(az);
112   ca = cos(az);
113   se = sin(el);
114   ce = cos(el);
115   sp = sin(phi);
116   cp = cos(phi);
117 
118   /*  HA,Dec as x,y,z */
119   x = -ca * ce * sp + se * cp;
120   y = -sa * ce;
121   z = ca * ce * cp + se * sp;
122 
123   /*  To HA,Dec */
124   r = sqrt(x * x + y * y);
125   if (r == 0.) {
126     *ha = 0.;
127   } else {
128     *ha = atan2(y, x);
129   }
130   *dec = atan2(z, r);
131 
132   return;
133 }
134