| 1 | #include "slalib.h"
|
|---|
| 2 | #include "slamac.h"
|
|---|
| 3 |
|
|---|
| 4 | static void atmt ( double, double, double, double, double, double,
|
|---|
| 5 | double, double, double, double, double, double,
|
|---|
| 6 | double*, double*, double* );
|
|---|
| 7 | static void atms ( double, double, double, double, double,
|
|---|
| 8 | double*, double* );
|
|---|
| 9 |
|
|---|
| 10 | void 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 |
|
|---|
| 369 | static 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 |
|
|---|
| 424 | static 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 | }
|
|---|
