File: palRefro.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 (437 lines) | stat: -rw-r--r-- 14,283 bytes parent folder | download
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
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
/*
*+
*  Name:
*     palRefro

*  Purpose:
*     Atmospheric refraction for radio and optical/IR wavelengths

*  Language:
*     Starlink ANSI C

*  Type of Module:
*     Library routine

*  Invocation:
*     void palRefro( double zobs, double hm, double tdk, double pmb,
*                    double rh, double wl, double phi, double tlr,
*                    double eps, double * ref ) {

*  Arguments:
*     zobs = double (Given)
*        Observed zenith distance of the source (radian)
*     hm = double (Given)
*        Height of the observer above sea level (metre)
*     tdk = double (Given)
*        Ambient temperature at the observer (K)
*     pmb = double (Given)
*        Pressure at the observer (millibar)
*     rh = double (Given)
*        Relative humidity at the observer (range 0-1)
*     wl = double (Given)
*        Effective wavelength of the source (micrometre)
*     phi = double (Given)
*        Latitude of the observer (radian, astronomical)
*     tlr = double (Given)
*        Temperature lapse rate in the troposphere (K/metre)
*     eps = double (Given)
*        Precision required to terminate iteration (radian)
*     ref = double * (Returned)
*        Refraction: in vacuao ZD minus observed ZD (radian)

*  Description:
*     Calculates the atmospheric refraction for radio and optical/IR
*     wavelengths.

*  Authors:
*     PTW: Patrick T. Wallace
*     TIMJ: Tim Jenness (JAC, Hawaii)
*     {enter_new_authors_here}

*  Notes:
*     - A suggested value for the TLR argument is 0.0065.  The
*     refraction is significantly affected by TLR, and if studies
*     of the local atmosphere have been carried out a better TLR
*     value may be available.  The sign of the supplied TLR value
*     is ignored.
*
*     - A suggested value for the EPS argument is 1E-8.  The result is
*     usually at least two orders of magnitude more computationally
*     precise than the supplied EPS value.
*
*     - The routine computes the refraction for zenith distances up
*     to and a little beyond 90 deg using the method of Hohenkerk
*     and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted
*     in the Explanatory Supplement, 1992 edition - see section 3.281).
*
*     - The code is a development of the optical/IR refraction subroutine
*     AREF of C.Hohenkerk (HMNAO, September 1984), with extensions to
*     support the radio case.  Apart from merely cosmetic changes, the
*     following modifications to the original HMNAO optical/IR refraction
*     code have been made:
*
*     .  The angle arguments have been changed to radians.
*
*     .  Any value of ZOBS is allowed (see note 6, below).
*
*     .  Other argument values have been limited to safe values.
*
*     .  Murray's values for the gas constants have been used
*        (Vectorial Astrometry, Adam Hilger, 1983).
*
*     .  The numerical integration phase has been rearranged for
*        extra clarity.
*
*     .  A better model for Ps(T) has been adopted (taken from
*        Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982).
*
*     .  More accurate expressions for Pwo have been adopted
*        (again from Gill 1982).
*
*     .  The formula for the water vapour pressure, given the
*        saturation pressure and the relative humidity, is from
*        Crane (1976), expression 2.5.5.

*     .  Provision for radio wavelengths has been added using
*        expressions devised by A.T.Sinclair, RGO (private
*        communication 1989).  The refractivity model currently
*        used is from J.M.Rueger, "Refractive Index Formulae for
*        Electronic Distance Measurement with Radio and Millimetre
*        Waves", in Unisurv Report S-68 (2002), School of Surveying
*        and Spatial Information Systems, University of New South
*        Wales, Sydney, Australia.
*
*     .  The optical refractivity for dry air is from Resolution 3 of
*        the International Association of Geodesy adopted at the XXIIth
*        General Assembly in Birmingham, UK, 1999.
*
*     .  Various small changes have been made to gain speed.
*
*     - The radio refraction is chosen by specifying WL > 100 micrometres.
*     Because the algorithm takes no account of the ionosphere, the
*     accuracy deteriorates at low frequencies, below about 30 MHz.
*
*     - Before use, the value of ZOBS is expressed in the range +/- pi.
*     If this ranged ZOBS is -ve, the result REF is computed from its
*     absolute value before being made -ve to match.  In addition, if
*     it has an absolute value greater than 93 deg, a fixed REF value
*     equal to the result for ZOBS = 93 deg is returned, appropriately
*     signed.
*
*     - As in the original Hohenkerk and Sinclair algorithm, fixed values
*     of the water vapour polytrope exponent, the height of the
*     tropopause, and the height at which refraction is negligible are
*     used.
*
*     - The radio refraction has been tested against work done by
*     Iain Coulson, JACH, (private communication 1995) for the
*     James Clerk Maxwell Telescope, Mauna Kea.  For typical conditions,
*     agreement at the 0.1 arcsec level is achieved for moderate ZD,
*     worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg.  At hot and
*     humid sea-level sites the accuracy will not be as good.
*
*     - It should be noted that the relative humidity RH is formally
*     defined in terms of "mixing ratio" rather than pressures or
*     densities as is often stated.  It is the mass of water per unit
*     mass of dry air divided by that for saturated air at the same
*     temperature and pressure (see Gill 1982).

*     - The algorithm is designed for observers in the troposphere. The
*     supplied temperature, pressure and lapse rate are assumed to be
*     for a point in the troposphere and are used to define a model
*     atmosphere with the tropopause at 11km altitude and a constant
*     temperature above that.  However, in practice, the refraction
*     values returned for stratospheric observers, at altitudes up to
*     25km, are quite usable.

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

*  Copyright:
*     Copyright (C) 2005 Patrick T. Wallace
*     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 "pal1.h"
#include "palmac.h"

void palRefro( double zobs, double hm, double tdk, double pmb,
               double rh, double wl, double phi, double tlr,
               double eps, double * ref ) {

  /*
   *  Fixed parameters
   */

  /*  93 degrees in radians */
  const double D93 = 1.623156204;
  /*  Universal gas constant */
  const double GCR = 8314.32;
  /*  Molecular weight of dry air */
  const double DMD = 28.9644;
  /*  Molecular weight of water vapour */
  const double DMW = 18.0152;
  /*  Mean Earth radius (metre) */
  const double S = 6378120.;
  /*  Exponent of temperature dependence of water vapour pressure */
  const double DELTA = 18.36;
  /*  Height of tropopause (metre) */
  const double HT = 11000.;
  /*  Upper limit for refractive effects (metre) */
  const double HS = 80000.;
  /*  Numerical integration: maximum number of strips. */
  const int ISMAX=16384l;

  /* Local variables */
  int is, k, n, i, j;

  int optic, loop; /* booleans */

  double zobs1,zobs2,hmok,tdkok,pmbok,rhok,wlok,alpha,
    tol,wlsq,gb,a,gamal,gamma,gamm2,delm2,
    tdc,psat,pwo,w,
    c1,c2,c3,c4,c5,c6,r0,tempo,dn0,rdndr0,sk0,f0,
    rt,tt,dnt,rdndrt,sine,zt,ft,dnts,rdndrp,zts,fts,
    rs,dns,rdndrs,zs,fs,refold,z0,zrange,fb,ff,fo,fe,
    h,r,sz,rg,dr,tg,dn,rdndr,t,f,refp,reft;

  /*  The refraction integrand */
#define refi(DN,RDNDR) RDNDR/(DN+RDNDR)

  /*  Transform ZOBS into the normal range. */
  zobs1 = palDrange(zobs);
  zobs2 = DMIN(fabs(zobs1),D93);

  /*  keep other arguments within safe bounds. */
  hmok = DMIN(DMAX(hm,-1e3),HS);
  tdkok = DMIN(DMAX(tdk,100.0),500.0);
  pmbok = DMIN(DMAX(pmb,0.0),10000.0);
  rhok = DMIN(DMAX(rh,0.0),1.0);
  wlok = DMAX(wl,0.1);
  alpha = DMIN(DMAX(fabs(tlr),0.001),0.01);

  /*  tolerance for iteration. */
  tol = DMIN(DMAX(fabs(eps),1e-12),0.1)/2.0;

  /*  decide whether optical/ir or radio case - switch at 100 microns. */
  optic = wlok < 100.0;

  /*  set up model atmosphere parameters defined at the observer. */
  wlsq = wlok*wlok;
  gb = 9.784*(1.0-0.0026*cos(phi+phi)-0.00000028*hmok);
  if (optic) {
    a = (287.6155+(1.62887+0.01360/wlsq)/wlsq) * 273.15e-6/1013.25;
  } else {
    a = 77.6890e-6;
  }
  gamal = (gb*DMD)/GCR;
  gamma = gamal/alpha;
  gamm2 = gamma-2.0;
  delm2 = DELTA-2.0;
  tdc = tdkok-273.15;
  psat = pow(10.0,(0.7859+0.03477*tdc)/(1.0+0.00412*tdc)) *
    (1.0+pmbok*(4.5e-6+6.0e-10*tdc*tdc));
  if (pmbok > 0.0) {
    pwo = rhok*psat/(1.0-(1.0-rhok)*psat/pmbok);
  } else {
    pwo = 0.0;
  }
  w = pwo*(1.0-DMW/DMD)*gamma/(DELTA-gamma);
  c1 = a*(pmbok+w)/tdkok;
  if (optic) {
    c2 = (a*w+11.2684e-6*pwo)/tdkok;
  } else {
    c2 = (a*w+6.3938e-6*pwo)/tdkok;
  }
  c3 = (gamma-1.0)*alpha*c1/tdkok;
  c4 = (DELTA-1.0)*alpha*c2/tdkok;
  if (optic) {
    c5 = 0.0;
    c6 = 0.0;
  } else {
    c5 = 375463e-6*pwo/tdkok;
    c6 = c5*delm2*alpha/(tdkok*tdkok);
  }

  /*  conditions at the observer. */
  r0 = S+hmok;
  pal1Atmt(r0,tdkok,alpha,gamm2,delm2,c1,c2,c3,c4,c5,c6,
           r0,&tempo,&dn0,&rdndr0);
  sk0 = dn0*r0*sin(zobs2);
  f0 = refi(dn0,rdndr0);

  /*  conditions in the troposphere at the tropopause. */
  rt = S+DMAX(HT,hmok);
  pal1Atmt(r0,tdkok,alpha,gamm2,delm2,c1,c2,c3,c4,c5,c6,
           rt,&tt,&dnt,&rdndrt);
  sine = sk0/(rt*dnt);
  zt = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
  ft = refi(dnt,rdndrt);

  /*  conditions in the stratosphere at the tropopause. */
  pal1Atms(rt,tt,dnt,gamal,rt,&dnts,&rdndrp);
  sine = sk0/(rt*dnts);
  zts = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
  fts = refi(dnts,rdndrp);

  /*  conditions at the stratosphere limit. */
  rs = S+HS;
  pal1Atms(rt,tt,dnt,gamal,rs,&dns,&rdndrs);
  sine = sk0/(rs*dns);
  zs = atan2(sine,sqrt(DMAX(1.0-sine*sine,0.0)));
  fs = refi(dns,rdndrs);

  /*  variable initialization to avoid compiler warning. */
  reft = 0.0;

  /*  integrate the refraction integral in two parts;  first in the
   *  troposphere (k=1), then in the stratosphere (k=2). */

  for (k=1; k<=2; k++) {

    /*     initialize previous refraction to ensure at least two iterations. */
    refold = 1.0;

    /*     start off with 8 strips. */
    is = 8;

    /*     start z, z range, and start and end values. */
    if (k==1) {
      z0 = zobs2;
      zrange = zt-z0;
      fb = f0;
      ff = ft;
    } else {
      z0 = zts;
      zrange = zs-z0;
      fb = fts;
      ff = fs;
    }

    /*     sums of odd and even values. */
    fo = 0.0;
    fe = 0.0;

    /*     first time through the loop we have to do every point. */
    n = 1;

    /*     start of iteration loop (terminates at specified precision). */
    loop = 1;
    while (loop) {

      /*        strip width. */
      h = zrange/((double)is);

      /*        initialize distance from earth centre for quadrature pass. */
      if (k == 1) {
        r = r0;
      } else {
        r = rt;
      }

      /*        one pass (no need to compute evens after first time). */
      for (i=1; i<is; i+=n) {

              /*           sine of observed zenith distance. */
        sz = sin(z0+h*(double)(i));

        /*           find r (to the nearest metre, maximum four iterations). */
        if (sz > 1e-20) {
          w = sk0/sz;
          rg = r;
          dr = 1.0e6;
          j = 0;
          while ( fabs(dr) > 1.0 && j < 4 ) {
            j++;
            if (k==1) {
              pal1Atmt(r0,tdkok,alpha,gamm2,delm2,
                       c1,c2,c3,c4,c5,c6,rg,&tg,&dn,&rdndr);
            } else {
              pal1Atms(rt,tt,dnt,gamal,rg,&dn,&rdndr);
            }
            dr = (rg*dn-w)/(dn+rdndr);
            rg = rg-dr;
          }
          r = rg;
        }

        /*           find the refractive index and integrand at r. */
        if (k==1) {
          pal1Atmt(r0,tdkok,alpha,gamm2,delm2,
                   c1,c2,c3,c4,c5,c6,r,&t,&dn,&rdndr);
        } else {
          pal1Atms(rt,tt,dnt,gamal,r,&dn,&rdndr);
        }
        f = refi(dn,rdndr);

        /*           accumulate odd and (first time only) even values. */
        if (n==1 && i%2 == 0) {
          fe += f;
        } else {
          fo += f;
        }
      }

      /*        evaluate the integrand using simpson's rule. */
      refp = h*(fb+4.0*fo+2.0*fe+ff)/3.0;

      /*        has the required precision been achieved (or can't be)? */
      if (fabs(refp-refold) > tol && is < ISMAX) {

        /*           no: prepare for next iteration.*/

        /*           save current value for convergence test. */
        refold = refp;

        /*           double the number of strips. */
        is += is;

        /*           sum of all current values = sum of next pass's even values. */
        fe += fo;

        /*           prepare for new odd values. */
        fo = 0.0;

        /*           skip even values next time. */
        n = 2;
      } else {

        /*           yes: save troposphere component and terminate the loop. */
        if (k==1) reft = refp;
        loop = 0;
      }
    }
  }

  /*  result. */
  *ref = reft+refp;
  if (zobs1 < 0.0) *ref = -(*ref);

}