From d54fe7c1f704a63824c5bfa0ece65245572e9b27 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 4 Mar 2015 21:21:30 -0500 Subject: Initial commit --- src/slalib/refcoq.f | 208 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 208 insertions(+) create mode 100644 src/slalib/refcoq.f (limited to 'src/slalib/refcoq.f') diff --git a/src/slalib/refcoq.f b/src/slalib/refcoq.f new file mode 100644 index 0000000..15e3af0 --- /dev/null +++ b/src/slalib/refcoq.f @@ -0,0 +1,208 @@ + SUBROUTINE sla_REFCOQ (TDK, PMB, RH, WL, REFA, REFB) +*+ +* - - - - - - - +* R E F C O Q +* - - - - - - - +* +* Determine the constants A and B in the atmospheric refraction +* model dZ = A tan Z + B tan**3 Z. This is a fast alternative +* to the sla_REFCO routine - see notes. +* +* Z is the "observed" zenith distance (i.e. affected by refraction) +* and dZ is what to add to Z to give the "topocentric" (i.e. in vacuo) +* zenith distance. +* +* Given: +* TDK d ambient temperature at the observer (deg K) +* PMB d pressure at the observer (millibar) +* RH d relative humidity at the observer (range 0-1) +* WL d effective wavelength of the source (micrometre) +* +* Returned: +* REFA d tan Z coefficient (radian) +* REFB d tan**3 Z coefficient (radian) +* +* The radio refraction is chosen by specifying WL > 100 micrometres. +* +* Notes: +* +* 1 The model is an approximation, for moderate zenith distances, +* to the predictions of the sla_REFRO routine. The approximation +* is maintained across a range of conditions, and applies to +* both optical/IR and radio. +* +* 2 The algorithm is a fast alternative to the sla_REFCO routine. +* The latter calls the sla_REFRO routine itself: this involves +* integrations through a model atmosphere, and is costly in +* processor time. However, the model which is produced is precisely +* correct for two zenith distance (45 degrees and about 76 degrees) +* and at other zenith distances is limited in accuracy only by the +* A tan Z + B tan**3 Z formulation itself. The present routine +* is not as accurate, though it satisfies most practical +* requirements. +* +* 3 The model omits the effects of (i) height above sea level (apart +* from the reduced pressure itself), (ii) latitude (i.e. the +* flattening of the Earth) and (iii) variations in tropospheric +* lapse rate. +* +* The model was tested using the following range of conditions: +* +* lapse rates 0.0055, 0.0065, 0.0075 deg/metre +* latitudes 0, 25, 50, 75 degrees +* heights 0, 2500, 5000 metres ASL +* pressures mean for height -10% to +5% in steps of 5% +* temperatures -10 deg to +20 deg with respect to 280 deg at SL +* relative humidity 0, 0.5, 1 +* wavelengths 0.4, 0.6, ... 2 micron, + radio +* zenith distances 15, 45, 75 degrees +* +* The accuracy with respect to direct use of the sla_REFRO routine +* was as follows: +* +* worst RMS +* +* optical/IR 62 mas 8 mas +* radio 319 mas 49 mas +* +* For this particular set of conditions: +* +* lapse rate 0.0065 degK/metre +* latitude 50 degrees +* sea level +* pressure 1005 mB +* temperature 280.15 degK +* humidity 80% +* wavelength 5740 Angstroms +* +* the results were as follows: +* +* ZD sla_REFRO sla_REFCOQ Saastamoinen +* +* 10 10.27 10.27 10.27 +* 20 21.19 21.20 21.19 +* 30 33.61 33.61 33.60 +* 40 48.82 48.83 48.81 +* 45 58.16 58.18 58.16 +* 50 69.28 69.30 69.27 +* 55 82.97 82.99 82.95 +* 60 100.51 100.54 100.50 +* 65 124.23 124.26 124.20 +* 70 158.63 158.68 158.61 +* 72 177.32 177.37 177.31 +* 74 200.35 200.38 200.32 +* 76 229.45 229.43 229.42 +* 78 267.44 267.29 267.41 +* 80 319.13 318.55 319.10 +* +* deg arcsec arcsec arcsec +* +* The values for Saastamoinen's formula (which includes terms +* up to tan^5) are taken from Hohenkerk and Sinclair (1985). +* +* The results from the much slower but more accurate sla_REFCO +* routine have not been included in the tabulation as they are +* identical to those in the sla_REFRO column to the 0.01 arcsec +* resolution used. +* +* 4 Outlandish input parameters are silently limited to mathematically +* safe values. Zero pressure is permissible, and causes zeroes to +* be returned. +* +* 5 The algorithm draws on several sources, as follows: +* +* a) The formula for the saturation vapour pressure of water as +* a function of temperature and temperature is taken from +* expressions A4.5-A4.7 of Gill (1982). +* +* b) The formula for the water vapour pressure, given the +* saturation pressure and the relative humidity, is from +* Crane (1976), expression 2.5.5. +* +* c) The refractivity of air is a function of temperature, +* total pressure, water-vapour pressure and, in the case +* of optical/IR but not radio, wavelength. The formulae +* for the two cases are developed from the Essen and Froome +* expressions adopted in Resolution 1 of the 12th International +* Geodesy Association General Assembly (1963). +* +* The above three items are as used in the sla_REFRO routine. +* +* d) The formula for beta, the ratio of the scale height of the +* atmosphere to the geocentric distance of the observer, is +* an adaption of expression 9 from Stone (1996). The +* adaptations, arrived at empirically, consist of (i) a +* small adjustment to the coefficient and (ii) a humidity +* term for the radio case only. +* +* e) The formulae for the refraction constants as a function of +* n-1 and beta are from Green (1987), expression 4.31. +* +* References: +* +* Crane, R.K., Meeks, M.L. (ed), "Refraction Effects in the Neutral +* Atmosphere", Methods of Experimental Physics: Astrophysics 12B, +* Academic Press, 1976. +* +* Gill, Adrian E., "Atmosphere-Ocean Dynamics", Academic Press, 1982. +* +* Hohenkerk, C.Y., & Sinclair, A.T., NAO Technical Note No. 63, 1985. +* +* International Geodesy Association General Assembly, Bulletin +* Geodesique 70 p390, 1963. +* +* Stone, Ronald C., P.A.S.P. 108 1051-1058, 1996. +* +* Green, R.M., "Spherical Astronomy", Cambridge University Press, 1987. +* +* P.T.Wallace Starlink 4 June 1997 +* +* Copyright (C) 1997 Rutherford Appleton Laboratory +*- + + IMPLICIT NONE + + DOUBLE PRECISION TDK,PMB,RH,WL,REFA,REFB + + LOGICAL OPTIC + DOUBLE PRECISION T,P,R,W,TDC,PS,PW,WLSQ,GAMMA,BETA + + + +* Decide whether optical/IR or radio case: switch at 100 microns. + OPTIC = WL.LE.100D0 + +* Restrict parameters to safe values. + T = MIN(MAX(TDK,100D0),500D0) + P = MIN(MAX(PMB,0D0),10000D0) + R = MIN(MAX(RH,0D0),1D0) + W = MIN(MAX(WL,0.1D0),1D6) + +* Water vapour pressure at the observer. + IF (P.GT.0D0) THEN + TDC = T-273.15D0 + PS = 10D0**((0.7859D0+0.03477D0*TDC)/(1D0+0.00412D0*TDC))* + : (1D0+P*(4.5D-6+6D-10*TDC*TDC)) + PW = R*PS/(1D0-(1D0-R)*PS/P) + ELSE + PW = 0D0 + END IF + +* Refractive index minus 1 at the observer. + IF (OPTIC) THEN + WLSQ = WL*WL + GAMMA = ((77.532D-6+(4.391D-7+3.57D-9/WLSQ)/WLSQ)*P + : -11.2684D-6*PW)/T + ELSE + GAMMA = (77.624D-6*P-(12.92D-6-0.371897D0/T)*PW)/T + END IF + +* Formula for beta adapted from Stone, with empirical adjustments. + BETA=4.4474D-6*T + IF (.NOT.OPTIC) BETA=BETA-0.0074D0*PW*BETA + +* Refraction constants from Green. + REFA = GAMMA*(1D0-BETA) + REFB = -GAMMA*(BETA-GAMMA/2D0) + + END -- cgit