source: trunk/FACT++/pal/palRefro.c

Last change on this file was 18347, checked in by tbretz, 9 years ago
File size: 13.9 KB
Line 
1/*
2*+
3* Name:
4* palRefro
5
6* Purpose:
7* Atmospheric refraction for radio and optical/IR wavelengths
8
9* Language:
10* Starlink ANSI C
11
12* Type of Module:
13* Library routine
14
15* Invocation:
16* void palRefro( double zobs, double hm, double tdk, double pmb,
17* double rh, double wl, double phi, double tlr,
18* double eps, double * ref ) {
19
20* Arguments:
21* zobs = double (Given)
22* Observed zenith distance of the source (radian)
23* hm = double (Given)
24* Height of the observer above sea level (metre)
25* tdk = double (Given)
26* Ambient temperature at the observer (K)
27* pmb = double (Given)
28* Pressure at the observer (millibar)
29* rh = double (Given)
30* Relative humidity at the observer (range 0-1)
31* wl = double (Given)
32* Effective wavelength of the source (micrometre)
33* phi = double (Given)
34* Latitude of the observer (radian, astronomical)
35* tlr = double (Given)
36* Temperature lapse rate in the troposphere (K/metre)
37* eps = double (Given)
38* Precision required to terminate iteration (radian)
39* ref = double * (Returned)
40* Refraction: in vacuao ZD minus observed ZD (radian)
41
42* Description:
43* Calculates the atmospheric refraction for radio and optical/IR
44* wavelengths.
45
46* Authors:
47* PTW: Patrick T. Wallace
48* TIMJ: Tim Jenness (JAC, Hawaii)
49* {enter_new_authors_here}
50
51* Notes:
52* - A suggested value for the TLR argument is 0.0065. The
53* refraction is significantly affected by TLR, and if studies
54* of the local atmosphere have been carried out a better TLR
55* value may be available. The sign of the supplied TLR value
56* is ignored.
57*
58* - A suggested value for the EPS argument is 1E-8. The result is
59* usually at least two orders of magnitude more computationally
60* precise than the supplied EPS value.
61*
62* - The routine computes the refraction for zenith distances up
63* to and a little beyond 90 deg using the method of Hohenkerk
64* and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted
65* in the Explanatory Supplement, 1992 edition - see section 3.281).
66*
67* - The code is a development of the optical/IR refraction subroutine
68* AREF of C.Hohenkerk (HMNAO, September 1984), with extensions to
69* support the radio case. Apart from merely cosmetic changes, the
70* following modifications to the original HMNAO optical/IR refraction
71* code have been made:
72*
73* . The angle arguments have been changed to radians.
74*
75* . Any value of ZOBS is allowed (see note 6, below).
76*
77* . Other argument values have been limited to safe values.
78*
79* . Murray's values for the gas constants have been used
80* (Vectorial Astrometry, Adam Hilger, 1983).
81*
82* . The numerical integration phase has been rearranged for
83* extra clarity.
84*
85* . A better model for Ps(T) has been adopted (taken from
86* Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982).
87*
88* . More accurate expressions for Pwo have been adopted
89* (again from Gill 1982).
90*
91* . The formula for the water vapour pressure, given the
92* saturation pressure and the relative humidity, is from
93* Crane (1976), expression 2.5.5.
94
95* . Provision for radio wavelengths has been added using
96* expressions devised by A.T.Sinclair, RGO (private
97* communication 1989). The refractivity model currently
98* used is from J.M.Rueger, "Refractive Index Formulae for
99* Electronic Distance Measurement with Radio and Millimetre
100* Waves", in Unisurv Report S-68 (2002), School of Surveying
101* and Spatial Information Systems, University of New South
102* Wales, Sydney, Australia.
103*
104* . The optical refractivity for dry air is from Resolution 3 of
105* the International Association of Geodesy adopted at the XXIIth
106* General Assembly in Birmingham, UK, 1999.
107*
108* . Various small changes have been made to gain speed.
109*
110* - The radio refraction is chosen by specifying WL > 100 micrometres.
111* Because the algorithm takes no account of the ionosphere, the
112* accuracy deteriorates at low frequencies, below about 30 MHz.
113*
114* - Before use, the value of ZOBS is expressed in the range +/- pi.
115* If this ranged ZOBS is -ve, the result REF is computed from its
116* absolute value before being made -ve to match. In addition, if
117* it has an absolute value greater than 93 deg, a fixed REF value
118* equal to the result for ZOBS = 93 deg is returned, appropriately
119* signed.
120*
121* - As in the original Hohenkerk and Sinclair algorithm, fixed values
122* of the water vapour polytrope exponent, the height of the
123* tropopause, and the height at which refraction is negligible are
124* used.
125*
126* - The radio refraction has been tested against work done by
127* Iain Coulson, JACH, (private communication 1995) for the
128* James Clerk Maxwell Telescope, Mauna Kea. For typical conditions,
129* agreement at the 0.1 arcsec level is achieved for moderate ZD,
130* worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg. At hot and
131* humid sea-level sites the accuracy will not be as good.
132*
133* - It should be noted that the relative humidity RH is formally
134* defined in terms of "mixing ratio" rather than pressures or
135* densities as is often stated. It is the mass of water per unit
136* mass of dry air divided by that for saturated air at the same
137* temperature and pressure (see Gill 1982).
138
139* - The algorithm is designed for observers in the troposphere. The
140* supplied temperature, pressure and lapse rate are assumed to be
141* for a point in the troposphere and are used to define a model
142* atmosphere with the tropopause at 11km altitude and a constant
143* temperature above that. However, in practice, the refraction
144* values returned for stratospheric observers, at altitudes up to
145* 25km, are quite usable.
146
147* History:
148* 2012-08-24 (TIMJ):
149* Initial version, direct port of SLA Fortran source.
150* Adapted with permission from the Fortran SLALIB library.
151* {enter_further_changes_here}
152
153* Copyright:
154* Copyright (C) 2005 Patrick T. Wallace
155* Copyright (C) 2012 Science and Technology Facilities Council.
156* All Rights Reserved.
157
158* Licence:
159* This program is free software; you can redistribute it and/or
160* modify it under the terms of the GNU General Public License as
161* published by the Free Software Foundation; either version 3 of
162* the License, or (at your option) any later version.
163*
164* This program is distributed in the hope that it will be
165* useful, but WITHOUT ANY WARRANTY; without even the implied
166* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
167* PURPOSE. See the GNU General Public License for more details.
168*
169* You should have received a copy of the GNU General Public License
170* along with this program; if not, write to the Free Software
171* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
172* MA 02110-1301, USA.
173
174* Bugs:
175* {note_any_bugs_here}
176*-
177*/
178
179#include <math.h>
180
181#include "pal.h"
182#include "pal1.h"
183#include "palmac.h"
184
185void palRefro( double zobs, double hm, double tdk, double pmb,
186 double rh, double wl, double phi, double tlr,
187 double eps, double * ref ) {
188
189 /*
190 * Fixed parameters
191 */
192
193 /* 93 degrees in radians */
194 const double D93 = 1.623156204;
195 /* Universal gas constant */
196 const double GCR = 8314.32;
197 /* Molecular weight of dry air */
198 const double DMD = 28.9644;
199 /* Molecular weight of water vapour */
200 const double DMW = 18.0152;
201 /* Mean Earth radius (metre) */
202 const double S = 6378120.;
203 /* Exponent of temperature dependence of water vapour pressure */
204 const double DELTA = 18.36;
205 /* Height of tropopause (metre) */
206 const double HT = 11000.;
207 /* Upper limit for refractive effects (metre) */
208 const double HS = 80000.;
209 /* Numerical integration: maximum number of strips. */
210 const int ISMAX=16384l;
211
212 /* Local variables */
213 int is, k, n, i, j;
214
215 int optic, loop; /* booleans */
216
217 double zobs1,zobs2,hmok,tdkok,pmbok,rhok,wlok,alpha,
218 tol,wlsq,gb,a,gamal,gamma,gamm2,delm2,
219 tdc,psat,pwo,w,
220 c1,c2,c3,c4,c5,c6,r0,tempo,dn0,rdndr0,sk0,f0,
221 rt,tt,dnt,rdndrt,sine,zt,ft,dnts,rdndrp,zts,fts,
222 rs,dns,rdndrs,zs,fs,refold,z0,zrange,fb,ff,fo,fe,
223 h,r,sz,rg,dr,tg,dn,rdndr,t,f,refp,reft;
224
225 /* The refraction integrand */
226#define refi(DN,RDNDR) RDNDR/(DN+RDNDR)
227
228 /* Transform ZOBS into the normal range. */
229 zobs1 = palDrange(zobs);
230 zobs2 = DMIN(fabs(zobs1),D93);
231
232 /* keep other arguments within safe bounds. */
233 hmok = DMIN(DMAX(hm,-1e3),HS);
234 tdkok = DMIN(DMAX(tdk,100.0),500.0);
235 pmbok = DMIN(DMAX(pmb,0.0),10000.0);
236 rhok = DMIN(DMAX(rh,0.0),1.0);
237 wlok = DMAX(wl,0.1);
238 alpha = DMIN(DMAX(fabs(tlr),0.001),0.01);
239
240 /* tolerance for iteration. */
241 tol = DMIN(DMAX(fabs(eps),1e-12),0.1)/2.0;
242
243 /* decide whether optical/ir or radio case - switch at 100 microns. */
244 optic = wlok < 100.0;
245
246 /* set up model atmosphere parameters defined at the observer. */
247 wlsq = wlok*wlok;
248 gb = 9.784*(1.0-0.0026*cos(phi+phi)-0.00000028*hmok);
249 if (optic) {
250 a = (287.6155+(1.62887+0.01360/wlsq)/wlsq) * 273.15e-6/1013.25;
251 } else {
252 a = 77.6890e-6;
253 }
254 gamal = (gb*DMD)/GCR;
255 gamma = gamal/alpha;
256 gamm2 = gamma-2.0;
257 delm2 = DELTA-2.0;
258 tdc = tdkok-273.15;
259 psat = pow(10.0,(0.7859+0.03477*tdc)/(1.0+0.00412*tdc)) *
260 (1.0+pmbok*(4.5e-6+6.0e-10*tdc*tdc));
261 if (pmbok > 0.0) {
262 pwo = rhok*psat/(1.0-(1.0-rhok)*psat/pmbok);
263 } else {
264 pwo = 0.0;
265 }
266 w = pwo*(1.0-DMW/DMD)*gamma/(DELTA-gamma);
267 c1 = a*(pmbok+w)/tdkok;
268 if (optic) {
269 c2 = (a*w+11.2684e-6*pwo)/tdkok;
270 } else {
271 c2 = (a*w+6.3938e-6*pwo)/tdkok;
272 }
273 c3 = (gamma-1.0)*alpha*c1/tdkok;
274 c4 = (DELTA-1.0)*alpha*c2/tdkok;
275 if (optic) {
276 c5 = 0.0;
277 c6 = 0.0;
278 } else {
279 c5 = 375463e-6*pwo/tdkok;
280 c6 = c5*delm2*alpha/(tdkok*tdkok);
281 }
282
283 /* conditions at the observer. */
284 r0 = S+hmok;
285 pal1Atmt(r0,tdkok,alpha,gamm2,delm2,c1,c2,c3,c4,c5,c6,
286 r0,&tempo,&dn0,&rdndr0);
287 sk0 = dn0*r0*sin(zobs2);
288 f0 = refi(dn0,rdndr0);
289
290 /* conditions in the troposphere at the tropopause. */
291 rt = S+DMAX(HT,hmok);
292 pal1Atmt(r0,tdkok,alpha,gamm2,delm2,c1,c2,c3,c4,c5,c6,
293 rt,&tt,&dnt,&rdndrt);
294 sine = sk0/(rt*dnt);
295 zt = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
296 ft = refi(dnt,rdndrt);
297
298 /* conditions in the stratosphere at the tropopause. */
299 pal1Atms(rt,tt,dnt,gamal,rt,&dnts,&rdndrp);
300 sine = sk0/(rt*dnts);
301 zts = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
302 fts = refi(dnts,rdndrp);
303
304 /* conditions at the stratosphere limit. */
305 rs = S+HS;
306 pal1Atms(rt,tt,dnt,gamal,rs,&dns,&rdndrs);
307 sine = sk0/(rs*dns);
308 zs = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
309 fs = refi(dns,rdndrs);
310
311 /* variable initialization to avoid compiler warning. */
312 reft = 0.0;
313
314 /* integrate the refraction integral in two parts; first in the
315 * troposphere (k=1), then in the stratosphere (k=2). */
316
317 for (k=1; k<=2; k++) {
318
319 /* initialize previous refraction to ensure at least two iterations. */
320 refold = 1.0;
321
322 /* start off with 8 strips. */
323 is = 8;
324
325 /* start z, z range, and start and end values. */
326 if (k==1) {
327 z0 = zobs2;
328 zrange = zt-z0;
329 fb = f0;
330 ff = ft;
331 } else {
332 z0 = zts;
333 zrange = zs-z0;
334 fb = fts;
335 ff = fs;
336 }
337
338 /* sums of odd and even values. */
339 fo = 0.0;
340 fe = 0.0;
341
342 /* first time through the loop we have to do every point. */
343 n = 1;
344
345 /* start of iteration loop (terminates at specified precision). */
346 loop = 1;
347 while (loop) {
348
349 /* strip width. */
350 h = zrange/((double)is);
351
352 /* initialize distance from earth centre for quadrature pass. */
353 if (k == 1) {
354 r = r0;
355 } else {
356 r = rt;
357 }
358
359 /* one pass (no need to compute evens after first time). */
360 for (i=1; i<is; i+=n) {
361
362 /* sine of observed zenith distance. */
363 sz = sin(z0+h*(double)(i));
364
365 /* find r (to the nearest metre, maximum four iterations). */
366 if (sz > 1e-20) {
367 w = sk0/sz;
368 rg = r;
369 dr = 1.0e6;
370 j = 0;
371 while ( fabs(dr) > 1.0 && j < 4 ) {
372 j++;
373 if (k==1) {
374 pal1Atmt(r0,tdkok,alpha,gamm2,delm2,
375 c1,c2,c3,c4,c5,c6,rg,&tg,&dn,&rdndr);
376 } else {
377 pal1Atms(rt,tt,dnt,gamal,rg,&dn,&rdndr);
378 }
379 dr = (rg*dn-w)/(dn+rdndr);
380 rg = rg-dr;
381 }
382 r = rg;
383 }
384
385 /* find the refractive index and integrand at r. */
386 if (k==1) {
387 pal1Atmt(r0,tdkok,alpha,gamm2,delm2,
388 c1,c2,c3,c4,c5,c6,r,&t,&dn,&rdndr);
389 } else {
390 pal1Atms(rt,tt,dnt,gamal,r,&dn,&rdndr);
391 }
392 f = refi(dn,rdndr);
393
394 /* accumulate odd and (first time only) even values. */
395 if (n==1 && i%2 == 0) {
396 fe += f;
397 } else {
398 fo += f;
399 }
400 }
401
402 /* evaluate the integrand using simpson's rule. */
403 refp = h*(fb+4.0*fo+2.0*fe+ff)/3.0;
404
405 /* has the required precision been achieved (or can't be)? */
406 if (fabs(refp-refold) > tol && is < ISMAX) {
407
408 /* no: prepare for next iteration.*/
409
410 /* save current value for convergence test. */
411 refold = refp;
412
413 /* double the number of strips. */
414 is += is;
415
416 /* sum of all current values = sum of next pass's even values. */
417 fe += fo;
418
419 /* prepare for new odd values. */
420 fo = 0.0;
421
422 /* skip even values next time. */
423 n = 2;
424 } else {
425
426 /* yes: save troposphere component and terminate the loop. */
427 if (k==1) reft = refp;
428 loop = 0;
429 }
430 }
431 }
432
433 /* result. */
434 *ref = reft+refp;
435 if (zobs1 < 0.0) *ref = -(*ref);
436
437}
Note: See TracBrowser for help on using the repository browser.