1 /*
2 *+
3 * Name:
4 * palAltaz
5
6 * Purpose:
7 * Positions, velocities and accelerations for an altazimuth telescope mount
8
9 * Language:
10 * Starlink ANSI C
11
12 * Type of Module:
13 * Library routine
14
15 * Invocation:
16 * palAltaz ( double ha, double dec, double phi,
17 * double *az, double *azd, double *azdd,
18 * double *el, double *eld, double *eldd,
19 * double *pa, double *pad, double *padd );
20
21 * Arguments:
22 * ha = double (Given)
23 * Hour angle (radians)
24 * dec = double (Given)
25 * Declination (radians)
26 * phi = double (Given)
27 * Observatory latitude (radians)
28 * az = double * (Returned)
29 * Azimuth (radians)
30 * azd = double * (Returned)
31 * Azimuth velocity (radians per radian of HA)
32 * azdd = double * (Returned)
33 * Azimuth acceleration (radians per radian of HA squared)
34 * el = double * (Returned)
35 * Elevation (radians)
36 * eld = double * (Returned)
37 * Elevation velocity (radians per radian of HA)
38 * eldd = double * (Returned)
39 * Elevation acceleration (radians per radian of HA squared)
40 * pa = double * (Returned)
41 * Parallactic angle (radians)
42 * pad = double * (Returned)
43 * Parallactic angle velocity (radians per radian of HA)
44 * padd = double * (Returned)
45 * Parallactic angle acceleration (radians per radian of HA squared)
46
47
48 * Description:
49 * Positions, velocities and accelerations for an altazimuth
50 * telescope mount.
51
52 * Authors:
53 * PTW: P. T. Wallace
54 * TIMJ: Tim Jenness (Cornell)
55 * {enter_new_authors_here}
56
57 * Notes:
58 * - Natural units are used throughout. HA, DEC, PHI, AZ, EL
59 * and ZD are in radians. The velocities and accelerations
60 * assume constant declination and constant rate of change of
61 * hour angle (as for tracking a star); the units of AZD, ELD
62 * and PAD are radians per radian of HA, while the units of AZDD,
63 * ELDD and PADD are radians per radian of HA squared. To
64 * convert into practical degree- and second-based units:
65 *
66 * angles * 360/2pi -> degrees
67 * velocities * (2pi/86400)*(360/2pi) -> degree/sec
68 * accelerations * ((2pi/86400)**2)*(360/2pi) -> degree/sec/sec
69 *
70 * Note that the seconds here are sidereal rather than SI. One
71 * sidereal second is about 0.99727 SI seconds.
72 *
73 * The velocity and acceleration factors assume the sidereal
74 * tracking case. Their respective numerical values are (exactly)
75 * 1/240 and (approximately) 1/3300236.9.
76 *
77 * - Azimuth is returned in the range 0-2pi; north is zero,
78 * and east is +pi/2. Elevation and parallactic angle are
79 * returned in the range +/-pi. Parallactic angle is +ve for
80 * a star west of the meridian and is the angle NP-star-zenith.
81 *
82 * - The latitude is geodetic as opposed to geocentric. The
83 * hour angle and declination are topocentric. Refraction and
84 * deficiencies in the telescope mounting are ignored. The
85 * purpose of the routine is to give the general form of the
86 * quantities. The details of a real telescope could profoundly
87 * change the results, especially close to the zenith.
88 *
89 * - No range checking of arguments is carried out.
90 *
91 * - In applications which involve many such calculations, rather
92 * than calling the present routine it will be more efficient to
93 * use inline code, having previously computed fixed terms such
94 * as sine and cosine of latitude, and (for tracking a star)
95 * sine and cosine of declination.
96
97 * History:
98 * 2014-09-30 (TIMJ):
99 * Initial version. Ported from Fortran SLA
100 * {enter_further_changes_here}
101
102 * Copyright:
103 * Copyright (C) 2004 P.T. Wallace
104 * Copyright (C) 2014 Cornell University
105 * All Rights Reserved.
106
107 * Licence:
108 * This program is free software; you can redistribute it and/or
109 * modify it under the terms of the GNU General Public License as
110 * published by the Free Software Foundation; either version 3 of
111 * the License, or (at your option) any later version.
112 *
113 * This program is distributed in the hope that it will be
114 * useful, but WITHOUT ANY WARRANTY; without even the implied
115 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
116 * PURPOSE. See the GNU General Public License for more details.
117 *
118 * You should have received a copy of the GNU General Public License
119 * along with this program. If not, see <http://www.gnu.org/licenses/>.
120
121 * Bugs:
122 * {note_any_bugs_here}
123 *-
124 */
125
126 #include <math.h>
127
128 #include "pal.h"
129 #include "palmac.h"
130
131 void
palAltaz(double ha,double dec,double phi,double * az,double * azd,double * azdd,double * el,double * eld,double * eldd,double * pa,double * pad,double * padd)132 palAltaz ( double ha, double dec, double phi,
133 double *az, double *azd, double *azdd,
134 double *el, double *eld, double *eldd,
135 double *pa, double *pad, double *padd ) {
136
137 const double TINY = 1E-30;
138
139 double sh,ch,sd,cd,sp,cp,chcd,sdcp,x,y,z,rsq,r,a,e,c,s,
140 q,qd,ad,ed,edr,add,edd,qdd;
141
142
143 /* Useful functions */
144 sh=sin(ha);
145 ch=cos(ha);
146 sd=sin(dec);
147 cd=cos(dec);
148 sp=sin(phi);
149 cp=cos(phi);
150 chcd=ch*cd;
151 sdcp=sd*cp;
152 x=-chcd*sp+sdcp;
153 y=-sh*cd;
154 z=chcd*cp+sd*sp;
155 rsq=x*x+y*y;
156 r=sqrt(rsq);
157
158 /* Azimuth and elevation */
159 if (rsq == 0.0) {
160 a=0.0;
161 } else {
162 a=atan2(y,x);
163 }
164 if (a < 0.0) a += PAL__D2PI;
165 e=atan2(z,r);
166
167 /* Parallactic angle */
168 c=cd*sp-ch*sdcp;
169 s=sh*cp;
170 if (c*c+s*s > 0) {
171 q=atan2(s,c);
172 } else {
173 q= PAL__DPI - ha;
174 }
175
176 /* Velocities and accelerations (clamped at zenith/nadir) */
177 if (rsq < TINY) {
178 rsq=TINY;
179 r=sqrt(rsq);
180 }
181 qd=-x*cp/rsq;
182 ad=sp+z*qd;
183 ed=cp*y/r;
184 edr=ed/r;
185 add=edr*(z*sp+(2.0-rsq)*qd);
186 edd=-r*qd*ad;
187 qdd=edr*(sp+2.0*z*qd);
188
189 /* Results */
190 *az=a;
191 *azd=ad;
192 *azdd=add;
193 *el=e;
194 *eld=ed;
195 *eldd=edd;
196 *pa=q;
197 *pad=qd;
198 *padd=qdd;
199
200 }
201