From c1f94e4d95ba158aafa3c3a784629072620b3be0 Mon Sep 17 00:00:00 2001 From: yosata Date: Wed, 5 Aug 2020 11:37:49 +0200 Subject: [PATCH] ported slalib.sla_refro, some refactoring --- AstroAtmosphere/__init__.py | 4 +- AstroAtmosphere/dispersionModels.py | 4 +- AstroAtmosphere/examples.py | 207 ------- AstroAtmosphere/misc.py | 364 ++++++++++++ AstroAtmosphere/neoslalib.py | 754 +++++++++++++++++++++++++ AstroAtmosphere/refractivityModels.py | 104 +++- README.md | 14 +- scripts/compare_refractivity_models.py | 20 +- setup.py | 4 +- 9 files changed, 1241 insertions(+), 234 deletions(-) delete mode 100644 AstroAtmosphere/examples.py create mode 100644 AstroAtmosphere/misc.py create mode 100644 AstroAtmosphere/neoslalib.py diff --git a/AstroAtmosphere/__init__.py b/AstroAtmosphere/__init__.py index ac46895..b6119c9 100644 --- a/AstroAtmosphere/__init__.py +++ b/AstroAtmosphere/__init__.py @@ -5,6 +5,6 @@ from .ciddorModel import * from .dispersionModels import * from .refractionModels import * from .refractionModels import * -from .examples import * +from .misc import * from .refractivityModels import * - +from .neoslalib import * diff --git a/AstroAtmosphere/dispersionModels.py b/AstroAtmosphere/dispersionModels.py index d61665a..f577013 100644 --- a/AstroAtmosphere/dispersionModels.py +++ b/AstroAtmosphere/dispersionModels.py @@ -716,8 +716,8 @@ class dispersion(Observatory): _h = heightData - _rList = [] - zenith = np.deg2rad(zenith) + _rList = [] + zenith = np.deg2rad(zenith) # Integration of the refraction integral over the atmosphere layers. for _n in [n1, n2]: diff --git a/AstroAtmosphere/examples.py b/AstroAtmosphere/examples.py deleted file mode 100644 index d13257a..0000000 --- a/AstroAtmosphere/examples.py +++ /dev/null @@ -1,207 +0,0 @@ -#!/usr/bin/env Python3 -from .ciddorModel import Observatory -from .refractionModels import refraction -from .dispersionModels import dispersion - -''' -@ Author: Joost van den Born -@ Contact: born@astron.nl -@ Description: Some functions for quick testing. -''' - - -def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325, - H=0.0, xc=450, lat=0, h=0): - ''' - Function to quickly calculate the refractive index - of atmospheric air at reference conditions, given - a wavelength l. - - Parameters - ---------- - l : float - Wavelength in microns - conditions: string - 'STANDARD' or 'CERRO_ARMAZONES', refers - to the standard atmospheric conditions. - T : float (optional) - Temperature in Kelvin - p : float (optional) - Atmospheric pressure in Pa - H : float (optional) - Relative humidity - xc : float (optional) - CO2 density in parts per million - lat : float (optional) - Latitude of the observer in degrees - h : float (optional) - Altitude of the observer in meters - - Returns - ------- - n : float - Refractive index for the given conditions. - ''' - if conditions == 'STANDARD': - # Parameters at Cerro Armazones - T = 288.15 # k - p = 101325 # Pa - H = 0.0 - xc = 450 # ppm - lat = 0 # degrees - h = 0 # m - - elif conditions == 'CERRO_ARMAZONES': - # Parameters at Cerro Armazones - T = 279.65 # k - p = 71200 # Pa - H = 0.22 - xc = 450 # ppm - lat = -24.5983 # degrees - h = 3064 # m - - # Initializing dispersion model - at = Observatory() - - # Calculating indices of refraction for l - n = at.n_tph(l=l, T=T, p=p, RH=H, xc=xc) - return n - - -def quick_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325, - H=0.0, xc=450, lat=0, h=0): - ''' - Function to quickly calculate the refraction - of atmospheric air at reference conditions, - given a wavelength l. - - Parameters - ---------- - l : float - Wavelength in microns - zenith : float - Angle of observation in degrees. - zenith = 0 means the observation is directly overhead. - conditions : string - 'STANDARD' or 'CERRO_ARMAZONES', refers - to the standard atmospheric conditions. - T : float (optional) - Temperature in Kelvin - p : float (optional) - Atmospheric pressure in Pa - H : float (optional) - Relative humidity - xc : float (optional) - CO2 density in parts per million - lat : float (optional) - Latitude of the observer in degrees - h : float (optional) - Altitude of the observer in meters - - Returns - ------- - R : float - Refraction in degrees - ''' - if conditions == 'STANDARD': - # Parameters at Cerro Armazones - T = 288.15 # k - p = 101325 # Pa - H = 0.0 - xc = 450 # ppm - lat = 0 # degrees - h = 0 # m - - elif conditions == 'CERRO_ARMAZONES': - # Parameters at Cerro Armazones - T = 279.65 # k - p = 71200 # Pa - H = 0.22 - xc = 450 # ppm - lat = -24.5983 # degrees - h = 3064 # m - - # Initializing dispersion model - at = Observatory() - - # Calculating indices of refraction for l - n = at.n_tph(l=l, T=T, p=p, RH=H, xc=xc) - - # Density of the atmosphere (following CIPM-81/91 equations) - rho = at.rho(p=p, T=T, RH=H, xc=xc) - - # Initializing refraction model and setting the reduced height - ref = refraction(lat, h) - ref.setReducedHeight(p, rho) - return ref.cassini(n, zenith) - - -def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, - p=101325, H=0.0, xc=450, lat=0, h=0): - ''' - Function to quickly calculate the dispersion - of atmospheric air at reference conditions, - given a wavelength l. - - Parameters - ---------- - l1 : float - Shortest wavelength of interest in micron. - l2 : float - Longest wavelength of interest in micron. - zenith : float - Angle of observation in degrees. - zenith = 0 means the observation is directly overhead. - conditions : string - 'STANDARD' or 'CERRO_ARMAZONES', refers - to the standard atmospheric conditions. - T : float (optional) - Temperature in Kelvin - p : float (optional) - Atmospheric pressure in Pa - H : float (optional) - Relative humidity - xc : float (optional) - CO2 density in parts per million - lat : float (optional) - Latitude of the observer in degrees - h : float (optional) - Altitude of the observer in meters - - Returns - ------- - dR : float - Atmospheric dispersion in degrees - ''' - if conditions == 'STANDARD': - # Parameters at Cerro Armazones - T = 288.15 # k - p = 101325 # Pa - H = 0.0 - xc = 450 # ppm - lat = 0 # degrees - h = 0 # m - - elif conditions == 'CERRO_ARMAZONES': - # Parameters at Cerro Armazones - T = 279.65 # k - p = 71200 # Pa - H = 0.22 - xc = 450 # ppm - lat = -24.5983 # degrees - h = 3064 # m - - # Initializing dispersion model - at = Observatory() - - # Calculating indices of refraction for l1 and l2 - n1 = at.n_tph(l=l1, T=T, p=p, RH=H, xc=xc) - n2 = at.n_tph(l=l2, T=T, p=p, RH=H, xc=xc) - - # Density of the atmosphere (following CIPM-81/91 equations) - rho = at.rho(p=p, T=T, RH=H, xc=xc) - - # Initializing refraction model and setting the reduced height - disp = dispersion(lat, h) - disp.setReducedHeight(p, rho) - return disp.cassini(n1, n2, zenith) \ No newline at end of file diff --git a/AstroAtmosphere/misc.py b/AstroAtmosphere/misc.py new file mode 100644 index 0000000..f71ff93 --- /dev/null +++ b/AstroAtmosphere/misc.py @@ -0,0 +1,364 @@ +#!/usr/bin/env Python3 +from .ciddorModel import Observatory +from .refractionModels import refraction +from .dispersionModels import dispersion +from .neoslalib import sla_refro + +''' +@ Author: Joost van den Born +@ Contact: born@astron.nl +@ Description: Some functions for quick testing. +''' + + +def quick_refractive_index(l, conditions='STANDARD', T=288.15, p=101325, + RH=0.0, xc=450, lat=0, h=0): + ''' + Function to quickly calculate the refractive index + of atmospheric air at reference conditions, given + a wavelength l. + + Parameters + ---------- + l : float + Wavelength in microns + conditions: string + 'STANDARD' or 'CERRO_ARMAZONES', refers + to the standard atmospheric conditions. + T : float (optional) + Temperature in Kelvin + p : float (optional) + Atmospheric pressure in Pa + RH : float (optional) + Relative humidity + xc : float (optional) + CO2 density in parts per million + lat : float (optional) + Latitude of the observer in degrees + h : float (optional) + Altitude of the observer in meters + + Returns + ------- + n : float + Refractive index for the given conditions. + ''' + if conditions == 'STANDARD': + # Parameters at Cerro Armazones + T = 288.15 # k + p = 101325 # Pa + RH = 0.0 + xc = 450 # ppm + lat = 0 # degrees + h = 0 # m + + elif conditions == 'CERRO_ARMAZONES': + # Parameters at Cerro Armazones + T = 279.65 # k + p = 71200 # Pa + RH = 0.22 + xc = 450 # ppm + lat = -24.5983 # degrees + h = 3064 # m + + # Initializing dispersion model + at = Observatory() + + # Calculating indices of refraction for l + n = at.n_tph(l=l, T=T, p=p, RH=RH, xc=xc) + return n + + +def quick_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325, + RH=0.0, xc=450, lat=0, h=0): + ''' + Function to quickly calculate the refraction + of atmospheric air at reference conditions, + given a wavelength l. + + Parameters + ---------- + l : float + Wavelength in microns + zenith : float + Angle of observation in degrees. + zenith = 0 means the observation is directly overhead. + conditions : string + 'STANDARD' or 'CERRO_ARMAZONES', refers + to the standard atmospheric conditions. + T : float (optional) + Temperature in Kelvin + p : float (optional) + Atmospheric pressure in Pa + RH : float (optional) + Relative humidity + xc : float (optional) + CO2 density in parts per million + lat : float (optional) + Latitude of the observer in degrees + h : float (optional) + Altitude of the observer in meters + + Returns + ------- + R : float + Refraction in degrees + ''' + if conditions == 'STANDARD': + # Parameters at Cerro Armazones + T = 288.15 # k + p = 101325 # Pa + RH = 0.0 + xc = 450 # ppm + lat = 0 # degrees + h = 0 # m + + elif conditions == 'CERRO_ARMAZONES': + # Parameters at Cerro Armazones + T = 279.65 # k + p = 71200 # Pa + RH = 0.22 + xc = 450 # ppm + lat = -24.5983 # degrees + h = 3064 # m + + # Initializing dispersion model + at = Observatory() + + # Calculating indices of refraction for l + n = at.n_tph(l=l, T=T, p=p, RH=RH, xc=xc) + + # Density of the atmosphere (following CIPM-81/91 equations) + rho = at.rho(p=p, T=T, RH=RH, xc=xc) + + # Initializing refraction model and setting the reduced height + ref = refraction(lat, h) + ref.setReducedHeight(p, rho) + return ref.cassini(n, zenith) + + +def quick_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, + p=101325, RH=0.0, xc=450, lat=0, h=0): + ''' + Function to quickly calculate the dispersion + of atmospheric air at reference conditions, + given a wavelength l. + + Parameters + ---------- + l1 : float + Shortest wavelength of interest in micron. + l2 : float + Longest wavelength of interest in micron. + zenith : float + Angle of observation in degrees. + zenith = 0 means the observation is directly overhead. + conditions : string + 'STANDARD' or 'CERRO_ARMAZONES', refers + to the standard atmospheric conditions. + T : float (optional) + Temperature in Kelvin + p : float (optional) + Atmospheric pressure in Pa + RH : float (optional) + Relative humidity + xc : float (optional) + CO2 density in parts per million + lat : float (optional) + Latitude of the observer in degrees + h : float (optional) + Altitude of the observer in meters + + Returns + ------- + dR : float + Atmospheric dispersion in degrees + ''' + if conditions == 'STANDARD': + # Parameters at Cerro Armazones + T = 288.15 # k + p = 101325 # Pa + RH = 0.0 + xc = 450 # ppm + lat = 0 # degrees + h = 0 # m + + elif conditions == 'CERRO_ARMAZONES': + # Parameters at Cerro Armazones + T = 279.65 # k + p = 71200 # Pa + RH = 0.22 + xc = 450 # ppm + lat = -24.5983 # degrees + h = 3064 # m + + # Initializing dispersion model + at = Observatory() + + # Calculating indices of refraction for l1 and l2 + n1 = at.n_tph(l=l1, T=T, p=p, RH=RH, xc=xc) + n2 = at.n_tph(l=l2, T=T, p=p, RH=RH, xc=xc) + + # Density of the atmosphere (following CIPM-81/91 equations) + rho = at.rho(p=p, T=T, RH=RH, xc=xc) + + # Initializing refraction model and setting the reduced height + disp = dispersion(lat, h) + disp.setReducedHeight(p, rho) + return disp.cassini(n1, n2, zenith) + + +def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325, + RH=0.0, xc=450, lat=0, h=0, tlr=0.0065, eps=1e-8, use_pyslalib=False): + ''' + Calculates the atmospheric refraction assuming the SLALIB + package. SLALIB uses the same atmospheric model as ZEMAX + based on the works of Hohenkerk & Sinclair and Seidelmann. + They take the B&S refraction formula as a basis. + + The pySLALIB FORTRAN has been ported, so pyslalib is not + needed, but can still be used. There are minor numerical + differences, likely due to precision or declaration issues. + + Parameters + ---------- + l : float + Wavelength in microns + zenith : float + Angle of observation in degrees. + zenith = 0 means the observation is directly overhead. + conditions : string + 'STANDARD' or 'CERRO_ARMAZONES', refers + to the standard atmospheric conditions. + T : float (optional) + Temperature in Kelvin + p : float (optional) + Atmospheric pressure in Pa + RH : float (optional) + Relative humidity + xc : float (optional) + CO2 density in parts per million + lat : float (optional) + Latitude of the observer in degrees (default=0) + h : float (optional) + Altitude of the observer in meters + tlr : float + Temperature lapse rate (default=0.0065) + eps : float + Precision at which the calculation is aborted (default=1e-10) + use_pyslalib : bool + Flag that indicates that the user wants to use the pySlalib + package. If False, the internal python port of the original + FORTRAN code will be used. (Default=False). + + Returns + ------- + R : float + Refraction in degrees + ''' + + if conditions == 'STANDARD': + # Parameters at Cerro Armazones + T = 288.15 # k + p = 101325 # Pa + RH = 0.0 + xc = 450 # ppm + lat = 0 # degrees + h = 0 # m + + elif conditions == 'CERRO_ARMAZONES': + # Parameters at Cerro Armazones + T = 279.65 # k + p = 71200 # Pa + RH = 0.22 + xc = 450 # ppm + lat = -24.5983 # degrees + h = 3064 # m + + # Converting units to the units that sla_refro() requires. + _zenith = np.deg2rad(zenith) # From degrees to radians + _P = p / 100 # From Pa to mbar + _phi = np.deg2rad(lat) # From degrees to radians + + R = [] + if type(l) == float or type(l) == int: + l = [l] + if use_pyslalib: + try: + from pyslalib import slalib + except ImportError as e: + print('pySLALIB package not found! Try installing it first.') + for wl in l: + # If using original slalib (e.g. PySLALIB) + Refrac = slalib.sla_refro(_zenith, h, T, _P, RH, wl, phi, tlr, eps) + R.append(np.degrees(Refrac)) + + # If using internal slalib port + else: + for wl in l: + Refrac = sla_refro(_zenith, h, T, _P, RH, wl, phi, tlr, eps) + R.append(np.degrees(Refrac)) + + if len(l) == 1: + return R[0] + else: + return np.asarray(R) + +def slalib_dispersion(l1, l2, zenith, conditions='STANDARD', T=288.15, p=101325, + RH=0.0, xc=450, lat=0, h=0, tlr=0.0065, eps=1e-8, use_pyslalib=False): + ''' + Calculates the atmospheric refraction assuming the SLALIB + package. SLALIB uses the same atmospheric model as ZEMAX + based on the works of Hohenkerk & Sinclair and Seidelmann. + They take the B&S refraction formula as a basis. + + The pySLALIB FORTRAN has been ported, so pyslalib is not + needed, but can still be used. There are minor numerical + differences, likely due to precision or declaration issues. + + Parameters + ---------- + l1 : float + Shortest wavelength of interest in micron. + l2 : float + Longest wavelength of interest in micron. + zenith : float + Angle of observation in degrees. + zenith = 0 means the observation is directly overhead. + conditions : string + 'STANDARD' or 'CERRO_ARMAZONES', refers + to the standard atmospheric conditions. + T : float (optional) + Temperature in Kelvin + p : float (optional) + Atmospheric pressure in Pa + RH : float (optional) + Relative humidity + xc : float (optional) + CO2 density in parts per million + lat : float (optional) + Latitude of the observer in degrees (default=0) + h : float (optional) + Altitude of the observer in meters + tlr : float + Temperature lapse rate (default=0.0065) + eps : float + Precision at which the calculation is aborted (default=1e-10) + use_pyslalib : bool + Flag that indicates that the user wants to use the pySlalib + package. If False, the internal python port of the original + FORTRAN code will be used. (Default=False). + + Returns + ------- + dR : float + Refraction in degrees + ''' + R1 = slalib_refraction(l1, zenith, conditions, T, p, RH, xc,\ + lat, h, tlr, eps, use_pyslalib) + + R2 = slalib_refraction(l2, zenith, conditions, T, p, RH, xc,\ + lat, h, tlr, eps, use_pyslalib) + + dR = R1 - R2 + return dR diff --git a/AstroAtmosphere/neoslalib.py b/AstroAtmosphere/neoslalib.py new file mode 100644 index 0000000..0488fdf --- /dev/null +++ b/AstroAtmosphere/neoslalib.py @@ -0,0 +1,754 @@ +#!/usr/bin/env Python3 +import numpy as np + +def sla_atms(RT, TT, DNT, GAMAL, R): + ''' + Internal routine used by sla_refro(), ported from the SLALIB + FORTRAN code (Copyright P.T.Wallace. All rights reserved.) + + The method is also described in "Explanatory Supplement + Astronomical Almanac" by P.K. Seidelmann, 1992, section 3.281. + + Refractive index and derivative with respect to height for the + stratosphere. + + Parameters + ---------- + RT : float + Height of tropopause from center of the Earth (meters) + TT : float + Temperature at the tropopause (K) + DNT : float + Refractive index at the tropopause + GAMAL : float + Constant of the atmospheric model = G*MD/R + R : float + Current distance from the centre of the Earth (metre) + + Returns + ------- + DN : float + Refractive index at R + RDNDR : float + R * rate the refractive index is changing at R + ''' + _B = GAMAL / TT + _W = (DNT - 1) * np.exp( -B * (R - RT) ) + DN = 1+_W + RDNDR = -R * _B * _W + return DN, RDNDR + + +def sla_atmt(R0, T0, ALPHA, GAMM2, DELM2, C1, C2, C3, C4, C5, C6, R): + ''' + Internal routine used by sla_refro(), ported from the SLALIB + FORTRAN code (Copyright P.T.Wallace. All rights reserved.) + + The method is also described in "Explanatory Supplement + Astronomical Almanac" by P.K. Seidelmann, 1992, section 3.281. + + Refractive index and derivative with respect to height for the + troposphere. + + Parameters + ---------- + R0 : float + Height of observer from center of the Earth (meter) + T0 : float + Temperature at the observer (K) + ALPHA : float + Alpha + GAMM2 : float + Gamma minus 2 - see HMNAO paper + DELM2 : float + Delta minus 2 + C1 : float + Useful term + C2 : float + Useful term + C3 : float + Useful term + C4 : float + Useful term + C5 : float + Useful term + C6 : float + Useful term + R : float + Current distance from the center of the Earth (meters) + + Returns + ------- + T : float + Temperature at R (K) + DN : float + Refractive index at R + RDNDR : float + R * rate the refractive index is changing at R + + Note that in the optical case C5 and C6 are zero. + ''' + T = max(min(T0 - ALPHA * (R - R0), 320), 100) + TT0 = T / T0 + TT0GM2 = TT0**GAMM2 + TT0DM2 = TT0**DELM2 + DN = 1 + ( C1 * TT0GM2 - (C2 - C5 / T) * TT0DM2) * TT0 + RDNDR = R * (-C3 * TT0GM2 + (C4 - C6 / TT0) * TT0DM2) + return T, DN, RDNDR + + +def sla_drange(a): + ''' + Helper function for sla_refro(). + ''' + drange = a % (2*np.pi) + if drange >= np.pi: + drange -= np.abs(2*np.pi)*np.sign(a) + return drange + + +def sla_refi(dn, rdndr): + ''' + Helper function for sla_refro(), that calculates + the refraction integrand. + ''' + return rdndr/(dn+rdndr) + +def sla_refro(ZOBS, HM, TDK, PMB, RH, WL, PHI, TLR, EPS): + ''' + This is a direct port from the original FORTRAN SLA_REFRO() + routine. + + Calculates Atmospheric refraction for radio and + optical/IR wavelengths, ported from the original SLALIB + FORTRAN code (Copyright P.T.Wallace. All rights reserved.). + + The method is also described in "Explanatory Supplement + Astronomical Almanac" by P.K. Seidelmann, 1992, section 3.281. + + Parameters + ---------- + ZOBS : float + Observed zenith distance of the source (radian) + HM : float + Height of the observer above sea level (meter) + TDK : float + Ambient temperature at the observer (K) + PMB : float + Pressure at the observer (millibar) + RH : float + Relative humidity at the observer (range 0-1) + WL : float + Effective wavelength of the source (micrometre) + PHI : float + Latitude of the observer (radian, astronomical) + TLR : float + Temperature lapse rate in the troposphere (K/metre) + EPS : float + Precision required to terminate iteration (radian) + + Returns + ------- + REF : float + 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. The sign of the supplied TLR value + is ignored. + 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). + . The formula for the water vapour pressure, given the + saturation pressure and the relative humidity, is from + Crane (1976), expression 2.5.5. + . Provision for radio wavelengths has been added using + expressions devised by A.T.Sinclair, RGO (private + communication 1989). The refractivity model currently + used is from J.M.Rueger, "Refractive Index Formulae for + Electronic Distance Measurement with Radio and Millimetre + Waves", in Unisurv Report S-68 (2002), School of Surveying + and Spatial Information Systems, University of New South + Wales, Sydney, Australia. + . The optical refractivity for dry air is from Resolution 3 of + the International Association of Geodesy adopted at the XXIIth + General Assembly in Birmingham, UK, 1999. + . Various small changes have been made to gain speed. + 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). + 10 The algorithm is designed for observers in the troposphere. The + supplied temperature, pressure and lapse rate are assumed to be + for a point in the troposphere and are used to define a model + atmosphere with the tropopause at 11km altitude and a constant + temperature above that. However, in practice, the refraction + values returned for stratospheric observers, at altitudes up to + 25km, are quite usable. + ''' + + # Fixed parameters + D93 = np.deg2rad(93)# 93 degrees in radians + GCR = 8314.32 # Universal gas constant + DMD = 28.9644 # Molecular weight of dry air + DMW = 18.0152 # Molecular weight of water vapour + S = 6378120 # Mean Earth radius (metre) + DELTA = 18.36 # Exponent of temperature dependence of water vapour pressure + HT = 11000 # Height of tropopause (metre) + HS = 80000 # Upper limit for refractive effects (metre) + ISMAX = 16384 # Numerical integration: maximum number of strips. + + # 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,-1E3),HS) + TDKOK = min(max(TDK,100),500) + PMBOK = min(max(PMB,0),10000) + RHOK = min(max(RH,0),1) + WLOK = max(WL,0.1) + ALPHA = min(max(abs(TLR),0.001),0.01) + + # Tolerance for iteration. + TOL = min(max(abs(EPS),1E-12),0.1)/2 + + # Decide whether optical/IR or radio case - switch at 100 microns. + OPTIC = (WLOK<=100) + + # Set up model atmosphere parameters defined at the observer. + WLSQ = WLOK*WLOK + + # Local gravity, can be made consistent with AstroAtmosphere + GB = 9.784*(1-0.0026*np.cos(PHI+PHI)-0.00000028*HMOK) + + if OPTIC: + # Refractivity from Barrell & Sears (1939)! + A = (287.6155+(1.62887+0.01360/WLSQ)/WLSQ)*273.15E-6/1013.25 + else: + A = 77.6890E-6 + + GAMAL = (GB*DMD)/GCR + GAMMA = GAMAL/ALPHA + GAMM2 = GAMMA-2 + DELM2 = DELTA-2 + TDC = TDKOK-273.15 + PSAT = 10**((0.7859+0.03477*TDC)/(1+0.00412*TDC))\ + *(1+PMBOK*(4.5E-6+6E-10*TDC*TDC)) + + if PMBOK > 0: + PWO = RHOK*PSAT/(1-(1-RHOK)*PSAT/PMBOK) + else: + PWO = 0 + + W = PWO*(1-DMW/DMD)*GAMMA/(DELTA-GAMMA) + C1 = A*(PMBOK+W)/TDKOK + if OPTIC: + C2 = (A*W+11.2684E-6*PWO)/TDKOK + else: + C2 = (A*W+6.3938E-6*PWO)/TDKOK + + C3 = (GAMMA-1)*ALPHA*C1/TDKOK + C4 = (DELTA-1)*ALPHA*C2/TDKOK + if OPTIC: + C5 = 0 + C6 = 0 + else: + C5 = 375463E-6*PWO/TDKOK + C6 = C5*DELM2*ALPHA/(TDKOK*TDKOK) + + # Conditions at the observer. + R0 = S+HMOK + TEMPO, DN0, RDNDR0 = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,R0) + SK0 = DN0*R0*np.sin(ZOBS2) + F0 = sla_refi(DN0,RDNDR0) + + # Conditions in the troposphere at the tropopause. + RT = S+max(HT,HMOK) + TT, DNT, RDNDRT = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,RT) + SINE = SK0/(RT*DNT) + ZT = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0))) + FT = sla_refi(DNT,RDNDRT) + + # Conditions in the stratosphere at the tropopause. + DNTS,RDNDRP = sla_atms(RT,TT,DNT,GAMAL,RT) + SINE = SK0/(RT*DNTS) + ZTS = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0))) + FTS = sla_refi(DNTS,RDNDRP) + + # Conditions at the stratosphere limit. + RS = S+HS + DNS,RDNDRS = sla_atms(RT,TT,DNT,GAMAL,RS) + SINE = SK0/(RS*DNS) + ZS = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0))) + FS = sla_refi(DNS,RDNDRS) + + # Variable initialization to avoid compiler warning. + REFT = 0 + + # Integrate the refraction integral in two parts; first in the + # troposphere (K=1), then in the stratosphere (K=2). + + for K in [1,2]: + # Initialize previous refraction to ensure at least two iterations. + REFOLD = 1 + + # Start off with 8 strips. + IS = 8 + + # Start Z, Z range, and start and end values. + if K==1: + Z0 = ZOBS2 + ZRANGE = ZT-Z0 + FB = F0 + FF = FT + else: + Z0 = ZTS + ZRANGE = ZS-Z0 + FB = FTS + FF = FS + + # Sums of odd and even values. + FO = 0 + FE = 0 + + # First time through the loop we have to do every point. + N = 1 + + # Start of iteration loop (terminates at specified precision). + LOOP = True + while LOOP: + + # Strip width. + H = ZRANGE/IS + + # Initialize distance from Earth centre for quadrature pass. + if K==1: + R = R0 + else: + R = RT + + # One pass (no need to compute evens after first time). + # DO I=1,IS-1,N + for I in np.arange(1, IS-1, N): + + # Sine of observed zenith distance. + SZ = np.sin(Z0+H*I) + + # Find R (to the nearest metre, maximum four iterations). + if (SZ>1E-20): + W = SK0/SZ + RG = R + DR = 1E6 + J = 0 + while (abs(DR)>1 and J<4): + J+=1 + if K==1: + TG,DN,RDNDR = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,RG) + else: + DN,RDNDR = sla_atms(RT,TT,DNT,GAMAL,RG) + DR = (RG*DN-W)/(DN+RDNDR) + RG = RG-DR + R = RG + + # Find the refractive index and integrand at R. + if K==1: + T,DN,RDNDR = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,R) + else: + DN,RDNDR = sla_atms(RT,TT,DNT,GAMAL,R) + + F = sla_refi(DN,RDNDR) + + # Accumulate odd and (first time only) even values. + if N==1 and (I%2)==0: + FE = FE+F + else: + FO = FO+F + + + # Evaluate the integrand using Simpson's Rule. + REFP = H*(FB+4*FO+2*FE+FF)/3 + + # Has the required precision been achieved (or can't be)? + if (abs(REFP-REFOLD)>TOL and IS 0: +# PWO = RHOK*PSAT/(1-(1-RHOK)*PSAT/PMBOK) +# else: +# PWO = 0 + +# W = PWO*(1-DMW/DMD)*GAMMA/(DELTA-GAMMA) +# C1 = A*(PMBOK+W)/TDKOK +# C2 = (A*W+11.2684E-6*PWO)/TDKOK + +# C3 = (GAMMA-1)*ALPHA*C1/TDKOK +# C4 = (DELTA-1)*ALPHA*C2/TDKOK +# C5 = 0 +# C6 = 0 + +# # Conditions at the observer. +# R0 = S+HMOK +# TEMPO, DN0, RDNDR0 = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,R0) +# SK0 = DN0*R0*np.sin(ZOBS2) +# F0 = sla_refi(DN0,RDNDR0) + +# # Conditions in the troposphere at the tropopause. +# RT = S+max(HT,HMOK) +# TT, DNT, RDNDRT = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,RT) +# SINE = SK0/(RT*DNT) +# ZT = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0))) +# FT = sla_refi(DNT,RDNDRT) + +# # Conditions in the stratosphere at the tropopause. +# DNTS,RDNDRP = sla_atms(RT,TT,DNT,GAMAL,RT) +# SINE = SK0/(RT*DNTS) +# ZTS = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0))) +# FTS = sla_refi(DNTS,RDNDRP) + +# # Conditions at the stratosphere limit. +# RS = S+HS +# DNS,RDNDRS = sla_atms(RT,TT,DNT,GAMAL,RS) +# SINE = SK0/(RS*DNS) +# ZS = np.arctan2(SINE,np.sqrt(max(1-SINE*SINE,0))) +# FS = sla_refi(DNS,RDNDRS) + +# # Variable initialization to avoid compiler warning. +# REFT = 0 + +# # Integrate the refraction integral in two parts; first in the +# # troposphere (K=1), then in the stratosphere (K=2). +# for K in [1,2]: +# # Initialize previous refraction to ensure at least two iterations. +# REFOLD = 1 + +# # Start off with 8 strips. +# IS = 8 + +# # Start Z, Z range, and start and end values. +# if K==1: +# Z0 = ZOBS2 +# ZRANGE = ZT-Z0 +# FB = F0 +# FF = FT +# else: +# Z0 = ZTS +# ZRANGE = ZS-Z0 +# FB = FTS +# FF = FS + +# # Sums of odd and even values. +# FO = 0 +# FE = 0 + +# # First time through the loop we have to do every point. +# N = 1 + +# # Start of iteration loop (terminates at specified precision). +# LOOP = True +# while LOOP: + +# # Strip width. +# H = ZRANGE/IS + +# # Initialize distance from Earth centre for quadrature pass. +# if K==1: +# R = R0 +# else: +# R = RT + +# # One pass (no need to compute evens after first time). +# # DO I=1,IS-1,N +# for I in np.arange(1, IS-1, N): + +# # Sine of observed zenith distance. +# SZ = np.sin(Z0+H*I) + +# # Find R (to the nearest metre, maximum four iterations). +# if (SZ>1E-20): +# W = SK0/SZ +# RG = R +# DR = 1E6 +# J = 0 +# while (abs(DR)>1 and J<4): +# J+=1 +# if K==1: +# TG,DN,RDNDR = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,RG) +# else: +# DN,RDNDR = sla_atms(RT,TT,DNT,GAMAL,RG) +# DR = (RG*DN-W)/(DN+RDNDR) +# RG = RG-DR +# R = RG + +# # Find the refractive index and integrand at R. +# if K==1: +# T,DN,RDNDR = sla_atmt(R0,TDKOK,ALPHA,GAMM2,DELM2,C1,C2,C3,C4,C5,C6,R) +# else: +# DN,RDNDR = sla_atms(RT,TT,DNT,GAMAL,R) + +# F = sla_refi(DN,RDNDR) + +# # Accumulate odd and (first time only) even values. +# if N==1 and (I%2)==0: +# FE = FE+F +# else: +# FO = FO+F + + +# # Evaluate the integrand using Simpson's Rule. +# REFP = H*(FB+4*FO+2*FE+FF)/3 + +# # Has the required precision been achieved (or can't be)? +# if (abs(REFP-REFOLD)>TOL and IS