| 1 | /*
|
|---|
| 2 | *+
|
|---|
| 3 | * Name:
|
|---|
| 4 | * palRefv
|
|---|
| 5 |
|
|---|
| 6 | * Purpose:
|
|---|
| 7 | * Adjust an unrefracted Cartesian vector to include the effect of atmospheric refraction
|
|---|
| 8 |
|
|---|
| 9 | * Language:
|
|---|
| 10 | * Starlink ANSI C
|
|---|
| 11 |
|
|---|
| 12 | * Type of Module:
|
|---|
| 13 | * Library routine
|
|---|
| 14 |
|
|---|
| 15 | * Invocation:
|
|---|
| 16 | * void palRefv ( double vu[3], double refa, double refb, double vr[3] );
|
|---|
| 17 |
|
|---|
| 18 | * Arguments:
|
|---|
| 19 | * vu[3] = double (Given)
|
|---|
| 20 | * Unrefracted position of the source (Az/El 3-vector)
|
|---|
| 21 | * refa = double (Given)
|
|---|
| 22 | * tan Z coefficient (radian)
|
|---|
| 23 | * refb = double (Given)
|
|---|
| 24 | * tan**3 Z coefficient (radian)
|
|---|
| 25 | * vr[3] = double (Returned)
|
|---|
| 26 | * Refracted position of the source (Az/El 3-vector)
|
|---|
| 27 |
|
|---|
| 28 | * Description:
|
|---|
| 29 | * Adjust an unrefracted Cartesian vector to include the effect of
|
|---|
| 30 | * atmospheric refraction, using the simple A tan Z + B tan**3 Z
|
|---|
| 31 | * model.
|
|---|
| 32 |
|
|---|
| 33 | * Authors:
|
|---|
| 34 | * TIMJ: Tim Jenness
|
|---|
| 35 | * PTW: Patrick Wallace
|
|---|
| 36 | * {enter_new_authors_here}
|
|---|
| 37 |
|
|---|
| 38 | * Notes:
|
|---|
| 39 | * - This routine applies the adjustment for refraction in the
|
|---|
| 40 | * opposite sense to the usual one - it takes an unrefracted
|
|---|
| 41 | * (in vacuo) position and produces an observed (refracted)
|
|---|
| 42 | * position, whereas the A tan Z + B tan**3 Z model strictly
|
|---|
| 43 | * applies to the case where an observed position is to have the
|
|---|
| 44 | * refraction removed. The unrefracted to refracted case is
|
|---|
| 45 | * harder, and requires an inverted form of the text-book
|
|---|
| 46 | * refraction models; the algorithm used here is equivalent to
|
|---|
| 47 | * one iteration of the Newton-Raphson method applied to the above
|
|---|
| 48 | * formula.
|
|---|
| 49 | *
|
|---|
| 50 | * - Though optimized for speed rather than precision, the present
|
|---|
| 51 | * routine achieves consistency with the refracted-to-unrefracted
|
|---|
| 52 | * A tan Z + B tan**3 Z model at better than 1 microarcsecond within
|
|---|
| 53 | * 30 degrees of the zenith and remains within 1 milliarcsecond to
|
|---|
| 54 | * beyond ZD 70 degrees. The inherent accuracy of the model is, of
|
|---|
| 55 | * course, far worse than this - see the documentation for palRefco
|
|---|
| 56 | * for more information.
|
|---|
| 57 | *
|
|---|
| 58 | * - At low elevations (below about 3 degrees) the refraction
|
|---|
| 59 | * correction is held back to prevent arithmetic problems and
|
|---|
| 60 | * wildly wrong results. For optical/IR wavelengths, over a wide
|
|---|
| 61 | * range of observer heights and corresponding temperatures and
|
|---|
| 62 | * pressures, the following levels of accuracy (arcsec, worst case)
|
|---|
| 63 | * are achieved, relative to numerical integration through a model
|
|---|
| 64 | * atmosphere:
|
|---|
| 65 | *
|
|---|
| 66 | * ZD error
|
|---|
| 67 | *
|
|---|
| 68 | * 80 0.7
|
|---|
| 69 | * 81 1.3
|
|---|
| 70 | * 82 2.5
|
|---|
| 71 | * 83 5
|
|---|
| 72 | * 84 10
|
|---|
| 73 | * 85 20
|
|---|
| 74 | * 86 55
|
|---|
| 75 | * 87 160
|
|---|
| 76 | * 88 360
|
|---|
| 77 | * 89 640
|
|---|
| 78 | * 90 1100
|
|---|
| 79 | * 91 1700 } relevant only to
|
|---|
| 80 | * 92 2600 } high-elevation sites
|
|---|
| 81 | *
|
|---|
| 82 | * The results for radio are slightly worse over most of the range,
|
|---|
| 83 | * becoming significantly worse below ZD=88 and unusable beyond
|
|---|
| 84 | * ZD=90.
|
|---|
| 85 | *
|
|---|
| 86 | * - See also the routine palRefz, which performs the adjustment to
|
|---|
| 87 | * the zenith distance rather than in Cartesian Az/El coordinates.
|
|---|
| 88 | * The present routine is faster than palRefz and, except very low down,
|
|---|
| 89 | * is equally accurate for all practical purposes. However, beyond
|
|---|
| 90 | * about ZD 84 degrees palRefz should be used, and for the utmost
|
|---|
| 91 | * accuracy iterative use of palRefro should be considered.
|
|---|
| 92 |
|
|---|
| 93 | * History:
|
|---|
| 94 | * 2014-07-15 (TIMJ):
|
|---|
| 95 | * Initial version. A direct copy of the Fortran SLA implementation.
|
|---|
| 96 | * Adapted with permission from the Fortran SLALIB library.
|
|---|
| 97 | * {enter_further_changes_here}
|
|---|
| 98 |
|
|---|
| 99 | * Copyright:
|
|---|
| 100 | * Copyright (C) 2014 Tim Jenness
|
|---|
| 101 | * Copyright (C) 2004 Patrick Wallace
|
|---|
| 102 | * All Rights Reserved.
|
|---|
| 103 |
|
|---|
| 104 | * Licence:
|
|---|
| 105 | * This program is free software; you can redistribute it and/or
|
|---|
| 106 | * modify it under the terms of the GNU General Public License as
|
|---|
| 107 | * published by the Free Software Foundation; either version 3 of
|
|---|
| 108 | * the License, or (at your option) any later version.
|
|---|
| 109 | *
|
|---|
| 110 | * This program is distributed in the hope that it will be
|
|---|
| 111 | * useful, but WITHOUT ANY WARRANTY; without even the implied
|
|---|
| 112 | * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|---|
| 113 | * PURPOSE. See the GNU General Public License for more details.
|
|---|
| 114 | *
|
|---|
| 115 | * You should have received a copy of the GNU General Public License
|
|---|
| 116 | * along with this program; if not, write to the Free Software
|
|---|
| 117 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|---|
| 118 | * MA 02110-1301, USA.
|
|---|
| 119 |
|
|---|
| 120 | * Bugs:
|
|---|
| 121 | * {note_any_bugs_here}
|
|---|
| 122 | *-
|
|---|
| 123 | */
|
|---|
| 124 |
|
|---|
| 125 | #include "pal.h"
|
|---|
| 126 | #include "palmac.h"
|
|---|
| 127 | #include <math.h>
|
|---|
| 128 |
|
|---|
| 129 | void palRefv ( double vu[3], double refa, double refb, double vr[3] ) {
|
|---|
| 130 |
|
|---|
| 131 | double x,y,z1,z,zsq,rsq,r,wb,wt,d,cd,f;
|
|---|
| 132 |
|
|---|
| 133 | /* Initial estimate = unrefracted vector */
|
|---|
| 134 | x = vu[0];
|
|---|
| 135 | y = vu[1];
|
|---|
| 136 | z1 = vu[2];
|
|---|
| 137 |
|
|---|
| 138 | /* Keep correction approximately constant below about 3 deg elevation */
|
|---|
| 139 | z = DMAX(z1,0.05);
|
|---|
| 140 |
|
|---|
| 141 | /* One Newton-Raphson iteration */
|
|---|
| 142 | zsq = z*z;
|
|---|
| 143 | rsq = x*x+y*y;
|
|---|
| 144 | r = sqrt(rsq);
|
|---|
| 145 | wb = refb*rsq/zsq;
|
|---|
| 146 | wt = (refa+wb)/(1.0+(refa+3.0*wb)*(zsq+rsq)/zsq);
|
|---|
| 147 | d = wt*r/z;
|
|---|
| 148 | cd = 1.0-d*d/2.0;
|
|---|
| 149 | f = cd*(1.0-wt);
|
|---|
| 150 |
|
|---|
| 151 | /* Post-refraction x,y,z */
|
|---|
| 152 | vr[0] = x*f;
|
|---|
| 153 | vr[1] = y*f;
|
|---|
| 154 | vr[2] = cd*(z+d*r)+(z1-z);
|
|---|
| 155 | }
|
|---|
