From ee0e7bccac5d3582fc79e46bc074a6ba2bac1f73 Mon Sep 17 00:00:00 2001
From: Joost van den Born <born@astron.nl>
Date: Thu, 21 Nov 2019 11:28:51 +0100
Subject: [PATCH] Bugfixes and removed old code

---
 atmospheric refraction example.html  |  46 +----
 atmospheric refraction example.ipynb |   7 +-
 atmosphericModel.py                  | 273 +--------------------------
 3 files changed, 11 insertions(+), 315 deletions(-)

diff --git a/atmospheric refraction example.html b/atmospheric refraction example.html
index 20491b1..65634ed 100644
--- a/atmospheric refraction example.html	
+++ b/atmospheric refraction example.html	
@@ -13017,45 +13017,6 @@ ul.typeahead-list  > li > a.pull-right {
 .highlight .vm { color: #19177C } /* Name.Variable.Magic */
 .highlight .il { color: #666666 } /* Literal.Number.Integer.Long */
     </style>
-<style type="text/css">
-    
-/* Temporary definitions which will become obsolete with Notebook release 5.0 */
-.ansi-black-fg { color: #3E424D; }
-.ansi-black-bg { background-color: #3E424D; }
-.ansi-black-intense-fg { color: #282C36; }
-.ansi-black-intense-bg { background-color: #282C36; }
-.ansi-red-fg { color: #E75C58; }
-.ansi-red-bg { background-color: #E75C58; }
-.ansi-red-intense-fg { color: #B22B31; }
-.ansi-red-intense-bg { background-color: #B22B31; }
-.ansi-green-fg { color: #00A250; }
-.ansi-green-bg { background-color: #00A250; }
-.ansi-green-intense-fg { color: #007427; }
-.ansi-green-intense-bg { background-color: #007427; }
-.ansi-yellow-fg { color: #DDB62B; }
-.ansi-yellow-bg { background-color: #DDB62B; }
-.ansi-yellow-intense-fg { color: #B27D12; }
-.ansi-yellow-intense-bg { background-color: #B27D12; }
-.ansi-blue-fg { color: #208FFB; }
-.ansi-blue-bg { background-color: #208FFB; }
-.ansi-blue-intense-fg { color: #0065CA; }
-.ansi-blue-intense-bg { background-color: #0065CA; }
-.ansi-magenta-fg { color: #D160C4; }
-.ansi-magenta-bg { background-color: #D160C4; }
-.ansi-magenta-intense-fg { color: #A03196; }
-.ansi-magenta-intense-bg { background-color: #A03196; }
-.ansi-cyan-fg { color: #60C6C8; }
-.ansi-cyan-bg { background-color: #60C6C8; }
-.ansi-cyan-intense-fg { color: #258F8F; }
-.ansi-cyan-intense-bg { background-color: #258F8F; }
-.ansi-white-fg { color: #C5C1B4; }
-.ansi-white-bg { background-color: #C5C1B4; }
-.ansi-white-intense-fg { color: #A1A6B2; }
-.ansi-white-intense-bg { background-color: #A1A6B2; }
-
-.ansi-bold { font-weight: bold; }
-
-    </style>
 
 
 <style type="text/css">
@@ -13194,8 +13155,7 @@ div#notebook {
 <div class="prompt input_prompt">In&nbsp;[3]:</div>
 <div class="inner_cell">
     <div class="input_area">
-<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># Calculating refraction for edge FoV</span>
-<span class="n">zenith_obs</span> <span class="o">=</span> <span class="mi">45</span>
+<div class=" highlight hl-ipython3"><pre><span></span><span class="n">zenith_obs</span> <span class="o">=</span> <span class="mi">45</span>
 <span class="n">refraction_obs</span> <span class="o">=</span> <span class="n">ref</span><span class="o">.</span><span class="n">cassini</span><span class="p">(</span><span class="n">n1</span><span class="p">,</span> <span class="n">zenith_obs</span><span class="p">)</span> <span class="o">*</span> <span class="mf">3.6e3</span>
 <span class="nb">print</span> <span class="p">(</span><span class="s1">&#39;The refraction is </span><span class="si">%.03f</span><span class="s1"> arc seconds&#39;</span> <span class="o">%</span><span class="p">(</span><span class="n">refraction_obs</span><span class="p">))</span>
 </pre></div>
@@ -13242,7 +13202,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[17]:</div>
+<div class="prompt input_prompt">In&nbsp;[4]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># Calculating the uncertainty in the refractive index</span>
@@ -13288,7 +13248,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[23]:</div>
+<div class="prompt input_prompt">In&nbsp;[5]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># Initializing dispersion object</span>
diff --git a/atmospheric refraction example.ipynb b/atmospheric refraction example.ipynb
index 59a1783..4d61886 100644
--- a/atmospheric refraction example.ipynb	
+++ b/atmospheric refraction example.ipynb	
@@ -81,7 +81,6 @@
     }
    ],
    "source": [
-    "# Calculating refraction for edge FoV\n",
     "zenith_obs = 45\n",
     "refraction_obs = ref.cassini(n1, zenith_obs) * 3.6e3\n",
     "print ('The refraction is %.03f arc seconds' %(refraction_obs))"
@@ -104,7 +103,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
@@ -136,7 +135,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
     {
@@ -178,7 +177,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.1"
+   "version": "3.7.4"
   }
  },
  "nbformat": 4,
diff --git a/atmosphericModel.py b/atmosphericModel.py
index ff3301b..dc651e5 100644
--- a/atmosphericModel.py
+++ b/atmosphericModel.py
@@ -178,7 +178,6 @@ class Observatory:
 		else:
 			# Partial pressure of water vapour is p_water = h*svp
 			_rho_w 	 = self.rho(p=RH*self.svp(T), T=T, RH=RH, xc=xc)			# Moist air component for the actual conditions.
-		# _rho_ws	 = self.rho(p=1333, T=293.15, RH=1/1.7564043320974863, xc=xc) 	# Moist air component at standard conditions.
 		_rho_ws	 = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc) 	# Moist air component at standard conditions.
 		
 		_Ltotal  = ( _rho_a / _rho_axs ) * _La + ( _rho_w / _rho_ws ) * _Lw
@@ -384,276 +383,14 @@ class Observatory:
 		# Final error:
 		_nError = np.sqrt( (_dndl * dl)**2 + (_dndT * dT)**2  + (_dndP * dP)**2 + (_dndRH * dRH)**2 + (_dndCO2 * dCO2)**2 )
 
-		return _nError
-
-###
-### THIS ERROR CLASS IS OUTDATED
-### Use Observatory.dn_tph(), dispersion.cassiniError() or refraction.cassiniError() instead.
-###
-class error(Observatory):
-	def __init__(self, dT=0.2, dP=20, dH=0.02, dCO2=20, dl=0.001):
-		Observatory.__init__(self)
-		self.dT 	= dT 	# in K
-		self.dP 	= dP 	# in Pa
-		self.dH 	= dH	# relative humidity between 0 and 1
-		self.dCO2 	= dCO2	# in ppm
-		self.dl 	= dl 	# in micron
-		print('THIS ERROR CLASS IS OUTDATED\nUse Observatory.dn_tph(), dispersion.cassiniError() or refraction.cassiniError() instead.')
-
-	def dn_as(self, l, dl=None):
-		'''
-			Formula that calculate the error on the refractive index for standard 
-			atmospheric conditions, i.e. 15 C, 101,325 Pa, 0% RH, 450 ppm CO2.
-			It returns dn_as, not d(n_as - 1) * 1e8!
-		'''
-		if not dl:
-			dl = self.dl
-
-		_k0 = 238.0185 		# um^-2
-		_k1 = 5792105  		# um^-2
-		_k2 = 57.362		# um^-2
-		_k3 = 167917		# um^-2	
-		_dn_as = np.sqrt((( - (2 * _k1 * l ) / ((_k0 * l**2 - 1)**2) - (2 * _k3 * l) / ((_k2 * l**2 - 1)**2) )*dl)**2)
-		return _dn_as / 1e8
-
-
-	def dn_axs(self, l, xc, dCO2=None, dl=None):
-		'''
-			Formula that calculates the error on the refractive index for standard 
-			atmospheric	conditions, i.e. 15 C, 101,325 Pa, 0% RH, but with xc ppm of CO2.
-			It returns dn_axs, not d(n_axs - 1)
-		'''
-		if not dCO2:
-			dCO2 = self.dCO2
-		if not dl:
-			dl = self.dl
-		_dn_as = self.dn_as(l, dl)
-		_n_as  = self.n_as(l)
-		_dn_axs = np.sqrt( ( (_n_as - 1) * 0.534e-6 * dCO2 )**2 + ( (0.534e-6 * (xc - 450) + 1) * _dn_as )**2 )
-		return _dn_axs
-
-
-	def dn_ws(self, l, dl=None):
-		'''
-			Formula that calculates the refractivity of air for standard moist conditions
-			as defined by B&S (1939), i.e. 20 C, 1333 Pa
-			Returns n_ws, not (n_ws - 1) * 1e8
-		'''
-		if not dl:
-			dl = self.dl
-		_cf = 1.022
-		_w0 = 295.235 	# um^-2
-		_w1 = 2.6422	# um^-2
-		_w2 = -0.032380	# um^-4
-		_w3 = 0.004028	# um^-6
-		
-		_dndl = _cf * (-1 * _w1 * 2/l**3 - _w2 * 4/l**5 - _w3 * 6/l**7)
-		return np.sqrt( (_dndl * dl)**2 ) / 1e8
-
-
-	def df(self, p, T, dP=None, dT=None):
-		'''
-			Uncertainty in enhancement factor of
-			water vapour in air.
-			t in degrees celsius
-		'''
-		if not dP:
-			dP = self.dP
-		if not dT:
-			dT = self.dT
-
-		_t = T - 273.15
-		_b = 3.14e-8 	# Pa^-1
-		_c = 5.6e-7		# deg C^-2
-		return np.sqrt( (_b*dP)**2 + (2*_c*_t*dT)**2 )
-
-
-	def dsvp(self, T, dT=None):
-		'''
-			Uncertainty in saturation vapour pressure in Pa.
-		'''
-		if not dT:
-			dT = self.dT
-		_A =  1.2378847e-5 # K^-2
-		_B = -1.9121316e-2 # K^-2
-		_C = 33.93711047
-		_D = -6.3431645e3  # K^-2
-		_p1 = 2*_A*T + _B - _D/(T**2)
-		_p2 = np.exp(_A*T**2 + _B*T + _C + _D/T)
-		return np.sqrt( (_p1 * _p2 * dT)**2 )
-
-
-	def dxw(self, p, T, RH, dP=None, dT=None, dRH=None):
-		'''
-			Uncertainty in molar fraction of water vapour in moist air.
-		'''
-		if not dP:
-			dP = self.dP
-		if not dT:
-			dT = self.dT
-		if not dRH:
-			dRH = self.dH
-
-		_f 		= self.f(p, T)
-		_svp 	= self.svp(T)
-
-		_df 	= self.df(p=p, T=T, dP=dP, dT=dT)
-		_dsvp	= self.dsvp(T=T, dT=dT)
-		_p1		= (RH * _svp / p) * _df
-		_p2		= (_f * _svp / p) * dRH
-		_p3		= (_f * RH  / p) * _dsvp
-		_p4		= (-1 * _f * RH * _svp / p**2) * dP
-		return np.sqrt( _p1**2 + _p2**2 + _p3**2 + _p4**2)
-
-
-	def dZ(self, p, T, RH, xw='', dP=None, dT=None, dRH=None):
-		'''
-			Uncertainty in compressibility of moist air
-		'''
-		if not dP:
-			dP = self.dP
-		if not dT:
-			dT = self.dT
-		if not dRH:
-			dRH = self.dH
-		_dt  = dT
-		_dxw = self.dxw(p=p, T=T, RH=RH, dP=dP, dT=dT, dRH=dRH)
-
-		_a0 =  1.58123e-6 # K Pa^-1
-		_a1 = -2.9331e-8  # Pa^-1
-		_a2 =  1.1043e-10 # K^-1 Pa^-1
-		_b0 =  5.707e-6   # K Pa^-1
-		_b1 = -2.051e-8   # Pa^-1
-		_c0 =  1.9898e-4  # K Pa^-1
-		_c1 = -2.376e-6   # Pa^-1
-		_d  =  1.83e-11   # K^2 Pa^-2
-		_e  = -0.765e-8   # K^2 Pa^-2
-		_t  = T - 273.15
-		if xw == '':
-			_xw = self.xw(p, T, RH)
-		else:
-			_xw = xw
-
-		_dZdp  = -1/T * ( _a0 + _a1*_t + _a2*_t**2 + (_b0 + _b1*_t)*_xw + (_c0 + _c1*_t)*_xw**2) + (2*p/T**2) * (_d + _e*_xw**2)
-		_dZdT  = p/T**2 * ( _a0 + _a1*_t + _a2*_t**2 + (_b0 + _b1*_t)*_xw + (_c0 + _c1*_t)*_xw**2) - 2*p**2/T**3 * (_d + _e*_xw**2)
-		_dZdt  = -1*p/T * (_a1 + 2*_a2*_t + _b1*_xw + _c1*_xw**2)
-		_dZdxw = -1*p/T * ( (_b0 + _b1*_t) + 2*(_c0 + _c1*_t)*_xw ) + 2 * (p/T)**2 * _e * _xw
-		_dZ    = np.sqrt( (_dZdp * dP)**2 + (_dZdT * dT)**2 + (_dZdt * _dt)**2 + (_dZdxw * _dxw)**2 )
-		return _dZ
-
-
-	def drho(self, p, T, RH, xc, xw='', dP=None, dT=None, dRH=None, dCO2=None):
-		'''
-			Calculates the uncertainty in the density of air for the given parameters.
-			T in Kelvin
-			p in Pascal
-			RH between 0 and 1
-			xc in ppm
-		'''
-		if not dP:
-			dP = self.dP
-		if not dT:
-			dT = self.dT
-		if not dRH:
-			dRH = self.dH
-		if not dCO2:
-			dCO2 = self.dCO2
-
-		_Mw  = 0.018015 										# kg/mol
-		_R   = 8.314510										# J mol^-1 K^-1, gas constant
-		_Ma  = 1e-3 * (28.9635 + 12.011e-6 * (xc - 400) )	# Molar mass in kg/mol of dry air containing xc ppm of CO2
-		if xw == '':
-			_xw = self.xw(p, T, RH)
-		else: 
-			_xw = xw
-		_Z   = self.Z(p=p, T=T, RH=RH, xw=_xw)
-
-		_dMa = np.sqrt( (12.011e-9*dCO2)**2 )
-		_dZ  = self.dZ(p=p, T=T, RH=RH, xw=_xw, dP=dP, dT=dT, dRH=dRH)
-		_dxw = self.dxw(p=p, T=T, RH=RH, dP=dP, dT=dT, dRH=dRH)
-
-		_drhodp  = _Ma / (_Z * _R * T) * (1 - _xw * (1 - _Mw / _Ma))
-		_drhodMa =   p / (_Z * _R * T) * (1 - _xw)
-		_drhodZ  = -1 * (p * _Ma) / (_Z**2 * _R * T) * (1 - _xw * (1 - _Mw / _Ma))
-		_drhodT  = -1 * (p * _Ma) / (_Z * _R * T**2) * (1 - _xw * (1 - _Mw / _Ma))
-		_drhodxw = p / (_Z * _R * T) * (_Mw - _Ma)
-
-		return np.sqrt( (_drhodp * dP)**2 + (_drhodMa * _dMa)**2 + (_drhodZ * _dZ)**2 + (_drhodT * dT)**2 + (_drhodxw * _dxw)**2 )
-
-
-	def dn_tph(self, l, T, p, RH, xc, dl=None, dT=None, dP=None, dRH=None, dCO2=None):
-		'''
-			[DEPRECATED]
-			Calculates the uncertainty in the refraction of the atmosphere for any conditions.
-		'''
-		if not dP:
-			dP = self.dP
-		if not dT:
-			dT = self.dT
-		if not dRH:
-			dRH = self.dH
-		if not dCO2:
-			dCO2 = self.dCO2
-
-		_n_a 	 = self.n_axs(l, xc)
-		_n_w 	 = self.n_ws(l)
-		_dn_a    = self.dn_axs(l=l, xc=xc, dl=dl, dCO2=dCO2)
-		_dn_w    = self.dn_ws(l=l, dl=dl)
-
-		_La 	 = (_n_a**2 - 1) / (_n_a**2 + 2)
-		_Lw 	 = (_n_w**2 - 1) / (_n_w**2 + 2)
-		_dLa	 = np.sqrt( ( (6*_n_a) / (_n_a**2+2)**2 * _dn_a)**2 )
-		_dLw	 = np.sqrt( ( (6*_n_w) / (_n_w**2+2)**2 * _dn_w)**2 )
-		
-		# Calculate the density of the dry air components, using the partial pressure of the dry air, i.e. p_dry = p_atm - p_water
-		_pw 	 = RH*self.svp(T)
-		_pa 	 = p - _pw
-		_rho_a 	 = self.rho(p=_pa, T=T, RH=0, xc=xc)				# Dry air component for the actual conditions.
-		_rho_axs = self.rho(p=101325, T=288.15, RH=0, xc=xc) 		# Dry air component at standard conditions.
-
-		_dpa 	 = np.sqrt( (dP)**2 + (-1*self.svp(T) * dRH)**2 + (-1 * RH * self.dsvp(T=T, dT=dT) )**2 )
-		_dpw 	 = np.sqrt( (RH * self.dsvp(T=T, dT=dT))**2 + (self.svp(T) * dRH)**2 )
-		_drho_a  = self.drho(p=_pa, T=T, RH=RH, xc=xc, dP=_dpa, dT=dT, dRH=dRH, dCO2=dCO2)
-		_drho_axs= 0
-		
-		if RH == 0:
-			# nan is returned for _rho_w if we put in a (partial) pressure of 0
-			_rho_w = 0
-			_drho_w = 0
-		else:
-			# Partial pressure of water vapour is p_water = h*svp
-			_rho_w 	 = self.rho(p=_pw, T=T, RH=RH, xc=xc)			# Moist air component for the actual conditions.
-			_drho_w  = self.drho(p=_pw, T=T, RH=RH, xc=xc, dP=_dpw, dT=dT, dRH=dRH, dCO2=dCO2)
-		_rho_ws	 = self.rho(p=1333, T=293.15, RH=1333/self.svp(293.15), xc=xc) 	# Moist air component at standard conditions.
-		_drho_ws = 0
-
-		
-		_dLdL_a 	= _rho_a / _rho_axs 
-		_dLdL_w 	= _rho_w / _rho_ws
-		_dLdrho_a 	= _La / _rho_axs
-		_dLdrho_axs = - (_rho_a * _La) / (_rho_axs**2)
-		_dLdrho_w 	= _Lw / _rho_ws
-		_dLdrho_ws 	= - (_rho_w * _Lw) / _rho_ws**2
-
-		_Ltotal  = ( _rho_a / _rho_axs ) * _La + ( _rho_w / _rho_ws ) * _Lw
-		_dLtotal = np.sqrt( (_dLdL_a * _dLa)**2 + (_dLdL_w * _dLw)**2 + (_dLdrho_a * _drho_a)**2 
-			+ (_dLdrho_axs * _drho_axs)**2 + (_dLdrho_w * _drho_w)**2 + (_dLdrho_ws * _drho_ws)**2 )
-
-		_nLL = np.sqrt( (1 + 2 * _Ltotal) / (1 - _Ltotal) )
-		_dnLL = np.sqrt( ( (3*np.sqrt(1 - _Ltotal )) / ( 2 * (1 - _Ltotal)**2 * np.sqrt(1 + 2*_Ltotal) ) * _dLtotal )**2 )
-		
-		return (_nLL, _dnLL)
-###
-### END OF OUTDATED ERROR CLASS
-###
-	
+		return _nError	
 
 
 class dispersion(Observatory):
 	def __init__(self, lat, h):
 		Observatory.__init__(self)
-		self.h  	= h 					# height above sea level in meters
-		self.lat 	= lat					# Latitude of the observer in degrees
+		self.h  	= h 						# height above sea level in meters
+		self.lat 	= lat						# Latitude of the observer in degrees
 		self.rc 	= self.set_rc() * 1000		# Radius of curvature of the earth at the observer
 
 	def setReducedHeight(self, P, rho):
@@ -759,7 +496,7 @@ class dispersion(Observatory):
 		if not R0:
 			R0 = self.rc
 		if useStandardAtmosphere:
-			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('./US1976 Standard Atmosphere/si1py_modified.txt', unpack=True, skip_header=6)
+			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.txt', unpack=True, skip_header=6)
 		else:
 			_Dens = rhoData
 			_h 	  = heightData 
@@ -1184,7 +921,7 @@ class refraction(Observatory):
 		if not R0:
 			R0 = self.rc
 		if useStandardAtmosphere:
-			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('./US1976 Standard Atmosphere/si1py_modified.txt', unpack=True, skip_header=6)
+			_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('1976USSA.txt', unpack=True, skip_header=6)
 		else:
 			_Dens = rhoData
 			_h 	  = heightData