source: trunk/MagicSoft/slalib/aopqk.c

Last change on this file was 731, checked in by tbretz, 23 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 8.8 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaAopqk ( double rap, double dap, double aoprms[14],
4 double *aob, double *zob, double *hob,
5 double *dob, double *rob )
6/*
7** - - - - - - - - -
8** s l a A o p q k
9** - - - - - - - - -
10**
11** Quick apparent to observed place (but see note 8, below, for
12** remarks about speed).
13**
14** Given:
15** rap double geocentric apparent right ascension
16** dap double geocentric apparent declination
17** aoprms double[14] star-independent apparent-to-observed parameters:
18**
19** (0) geodetic latitude (radians)
20** (1,2) sine and cosine of geodetic latitude
21** (3) magnitude of diurnal aberration vector
22** (4) height (hm)
23** (5) ambient temperature (t)
24** (6) pressure (p)
25** (7) relative humidity (rh)
26** (8) wavelength (wl)
27** (9) lapse rate (tlr)
28** (10,11) refraction constants A and B (radians)
29** (12) longitude + eqn of equinoxes + sidereal DUT (radians)
30** (13) local apparent sidereal time (radians)
31**
32** Returned:
33** *aob double observed azimuth (radians: N=0,E=90)
34** *zob double observed zenith distance (radians)
35** *hob double observed hour angle (radians)
36** *dob double observed declination (radians)
37** *rob double observed right ascension (radians)
38**
39** Notes:
40**
41** 1) This routine returns zenith distance rather than elevation
42** in order to reflect the fact that no allowance is made for
43** depression of the horizon.
44**
45** 2) The accuracy of the result is limited by the corrections for
46** refraction. Providing the meteorological parameters are
47** known accurately and there are no gross local effects, the
48** observed RA,Dec predicted by this routine should be within
49** about 0.1 arcsec for a zenith distance of less than 70 degrees.
50** Even at a topocentric zenith distance of 90 degrees, the
51** accuracy in elevation should be better than 1 arcmin; useful
52** results are available for a further 3 degrees, beyond which
53** the slaRefro routine returns a fixed value of the refraction.
54** The complementary routines slaAop (or slaAopqk) and slaOap
55** (or slaOapqk) are self-consistent to better than 1 micro-
56** arcsecond all over the celestial sphere.
57**
58** 3) It is advisable to take great care with units, as even
59** unlikely values of the input parameters are accepted and
60** processed in accordance with the models used.
61**
62** 4) "Apparent" place means the geocentric apparent right ascension
63** and declination, which is obtained from a catalogue mean place
64** by allowing for space motion, parallax, precession, nutation,
65** annual aberration, and the Sun's gravitational lens effect. For
66** star positions in the FK5 system (i.e. J2000), these effects can
67** be applied by means of the slaMap etc routines. Starting from
68** other mean place systems, additional transformations will be
69** needed; for example, FK4 (i.e. B1950) mean places would first
70** have to be converted to FK5, which can be done with the
71** slaFk425 etc routines.
72**
73** 5) "Observed" Az,El means the position that would be seen by a
74** perfect theodolite located at the observer. This is obtained
75** from the geocentric apparent RA,Dec by allowing for Earth
76** orientation and diurnal aberration, rotating from equator
77** to horizon coordinates, and then adjusting for refraction.
78** The HA,Dec is obtained by rotating back into equatorial
79** coordinates, using the geodetic latitude corrected for polar
80** motion, and is the position that would be seen by a perfect
81** equatorial located at the observer and with its polar axis
82** aligned to the Earth's axis of rotation (n.b. not to the
83** refracted pole). Finally, the RA is obtained by subtracting
84** the HA from the local apparent ST.
85**
86** 6) To predict the required setting of a real telescope, the
87** observed place produced by this routine would have to be
88** adjusted for the tilt of the azimuth or polar axis of the
89** mounting (with appropriate corrections for mount flexures),
90** for non-perpendicularity between the mounting axes, for the
91** position of the rotator axis and the pointing axis relative
92** to it, for tube flexure, for gear and encoder errors, and
93** finally for encoder zero points. Some telescopes would, of
94** course, exhibit other properties which would need to be
95** accounted for at the appropriate point in the sequence.
96**
97** 7) The star-independent apparent-to-observed-place parameters
98** in aoprms may be computed by means of the slaAoppa routine.
99** If nothing has changed significantly except the time, the
100** slaAoppat routine may be used to perform the requisite
101** partial recomputation of aoprms.
102**
103** 8) At zenith distances beyond about 76 degrees, the need for
104** special care with the corrections for refraction causes a
105** marked increase in execution time. Moreover, the effect
106** gets worse with increasing zenith distance. Adroit
107** programming in the calling application may allow the
108** problem to be reduced. Prepare an alternative aoprms array,
109** computed for zero air-pressure; this will disable the
110** refraction corrections and cause rapid execution. Using
111** this aoprms array, a preliminary call to the present routine
112** will, depending on the application, produce a rough position
113** which may be enough to establish whether the full, slow
114** calculation (using the real aoprms array) is worthwhile.
115** For example, there would be no need for the full calculation
116** if the preliminary call had already established that the
117** source was well below the elevation limits for a particular
118** telescope.
119**
120** 9) The azimuths etc produced by the present routine are with
121** respect to the celestial pole. Corrections to the terrestrial
122** pole can be computed using slaPolmo.
123**
124** Called: slaDcs2c, slaRefz, slaRefro, slaDcc2s, slaDranrm
125**
126** Last revision: 22 February 1996
127**
128** Copyright P.T.Wallace. All rights reserved.
129*/
130{
131/*
132** Breakpoint for fast/slow refraction algorithm:
133** ZD greater than arctan(4), (see slaRefco routine)
134** or vector z less than cosine(arctan(z)) = 1/sqrt(17)
135*/
136 static double zbreak = 0.242535625;
137
138 int i;
139
140 double sphi, cphi, st, v[3], xhd, yhd, zhd, diurab, f,
141 xhdt, yhdt, zhdt, xaet, yaet, zaet, azobs,
142 zdt, refa, refb, zdobs, dzd, dref, ce,
143 xaeo, yaeo, zaeo, hmobs, dcobs, raobs;
144
145/* Sin, cos of latitude */
146 sphi = aoprms[1];
147 cphi = aoprms[2];
148
149/* Local apparent sidereal time */
150 st = aoprms[13];
151
152/* Apparent RA,Dec to Cartesian -HA,Dec */
153 slaDcs2c(rap - st, dap, v);
154 xhd = v[0];
155 yhd = v[1];
156 zhd = v[2];
157
158/* Diurnal aberration */
159 diurab = aoprms[3];
160 f = 1.0 - diurab * yhd;
161 xhdt = f * xhd;
162 yhdt = f * ( yhd + diurab );
163 zhdt = f * zhd;
164
165/* Cartesian -HA,Dec to Cartesian az,el (S=0,E=90) */
166 xaet = sphi * xhdt - cphi * zhdt;
167 yaet = yhdt;
168 zaet = cphi * xhdt + sphi * zhdt;
169
170/* Azimuth (N=0,E=90) */
171 azobs = ( (xaet == 0.0) && (yaet == 0.0) ) ?
172 0.0 : atan2 ( yaet, -xaet );
173
174/* Topocentric zenith distance */
175 zdt = atan2 ( sqrt ( xaet * xaet + yaet * yaet ), zaet );
176
177/*
178** Refraction
179** ----------
180*/
181
182/* Fast algorithm using two constant model */
183 refa = aoprms[10];
184 refb = aoprms[11];
185 slaRefz ( zdt, refa, refb, &zdobs );
186
187/* Large zenith distance? */
188 if ( cos ( zdobs ) < zbreak ) {
189
190 /* Yes: use rigorous algorithm */
191
192 /* Initialize loop (maximum of 10 iterations) */
193 i = 1;
194 do {
195
196 /* Compute refraction using current estimate of observed ZD */
197 slaRefro ( zdobs, aoprms[4], aoprms[5], aoprms[6],
198 aoprms[7], aoprms[8], aoprms[0],
199 aoprms[9], 1e-8, &dref );
200
201 /* Remaining discrepancy */
202 dzd = zdobs + dref - zdt;
203
204 /* Update the estimate */
205 zdobs -= dzd;
206
207 /* Increment the iteration counter */
208 i++;
209
210 } while ( fabs ( dzd ) > 1e-10 && i <= 10 );
211 }
212
213/* To Cartesian az/ZD */
214 ce = sin ( zdobs );
215 xaeo = -cos ( azobs ) * ce;
216 yaeo = sin ( azobs ) * ce;
217 zaeo = cos ( zdobs );
218
219/* Cartesian az/ZD to Cartesian -HA,Dec */
220 v[0] = sphi * xaeo + cphi * zaeo;
221 v[1] = yaeo;
222 v[2] = -cphi * xaeo + sphi * zaeo;
223
224/* To spherical -HA,dec */
225 slaDcc2s ( v, &hmobs, &dcobs );
226
227/* Right ascension */
228 raobs = slaDranrm ( st + hmobs );
229
230/* Return the results */
231 *aob = azobs;
232 *zob = zdobs;
233 *hob = -hmobs;
234 *dob = dcobs;
235 *rob = raobs;
236}
Note: See TracBrowser for help on using the repository browser.