File Coverage

palsrc/palRefv.c
Criterion Covered Total %
statement 17 17 100.0
branch 2 2 100.0
condition n/a
subroutine n/a
pod n/a
total 19 19 100.0


line stmt bran cond sub pod time code
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
128              
129 2           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 2           x = vu[0];
135 2           y = vu[1];
136 2           z1 = vu[2];
137              
138             /* Keep correction approximately constant below about 3 deg elevation */
139 2 100         z = DMAX(z1,0.05);
140              
141             /* One Newton-Raphson iteration */
142 2           zsq = z*z;
143 2           rsq = x*x+y*y;
144 2           r = sqrt(rsq);
145 2           wb = refb*rsq/zsq;
146 2           wt = (refa+wb)/(1.0+(refa+3.0*wb)*(zsq+rsq)/zsq);
147 2           d = wt*r/z;
148 2           cd = 1.0-d*d/2.0;
149 2           f = cd*(1.0-wt);
150              
151             /* Post-refraction x,y,z */
152 2           vr[0] = x*f;
153 2           vr[1] = y*f;
154 2           vr[2] = cd*(z+d*r)+(z1-z);
155 2           }