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