diff options
Diffstat (limited to 'src/slalib/refro.f')
| -rw-r--r-- | src/slalib/refro.f | 374 |
1 files changed, 374 insertions, 0 deletions
diff --git a/src/slalib/refro.f b/src/slalib/refro.f new file mode 100644 index 0000000..8d11184 --- /dev/null +++ b/src/slalib/refro.f @@ -0,0 +1,374 @@ + SUBROUTINE sla_REFRO (ZOBS, HM, TDK, PMB, RH, WL, PHI, TLR, + : EPS, REF) +*+ +* - - - - - - +* R E F R O +* - - - - - - +* +* Atmospheric refraction for radio and optical/IR wavelengths. +* +* Given: +* ZOBS d observed zenith distance of the source (radian) +* HM d height of the observer above sea level (metre) +* 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) +* PHI d latitude of the observer (radian, astronomical) +* TLR d temperature lapse rate in the troposphere (degK/metre) +* EPS d precision required to terminate iteration (radian) +* +* Returned: +* REF d refraction: in vacuo ZD minus observed ZD (radian) +* +* Notes: +* +* 1 A suggested value for the TLR argument is 0.0065D0. 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. +* +* 2 A suggested value for the EPS argument is 1D-8. The result is +* usually at least two orders of magnitude more computationally +* precise than the supplied EPS value. +* +* 3 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). +* +* 4 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). +* +* . Provision for radio wavelengths has been added using +* expressions devised by A.T.Sinclair, RGO (private +* communication 1989), based on the Essen & Froome +* refractivity formula adopted in Resolution 1 of the +* 13th International Geodesy Association General Assembly +* (Bulletin Geodesique 70 p390, 1963). +* +* . Various small changes have been made to gain speed. +* +* None of the changes significantly affects the optical/IR results +* with respect to the algorithm given in the 1992 Explanatory +* Supplement. For example, at 70 deg zenith distance the present +* routine agrees with the ES algorithm to better than 0.05 arcsec +* for any reasonable combination of parameters. However, the +* improved water-vapour expressions do make a significant difference +* in the radio band, at 70 deg zenith distance reaching almost +* 4 arcsec for a hot, humid, low-altitude site during a period of +* low pressure. +* +* 5 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. +* +* 6 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. +* +* 7 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. +* +* 8 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. +* +* 9 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). +* +* Called: sla_DRANGE, sla__ATMT, sla__ATMS +* +* P.T.Wallace Starlink 3 June 1997 +* +* Copyright (C) 1997 Rutherford Appleton Laboratory +*- + + IMPLICIT NONE + + DOUBLE PRECISION ZOBS,HM,TDK,PMB,RH,WL,PHI,TLR,EPS,REF + +* +* Fixed parameters +* + DOUBLE PRECISION D93,GCR,DMD,DMW,S,DELTA,HT,HS +* 93 degrees in radians + PARAMETER (D93=1.623156204D0) +* Universal gas constant + PARAMETER (GCR=8314.32D0) +* Molecular weight of dry air + PARAMETER (DMD=28.9644D0) +* Molecular weight of water vapour + PARAMETER (DMW=18.0152D0) +* Mean Earth radius (metre) + PARAMETER (S=6378120D0) +* Exponent of temperature dependence of water vapour pressure + PARAMETER (DELTA=18.36D0) +* Height of tropopause (metre) + PARAMETER (HT=11000D0) +* Upper limit for refractive effects (metre) + PARAMETER (HS=80000D0) + + INTEGER IS,K,N,I,J + LOGICAL OPTIC,LOOP + DOUBLE PRECISION 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 + + DOUBLE PRECISION sla_DRANGE + +* The refraction integrand + DOUBLE PRECISION REFI + REFI(R,DN,RDNDR) = RDNDR/(DN+RDNDR) + + + +* Transform ZOBS into the normal range. + ZOBS1 = sla_DRANGE(ZOBS) + ZOBS2 = MIN(ABS(ZOBS1),D93) + +* Keep other arguments within safe bounds. + HMOK = MIN(MAX(HM,-1D3),10D3) + TDKOK = MIN(MAX(TDK,100D0),500D0) + PMBOK = MIN(MAX(PMB,0D0),10000D0) + RHOK = MIN(MAX(RH,0D0),1D0) + WLOK = MAX(WL,0.1D0) + ALPHA = MIN(MAX(ABS(TLR),0.001D0),0.01D0) + +* Tolerance for iteration. + TOL = MIN(MAX(ABS(EPS),1D-12),0.1D0)/2D0 + +* Decide whether optical/IR or radio case - switch at 100 microns. + OPTIC = WLOK.LE.100D0 + +* Set up model atmosphere parameters defined at the observer. + WLSQ = WLOK*WLOK + GB = 9.784D0*(1D0-0.0026D0*COS(PHI+PHI)-0.00000028D0*HMOK) + IF (OPTIC) THEN + A = (287.604D0+(1.6288D0+0.0136D0/WLSQ)/WLSQ) + : *273.15D-6/1013.25D0 + ELSE + A = 77.624D-6 + END IF + GAMAL = (GB*DMD)/GCR + GAMMA = GAMAL/ALPHA + GAMM2 = GAMMA-2D0 + DELM2 = DELTA-2D0 + TDC = TDKOK-273.15D0 + PSAT = 10D0**((0.7859D0+0.03477D0*TDC)/(1D0+0.00412D0*TDC))* + : (1D0+PMBOK*(4.5D-6+6D-10*TDC*TDC)) + IF (PMBOK.GT.0D0) THEN + PWO = RHOK*PSAT/(1D0-(1D0-RHOK)*PSAT/PMBOK) + ELSE + PWO = 0D0 + END IF + W = PWO*(1D0-DMW/DMD)*GAMMA/(DELTA-GAMMA) + C1 = A*(PMBOK+W)/TDKOK + IF (OPTIC) THEN + C2 = (A*W+11.2684D-6*PWO)/TDKOK + ELSE + C2 = (A*W+12.92D-6*PWO)/TDKOK + END IF + C3 = (GAMMA-1D0)*ALPHA*C1/TDKOK + C4 = (DELTA-1D0)*ALPHA*C2/TDKOK + IF (OPTIC) THEN + C5 = 0D0 + C6 = 0D0 + ELSE + C5 = 371897D-6*PWO/TDKOK + C6 = C5*DELM2*ALPHA/(TDKOK*TDKOK) + END IF + +* Conditions at the observer. + R0 = S+HMOK + CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6, + : R0,TEMPO,DN0,RDNDR0) + SK0 = DN0*R0*SIN(ZOBS2) + F0 = REFI(R0,DN0,RDNDR0) + +* Conditions in the troposphere at the tropopause. + RT = S+HT + CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6, + : RT,TT,DNT,RDNDRT) + SINE = SK0/(RT*DNT) + ZT = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) + FT = REFI(RT,DNT,RDNDRT) + +* Conditions in the stratosphere at the tropopause. + CALL sla__ATMS(RT,TT,DNT,GAMAL,RT,DNTS,RDNDRP) + SINE = SK0/(RT*DNTS) + ZTS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) + FTS = REFI(RT,DNTS,RDNDRP) + +* Conditions at the stratosphere limit. + RS = S+HS + CALL sla__ATMS(RT,TT,DNT,GAMAL,RS,DNS,RDNDRS) + SINE = SK0/(RS*DNS) + ZS = ATAN2(SINE,SQRT(MAX(1D0-SINE*SINE,0D0))) + FS = REFI(RS,DNS,RDNDRS) + +* +* Integrate the refraction integral in two parts; first in the +* troposphere (K=1), then in the stratosphere (K=2). +* + +* Initialize previous refraction to ensure at least two iterations. + REFOLD = 1D6 + +* Start off with 8 strips for the troposphere integration, and then +* use the final troposphere value for the stratosphere integration, +* which tends to need more strips. + IS = 8 + +* Troposphere then stratosphere. + DO K = 1,2 + +* Start Z, Z range, and start and end values. + IF (K.EQ.1) THEN + Z0 = ZOBS2 + ZRANGE = ZT-Z0 + FB = F0 + FF = FT + ELSE + Z0 = ZTS + ZRANGE = ZS-Z0 + FB = FTS + FF = FS + END IF + +* Sums of odd and even values. + FO = 0D0 + FE = 0D0 + +* First time through the loop we have to do every point. + N = 1 + +* Start of iteration loop (terminates at specified precision). + LOOP = .TRUE. + DO WHILE (LOOP) + +* Strip width. + H = ZRANGE/DBLE(IS) + +* Initialize distance from Earth centre for quadrature pass. + IF (K.EQ.1) THEN + R = R0 + ELSE + R = RT + END IF + +* One pass (no need to compute evens after first time). + DO I=1,IS-1,N + +* Sine of observed zenith distance. + SZ = SIN(Z0+H*DBLE(I)) + +* Find R (to the nearest metre, maximum four iterations). + IF (SZ.GT.1D-20) THEN + W = SK0/SZ + RG = R + DR = 1D6 + J = 0 + DO WHILE (ABS(DR).GT.1D0.AND.J.LT.4) + J=J+1 + IF (K.EQ.1) THEN + CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2, + : C1,C2,C3,C4,C5,C6,RG,TG,DN,RDNDR) + ELSE + CALL sla__ATMS(RT,TT,DNT,GAMAL,RG,DN,RDNDR) + END IF + DR = (RG*DN-W)/(DN+RDNDR) + RG = RG-DR + END DO + R = RG + END IF + +* Find the refractive index and integrand at R. + IF (K.EQ.1) THEN + CALL sla__ATMT(R0,TDKOK,ALPHA,GAMM2,DELM2, + : C1,C2,C3,C4,C5,C6,R,T,DN,RDNDR) + ELSE + CALL sla__ATMS(RT,TT,DNT,GAMAL,R,DN,RDNDR) + END IF + F = REFI(R,DN,RDNDR) + +* Accumulate odd and (first time only) even values. + IF (N.EQ.1.AND.MOD(I,2).EQ.0) THEN + FE = FE+F + ELSE + FO = FO+F + END IF + END DO + +* Evaluate the integrand using Simpson's Rule. + REFP = H*(FB+4D0*FO+2D0*FE+FF)/3D0 + +* Has the required precision been achieved? + IF (ABS(REFP-REFOLD).GT.TOL) THEN + +* No: prepare for next iteration. + +* Save current value for convergence test. + REFOLD = REFP + +* Double the number of strips. + IS = IS+IS + +* Sum of all current values = sum of next pass's even values. + FE = FE+FO + +* Prepare for new odd values. + FO = 0D0 + +* Skip even values next time. + N = 2 + ELSE + +* Yes: save troposphere component and terminate the loop. + IF (K.EQ.1) REFT = REFP + LOOP = .FALSE. + END IF + END DO + END DO + +* Result. + REF = REFT+REFP + IF (ZOBS1.LT.0D0) REF = -REF + + END |