| 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 |
|
|---|
| 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 | }
|
|---|
