source: trunk/MagicSoft/slalib/refz.c@ 13961

Last change on this file since 13961 was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 4.2 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaRefz ( double zu, double refa, double refb, double *zr )
4/*
5** - - - - - - - -
6** s l a R e f z
7** - - - - - - - -
8**
9** Adjust an unrefracted zenith distance to include the effect of
10** atmospheric refraction, using the simple A tan z + B tan^3 z
11** model.
12**
13** Given:
14** zu double unrefracted zenith distance of the source (radian)
15** refa double A: tan z coefficient (radian)
16** refb double B: tan^3 z coefficient (radian)
17**
18** Returned:
19** *zr double refracted zenith distance (radian)
20**
21** Notes:
22**
23** 1 This routine applies the adjustment for refraction in the
24** opposite sense to the usual one - it takes an unrefracted
25** (in vacuo) position and produces an observed (refracted)
26** position, whereas the A tan Z + B tan^3 Z model strictly
27** applies to the case where an observed position is to have the
28** refraction removed. The unrefracted to refracted case is
29** harder, and requires an inverted form of the text-book
30** refraction models; the formula used here is based on the
31** Newton-Raphson method. For the utmost numerical consistency
32** with the refracted to unrefracted model, two iterations are
33** carried out, achieving agreement at the 1D-11 arcseconds level
34** for a ZD of 80 degrees. The inherent accuracy of the model
35** is, of course, far worse than this - see the documentation for
36** slaRefco for more information.
37**
38** 2 At ZD 83 degrees, the rapidly-worsening A tan Z + B tan^3 Z
39** model is abandoned and an empirical formula takes over. Over a
40** wide range of observer heights and corresponding temperatures and
41** pressures, the following levels of accuracy (arcsec) are achieved,
42** relative to numerical integration through a model atmosphere:
43**
44** ZR error
45**
46** 80 0.4
47** 81 0.8
48** 82 1.5
49** 83 3.2
50** 84 4.9
51** 85 5.8
52** 86 6.1
53** 87 7.1
54** 88 10
55** 89 20
56** 90 40
57** 91 100 } relevant only to
58** 92 200 } high-elevation sites
59**
60** The high-ZD model is scaled to match the normal model at the
61** transition point; there is no glitch.
62**
63**
64** 3 Beyond 93 deg zenith distance, the refraction is held at its
65** 93 deg value.
66**
67** 4 See also the routine slaRefv, which performs the adjustment in
68** Cartesian Az/El coordinates, and with the emphasis on speed
69** rather than numerical accuracy.
70**
71** Defined in slamac.h: DR2D
72**
73** Last revision: 4 June 1997
74**
75** Copyright P.T.Wallace. All rights reserved.
76*/
77{
78 double zu1, zl, s, c, t, tsq, tcu, ref, e, e2;
79
80/* Coefficients for high ZD model (used beyond ZD 83 deg */
81 const double c1 = 0.55445,
82 c2 = -0.01133,
83 c3 = 0.00202,
84 c4 = 0.28385,
85 c5 = 0.02390;
86
87/* Largest usable ZD (deg) */
88 const double z93 = 93.0;
89
90/* ZD at which one model hands over to the other (radians) */
91 const double z83 = 83.0 / DR2D;
92
93/* High-ZD-model prediction (deg) for that point */
94 const double ref83 = ( c1 + c2 * 7.0 + c3 * 49.0 ) /
95 ( 1.0 + c4 * 7.0 + c5 * 49.0 );
96
97
98/* Perform calculations for zu or 83 deg, whichever is smaller */
99 zu1 = gmin ( zu, z83 );
100
101/* Functions of ZD */
102 zl = zu1;
103 s = sin ( zl );
104 c = cos ( zl );
105 t = s / c;
106 tsq = t * t;
107 tcu = t * tsq;
108
109/* Refracted ZD (mathematically to better than 1mas at 70 deg) */
110 zl -= ( refa * t + refb * tcu )
111 / ( 1.0 + ( refa + 3.0 * refb * tsq ) / ( c * c ) );
112
113/* Further iteration */
114 s = sin ( zl );
115 c = cos ( zl );
116 t = s / c;
117 tsq = t * t;
118 tcu = t * tsq;
119 ref = zu1 - zl + ( zl - zu1 + refa * t + refb * tcu )
120 / ( 1.0 + ( refa + 3.0 * refb * tsq ) / ( c * c ) );
121
122/* Special handling for large zu */
123 if ( zu > zu1 ) {
124 e = 90.0 - gmin ( z93, zu * DR2D );
125 e2 = e * e;
126 ref = ( ref / ref83 ) * ( c1 + c2 * e + c3 * e2 ) /
127 ( 1.0 + c4 * e + c5 * e2 );
128 }
129
130/* Refracted ZD */
131 *zr = zu - ref;
132}
Note: See TracBrowser for help on using the repository browser.