File: palRefz.c

package info (click to toggle)
starlink-pal 0.9.8-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 1,808 kB
  • sloc: ansic: 6,689; makefile: 128; sh: 81
file content (191 lines) | stat: -rw-r--r-- 6,141 bytes parent folder | download | duplicates (3)
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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
/*
*+
*  Name:
*     palRefz

*  Purpose:
*     Adjust unrefracted zenith distance

*  Language:
*     Starlink ANSI C

*  Type of Module:
*     Library routine

*  Invocation:
*     void palRefz ( double zu, double refa, double refb, double *zr );

*  Arguments:
*     zu = double (Given)
*         Unrefracted zenith distance of the source (radians)
*     refa = double (Given)
*         tan Z coefficient (radians)
*     refb = double (Given)
*         tan**3 Z coefficient (radian)
*     zr = double * (Returned)
*         Refracted zenith distance (radians)

*  Description:
*     Adjust an unrefracted zenith distance to include the effect of
*     atmospheric refraction, using the simple A tan Z + B tan**3 Z
*     model (plus special handling for large ZDs).

*  Authors:
*     PTW: Patrick T. Wallace
*     TIMJ: Tim Jenness (JAC, Hawaii)
*     {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 text-book
*     refraction models;  the formula used here is based on the
*     Newton-Raphson method.  For the utmost numerical consistency
*     with the refracted to unrefracted model, two iterations are
*     carried out, achieving agreement at the 1D-11 arcseconds level
*     for a ZD of 80 degrees.  The inherent accuracy of the model
*     is, of course, far worse than this - see the documentation for
*     palRefco for more information.
*
*     - At ZD 83 degrees, the rapidly-worsening A tan Z + B tan^3 Z
*     model is abandoned and an empirical formula takes over.  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:
*
*              ZR    error
*
*              80      0.7
*              81      1.3
*              82      2.4
*              83      4.7
*              84      6.2
*              85      6.4
*              86      8
*              87     10
*              88     15
*              89     30
*              90     60
*              91    150         } relevant only to
*              92    400         } high-elevation sites
*
*     For radio wavelengths the errors are typically 50% larger than
*     the optical figures and by ZD 85 deg are twice as bad, worsening
*     rapidly below that.  To maintain 1 arcsec accuracy down to ZD=85
*     at the Green Bank site, Condon (2004) has suggested amplifying
*     the amount of refraction predicted by palRefz below 10.8 deg
*     elevation by the factor (1+0.00195*(10.8-E_t)), where E_t is the
*     unrefracted elevation in degrees.
*
*     The high-ZD model is scaled to match the normal model at the
*     transition point;  there is no glitch.
*
*     - Beyond 93 deg zenith distance, the refraction is held at its
*     93 deg value.
*
*     - See also the routine palRefv, which performs the adjustment in
*     Cartesian Az/El coordinates, and with the emphasis on speed
*     rather than numerical accuracy.

*  References:
*     Condon,J.J., Refraction Corrections for the GBT, PTCS/PN/35.2,
*     NRAO Green Bank, 2004.

*  History:
*     2012-08-24 (TIMJ):
*        Initial version, ported directly from Fortran SLA
*        Adapted with permission from the Fortran SLALIB library.
*     {enter_further_changes_here}

*  Copyright:
*     Copyright (C) 2004 Rutherford Appleton Laboratory
*     Copyright (C) 2012 Science and Technology Facilities Council.
*     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 02110-1301, USA.

*  Bugs:
*     {note_any_bugs_here}
*-
*/

#include <math.h>

#include "pal.h"
#include "palmac.h"

void palRefz ( double zu, double refa, double refb, double *zr ) {

  /* Constants */

  /* Largest usable ZD (deg) */
  const double D93 = 93.0;

  /* ZD at which one model hands over to the other (radians) */
  const double Z83 = 83.0 * PAL__DD2R;

  /* coefficients for high ZD model (used beyond ZD 83 deg) */
  const double C1 = +0.55445;
  const double C2 = -0.01133;
  const double C3 = +0.00202;
  const double C4 = +0.28385;
  const double C5 = +0.02390;

  /* High-ZD-model prefiction (deg) for that point */
  const double REF83 = (C1+C2*7.0+C3*49.0)/(1.0+C4*7.0+C5*49.0);

  double zu1,zl,s,c,t,tsq,tcu,ref,e,e2;

  /*  perform calculations for zu or 83 deg, whichever is smaller */
  zu1 = DMIN(zu,Z83);

  /*  functions of ZD */
  zl = zu1;
  s = sin(zl);
  c = cos(zl);
  t = s/c;
  tsq = t*t;
  tcu = t*tsq;

  /*  refracted zd (mathematically to better than 1 mas at 70 deg) */
  zl = zl-(refa*t+refb*tcu)/(1.0+(refa+3.0*refb*tsq)/(c*c));

  /*  further iteration */
  s = sin(zl);
  c = cos(zl);
  t = s/c;
  tsq = t*t;
  tcu = t*tsq;
  ref = zu1-zl+
    (zl-zu1+refa*t+refb*tcu)/(1.0+(refa+3.0*refb*tsq)/(c*c));

  /*  special handling for large zu */
  if (zu > zu1) {
    e = 90.0-DMIN(D93,zu*PAL__DR2D);
    e2 = e*e;
    ref = (ref/REF83)*(C1+C2*e+C3*e2)/(1.0+C4*e+C5*e2);
  }

  /*  return refracted zd */
  *zr = zu-ref;

}