FACT++  1.0
palOapqk.c
Go to the documentation of this file.
1 /*
2 *+
3 * Name:
4 * palOapqk
5 
6 * Purpose:
7 * Quick observed to apparent place
8 
9 * Language:
10 * Starlink ANSI C
11 
12 * Type of Module:
13 * Library routine
14 
15 * Invocation:
16 * void palOapqk ( const char *type, double ob1, double ob2,
17 * const double aoprms[14], double *rap, double *dap );
18 
19 * Arguments:
20 * Quick observed to apparent place.
21 
22 * Description:
23 * type = const char * (Given)
24 * Type of coordinates - 'R', 'H' or 'A' (see below)
25 * ob1 = double (Given)
26 * Observed Az, HA or RA (radians; Az is N=0;E=90)
27 * ob2 = double (Given)
28 * Observed ZD or Dec (radians)
29 * aoprms = const double [14] (Given)
30 * Star-independent apparent-to-observed parameters.
31 * See palAopqk for details.
32 * rap = double * (Given)
33 * Geocentric apparent right ascension
34 * dap = double * (Given)
35 * Geocentric apparent declination
36 
37 * Authors:
38 * PTW: Patrick T. Wallace
39 * TIMJ: Tim Jenness (JAC, Hawaii)
40 * {enter_new_authors_here}
41 
42 * Notes:
43 * - Only the first character of the TYPE argument is significant.
44 * 'R' or 'r' indicates that OBS1 and OBS2 are the observed right
45 * ascension and declination; 'H' or 'h' indicates that they are
46 * hour angle (west +ve) and declination; anything else ('A' or
47 * 'a' is recommended) indicates that OBS1 and OBS2 are azimuth
48 * (north zero, east 90 deg) and zenith distance. (Zenith distance
49 * is used rather than elevation in order to reflect the fact that
50 * no allowance is made for depression of the horizon.)
51 *
52 * - The accuracy of the result is limited by the corrections for
53 * refraction. Providing the meteorological parameters are
54 * known accurately and there are no gross local effects, the
55 * predicted apparent RA,Dec should be within about 0.1 arcsec
56 * for a zenith distance of less than 70 degrees. Even at a
57 * topocentric zenith distance of 90 degrees, the accuracy in
58 * elevation should be better than 1 arcmin; useful results
59 * are available for a further 3 degrees, beyond which the
60 * palREFRO routine returns a fixed value of the refraction.
61 * The complementary routines palAop (or palAopqk) and palOap
62 * (or palOapqk) are self-consistent to better than 1 micro-
63 * arcsecond all over the celestial sphere.
64 *
65 * - It is advisable to take great care with units, as even
66 * unlikely values of the input parameters are accepted and
67 * processed in accordance with the models used.
68 *
69 * - "Observed" Az,El means the position that would be seen by a
70 * perfect theodolite located at the observer. This is
71 * related to the observed HA,Dec via the standard rotation, using
72 * the geodetic latitude (corrected for polar motion), while the
73 * observed HA and RA are related simply through the local
74 * apparent ST. "Observed" RA,Dec or HA,Dec thus means the
75 * position that would be seen by a perfect equatorial located
76 * at the observer and with its polar axis aligned to the
77 * Earth's axis of rotation (n.b. not to the refracted pole).
78 * By removing from the observed place the effects of
79 * atmospheric refraction and diurnal aberration, the
80 * geocentric apparent RA,Dec is obtained.
81 *
82 * - Frequently, mean rather than apparent RA,Dec will be required,
83 * in which case further transformations will be necessary. The
84 * palAmp etc routines will convert the apparent RA,Dec produced
85 * by the present routine into an "FK5" (J2000) mean place, by
86 * allowing for the Sun's gravitational lens effect, annual
87 * aberration, nutation and precession. Should "FK4" (1950)
88 * coordinates be needed, the routines palFk524 etc will also
89 * need to be applied.
90 *
91 * - To convert to apparent RA,Dec the coordinates read from a
92 * real telescope, corrections would have to be applied for
93 * encoder zero points, gear and encoder errors, tube flexure,
94 * the position of the rotator axis and the pointing axis
95 * relative to it, non-perpendicularity between the mounting
96 * axes, and finally for the tilt of the azimuth or polar axis
97 * of the mounting (with appropriate corrections for mount
98 * flexures). Some telescopes would, of course, exhibit other
99 * properties which would need to be accounted for at the
100 * appropriate point in the sequence.
101 *
102 * - The star-independent apparent-to-observed-place parameters
103 * in AOPRMS may be computed by means of the palAoppa routine.
104 * If nothing has changed significantly except the time, the
105 * palAoppat routine may be used to perform the requisite
106 * partial recomputation of AOPRMS.
107 *
108 * - The azimuths etc used by the present routine are with respect
109 * to the celestial pole. Corrections from the terrestrial pole
110 * can be computed using palPolmo.
111 
112 
113 * History:
114 * 2012-08-27 (TIMJ):
115 * Initial version, direct copy of Fortran SLA
116 * Adapted with permission from the Fortran SLALIB library.
117 * {enter_further_changes_here}
118 
119 * Copyright:
120 * Copyright (C) 2004 Patrick T. Wallace
121 * Copyright (C) 2012 Science and Technology Facilities Council.
122 * All Rights Reserved.
123 
124 * Licence:
125 * This program is free software; you can redistribute it and/or
126 * modify it under the terms of the GNU General Public License as
127 * published by the Free Software Foundation; either version 3 of
128 * the License, or (at your option) any later version.
129 *
130 * This program is distributed in the hope that it will be
131 * useful, but WITHOUT ANY WARRANTY; without even the implied
132 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
133 * PURPOSE. See the GNU General Public License for more details.
134 *
135 * You should have received a copy of the GNU General Public License
136 * along with this program; if not, write to the Free Software
137 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
138 * MA 02110-1301, USA.
139 
140 * Bugs:
141 * {note_any_bugs_here}
142 *-
143 */
144 
145 #include <math.h>
146 
147 #include "pal.h"
148 #include "palmac.h"
149 
150 void palOapqk ( const char *type, double ob1, double ob2, const double aoprms[14],
151  double *rap, double *dap ) {
152 
153  /* breakpoint for fast/slow refraction algorithm:
154  * zd greater than arctan(4), (see palRefco routine)
155  * or vector z less than cosine(arctan(z)) = 1/sqrt(17) */
156  const double zbreak = 0.242535625;
157 
158  char c;
159  double c1,c2,sphi,cphi,st,ce,xaeo,yaeo,zaeo,v[3],
160  xmhdo,ymhdo,zmhdo,az,sz,zdo,tz,dref,zdt,
161  xaet,yaet,zaet,xmhda,ymhda,zmhda,diurab,f,hma;
162 
163  /* coordinate type */
164  c = type[0];
165 
166  /* coordinates */
167  c1 = ob1;
168  c2 = ob2;
169 
170  /* sin, cos of latitude */
171  sphi = aoprms[1];
172  cphi = aoprms[2];
173 
174  /* local apparent sidereal time */
175  st = aoprms[13];
176 
177  /* standardise coordinate type */
178  if (c == 'r' || c == 'R') {
179  c = 'r';
180  } else if (c == 'h' || c == 'H') {
181  c = 'h';
182  } else {
183  c = 'a';
184  }
185 
186  /* if az,zd convert to cartesian (s=0,e=90) */
187  if (c == 'a') {
188  ce = sin(c2);
189  xaeo = -cos(c1)*ce;
190  yaeo = sin(c1)*ce;
191  zaeo = cos(c2);
192  } else {
193 
194  /* if ra,dec convert to ha,dec */
195  if (c == 'r') {
196  c1 = st-c1;
197  }
198 
199  /* to cartesian -ha,dec */
200  palDcs2c( -c1, c2, v );
201  xmhdo = v[0];
202  ymhdo = v[1];
203  zmhdo = v[2];
204 
205  /* to cartesian az,el (s=0,e=90) */
206  xaeo = sphi*xmhdo-cphi*zmhdo;
207  yaeo = ymhdo;
208  zaeo = cphi*xmhdo+sphi*zmhdo;
209  }
210 
211  /* azimuth (s=0,e=90) */
212  if (xaeo != 0.0 || yaeo != 0.0) {
213  az = atan2(yaeo,xaeo);
214  } else {
215  az = 0.0;
216  }
217 
218  /* sine of observed zd, and observed zd */
219  sz = sqrt(xaeo*xaeo+yaeo*yaeo);
220  zdo = atan2(sz,zaeo);
221 
222  /*
223  * refraction
224  * ---------- */
225 
226  /* large zenith distance? */
227  if (zaeo >= zbreak) {
228 
229  /* fast algorithm using two constant model */
230  tz = sz/zaeo;
231  dref = (aoprms[10]+aoprms[11]*tz*tz)*tz;
232 
233  } else {
234 
235  /* rigorous algorithm for large zd */
236  palRefro(zdo,aoprms[4],aoprms[5],aoprms[6],aoprms[7],
237  aoprms[8],aoprms[0],aoprms[9],1e-8,&dref);
238  }
239 
240  zdt = zdo+dref;
241 
242  /* to cartesian az,zd */
243  ce = sin(zdt);
244  xaet = cos(az)*ce;
245  yaet = sin(az)*ce;
246  zaet = cos(zdt);
247 
248  /* cartesian az,zd to cartesian -ha,dec */
249  xmhda = sphi*xaet+cphi*zaet;
250  ymhda = yaet;
251  zmhda = -cphi*xaet+sphi*zaet;
252 
253  /* diurnal aberration */
254  diurab = -aoprms[3];
255  f = (1.0-diurab*ymhda);
256  v[0] = f*xmhda;
257  v[1] = f*(ymhda+diurab);
258  v[2] = f*zmhda;
259 
260  /* to spherical -ha,dec */
261  palDcc2s(v,&hma,dap);
262 
263  /* Right Ascension */
264  *rap = palDranrm(st+hma);
265 
266 }
void palOapqk(const char *type, double ob1, double ob2, const double aoprms[14], double *rap, double *dap)
Definition: palOapqk.c:150
void palRefro(double zobs, double hm, double tdk, double pmb, double rh, double wl, double phi, double tlr, double eps, double *ref)
Definition: palRefro.c:185
void palDcc2s(double v[3], double *a, double *b)
Definition: palOne2One.c:327
double palDranrm(double angle)
Definition: palOne2One.c:766
int type
void palDcs2c(double a, double b, double v[3])
Definition: palOne2One.c:368