source: trunk/MagicSoft/slalib/refro.c@ 7489

Last change on this file since 7489 was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 16.4 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3
4static void atmt ( double, double, double, double, double, double,
5 double, double, double, double, double, double,
6 double*, double*, double* );
7static void atms ( double, double, double, double, double,
8 double*, double* );
9
10void slaRefro ( double zobs, double hm, double tdk, double pmb,
11 double rh, double wl, double phi, double tlr,
12 double eps, double *ref )
13/*
14** - - - - - - - - -
15** s l a R e f r o
16** - - - - - - - - -
17**
18** Atmospheric refraction for radio and optical/IR wavelengths.
19**
20** Given:
21** zobs double observed zenith distance of the source (radian)
22** hm double height of the observer above sea level (metre)
23** tdk double ambient temperature at the observer (deg K)
24** pmb double pressure at the observer (millibar)
25** rh double relative humidity at the observer (range 0-1)
26** wl double effective wavelength of the source (micrometre)
27** phi double latitude of the observer (radian, astronomical)
28** tlr double tropospheric lapse rate (degK/metre)
29** eps double precision required to terminate iteration (radian)
30**
31** Returned:
32** ref double refraction: in vacuo ZD minus observed ZD (radian)
33**
34** Notes:
35**
36** 1 A suggested value for the tlr argument is 0.0065. The
37** refraction is significantly affected by tlr, and if studies
38** of the local atmosphere have been carried out a better tlr
39** value may be available.
40**
41** 2 A suggested value for the eps argument is 1e-8. The result is
42** usually at least two orders of magnitude more computationally
43** precise than the supplied eps value.
44**
45** 3 The routine computes the refraction for zenith distances up
46** to and a little beyond 90 deg using the method of Hohenkerk
47** and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted
48** in the Explanatory Supplement, 1992 edition - see section 3.281).
49**
50** 4 The C code is an adaptation of the Fortran optical/IR refraction
51** subroutine AREF of C.Hohenkerk (HMNAO, September 1984), with
52** extensions to support the radio case. The following modifications
53** to the original HMNAO optical/IR refraction algorithm have been made:
54**
55** . The angle arguments have been changed to radians.
56**
57** . Any value of zobs is allowed (see note 6, below).
58**
59** . Other argument values have been limited to safe values.
60**
61** . Murray's values for the gas constants have been used
62** (Vectorial Astrometry, Adam Hilger, 1983).
63**
64** . The numerical integration phase has been rearranged for
65** extra clarity.
66**
67** . A better model for Ps(T) has been adopted (taken from
68** Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982).
69**
70** . More accurate expressions for Pwo have been adopted
71** (again from Gill 1982).
72**
73** . Provision for radio wavelengths has been added using
74** expressions devised by A.T.Sinclair, RGO (private
75** communication 1989), based on the Essen & Froome
76** refractivity formula adopted in Resolution 1 of the
77** 13th International Geodesy Association General Assembly
78** (Bulletin Geodesique 70 p390, 1963).
79**
80** . Various small changes have been made to gain speed.
81**
82** None of the changes significantly affects the optical/IR results
83** with respect to the algorithm given in the 1992 Explanatory
84** Supplement. For example, at 70 deg zenith distance the present
85** routine agrees with the ES algorithm to better than 0.05 arcsec
86** for any reasonable combination of parameters. However, the
87** improved water-vapour expressions do make a significant difference
88** in the radio band, at 70 deg zenith distance reaching almost
89** 4 arcsec for a hot, humid, low-altitude site during a period of
90** low pressure.
91**
92** 5 The radio refraction is chosen by specifying wl > 100 micrometres.
93** Because the algorithm takes no account of the ionosphere, the
94** accuracy deteriorates at low frequencies, below about 30 MHz.
95**
96** 6 Before use, the value of zobs is expressed in the range +/- pi.
97** If this ranged zobs is -ve, the result ref is computed from its
98** absolute value before being made -ve to match. In addition, if
99** it has an absolute value greater than 93 deg, a fixed ref value
100** equal to the result for zobs = 93 deg is returned, appropriately
101** signed.
102**
103** 7 As in the original Hohenkerk and Sinclair algorithm, fixed
104** values of the water vapour polytrope exponent, the height of the
105** tropopause, and the height at which refraction is negligible are
106** used.
107**
108** 8 The radio refraction has been tested against work done by
109** Iain Coulson, JACH, (private communication 1995) for the
110** James Clerk Maxwell Telescope, Mauna Kea. For typical conditions,
111** agreement at the 0.1 arcsec level is achieved for moderate ZD,
112** worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg. At hot and
113** humid sea-level sites the accuracy will not be as good.
114**
115** 9 It should be noted that the relative humidity rh is formally
116** defined in terms of "mixing ratio" rather than pressures or
117** densities as is often stated. It is the mass of water per unit
118** mass of dry air divided by that for saturated air at the same
119** temperature and pressure (see Gill 1982).
120**
121** Called: slaDrange, atmt, atms
122**
123** Defined in slamac.h: TRUE, FALSE
124**
125** Last revision: 3 June 1997
126**
127** Copyright P.T.Wallace. All rights reserved.
128*/
129{
130/* Fixed parameters */
131
132 static double d93 = 1.623156204; /* 93 degrees in radians */
133 static double gcr = 8314.32; /* Universal gas constant */
134 static double dmd = 28.9644; /* Molecular weight of dry air */
135 static double dmw = 18.0152; /* Molecular weight of water
136 vapour */
137 static double s = 6378120.0; /* Mean Earth radius (metre) */
138 static double delta = 18.36; /* Exponent of temperature
139 dependence of water vapour
140 pressure */
141 static double ht = 11000.0; /* Height of tropopause (metre) */
142 static double hs = 80000.0; /* Upper limit for refractive
143 effects (metre) */
144
145/* Variables used when calling the internal routine atmt */
146 double robs; /* height of observer from centre of Earth (metre) */
147 double tdkok; /* temperature at the observer (deg K) */
148 double alpha; /* alpha | */
149 double gamm2; /* gamma minus 2 | see ES */
150 double delm2; /* delta minus 2 | */
151 double c1,c2,c3,c4,c5,c6; /* various */
152
153/* Variables used when calling the internal routine atms */
154 double rt; /* height of tropopause from centre of Earth (metre) */
155 double tt; /* temperature at the tropopause (deg k) */
156 double dnt; /* refractive index at the tropopause */
157 double gamal; /* constant of the atmospheric model = g*md/r */
158
159 int is, k, n, i, j, optic;
160 double zobs1, zobs2, hmok, pmbok, rhok, wlok, tol, wlsq, gb,
161 a, gamma, tdc, psat, pwo, w, tempo, dn0, rdndr0, sk0,
162 f0, rdndrt, zt, ft, dnts, rdndrp, zts, fts, rs,
163 dns, rdndrs, zs, fs, refold, z0, zrange, fb, ff, fo,
164 fe, h, r, sz, rg, dr, tg, dn, rdndr, t, f, refp, reft;
165
166/* The refraction integrand */
167#define refi(R,DN,RDNDR) ((RDNDR)/(DN+RDNDR));
168
169
170
171/* Transform zobs into the normal range. */
172 zobs1 = slaDrange ( zobs );
173 zobs2 = fabs ( zobs1 );
174 zobs2 = gmin ( zobs2, d93 );
175
176/* Keep other arguments within safe bounds. */
177 hmok = gmax ( hm, -1000.0 );
178 hmok = gmin ( hmok, 10000.0 );
179 tdkok = gmax ( tdk, 100.0 );
180 tdkok = gmin ( tdkok, 500.0 );
181 pmbok = gmax ( pmb, 0.0 );
182 pmbok = gmin ( pmbok, 10000.0 );
183 rhok = gmax ( rh, 0.0 );
184 rhok = gmin ( rhok, 1.0 );
185 wlok = gmax ( wl, 0.1 );
186 alpha = fabs ( tlr );
187 alpha = gmax ( alpha, 0.001 );
188 alpha = gmin ( alpha, 0.01 );
189
190/* Tolerance for iteration. */
191 w = fabs ( eps );
192 w = gmax ( w, 1e-12 );
193 tol = gmin ( w, 0.1 ) / 2.0;
194
195/* Decide whether optical/IR or radio case - switch at 100 microns. */
196 optic = ( wlok <= 100.0 );
197
198/* Set up model atmosphere parameters defined at the observer. */
199 wlsq = wlok * wlok;
200 gb = 9.784 * ( 1.0 - 0.0026 * cos ( 2.0 * phi ) - 2.8e-7 * hmok );
201 a = ( optic ) ?
202 ( ( 287.604 + 1.6288 / wlsq + 0.0136 / ( wlsq * wlsq ) )
203 * 273.15 / 1013.25 ) * 1e-6
204 :
205 77.624e-6;
206 gamal = gb * dmd / gcr;
207 gamma = gamal / alpha;
208 gamm2 = gamma - 2.0;
209 delm2 = delta - 2.0;
210 tdc = tdkok - 273.15;
211 psat = pow ( 10.0, ( 0.7859 + 0.03477 * tdc ) /
212 ( 1.0 + 0.00412 * tdc ) ) *
213 ( 1.0 + pmbok * ( 4.5e-6 + 6e-10 * tdc * tdc ) );
214 pwo = ( pmbok > 0.0 ) ?
215 rhok * psat / ( 1.0 - ( 1.0 - rhok ) * psat / pmbok ) :
216 0.0;
217 w = pwo * ( 1.0 - dmw / dmd ) * gamma / ( delta - gamma );
218 c1 = a * ( pmbok + w ) / tdkok;
219 c2 = ( a * w + ( optic ? 11.2684e-6 : 12.92e-6 ) * pwo ) / tdkok;
220 c3 = ( gamma - 1.0 ) * alpha * c1 / tdkok;
221 c4 = ( delta - 1.0 ) * alpha * c2 / tdkok;
222 c5 = optic ? 0.0 : 371897e-6 * pwo / tdkok;
223 c6 = c5 * delm2 * alpha / ( tdkok * tdkok );
224
225/* Conditions at the observer. */
226 robs = s + hmok;
227 atmt ( robs, tdkok, alpha, gamm2, delm2, c1, c2, c3, c4, c5, c6, robs,
228 &tempo, &dn0, &rdndr0 );
229 sk0 = dn0 * robs * sin ( zobs2 );
230 f0 = refi ( robs, dn0, rdndr0 );
231
232/* Conditions at the tropopause in the troposphere. */
233 rt = s + ht;
234 atmt ( robs, tdkok, alpha, gamm2, delm2, c1, c2, c3, c4, c5, c6, rt,
235 &tt, &dnt, &rdndrt );
236 zt = asin ( sk0 / ( rt * dnt ) );
237 ft = refi ( rt, dnt, rdndrt );
238
239/* Conditions at the tropopause in the stratosphere. */
240 atms ( rt, tt, dnt, gamal, rt, &dnts, &rdndrp );
241 zts = asin ( sk0 / ( rt * dnts ) );
242 fts = refi ( rt, dnts, rdndrp );
243
244/* Conditions at the stratosphere limit. */
245 rs = s + hs;
246 atms ( rt, tt, dnt, gamal, rs, &dns, &rdndrs );
247 zs = asin ( sk0 / ( rs * dns ) );
248 fs = refi ( rs, dns, rdndrs );
249
250/*
251** Integrate the refraction integral in two parts; first in the
252** troposphere (k=1), then in the stratosphere (k=2).
253*/
254
255/* Initialize previous refraction to ensure at least two iterations. */
256 refold = 1e6;
257
258/*
259** Start off with 8 strips for the troposphere integration, and then
260** use the final troposphere value for the stratosphere integration,
261** which tends to need more strips.
262*/
263 is = 8;
264
265/* Troposphere then stratosphere. */
266 for ( k = 1; k <= 2; k++ ) {
267
268 /* Start z, z range, and start and end values. */
269 if ( k == 1 ) {
270 z0 = zobs2;
271 zrange = zt - z0;
272 fb = f0;
273 ff = ft;
274 } else {
275 z0 = zts;
276 zrange = zs - z0;
277 fb = fts;
278 ff = fs;
279 }
280
281 /* Sums of odd and even values. */
282 fo = 0.0;
283 fe = 0.0;
284
285 /* First time through the loop we have to do every point. */
286 n = 1;
287
288 /* Start of iteration loop (terminates at specified precision). */
289 for ( ; ; ) {
290
291 /* Strip width */
292 h = zrange / (double) is;
293
294 /* Initialize distance from Earth centre for quadrature pass. */
295 r = ( k == 1 ) ? robs : rt;
296
297 /* One pass (no need to compute evens after first time). */
298 for ( i = 1; i < is; i += n ) {
299
300 /* Sine of observed zenith distance. */
301 sz = sin ( z0 + h * (double) i );
302
303 /* Find r (to nearest metre, maximum four iterations). */
304 if ( sz > 1e-20 ) {
305 w = sk0 / sz;
306 rg = r;
307 j = 0;
308 do {
309 if ( k == 1 ) {
310 atmt ( robs, tdkok, alpha, gamm2, delm2,
311 c1, c2, c3, c4, c5, c6, rg,
312 &tg, &dn, &rdndr );
313 } else {
314 atms ( rt, tt, dnt, gamal, rg, &dn, &rdndr );
315 }
316 dr = ( rg * dn - w ) / ( dn + rdndr );
317 rg -= dr;
318 } while ( fabs ( dr ) > 1.0 && j++ <= 4 );
319 r = rg;
320 }
321
322 /* Find refractive index and integrand at r. */
323 if ( k == 1 ) {
324 atmt ( robs, tdkok, alpha, gamm2, delm2,
325 c1, c2, c3, c4, c5, c6, r,
326 &t, &dn, &rdndr );
327 } else {
328 atms ( rt, tt, dnt, gamal, r, &dn, &rdndr );
329 }
330 f = refi ( r, dn, rdndr );
331
332 /* Accumulate odd and (first time only) even values. */
333 if ( n == 1 && i % 2 == 0 ) {
334 fe += f;
335 } else {
336 fo += f;
337 }
338 }
339
340 /* Evaluate the integrand using Simpson's Rule. */
341 refp = h * ( fb + 4.0 * fo + 2.0 * fe + ff ) / 3.0;
342
343 /* Has the required precision been reached? */
344 if ( fabs ( refp - refold ) > tol ) {
345
346 /* No: prepare for next iteration. */
347 refold = refp; /* Save current value for convergence test */
348 is += is; /* Double the number of strips */
349 fe += fo; /* Sum of all = sum of evens next time */
350 fo = 0.0; /* Reset odds accumulator */
351 n = 2; /* Skip even values next time */
352
353 } else {
354
355 /* Yes: save troposphere component and terminate loop. */
356 if ( k == 1 ) reft = refp;
357 break;
358 }
359 }
360 }
361
362/* Result. */
363 *ref = reft + refp;
364 if ( zobs1 < 0.0 ) *ref = - ( *ref );
365}
366
367/*--------------------------------------------------------------------------*/
368
369static void atmt ( double robs, double tdkok, double alpha, double gamm2,
370 double delm2, double c1, double c2, double c3,
371 double c4, double c5, double c6, double r,
372 double *t, double *dn, double *rdndr )
373/*
374** - - - - -
375** a t m t
376** - - - - -
377**
378** Internal routine used by slaRefro:
379**
380** refractive index and derivative with respect to height for the
381** troposphere.
382**
383** Given:
384** robs double height of observer from centre of the Earth (metre)
385** tdkok double temperature at the observer (deg K)
386** alpha double alpha )
387** gamm2 double gamma minus 2 ) see ES
388** delm2 double delta minus 2 )
389** c1 double useful term )
390** c2 double useful term )
391** c3 double useful term ) see source of
392** c4 double useful term ) slaRefro main routine
393** c5 double useful term )
394** c6 double useful term )
395** r double current distance from the centre of the Earth (metre)
396**
397** Returned:
398** *t double temperature at r (deg K)
399** *dn double refractive index at r
400** *rdndr double r * rate the refractive index is changing at r
401**
402** This routine is derived from the ATMOSTRO routine (C.Hohenkerk,
403** HMNAO), with enhancements specified by A.T.Sinclair (RGO) to
404** handle the radio case.
405**
406** Note that in the optical case c5 and c6 are zero.
407*/
408{
409 double w, tt0, tt0gm2, tt0dm2;
410
411 w = tdkok - alpha * ( r - robs );
412 w = gmin ( w, 320.0 );
413 w = gmax ( w, 100.0 );
414 tt0 = w / tdkok;
415 tt0gm2 = pow ( tt0, gamm2 );
416 tt0dm2 = pow ( tt0, delm2 );
417 *t = w;
418 *dn = 1.0 + ( c1 * tt0gm2 - ( c2 - c5 / w ) * tt0dm2 ) * tt0;
419 *rdndr = r * ( - c3 * tt0gm2 + ( c4 - c6 / tt0 ) * tt0dm2 );
420}
421
422/*--------------------------------------------------------------------------*/
423
424static void atms ( double rt, double tt, double dnt, double gamal, double r,
425 double *dn, double *rdndr )
426/*
427** - - - - -
428** a t m s
429** - - - - -
430**
431** Internal routine used by slaRefro:
432**
433** refractive index and derivative with respect to height for the
434** stratosphere.
435**
436** Given:
437** rt double height of tropopause from centre of the Earth (metre)
438** tt double temperature at the tropopause (deg k)
439** dnt double refractive index at the tropopause
440** gamal double constant of the atmospheric model = g*md/r
441** r double current distance from the centre of the Earth (metre)
442**
443** Returned:
444** *dn double refractive index at r
445** *rdndr double r * rate the refractive index is changing at r
446**
447** This routine is derived from the ATMOSSTR routine (C.Hohenkerk, HMNAO).
448*/
449{
450 double b, w;
451
452 b = gamal / tt;
453 w = ( dnt - 1.0 ) * exp ( - b * ( r - rt ) );
454 *dn = 1.0 + w;
455 *rdndr = - r * b * w;
456}
Note: See TracBrowser for help on using the repository browser.