source: trunk/FACT++/pal/palAopqk.c

Last change on this file was 18347, checked in by tbretz, 9 years ago
File size: 10.0 KB
Line 
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
179void 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}
Note: See TracBrowser for help on using the repository browser.