From 7771c59268287b5351f6ea3078ca4d094cc8317c Mon Sep 17 00:00:00 2001
From: yosata <javandenborn@gmail.com>
Date: Wed, 10 Mar 2021 11:50:00 +0100
Subject: [PATCH] Added TPH dependency to Edlen and BonschPotulski functions.

---
 AstroAtmosphere/refractivityModels.py  | 92 +++++++++++++++++++++-----
 scripts/compare_refractivity_models.py | 21 +++++-
 2 files changed, 97 insertions(+), 16 deletions(-)

diff --git a/AstroAtmosphere/refractivityModels.py b/AstroAtmosphere/refractivityModels.py
index c8577b0..3787ea0 100644
--- a/AstroAtmosphere/refractivityModels.py
+++ b/AstroAtmosphere/refractivityModels.py
@@ -36,9 +36,9 @@ def BarrellAndSears(l, T=288.15, p=101325, RH=0, f=None):
 	'''
 	if f is None:
 		f = RH * svp(T)
-	f_mmHg = hPA2mmHg(f/100)
+	f_mmHg = hPa2mmHg(f/100)
 	T = K2C(T)
-	P = hPA2mmHg(p/100)
+	P = hPa2mmHg(p/100)
 	R = (0.378125 + 0.0021414 / l**2 + 0.00001793 / l**4) * (P*(1 + (1.049 - 0.0157*T) * P * 1e-6)/(1 + 0.003661*T)) - (0.0624 - 0.000680 / l**2) * f_mmHg / (1 + 0.003661*T)
 	n = R/1e6 + 1
 	return n
@@ -55,6 +55,8 @@ def Ciddor(l, CO2=450):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		CO2 : 	float
+			CO2 density in parts per million
 
 		Returns
 		-------
@@ -86,6 +88,8 @@ def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
 		f 	: 	float
 			Water vapor pressure in Pascal
 
@@ -98,11 +102,11 @@ def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 	R = 6432.8 + 2949810 / (146 - s**2) + 25540/(41-s**2)
 	n_standard = R/1e8 + 1
 
-	pHg = hPA2mmHg(p/100)
+	pHg = hPa2mmHg(p/100)
 	Tc  = K2C(T)
 	if f is None:
 		f = RH * svp(T)
-	f_mmHg = hPA2mmHg(f/100)
+	f_mmHg = hPa2mmHg(f/100)
 
 	# Changes in temperature and pressure, from Barrell (1951)
 	n_tp = R/1e8 * (pHg * (1 + (1.049 - 0.0157*Tc) * 1e-6 * pHg ) ) / ( 720.883 * (1 + 0.003661 * Tc)) + 1
@@ -113,10 +117,10 @@ def Edlen1953(l, T=288.15, p=101325, RH=0, f=None):
 	return n
 
 
-def Edlen1966(l):
+def Edlen1966(l, T=288.15, p=101325, RH=0, CO2=300, f=None):
 	'''
 		From Edlen, 1966
-		Equation 1
+		Equations 1, 12, 17, 22
 
 		At standard conditions only
 
@@ -124,6 +128,16 @@ def Edlen1966(l):
 		---------- 
 		l 	: 	float
 			Wavelength in microns
+		T 	: 	float
+			Temperature in Kelvin
+		p 	: 	float
+			Atmospheric pressure in Pascal
+		RH 	: 	float
+			Relative humidity (0=<H=<1)
+		CO2 : 	float
+			CO2 density in parts per million
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -133,6 +147,27 @@ def Edlen1966(l):
 	s = 1 / l
 	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2)
 	n = R/1e8 + 1
+
+	# Influence of temperature and pressure, eq. 12
+	Tc = K2C(T)
+	Ptorr = Pa2Torr(p)
+	R_tp = ( Ptorr * R/1e8 ) / 720.775 * (1 + Ptorr * (0.817 - 0.0133 * Tc) * 1e-6) / ( 1 + 0.0036610 * Tc)
+
+	# Allow for a CO2 fraction, eq. 17, but we include T and P dependence.
+	# This is okay from "Asthe refractivity of CO, is about 50% higher than
+	# that of air this CO,-content increases the refractivity over that of 
+	# CO2-free air by approximately 0.5 x 0.0003 x (n-1)s" and "As the 
+	# density factor of CO, is not much different from that of air this ratio r
+	# can be taken as independent of t and p ."
+	R_x = R_tp * (1 + 0.540 * (CO2/1e6 - 0.0003))
+	n_x = R_x + 1
+
+	# Correction for water vapour, eq. 22
+	if f is None:
+		f = RH * svp(T)
+	ftorr = Pa2Torr(f)
+
+	n = n_x - ftorr * (5.7224 - 0.0457 * s**2) * 1e-8
 	return n
 
 
@@ -159,7 +194,7 @@ def Owens(l, T=288.15, p=101325, RH=0):
 		n 	: 	float
 			Refractive index of atmospheric air at the given conditions
 	'''
-	P = hPA2mmHg(p/100)
+	P = hPa2mmHg(p/100)
 	T = K2C(T)
 	s = 1 / l
 	R = 8342.13 + 2406030 / (130 - s**2) + 15997 / (38.9 - s**2) * (P/720.775) * ((1 + P * (0.817 - 0.0133) * 1e-6)/(1 + 0.0036610*T)) - RH * (5.722 - 0.0457*s**2)
@@ -167,10 +202,10 @@ def Owens(l, T=288.15, p=101325, RH=0):
 	return n
 
 
-def BonschPotulski(l, T=288.15, p=101325, CO2=450):
+def BonschPotulski(l, T=288.15, p=101325, RH=0, CO2=400, f=None):
 	'''
 		From Bonsch & Potulski, 1998
-		Equations 6, 7, 8
+		Equations 6, 7, 8 and 9a
 		
 		Parameters
 		---------- 
@@ -180,8 +215,12 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 			Temperature in Kelvin
 		p 	: 	float
 			Atmospheric pressure in Pa
+		RH 	: 	float
+			Relative humidity (0<=H<=1)
 		CO2 : 	float
 			Particles per million of CO2
+		f 	: 	float
+			Water vapor pressure in Pascal
 
 		Returns
 		-------
@@ -189,14 +228,19 @@ def BonschPotulski(l, T=288.15, p=101325, CO2=450):
 			Refractive index of atmospheric air at the given conditions
 	'''
 	x = CO2 / 1e6
-	T = K2C(T)
+	Tc = K2C(T) 
 	s = 1 / l
-	# Rs = 8092.33 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
+
 	Rs = 8091.37 + 2333983 / (130 - s**2) + 15518 / (38.9 - s**2)
 	Rs = Rs / 1e8
 	Rx = Rs * (1 + 0.5327 * (x - 0.0004))
-	Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * T)*p) / (1 + 0.0036610 * T)
-	n = Rtp + 1
+	Rtp = Rx * p / 93214.60 * (1 + 1e-8 * (0.5953 - 0.009876 * Tc)*p) / (1 + 0.0036610 * Tc)
+
+	if f is None:
+		f = RH * svp(T)
+	
+	n_tp = Rtp + 1
+	n = n_tp - f * (3.8020 - 0.0384 * s**2) * 1e-10
 	return n
 
 
@@ -490,7 +534,7 @@ def Mathar(l, T=288.15, p=101325, RH=0):
 		σref = 1e4/20      # cm^−1, reference wavenumber
 		Tref = 273.15+17.5 # K, reference temperature
 		pref = 75000       # Pa, reference pressure
-		Href = 10          #%, reference humdity in percent
+		Href = 10          # %, reference humdity in percent
 		########################################################################
 
 	s = 1e4/l # cm^-1
@@ -525,7 +569,7 @@ def K2C(K):
 	return K - 273.15
 
 
-def hPA2mmHg(hPa):
+def hPa2mmHg(hPa):
 	'''
 		hPa to mmHg conversion
 
@@ -542,6 +586,24 @@ def hPA2mmHg(hPa):
 	return 0.7500616827 * hPa
 
 
+def Pa2Torr(Pa):
+	'''
+		Pa to Torr conversion
+
+		Parameters
+		----------
+		Pa 	:	float
+			Atmospheric pressure in Pa (Pascal)
+
+		Returns
+		-------
+		Torr 	: 	float
+			Atmospheric pressure in Torr
+	'''
+
+	return Pa * 760/101325
+
+
 def svp(T):
 	'''
 		Saturation vapour pressure in Pa.
diff --git a/scripts/compare_refractivity_models.py b/scripts/compare_refractivity_models.py
index fa41cc6..65c1ff1 100644
--- a/scripts/compare_refractivity_models.py
+++ b/scripts/compare_refractivity_models.py
@@ -81,7 +81,26 @@ class atmosphericModel:
 	def slalib(self, l=None):
 		if not l:
 			l = self.lrange
-		return(AstroAtmosphere.slalib(l, T=self.T, p=self.P, RH=self.H)-1)
+		try:
+			from pyslalib import slalib
+			return(AstroAtmosphere.slalib(l, T=self.T, p=self.P, RH=self.H)-1)
+		except ImportError as e:
+			print('Could not import pyslalib.slalib. Using slower Python native function sla_refro().')
+			h=0
+			lat=-24
+			tlr=0.0065
+			eps=1e-10
+			zObs = np.deg2rad(45.07182)
+			p = self.P/100
+			R = []
+			if type(l) == float:
+				l = [l]
+			for wl in l:
+				R.append(AstroAtmosphere.sla_refro(zObs, h, self.T, p, self.H, wl, lat, tlr, eps))
+			if len(l) == 1:
+				return R[0]
+			else:
+				return np.asarray(R)
 
 
 if __name__ == '__main__':