1 /*
2 *+
3 * Name:
4 * palAopqk
5
6 * Purpose:
7 * Quick apparent to observed place
8
9 * Language:
10 * Starlink ANSI C
11
12 * Type of Module:
13 * Library routine
14
15 * Invocation:
16 * void palAopqk ( double rap, double dap, const double aoprms[14],
17 * double *aob, double *zob, double *hob,
18 * double *dob, double *rob );
19
20 * Arguments:
21 * rap = double (Given)
22 * Geocentric apparent right ascension
23 * dap = double (Given)
24 * Geocentric apparent declination
25 * aoprms = const double [14] (Given)
26 * Star-independent apparent-to-observed parameters.
27 *
28 * [0] geodetic latitude (radians)
29 * [1,2] sine and cosine of geodetic latitude
30 * [3] magnitude of diurnal aberration vector
31 * [4] height (HM)
32 * [5] ambient temperature (T)
33 * [6] pressure (P)
34 * [7] relative humidity (RH)
35 * [8] wavelength (WL)
36 * [9] lapse rate (TLR)
37 * [10,11] refraction constants A and B (radians)
38 * [12] longitude + eqn of equinoxes + sidereal DUT (radians)
39 * [13] local apparent sidereal time (radians)
40 * aob = double * (Returned)
41 * Observed azimuth (radians: N=0,E=90)
42 * zob = double * (Returned)
43 * Observed zenith distance (radians)
44 * hob = double * (Returned)
45 * Observed Hour Angle (radians)
46 * dob = double * (Returned)
47 * Observed Declination (radians)
48 * rob = double * (Returned)
49 * Observed Right Ascension (radians)
50
51 * Description:
52 * Quick apparent to observed place.
53
54 * Authors:
55 * TIMJ: Tim Jenness (JAC, Hawaii)
56 * PTW: Patrick T. Wallace
57 * {enter_new_authors_here}
58
59 * Notes:
60 * - This routine returns zenith distance rather than elevation
61 * in order to reflect the fact that no allowance is made for
62 * depression of the horizon.
63 *
64 * - The accuracy of the result is limited by the corrections for
65 * refraction. Providing the meteorological parameters are
66 * known accurately and there are no gross local effects, the
67 * observed RA,Dec predicted by this routine should be within
68 * about 0.1 arcsec for a zenith distance of less than 70 degrees.
69 * Even at a topocentric zenith distance of 90 degrees, the
70 * accuracy in elevation should be better than 1 arcmin; useful
71 * results are available for a further 3 degrees, beyond which
72 * the palRefro routine returns a fixed value of the refraction.
73 * The complementary routines palAop (or palAopqk) and palOap
74 * (or palOapqk) are self-consistent to better than 1 micro-
75 * arcsecond all over the celestial sphere.
76 *
77 * - It is advisable to take great care with units, as even
78 * unlikely values of the input parameters are accepted and
79 * processed in accordance with the models used.
80 *
81 * - "Apparent" place means the geocentric apparent right ascension
82 * and declination, which is obtained from a catalogue mean place
83 * by allowing for space motion, parallax, precession, nutation,
84 * annual aberration, and the Sun's gravitational lens effect. For
85 * star positions in the FK5 system (i.e. J2000), these effects can
86 * be applied by means of the palMap etc routines. Starting from
87 * other mean place systems, additional transformations will be
88 * needed; for example, FK4 (i.e. B1950) mean places would first
89 * have to be converted to FK5, which can be done with the
90 * palFk425 etc routines.
91 *
92 * - "Observed" Az,El means the position that would be seen by a
93 * perfect theodolite located at the observer. This is obtained
94 * from the geocentric apparent RA,Dec by allowing for Earth
95 * orientation and diurnal aberration, rotating from equator
96 * to horizon coordinates, and then adjusting for refraction.
97 * The HA,Dec is obtained by rotating back into equatorial
98 * coordinates, using the geodetic latitude corrected for polar
99 * motion, and is the position that would be seen by a perfect
100 * equatorial located at the observer and with its polar axis
101 * aligned to the Earth's axis of rotation (n.b. not to the
102 * refracted pole). Finally, the RA is obtained by subtracting
103 * the HA from the local apparent ST.
104 *
105 * - To predict the required setting of a real telescope, the
106 * observed place produced by this routine would have to be
107 * adjusted for the tilt of the azimuth or polar axis of the
108 * mounting (with appropriate corrections for mount flexures),
109 * for non-perpendicularity between the mounting axes, for the
110 * position of the rotator axis and the pointing axis relative
111 * to it, for tube flexure, for gear and encoder errors, and
112 * finally for encoder zero points. Some telescopes would, of
113 * course, exhibit other properties which would need to be
114 * accounted for at the appropriate point in the sequence.
115 *
116 * - The star-independent apparent-to-observed-place parameters
117 * in AOPRMS may be computed by means of the palAoppa routine.
118 * If nothing has changed significantly except the time, the
119 * palAoppat routine may be used to perform the requisite
120 * partial recomputation of AOPRMS.
121 *
122 * - At zenith distances beyond about 76 degrees, the need for
123 * special care with the corrections for refraction causes a
124 * marked increase in execution time. Moreover, the effect
125 * gets worse with increasing zenith distance. Adroit
126 * programming in the calling application may allow the
127 * problem to be reduced. Prepare an alternative AOPRMS array,
128 * computed for zero air-pressure; this will disable the
129 * refraction corrections and cause rapid execution. Using
130 * this AOPRMS array, a preliminary call to the present routine
131 * will, depending on the application, produce a rough position
132 * which may be enough to establish whether the full, slow
133 * calculation (using the real AOPRMS array) is worthwhile.
134 * For example, there would be no need for the full calculation
135 * if the preliminary call had already established that the
136 * source was well below the elevation limits for a particular
137 * telescope.
138 *
139 * - The azimuths etc produced by the present routine are with
140 * respect to the celestial pole. Corrections to the terrestrial
141 * pole can be computed using palPolmo.
142
143 * History:
144 * 2012-08-25 (TIMJ):
145 * Initial version, copied from Fortran SLA
146 * Adapted with permission from the Fortran SLALIB library.
147 * {enter_further_changes_here}
148
149 * Copyright:
150 * Copyright (C) 2003 Rutherford Appleton Laboratory
151 * Copyright (C) 2012 Science and Technology Facilities Council.
152 * All Rights Reserved.
153
154 * Licence:
155 * This program is free software; you can redistribute it and/or
156 * modify it under the terms of the GNU General Public License as
157 * published by the Free Software Foundation; either version 3 of
158 * the License, or (at your option) any later version.
159 *
160 * This program is distributed in the hope that it will be
161 * useful, but WITHOUT ANY WARRANTY; without even the implied
162 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
163 * PURPOSE. See the GNU General Public License for more details.
164 *
165 * You should have received a copy of the GNU General Public License
166 * along with this program; if not, write to the Free Software
167 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
168 * MA 02110-1301, USA.
169
170 * Bugs:
171 * {note_any_bugs_here}
172 *-
173 */
174
175 #include <math.h>
176
177 #include "pal.h"
178
palAopqk(double rap,double dap,const double aoprms[14],double * aob,double * zob,double * hob,double * dob,double * rob)179 void palAopqk ( double rap, double dap, const double aoprms[14],
180 double *aob, double *zob, double *hob,
181 double *dob, double *rob ) {
182
183 /* Breakpoint for fast/slow refraction algorithm:
184 * ZD greater than arctan(4), (see palRefco routine)
185 * or vector Z less than cosine(arctan(Z)) = 1/sqrt(17) */
186 const double zbreak = 0.242535625;
187 int i;
188
189 double sphi,cphi,st,v[3],xhd,yhd,zhd,diurab,f,
190 xhdt,yhdt,zhdt,xaet,yaet,zaet,azobs,
191 zdt,refa,refb,zdobs,dzd,dref,ce,
192 xaeo,yaeo,zaeo,hmobs,dcobs,raobs;
193
194 /* sin, cos of latitude */
195 sphi = aoprms[1];
196 cphi = aoprms[2];
197
198 /* local apparent sidereal time */
199 st = aoprms[13];
200
201 /* apparent ra,dec to cartesian -ha,dec */
202 palDcs2c( rap-st, dap, v );
203 xhd = v[0];
204 yhd = v[1];
205 zhd = v[2];
206
207 /* diurnal aberration */
208 diurab = aoprms[3];
209 f = (1.0-diurab*yhd);
210 xhdt = f*xhd;
211 yhdt = f*(yhd+diurab);
212 zhdt = f*zhd;
213
214 /* cartesian -ha,dec to cartesian az,el (s=0,e=90) */
215 xaet = sphi*xhdt-cphi*zhdt;
216 yaet = yhdt;
217 zaet = cphi*xhdt+sphi*zhdt;
218
219 /* azimuth (n=0,e=90) */
220 if (xaet == 0.0 && yaet == 0.0) {
221 azobs = 0.0;
222 } else {
223 azobs = atan2(yaet,-xaet);
224 }
225
226 /* topocentric zenith distance */
227 zdt = atan2(sqrt(xaet*xaet+yaet*yaet),zaet);
228
229 /*
230 * refraction
231 * ---------- */
232
233 /* fast algorithm using two constant model */
234 refa = aoprms[10];
235 refb = aoprms[11];
236 palRefz(zdt,refa,refb,&zdobs);
237
238 /* large zenith distance? */
239 if (cos(zdobs) < zbreak) {
240
241 /* yes: use rigorous algorithm */
242
243 /* initialize loop (maximum of 10 iterations) */
244 i = 1;
245 dzd = 1.0e1;
246 while (fabs(dzd) > 1e-10 && i <= 10) {
247
248 /* compute refraction using current estimate of observed zd */
249 palRefro(zdobs,aoprms[4],aoprms[5],aoprms[6],
250 aoprms[7],aoprms[8],aoprms[0],
251 aoprms[9],1e-8,&dref);
252
253 /* remaining discrepancy */
254 dzd = zdobs+dref-zdt;
255
256 /* update the estimate */
257 zdobs = zdobs-dzd;
258
259 /* increment the iteration counter */
260 i++;
261 }
262 }
263
264 /* to cartesian az/zd */
265 ce = sin(zdobs);
266 xaeo = -cos(azobs)*ce;
267 yaeo = sin(azobs)*ce;
268 zaeo = cos(zdobs);
269
270 /* cartesian az/zd to cartesian -ha,dec */
271 v[0] = sphi*xaeo+cphi*zaeo;
272 v[1] = yaeo;
273 v[2] = -cphi*xaeo+sphi*zaeo;
274
275 /* to spherical -ha,dec */
276 palDcc2s(v,&hmobs,&dcobs);
277
278 /* right ascension */
279 raobs = palDranrm(st+hmobs);
280
281 /* return the results */
282 *aob = azobs;
283 *zob = zdobs;
284 *hob = -hmobs;
285 *dob = dcobs;
286 *rob = raobs;
287
288 }
289