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
palRefv(double vu[3],double refa,double refb,double vr[3])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 }
156