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