| 1 | #include "slalib.h" | 
|---|
| 2 | #include "slamac.h" | 
|---|
| 3 | void slaRefcoq ( double tdk, double pmb, double rh, double wl, | 
|---|
| 4 | double *refa, double *refb ) | 
|---|
| 5 | /* | 
|---|
| 6 | **  - - - - - - - - - - | 
|---|
| 7 | **   s l a R e f c o q | 
|---|
| 8 | **  - - - - - - - - - - | 
|---|
| 9 | ** | 
|---|
| 10 | **  Determine the constants A and B in the atmospheric refraction | 
|---|
| 11 | **  model dZ = A tan Z + B tan^3 Z.  This is a fast alternative | 
|---|
| 12 | **  to the slaRefco routine - see notes. | 
|---|
| 13 | ** | 
|---|
| 14 | **  Z is the "observed" zenith distance (i.e. affected by refraction) | 
|---|
| 15 | **  and dZ is what to add to Z to give the "topocentric" (i.e. in vacuo) | 
|---|
| 16 | **  zenith distance. | 
|---|
| 17 | ** | 
|---|
| 18 | **  Given: | 
|---|
| 19 | **    tdk    double    ambient temperature at the observer (deg K) | 
|---|
| 20 | **    pmb    double    pressure at the observer (millibar) | 
|---|
| 21 | **    rh     double    relative humidity at the observer (range 0-1) | 
|---|
| 22 | **    wl     double    effective wavelength of the source (micrometre) | 
|---|
| 23 | ** | 
|---|
| 24 | **  Returned: | 
|---|
| 25 | **    refa   double*   tan Z coefficient (radian) | 
|---|
| 26 | **    refb   double*   tan^3 Z coefficient (radian) | 
|---|
| 27 | ** | 
|---|
| 28 | **  The radio refraction is chosen by specifying WL > 100 micrometres. | 
|---|
| 29 | ** | 
|---|
| 30 | **  Notes: | 
|---|
| 31 | ** | 
|---|
| 32 | **  1  The model is an approximation, for moderate zenith distances, | 
|---|
| 33 | **     to the predictions of the slaRefro routine.  The approximation | 
|---|
| 34 | **     is maintained across a range of conditions, and applies to | 
|---|
| 35 | **     both optical/IR and radio. | 
|---|
| 36 | ** | 
|---|
| 37 | **  2  The algorithm is a fast alternative to the slaRefco routine. | 
|---|
| 38 | **     The latter calls the slaRefro routine itself:  this involves | 
|---|
| 39 | **     integrations through a model atmosphere, and is costly in | 
|---|
| 40 | **     processor time.  However, the model which is produced is precisely | 
|---|
| 41 | **     correct for two zenith distance (45 degrees and about 76 degrees) | 
|---|
| 42 | **     and at other zenith distances is limited in accuracy only by the | 
|---|
| 43 | **     A tan Z + B tan^3 Z formulation itself.  The present routine | 
|---|
| 44 | **     is not as accurate, though it satisfies most practical | 
|---|
| 45 | **     requirements. | 
|---|
| 46 | ** | 
|---|
| 47 | **  3  The model omits the effects of (i) height above sea level (apart | 
|---|
| 48 | **     from the reduced pressure itself), (ii) latitude (i.e. the | 
|---|
| 49 | **     flattening of the Earth) and (iii) variations in tropospheric | 
|---|
| 50 | **     lapse rate. | 
|---|
| 51 | ** | 
|---|
| 52 | **     The model was tested using the following range of conditions: | 
|---|
| 53 | ** | 
|---|
| 54 | **       lapse rates 0.0055, 0.0065, 0.0075 deg/metre | 
|---|
| 55 | **       latitudes 0, 25, 50, 75 degrees | 
|---|
| 56 | **       heights 0, 2500, 5000 metres ASL | 
|---|
| 57 | **       pressures mean for height -10% to +5% in steps of 5% | 
|---|
| 58 | **       temperatures -10 deg to +20 deg with respect to 280 deg at SL | 
|---|
| 59 | **       relative humidity 0, 0.5, 1 | 
|---|
| 60 | **       wavelengths 0.4, 0.6, ... 2 micron, + radio | 
|---|
| 61 | **       zenith distances 15, 45, 75 degrees | 
|---|
| 62 | ** | 
|---|
| 63 | **     The accuracy with respect to direct use of the slaRefro routine | 
|---|
| 64 | **     was as follows: | 
|---|
| 65 | ** | 
|---|
| 66 | **                            worst         RMS | 
|---|
| 67 | ** | 
|---|
| 68 | **       optical/IR           62 mas       8 mas | 
|---|
| 69 | **       radio               319 mas      49 mas | 
|---|
| 70 | ** | 
|---|
| 71 | **     For this particular set of conditions: | 
|---|
| 72 | ** | 
|---|
| 73 | **       lapse rate 0.0065 degK/metre | 
|---|
| 74 | **       latitude 50 degrees | 
|---|
| 75 | **       sea level | 
|---|
| 76 | **       pressure 1005 mB | 
|---|
| 77 | **       temperature 280.15 degK | 
|---|
| 78 | **       humidity 80% | 
|---|
| 79 | **       wavelength 5740 Angstroms | 
|---|
| 80 | ** | 
|---|
| 81 | **     the results were as follows: | 
|---|
| 82 | ** | 
|---|
| 83 | **       ZD        slaRefro    slaRefcoq   Saastamoinen | 
|---|
| 84 | ** | 
|---|
| 85 | **       10         10.27        10.27        10.27 | 
|---|
| 86 | **       20         21.19        21.20        21.19 | 
|---|
| 87 | **       30         33.61        33.61        33.60 | 
|---|
| 88 | **       40         48.82        48.83        48.81 | 
|---|
| 89 | **       45         58.16        58.18        58.16 | 
|---|
| 90 | **       50         69.28        69.30        69.27 | 
|---|
| 91 | **       55         82.97        82.99        82.95 | 
|---|
| 92 | **       60        100.51       100.54       100.50 | 
|---|
| 93 | **       65        124.23       124.26       124.20 | 
|---|
| 94 | **       70        158.63       158.68       158.61 | 
|---|
| 95 | **       72        177.32       177.37       177.31 | 
|---|
| 96 | **       74        200.35       200.38       200.32 | 
|---|
| 97 | **       76        229.45       229.43       229.42 | 
|---|
| 98 | **       78        267.44       267.29       267.41 | 
|---|
| 99 | **       80        319.13       318.55       319.10 | 
|---|
| 100 | ** | 
|---|
| 101 | **      deg        arcsec       arcsec       arcsec | 
|---|
| 102 | ** | 
|---|
| 103 | **     The values for Saastamoinen's formula (which includes terms | 
|---|
| 104 | **     up to tan^5) are taken from Hohenkerk and Sinclair (1985). | 
|---|
| 105 | ** | 
|---|
| 106 | **     The results from the much slower but more accurate slaRefco | 
|---|
| 107 | **     routine have not been included in the tabulation as they are | 
|---|
| 108 | **     identical to those in the slaRefro column to the 0.01 arcsec | 
|---|
| 109 | **     resolution used. | 
|---|
| 110 | ** | 
|---|
| 111 | **  4  Outlandish input parameters are silently limited to mathematically | 
|---|
| 112 | **     safe values.  Zero pressure is permissible, and causes zeroes to | 
|---|
| 113 | **     be returned. | 
|---|
| 114 | ** | 
|---|
| 115 | **  5  The algorithm draws on several sources, as follows: | 
|---|
| 116 | ** | 
|---|
| 117 | **     a) The formula for the saturation vapour pressure of water as | 
|---|
| 118 | **        a function of temperature and temperature is taken from | 
|---|
| 119 | **        expressions A4.5-A4.7 of Gill (1982). | 
|---|
| 120 | ** | 
|---|
| 121 | **     b) The formula for the water vapour pressure, Given the | 
|---|
| 122 | **        saturation pressure and the relative humidity, is from | 
|---|
| 123 | **        Crane (1976), expression 2.5.5. | 
|---|
| 124 | ** | 
|---|
| 125 | **     c) The refractivity of air is a function of temperature, | 
|---|
| 126 | **        total pressure, water-vapour pressure and, in the case | 
|---|
| 127 | **        of optical/IR but not radio, wavelength.  The formulae | 
|---|
| 128 | **        for the two cases are developed from the Essen and Froome | 
|---|
| 129 | **        expressions adopted in Resolution 1 of the 12th International | 
|---|
| 130 | **        Geodesy Association General Assembly (1963). | 
|---|
| 131 | ** | 
|---|
| 132 | **     The above three items are as used in the slaRefro routine. | 
|---|
| 133 | ** | 
|---|
| 134 | **     d) The formula for beta, the ratio of the scale height of the | 
|---|
| 135 | **        atmosphere to the geocentric distance of the observer, is | 
|---|
| 136 | **        an adaption of expression 9 from Stone (1996).  The | 
|---|
| 137 | **        adaptations, arrived at empirically, consist of (i) a | 
|---|
| 138 | **        small adjustment to the coefficient and (ii) a humidity | 
|---|
| 139 | **        term for the radio case only. | 
|---|
| 140 | ** | 
|---|
| 141 | **     e) The formulae for the refraction constants as a function of | 
|---|
| 142 | **        n-1 and beta are from Green (1987), expression 4.31. | 
|---|
| 143 | ** | 
|---|
| 144 | **  References: | 
|---|
| 145 | ** | 
|---|
| 146 | **     Crane, R.K., Meeks, M.L. (ed), "Refraction Effects in the Neutral | 
|---|
| 147 | **     Atmosphere", Methods of Experimental Physics: Astrophysics 12B, | 
|---|
| 148 | **     Academic Press, 1976. | 
|---|
| 149 | ** | 
|---|
| 150 | **     Gill, Adrian E., "Atmosphere-Ocean Dynamics", Academic Press, 1982. | 
|---|
| 151 | ** | 
|---|
| 152 | **     Hohenkerk, C.Y., & Sinclair, A.T., NAO Technical Note No. 63, 1985. | 
|---|
| 153 | ** | 
|---|
| 154 | **     International Geodesy Association General Assembly, Bulletin | 
|---|
| 155 | **     Geodesique 70 p390, 1963. | 
|---|
| 156 | ** | 
|---|
| 157 | **     Stone, Ronald C., P.A.S.P. 108 1051-1058, 1996. | 
|---|
| 158 | ** | 
|---|
| 159 | **     Green, R.M., "Spherical Astronomy", Cambridge University Press, 1987. | 
|---|
| 160 | ** | 
|---|
| 161 | **  Last revision:   17 March 1999 | 
|---|
| 162 | ** | 
|---|
| 163 | **  Copyright P.T.Wallace.  All rights reserved. | 
|---|
| 164 | */ | 
|---|
| 165 | { | 
|---|
| 166 | int optic; | 
|---|
| 167 | double t, p, r,w, tdc, ps, pw, wlsq, gamma, beta; | 
|---|
| 168 |  | 
|---|
| 169 |  | 
|---|
| 170 | /* Decide whether optical/IR or radio case:  switch at 100 microns. */ | 
|---|
| 171 | optic = ( wl <= 100.0 ); | 
|---|
| 172 |  | 
|---|
| 173 | /* Restrict parameters to safe values. */ | 
|---|
| 174 | t = gmax ( tdk, 100.0 ); | 
|---|
| 175 | t = gmin ( t, 500.0 ); | 
|---|
| 176 | p = gmax ( pmb, 0.0 ); | 
|---|
| 177 | p = gmin ( p, 10000.0 ); | 
|---|
| 178 | r = gmax ( rh, 0.0 ); | 
|---|
| 179 | r = gmin ( r, 1.0 ); | 
|---|
| 180 | w = gmax ( wl, 0.1 ); | 
|---|
| 181 | w = gmin ( w, 1e6 ); | 
|---|
| 182 |  | 
|---|
| 183 | /* Water vapour pressure at the observer. */ | 
|---|
| 184 | if ( p > 0.0 ) { | 
|---|
| 185 | tdc = t - 273.15; | 
|---|
| 186 | ps = pow ( 10.0, ( 0.7859 + 0.03477 * tdc ) / | 
|---|
| 187 | ( 1.0 + 0.00412 * tdc ) ) * | 
|---|
| 188 | ( 1.0 + p * ( 4.5e-6 + 6e-10 * tdc * tdc )  ); | 
|---|
| 189 | pw = r * ps / ( 1.0 - ( 1.0 - r ) * ps / p ); | 
|---|
| 190 | } else { | 
|---|
| 191 | pw = 0.0; | 
|---|
| 192 | } | 
|---|
| 193 |  | 
|---|
| 194 | /* Refractive index minus 1 at the observer. */ | 
|---|
| 195 | if ( optic ) { | 
|---|
| 196 | wlsq = wl * wl; | 
|---|
| 197 | gamma = ( ( 77.532e-6 + ( 4.391e-7 + 3.57e-9 / wlsq ) / wlsq ) * p | 
|---|
| 198 | - 11.2684e-6 * pw ) / t; | 
|---|
| 199 | } else { | 
|---|
| 200 | gamma = ( 77.624e-6 * p - ( 12.92e-6 - 0.371897 / t ) * pw ) / t; | 
|---|
| 201 | } | 
|---|
| 202 |  | 
|---|
| 203 | /* Formula for beta from Stone, with empirical adjustments. */ | 
|---|
| 204 | beta = 4.4474e-6 * t; | 
|---|
| 205 | if ( !optic ) beta -= 0.0074 * pw * beta; | 
|---|
| 206 |  | 
|---|
| 207 | /* Refraction constants from Green. */ | 
|---|
| 208 | *refa = gamma * ( 1.0 - beta ); | 
|---|
| 209 | *refb = - gamma * ( beta - gamma / 2.0 ); | 
|---|
| 210 | } | 
|---|