| 1 | /* | 
|---|
| 2 | *+ | 
|---|
| 3 | *  Name: | 
|---|
| 4 | *     palOap | 
|---|
| 5 |  | 
|---|
| 6 | *  Purpose: | 
|---|
| 7 | *     Observed to apparent place | 
|---|
| 8 |  | 
|---|
| 9 | *  Language: | 
|---|
| 10 | *     Starlink ANSI C | 
|---|
| 11 |  | 
|---|
| 12 | *  Type of Module: | 
|---|
| 13 | *     Library routine | 
|---|
| 14 |  | 
|---|
| 15 | *  Invocation: | 
|---|
| 16 | *     void palOap ( const char *type, double ob1, double ob2, double date, | 
|---|
| 17 | *                   double dut, double elongm, double phim, double hm, | 
|---|
| 18 | *                   double xp, double yp, double tdk, double pmb, | 
|---|
| 19 | *                   double rh, double wl, double tlr, | 
|---|
| 20 | *                   double *rap, double *dap ); | 
|---|
| 21 |  | 
|---|
| 22 | *  Arguments: | 
|---|
| 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 | *     date = double (Given) | 
|---|
| 30 | *        UTC date/time (Modified Julian Date, JD-2400000.5) | 
|---|
| 31 | *     dut = double (Given) | 
|---|
| 32 | *        delta UT: UT1-UTC (UTC seconds) | 
|---|
| 33 | *     elongm = double (Given) | 
|---|
| 34 | *        Mean longitude of the observer (radians, east +ve) | 
|---|
| 35 | *     phim = double (Given) | 
|---|
| 36 | *        Mean geodetic latitude of the observer (radians) | 
|---|
| 37 | *     hm = double (Given) | 
|---|
| 38 | *        Observer's height above sea level (metres) | 
|---|
| 39 | *     xp = double (Given) | 
|---|
| 40 | *        Polar motion x-coordinates (radians) | 
|---|
| 41 | *     yp = double (Given) | 
|---|
| 42 | *        Polar motion y-coordinates (radians) | 
|---|
| 43 | *     tdk = double (Given) | 
|---|
| 44 | *        Local ambient temperature (K; std=273.15) | 
|---|
| 45 | *     pmb = double (Given) | 
|---|
| 46 | *        Local atmospheric pressure (mb; std=1013.25) | 
|---|
| 47 | *     rh = double (Given) | 
|---|
| 48 | *        Local relative humidity (in the range 0.0-1.0) | 
|---|
| 49 | *     wl = double (Given) | 
|---|
| 50 | *        Effective wavelength (micron, e.g. 0.55) | 
|---|
| 51 | *     tlr = double (Given) | 
|---|
| 52 | *        Tropospheric laps rate (K/metre, e.g. 0.0065) | 
|---|
| 53 | *     rap = double * (Given) | 
|---|
| 54 | *        Geocentric apparent right ascension | 
|---|
| 55 | *     dap = double * (Given) | 
|---|
| 56 | *        Geocentric apparent declination | 
|---|
| 57 |  | 
|---|
| 58 | *  Description: | 
|---|
| 59 | *     Observed to apparent place. | 
|---|
| 60 |  | 
|---|
| 61 | *  Authors: | 
|---|
| 62 | *     PTW: Patrick T. Wallace | 
|---|
| 63 | *     TIMJ: Tim Jenness (JAC, Hawaii) | 
|---|
| 64 | *     {enter_new_authors_here} | 
|---|
| 65 |  | 
|---|
| 66 | *  Notes: | 
|---|
| 67 | *     - Only the first character of the TYPE argument is significant. | 
|---|
| 68 | *     'R' or 'r' indicates that OBS1 and OBS2 are the observed right | 
|---|
| 69 | *     ascension and declination;  'H' or 'h' indicates that they are | 
|---|
| 70 | *     hour angle (west +ve) and declination;  anything else ('A' or | 
|---|
| 71 | *     'a' is recommended) indicates that OBS1 and OBS2 are azimuth | 
|---|
| 72 | *     (north zero, east 90 deg) and zenith distance.  (Zenith | 
|---|
| 73 | *     distance is used rather than elevation in order to reflect the | 
|---|
| 74 | *     fact that no allowance is made for depression of the horizon.) | 
|---|
| 75 | * | 
|---|
| 76 | *     - The accuracy of the result is limited by the corrections for | 
|---|
| 77 | *     refraction.  Providing the meteorological parameters are | 
|---|
| 78 | *     known accurately and there are no gross local effects, the | 
|---|
| 79 | *     predicted apparent RA,Dec should be within about 0.1 arcsec | 
|---|
| 80 | *     for a zenith distance of less than 70 degrees.  Even at a | 
|---|
| 81 | *     topocentric zenith distance of 90 degrees, the accuracy in | 
|---|
| 82 | *     elevation should be better than 1 arcmin;  useful results | 
|---|
| 83 | *     are available for a further 3 degrees, beyond which the | 
|---|
| 84 | *     palRefro routine returns a fixed value of the refraction. | 
|---|
| 85 | *     The complementary routines palAop (or palAopqk) and palOap | 
|---|
| 86 | *     (or palOapqk) are self-consistent to better than 1 micro- | 
|---|
| 87 | *     arcsecond all over the celestial sphere. | 
|---|
| 88 | * | 
|---|
| 89 | *     - It is advisable to take great care with units, as even | 
|---|
| 90 | *     unlikely values of the input parameters are accepted and | 
|---|
| 91 | *     processed in accordance with the models used. | 
|---|
| 92 | * | 
|---|
| 93 | *     - "Observed" Az,El means the position that would be seen by a | 
|---|
| 94 | *     perfect theodolite located at the observer.  This is | 
|---|
| 95 | *     related to the observed HA,Dec via the standard rotation, using | 
|---|
| 96 | *     the geodetic latitude (corrected for polar motion), while the | 
|---|
| 97 | *     observed HA and RA are related simply through the local | 
|---|
| 98 | *     apparent ST.  "Observed" RA,Dec or HA,Dec thus means the | 
|---|
| 99 | *     position that would be seen by a perfect equatorial located | 
|---|
| 100 | *     at the observer and with its polar axis aligned to the | 
|---|
| 101 | *     Earth's axis of rotation (n.b. not to the refracted pole). | 
|---|
| 102 | *     By removing from the observed place the effects of | 
|---|
| 103 | *     atmospheric refraction and diurnal aberration, the | 
|---|
| 104 | *     geocentric apparent RA,Dec is obtained. | 
|---|
| 105 | * | 
|---|
| 106 | *     - Frequently, mean rather than apparent RA,Dec will be required, | 
|---|
| 107 | *     in which case further transformations will be necessary.  The | 
|---|
| 108 | *     palAmp etc routines will convert the apparent RA,Dec produced | 
|---|
| 109 | *     by the present routine into an "FK5" (J2000) mean place, by | 
|---|
| 110 | *     allowing for the Sun's gravitational lens effect, annual | 
|---|
| 111 | *     aberration, nutation and precession.  Should "FK4" (1950) | 
|---|
| 112 | *     coordinates be needed, the routines palFk524 etc will also | 
|---|
| 113 | *     need to be applied. | 
|---|
| 114 | * | 
|---|
| 115 | *     - To convert to apparent RA,Dec the coordinates read from a | 
|---|
| 116 | *     real telescope, corrections would have to be applied for | 
|---|
| 117 | *     encoder zero points, gear and encoder errors, tube flexure, | 
|---|
| 118 | *     the position of the rotator axis and the pointing axis | 
|---|
| 119 | *     relative to it, non-perpendicularity between the mounting | 
|---|
| 120 | *     axes, and finally for the tilt of the azimuth or polar axis | 
|---|
| 121 | *     of the mounting (with appropriate corrections for mount | 
|---|
| 122 | *     flexures).  Some telescopes would, of course, exhibit other | 
|---|
| 123 | *     properties which would need to be accounted for at the | 
|---|
| 124 | *     appropriate point in the sequence. | 
|---|
| 125 | * | 
|---|
| 126 | *     - This routine takes time to execute, due mainly to the rigorous | 
|---|
| 127 | *     integration used to evaluate the refraction.  For processing | 
|---|
| 128 | *     multiple stars for one location and time, call palAoppa once | 
|---|
| 129 | *     followed by one call per star to palOapqk.  Where a range of | 
|---|
| 130 | *     times within a limited period of a few hours is involved, and the | 
|---|
| 131 | *     highest precision is not required, call palAoppa once, followed | 
|---|
| 132 | *     by a call to palAoppat each time the time changes, followed by | 
|---|
| 133 | *     one call per star to palOapqk. | 
|---|
| 134 | * | 
|---|
| 135 | *     - The DATE argument is UTC expressed as an MJD.  This is, strictly | 
|---|
| 136 | *     speaking, wrong, because of leap seconds.  However, as long as | 
|---|
| 137 | *     the delta UT and the UTC are consistent there are no | 
|---|
| 138 | *     difficulties, except during a leap second.  In this case, the | 
|---|
| 139 | *     start of the 61st second of the final minute should begin a new | 
|---|
| 140 | *     MJD day and the old pre-leap delta UT should continue to be used. | 
|---|
| 141 | *     As the 61st second completes, the MJD should revert to the start | 
|---|
| 142 | *     of the day as, simultaneously, the delta UTC changes by one | 
|---|
| 143 | *     second to its post-leap new value. | 
|---|
| 144 | * | 
|---|
| 145 | *     - The delta UT (UT1-UTC) is tabulated in IERS circulars and | 
|---|
| 146 | *     elsewhere.  It increases by exactly one second at the end of | 
|---|
| 147 | *     each UTC leap second, introduced in order to keep delta UT | 
|---|
| 148 | *     within +/- 0.9 seconds. | 
|---|
| 149 | * | 
|---|
| 150 | *     - IMPORTANT -- TAKE CARE WITH THE LONGITUDE SIGN CONVENTION. | 
|---|
| 151 | *     The longitude required by the present routine is east-positive, | 
|---|
| 152 | *     in accordance with geographical convention (and right-handed). | 
|---|
| 153 | *     In particular, note that the longitudes returned by the | 
|---|
| 154 | *     palOBS routine are west-positive, following astronomical | 
|---|
| 155 | *     usage, and must be reversed in sign before use in the present | 
|---|
| 156 | *     routine. | 
|---|
| 157 | * | 
|---|
| 158 | *     - The polar coordinates XP,YP can be obtained from IERS | 
|---|
| 159 | *     circulars and equivalent publications.  The maximum amplitude | 
|---|
| 160 | *     is about 0.3 arcseconds.  If XP,YP values are unavailable, | 
|---|
| 161 | *     use XP=YP=0D0.  See page B60 of the 1988 Astronomical Almanac | 
|---|
| 162 | *     for a definition of the two angles. | 
|---|
| 163 | * | 
|---|
| 164 | *     - The height above sea level of the observing station, HM, | 
|---|
| 165 | *     can be obtained from the Astronomical Almanac (Section J | 
|---|
| 166 | *     in the 1988 edition), or via the routine palOBS.  If P, | 
|---|
| 167 | *     the pressure in millibars, is available, an adequate | 
|---|
| 168 | *     estimate of HM can be obtained from the expression | 
|---|
| 169 | * | 
|---|
| 170 | *            HM ~ -29.3*TSL*LOG(P/1013.25). | 
|---|
| 171 | * | 
|---|
| 172 | *     where TSL is the approximate sea-level air temperature in K | 
|---|
| 173 | *     (see Astrophysical Quantities, C.W.Allen, 3rd edition, | 
|---|
| 174 | *     section 52).  Similarly, if the pressure P is not known, | 
|---|
| 175 | *     it can be estimated from the height of the observing | 
|---|
| 176 | *     station, HM, as follows: | 
|---|
| 177 | * | 
|---|
| 178 | *            P ~ 1013.25*EXP(-HM/(29.3*TSL)). | 
|---|
| 179 | * | 
|---|
| 180 | *     Note, however, that the refraction is nearly proportional to the | 
|---|
| 181 | *     pressure and that an accurate P value is important for precise | 
|---|
| 182 | *     work. | 
|---|
| 183 | * | 
|---|
| 184 | *     - The azimuths etc. used by the present routine are with respect | 
|---|
| 185 | *     to the celestial pole.  Corrections from the terrestrial pole | 
|---|
| 186 | *     can be computed using palPolmo. | 
|---|
| 187 |  | 
|---|
| 188 | *  History: | 
|---|
| 189 | *     2012-08-27 (TIMJ): | 
|---|
| 190 | *        Initial version, copied from Fortran SLA | 
|---|
| 191 | *        Adapted with permission from the Fortran SLALIB library. | 
|---|
| 192 | *     {enter_further_changes_here} | 
|---|
| 193 |  | 
|---|
| 194 | *  Copyright: | 
|---|
| 195 | *     Copyright (C) 2005 Patrick T. Wallace | 
|---|
| 196 | *     Copyright (C) 2012 Science and Technology Facilities Council. | 
|---|
| 197 | *     All Rights Reserved. | 
|---|
| 198 |  | 
|---|
| 199 | *  Licence: | 
|---|
| 200 | *     This program is free software; you can redistribute it and/or | 
|---|
| 201 | *     modify it under the terms of the GNU General Public License as | 
|---|
| 202 | *     published by the Free Software Foundation; either version 3 of | 
|---|
| 203 | *     the License, or (at your option) any later version. | 
|---|
| 204 | * | 
|---|
| 205 | *     This program is distributed in the hope that it will be | 
|---|
| 206 | *     useful, but WITHOUT ANY WARRANTY; without even the implied | 
|---|
| 207 | *     warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR | 
|---|
| 208 | *     PURPOSE. See the GNU General Public License for more details. | 
|---|
| 209 | * | 
|---|
| 210 | *     You should have received a copy of the GNU General Public License | 
|---|
| 211 | *     along with this program; if not, write to the Free Software | 
|---|
| 212 | *     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, | 
|---|
| 213 | *     MA 02110-1301, USA. | 
|---|
| 214 |  | 
|---|
| 215 | *  Bugs: | 
|---|
| 216 | *     {note_any_bugs_here} | 
|---|
| 217 | *- | 
|---|
| 218 | */ | 
|---|
| 219 |  | 
|---|
| 220 | #include "pal.h" | 
|---|
| 221 |  | 
|---|
| 222 | void palOap ( const char *type, double ob1, double ob2, double date, | 
|---|
| 223 | double dut, double elongm, double phim, double hm, | 
|---|
| 224 | double xp, double yp, double tdk, double pmb, | 
|---|
| 225 | double rh, double wl, double tlr, | 
|---|
| 226 | double *rap, double *dap ) { | 
|---|
| 227 |  | 
|---|
| 228 | double aoprms[14]; | 
|---|
| 229 |  | 
|---|
| 230 | palAoppa(date,dut,elongm,phim,hm,xp,yp,tdk,pmb,rh,wl,tlr, | 
|---|
| 231 | aoprms); | 
|---|
| 232 | palOapqk(type,ob1,ob2,aoprms,rap,dap); | 
|---|
| 233 |  | 
|---|
| 234 | } | 
|---|