source: trunk/MagicSoft/slalib/oapqk.c

Last change on this file was 1953, checked in by tbretz, 22 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 7.7 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaOapqk ( char *type, double ob1, double ob2,
4 double aoprms[14], double *rap, double *dap )
5/*
6** - - - - - - - - -
7** s l a O a p q k
8** - - - - - - - - -
9**
10** Quick observed to apparent place.
11**
12** Given:
13** type char type of coordinates - 'r', 'h' or 'a' (see below)
14** ob1 double observed az, HA or RA (radians; az is n=0,e=90)
15** ob2 double observed ZD or Dec (radians)
16** aoprms double[14] star-independent apparent-to-observed parameters:
17**
18** (0) geodetic latitude (radians)
19** (1,2) sine and cosine of geodetic latitude
20** (3) magnitude of diurnal aberration vector
21** (4) height (hm)
22** (5) ambient temperature (t)
23** (6) pressure (p)
24** (7) relative humidity (rh)
25** (8) wavelength (wl)
26** (9) lapse rate (tlr)
27** (10,11) refraction constants a and b (radians)
28** (12) longitude + eqn of equinoxes + sidereal DUT (radians)
29** (13) local apparent sidereal time (radians)
30**
31** Returned:
32** *rap double geocentric apparent right ascension
33** *dap double geocentric apparent declination
34**
35** Notes:
36**
37** 1) Only the first character of the type argument is significant.
38** 'R' or 'r' indicates that obs1 and obs2 are the observed right
39** ascension and declination; 'H' or 'h' indicates that they are
40** hour angle (west +ve) and declination; anything else ('A' or
41** 'a' is recommended) indicates that obs1 and obs2 are azimuth
42** (north zero, east is 90 deg) and zenith distance. (Zenith
43** distance is used rather than elevation in order to reflect the
44** fact that no allowance is made for depression of the horizon.)
45**
46** 2) The accuracy of the result is limited by the corrections for
47** refraction. Providing the meteorological parameters are
48** known accurately and there are no gross local effects, the
49** predicted apparent RA,Dec should be within about 0.1 arcsec.
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** 5) "Observed" az,el means the position that would be seen by a
63** perfect theodolite located at the observer. This is
64** related to the observed HA,Dec via the standard rotation, using
65** the geodetic latitude (corrected for polar motion), while the
66** observed HA and RA are related simply through the local
67** apparent ST. "Observed" RA,Dec or HA,Dec thus means the
68** position that would be seen by a perfect equatorial located
69** at the observer and with its polar axis aligned to the
70** Earth's axis of rotation (n.b. not to the refracted pole).
71** by removing from the observed place the effects of
72** atmospheric refraction and diurnal aberration, the
73** geocentric apparent RA,Dec is obtained.
74**
75** 5) Frequently, mean rather than apparent RA,Dec will be required,
76** in which case further transformations will be necessary. The
77** slaAmp etc routines will convert the apparent RA,Dec produced
78** by the present routine into an "FK5" (J2000) mean place, by
79** allowing for the Sun's gravitational lens effect, annual
80** aberration, nutation and precession. Should "FK4" (1950)
81** coordinates be needed, the routines slaFk524 etc will also
82** need to be applied.
83**
84** 6) To convert to apparent RA,Dec the coordinates read from a
85** real telescope, corrections would have to be applied for
86** encoder zero points, gear and encoder errors, tube flexure,
87** the position of the rotator axis and the pointing axis
88** relative to it, non-perpendicularity between the mounting
89** axes, and finally for the tilt of the azimuth or polar axis
90** of the mounting (with appropriate corrections for mount
91** flexures). Some telescopes would, of course, exhibit other
92** properties which would need to be accounted for at the
93** appropriate point in the sequence.
94**
95** 7) The star-independent apparent-to-observed-place parameters
96** in aoprms may be computed by means of the slaAoppa routine.
97** If nothing has changed significantly except the time, the
98** slaAoppat routine may be used to perform the requisite
99** partial recomputation of aoprms.
100**
101** 8) The azimuths etc used by the present routine are with respect
102** to the celestial pole. Corrections from the terrestrial pole
103** can be computed using slaPolmo.
104**
105** Called: slaDcs2c, slaDcc2s, slaRefro, slaDranrm
106**
107** Last revision: 3 February 2000
108**
109** Copyright P.T.Wallace. All rights reserved.
110*/
111{
112
113/*
114** Breakpoint for fast/slow refraction algorithm:
115** ZD greater than arctan(4), (see slaRefco routine)
116** or vector z less than cosine(arctan(z)) = 1/sqrt(17)
117*/
118 static double zbreak = 0.242535625;
119
120 char c;
121 double c1, c2, sphi, cphi, st, ce, xaeo, yaeo, zaeo, v[3],
122 xmhdo, ymhdo, zmhdo, az, sz, zdo, tz, dref, zdt,
123 xaet, yaet, zaet, xmhda, ymhda, zmhda, diurab, f, hma;
124
125
126/* Coordinate type */
127 c = *type;
128
129/* Coordinates */
130 c1 = ob1;
131 c2 = ob2;
132
133/* Sin, cos of latitude */
134 sphi = aoprms[1];
135 cphi = aoprms[2];
136
137/* Local apparent sidereal time */
138 st = aoprms[13];
139
140/* Standardize coordinate type */
141 if ( c == 'r' || c == 'R' ) {
142 c = 'R';
143 } else if ( c == 'h' || c == 'H' ) {
144 c = 'H';
145 } else {
146 c = 'A';
147 }
148
149/* If az,ZD convert to Cartesian (S=0,E=90) */
150 if ( c == 'A' ) {
151 ce = sin ( c2 );
152 xaeo = - cos ( c1 ) * ce;
153 yaeo = sin ( c1 ) * ce;
154 zaeo = cos ( c2 );
155 } else {
156
157 /* If RA,Dec convert to HA,Dec */
158 if ( c == 'R' ) {
159 c1 = st - c1;
160 }
161
162 /* To Cartesian -HA,Dec */
163 slaDcs2c ( -c1, c2, v );
164 xmhdo = v[0];
165 ymhdo = v[1];
166 zmhdo = v[2];
167
168 /* To Cartesian az,el (S=0,E=90) */
169 xaeo = sphi * xmhdo - cphi * zmhdo;
170 yaeo = ymhdo;
171 zaeo = cphi * xmhdo + sphi * zmhdo;
172 }
173
174/* Azimuth (S=0,E=90) */
175 //*TB* az = xaeo != 0.0 && yaeo != 0.0 ? atan2 ( yaeo, xaeo ) : 0.0;
176 az = atan2 ( yaeo, xaeo );
177
178/* Sine of observed ZD, and observed ZD */
179 sz = sqrt ( xaeo * xaeo + yaeo * yaeo );
180 zdo = atan2 ( sz, zaeo );
181
182/*
183** Refraction
184** ----------
185*/
186
187/* Large zenith distance? */
188 if ( zaeo >= zbreak ) {
189
190 /* Fast algorithm using two constant model */
191 tz = sz / zaeo;
192 dref = ( aoprms[10] + aoprms[11] * tz * tz ) * tz;
193 } else {
194
195 /* Rigorous algorithm for large ZD */
196 slaRefro ( zdo, aoprms[4], aoprms[5], aoprms[6], aoprms[7],
197 aoprms[8], aoprms[0], aoprms[9], 1e-8, &dref );
198 }
199 zdt = zdo + dref;
200
201/* To Cartesian az,ZD */
202 ce = sin ( zdt );
203 xaet = cos ( az ) * ce;
204 yaet = sin ( az ) * ce;
205 zaet = cos ( zdt );
206
207/* Cartesian az,ZD to Cartesian -HA,Dec */
208 xmhda = sphi * xaet + cphi * zaet;
209 ymhda = yaet;
210 zmhda = - cphi * xaet + sphi * zaet;
211
212/* Diurnal aberration */
213 diurab = -aoprms[3];
214 f = 1.0 - diurab * ymhda;
215 v[0] = f * xmhda;
216 v[1] = f * ( ymhda + diurab );
217 v[2] = f * zmhda;
218
219/* To spherical -HA,Dec */
220 slaDcc2s ( v, &hma, dap );
221
222/* Right ascension */
223 *rap = slaDranrm ( st + hma );
224}
Note: See TracBrowser for help on using the repository browser.