Added uncertainty calculation independent of refractivity model and small bugfixes

AstroAtmosphere/dispersionModels.py

@@ -387,6 +387,73 @@ class dispersion(Observatory):
 		return np.degrees(_dispersionError)
 
 
+	def independentCassiniError(self, n1, n2, zenith, dn1=0, dn2=0, dz=0, 
+									lat=None, h=None):
+		'''
+			Calculates the uncertainty in the atmospheric dispersion, 
+			following the Cassini model, without assumptions about the
+			refractivity model.
+			This involves a much shorter use the chain rule. 
+
+			Not taken into account is the errror propagation in the scale height
+			of the atmosphere, from my results this should give a small or even 
+			negligible difference.
+
+			Parameters
+			----------
+			n1 	: 	float
+				Refractive index at the first wavelength of interest at the 
+				observer
+			n2 	: 	float
+				Refractive index at the second wavelength of interest at the 
+				observer
+			zenith : float
+				Angle of observation in degrees.
+				zenith = 0 means the observation is directly overhead.
+			dn1 : float
+				Uncertainty in the provided value of the refractive index
+				at the first wavelength
+			dn2 : float
+				Uncertainty in the provided value of the refractive index
+				at the second wavelength
+			dz 	: float (default: 0)
+				Uncertainty in the zenith angle, in degrees.
+			lat : float (default: None)
+				If entered, overrides the latitude attribute of the class.
+			h 	: float (default=None)
+				If entered, overrides the altitude attribute of the class.
+
+			Returns
+			-------
+			dRR	: float
+				Uncertainty in the atmospheric dispersion in degrees.
+		'''
+		if lat == None:
+			lat = self.lat
+		else:
+			self.lat = lat
+		if h == None:
+			h 	= self.h
+		else:
+			self.h = h
+
+		# Some precalculation and variable setting
+		_zenith = np.deg2rad(zenith)
+		_dz   	= np.deg2rad(dz)
+		_R 		= self.rc
+		_H 		= self.H
+
+		# Partial derivatives for derivation to n1, n2, and z
+		_dddn1  = - _R * np.sin(_zenith) / ( (_R + _H) * np.sqrt(1 - (n1 * _R * np.sin(_zenith) / (_R + _H))**2 ) )
+		_dddn2  =   _R * np.sin(_zenith) / ( (_R + _H) * np.sqrt(1 - (n2 * _R * np.sin(_zenith) / (_R + _H))**2 ) )
+		_dddz 	= (_R * np.cos(_zenith)) / (_R + _H) * ( n1 / np.sqrt(1-(n1**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) - n2 / np.sqrt(1-(n2**2 * _R**2 * np.sin(_zenith)**2)/(_R + _H)**2) )		
+		
+		# Combining the partial derivatives
+		_dispersionError = np.sqrt( (_dddn1 * dn1)**2 + (_dddn2 * dn2)**2 + (_dddz * dz)**2)
+
+		return np.degrees(_dispersionError)
+
+
 	def corbard(self, n1, n2, zenith, reduced_height=None):
 		'''
 			Corbard et al., 2019, mention an additional formula based on

AstroAtmosphere/misc.py

@@ -336,7 +336,7 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 	_phi  	= np.deg2rad(lat) 		# From degrees to radians
 
 	R = []
-	if type(l) == float or type(l) == int:
+	if isinstance(l, float) or isinstance(l, int): 
 		l = [l]
 	if use_pyslalib:
 		try:
@@ -345,13 +345,13 @@ def slalib_refraction(l, zenith, conditions='STANDARD', T=288.15, p=101325,
 		    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)
+			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)
+			Refrac = sla_refro(_zenith, h, T, _P, RH, wl, _phi, tlr, eps)
 			R.append(np.degrees(Refrac))
 
 	if len(l) == 1: