| 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 |
|
|---|
| 185 | void 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 | }
|
|---|
