1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155

/*
*+
* Name:
* palRefv
* Purpose:
* Adjust an unrefracted Cartesian vector to include the effect of atmospheric refraction
* Language:
* Starlink ANSI C
* Type of Module:
* Library routine
* Invocation:
* void palRefv ( double vu[3], double refa, double refb, double vr[3] );
* Arguments:
* vu[3] = double (Given)
* Unrefracted position of the source (Az/El 3vector)
* refa = double (Given)
* tan Z coefficient (radian)
* refb = double (Given)
* tan**3 Z coefficient (radian)
* vr[3] = double (Returned)
* Refracted position of the source (Az/El 3vector)
* Description:
* Adjust an unrefracted Cartesian vector to include the effect of
* atmospheric refraction, using the simple A tan Z + B tan**3 Z
* model.
* Authors:
* TIMJ: Tim Jenness
* PTW: Patrick Wallace
* {enter_new_authors_here}
* Notes:
*  This routine applies the adjustment for refraction in the
* opposite sense to the usual one  it takes an unrefracted
* (in vacuo) position and produces an observed (refracted)
* position, whereas the A tan Z + B tan**3 Z model strictly
* applies to the case where an observed position is to have the
* refraction removed. The unrefracted to refracted case is
* harder, and requires an inverted form of the textbook
* refraction models; the algorithm used here is equivalent to
* one iteration of the NewtonRaphson method applied to the above
* formula.
*
*  Though optimized for speed rather than precision, the present
* routine achieves consistency with the refractedtounrefracted
* A tan Z + B tan**3 Z model at better than 1 microarcsecond within
* 30 degrees of the zenith and remains within 1 milliarcsecond to
* beyond ZD 70 degrees. The inherent accuracy of the model is, of
* course, far worse than this  see the documentation for palRefco
* for more information.
*
*  At low elevations (below about 3 degrees) the refraction
* correction is held back to prevent arithmetic problems and
* wildly wrong results. For optical/IR wavelengths, over a wide
* range of observer heights and corresponding temperatures and
* pressures, the following levels of accuracy (arcsec, worst case)
* are achieved, relative to numerical integration through a model
* atmosphere:
*
* ZD error
*
* 80 0.7
* 81 1.3
* 82 2.5
* 83 5
* 84 10
* 85 20
* 86 55
* 87 160
* 88 360
* 89 640
* 90 1100
* 91 1700 } relevant only to
* 92 2600 } highelevation sites
*
* The results for radio are slightly worse over most of the range,
* becoming significantly worse below ZD=88 and unusable beyond
* ZD=90.
*
*  See also the routine palRefz, which performs the adjustment to
* the zenith distance rather than in Cartesian Az/El coordinates.
* The present routine is faster than palRefz and, except very low down,
* is equally accurate for all practical purposes. However, beyond
* about ZD 84 degrees palRefz should be used, and for the utmost
* accuracy iterative use of palRefro should be considered.
* History:
* 20140715 (TIMJ):
* Initial version. A direct copy of the Fortran SLA implementation.
* Adapted with permission from the Fortran SLALIB library.
* {enter_further_changes_here}
* Copyright:
* Copyright (C) 2014 Tim Jenness
* Copyright (C) 2004 Patrick Wallace
* All Rights Reserved.
* Licence:
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation; either version 3 of
* the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be
* useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
* PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
* MA 021101301, USA.
* Bugs:
* {note_any_bugs_here}
*
*/
#include "pal.h"
#include "palmac.h"
#include <math.h>
void palRefv ( double vu[3], double refa, double refb, double vr[3] ) {
double x,y,z1,z,zsq,rsq,r,wb,wt,d,cd,f;
/* Initial estimate = unrefracted vector */
x = vu[0];
y = vu[1];
z1 = vu[2];
/* Keep correction approximately constant below about 3 deg elevation */
z = DMAX(z1,0.05);
/* One NewtonRaphson iteration */
zsq = z*z;
rsq = x*x+y*y;
r = sqrt(rsq);
wb = refb*rsq/zsq;
wt = (refa+wb)/(1.0+(refa+3.0*wb)*(zsq+rsq)/zsq);
d = wt*r/z;
cd = 1.0d*d/2.0;
f = cd*(1.0wt);
/* Postrefraction x,y,z */
vr[0] = x*f;
vr[1] = y*f;
vr[2] = cd*(z+d*r)+(z1z);
}
