From 3468a9159c15bc537df1ceda4dc471559cdd7550 Mon Sep 17 00:00:00 2001 From: Joost van den Born Date: Tue, 4 Jun 2019 14:33:46 +0200 Subject: [PATCH] Fixed the full integration methods --- atmosphericModel.py | 60 +++++++++++++++++++++++++++++---------------- 1 file changed, 39 insertions(+), 21 deletions(-) diff --git a/atmosphericModel.py b/atmosphericModel.py index fcb148b..2708dcb 100644 --- a/atmosphericModel.py +++ b/atmosphericModel.py @@ -1,7 +1,7 @@ #!/usr/bin/env Python3 import numpy as np from scipy import integrate - +from pdb import set_trace ''' @ Author: Joost van den Born @@ -755,6 +755,7 @@ class dispersion(Observatory): different atmosphere layers. Recommended use with the US Standard Atmosphere from 1976. USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html + The method described in eq. 3 in Auer & Standish (2000) is used for the calculation. ''' if not R0: R0 = self.rc @@ -770,19 +771,26 @@ class dispersion(Observatory): # Integration of the refraction integral over the atmosphere layers. for _n in [n1, n2]: + _Dens = np.asarray(_Dens) + _h = np.asarray(_h) + _R = R0 + _h * 1000 - _refractionIntegralList = [] - _nList = [] + # Gladstone-Dale relation between refractive index and density + _n_h = 1 + (_n - 1) / _Dens[0] * _Dens - for _i in range(len(_h)): - _R = R0 + _h[_i]*1000 - _n_h = 1 + (_n - 1) / _Dens[0] * _Dens[_i] # Gladstone-Dale relation between refractive index and density - _nList.append(_n_h) - _refractionIntegralList.append( (_n * R0 * np.sin(z0))/(_n_h*np.sqrt((_n_h*_R)**2 - (_n * R0 * np.sin(z0))**2 )) ) - _dn = (_nList[0] - _nList[-1]) / ((len(_h)-1)) - - _rList.append(np.sum(np.degrees(np.asarray(_refractionIntegralList) * _dn))) + # Calculation of the zenith angle for given height and refractive index, from the refractive invariant + _z = np.arcsin( (_n * R0) / (_n_h * _R) * np.sin(z0) ) + + # Calculation of log space derivative of n and r + _lnn = np.log(_n_h) + _lnr = np.log(_R) + _dlnn = np.asarray([_lnn[i+1] - _lnn[i] for i in range(len(_lnn)-1)]) + _dlnr = np.asarray([_lnr[i+1] - _lnr[i] for i in range(len(_lnr)-1)]) + _dz = np.asarray([(_z[i] - _z[i+1]) for i in range(len(_z)-1)]) + # Calculation of the integrand in eq. 3 in Auer & Standish + _Integrand = -1*(_dlnn/_dlnr) / (1 + _dlnn/_dlnr) * _dz + _rList.append(np.degrees(np.sum(_Integrand))) return _rList[0] - _rList[1] @@ -1182,19 +1190,29 @@ class refraction(Observatory): _Dens = rhoData _h = heightData - _refractionIntegralList = [] - _nList = [] z0 = np.deg2rad(z0) + _Dens = np.asarray(_Dens) + _h = np.asarray(_h) + _R = R0 + _h * 1000 - # Integration of the refraction integral over the atmosphere layers. - for _i in range(len(_h)): - _R = R0 + _h[_i]*1000 - _n_h = 1 + (n0 - 1) / _Dens[0] * _Dens[_i] # Gladstone-Dale relation between refractive index and density - _nList.append(_n_h) - _refractionIntegralList.append( (n0 * R0 * np.sin(z0))/(_n_h*np.sqrt((_n_h*_R)**2 - (n0 * R0 * np.sin(z0))**2 )) ) - _dn = (_nList[0] - _nList[-1]) / ((len(_h)-1)) + # Gladstone-Dale relation between refractive index and density + _n_h = 1 + (n0 - 1) / _Dens[0] * _Dens + + # Calculation of the zenith angle for given height and refractive index, from the refractive invariant + _z = np.arcsin( (n0 * R0) / (_n_h * _R) * np.sin(z0) ) - return np.sum(np.degrees(np.asarray(_refractionIntegralList) * _dn)) + # Calculation of log space derivative of n and r + _lnn = np.log(_n_h) + _lnr = np.log(_R) + _dlnn = np.asarray([_lnn[i+1] - _lnn[i] for i in range(len(_lnn)-1)]) + _dlnr = np.asarray([_lnr[i+1] - _lnr[i] for i in range(len(_lnr)-1)]) + _dz = np.asarray([(_z[i] - _z[i+1]) for i in range(len(_z)-1)]) + + # Calculation of the integrand in eq. 3 in Auer & Standish + _Integrand = -1*(_dlnn/_dlnr) / (1 + _dlnn/_dlnr) * _dz + + + return np.degrees(np.sum(_Integrand))