Package ivs :: Package observations :: Module airmass
[hide private]
[frames] | no frames]

Source Code for Module ivs.observations.airmass

  1  """ 
  2  Calculate airmass according to various models. 
  3  """ 
  4   
  5  import numpy as np 
  6   
7 -def _deg2rad(degrees):
8 return degrees*np.pi/180.
9
10 -def _rad2deg(radians):
11 return radians*180./np.pi
12
13 -def _sec(radians):
14 return 1./np.cos(radians)
15
16 -def get_modelnames():
17 return(np.array(['YoungAndIrvine1967', 'Hardie1962', 'Rozenberg1966', 'KastenAndYoung1989', 'Young1994', 'Pickering2002']))
18
19 -def airmass(zz, model='Pickering2002'):
20 """ 21 Calculate the airmass from an interpolative model 22 23 Input: 24 @param z: array with zenith angles (90-altitude in degrees) 25 @type z: ndarray 26 @param model: name of the model 27 @type model: string 28 29 Output: 30 @return: airmass 31 @rtype: ndarray 32 33 Many formulas have been developed to fit tabular values of air mass: 34 + The first one by Young and Irvine (1967) included a simple corrective term: 35 36 X = sec(z)*(1.-0.0012*(sec(z)**2. - 1.)) 37 38 where zt is the true zenith angle. This gives usable results up to about 39 80 degrees, but the accuracy degrades rapidly at greater zenith angles. The 40 calculated air mass reaches a maximum of 11.13 at 86.6 degrees, becomes zero 41 at 88 degrees, and approaches negative infinity at the horizon. The plot of this 42 formula on the accompanying graph includes a correction for atmospheric 43 refraction so that the calculated air mass is for apparent rather than true 44 zenith angle. 45 46 + Hardie (1962) introduced a polynomial in: 47 48 X = sec(z) - 0.0018167*(sec(z) - 1.) - 0.002875*(sec(z) - 1.)**2. - 0.0008083*(sec(z)-1.)**3. 49 50 which gives usable results for zenith angles of up to perhaps 85 degrees. 51 As with the previous formula, the calculated air mass reaches a maximum, 52 and then approaches negative infinity at the horizon. 53 54 + Rozenberg (1966) suggested: 55 56 X = 1./(cos(z) + 0.025*exp(-11.*cos(z))) 57 58 which gives reasonable results for high zenith angles, with a horizon 59 airmass of 40. 60 61 + Kasten and Young (1989) developed: 62 63 X = 1./(cos(z) + 0.50572*(96.07995 - zd)**(-1.6364)) 64 65 which gives reasonable results for zenith angles of up to 90 degrees, with 66 an airmass of approximately 38 at the horizon. Here the zd term is z in 67 degrees. 68 69 + Young (1994) developed: 70 71 X = (1.002432*cos(z)**2. + 0.148386*cos(z) + 0.0096467)/(cos(z)**3. + 0.149864*cos(z)**2. + 0.0102963*cos(z) + 0.000303978) 72 73 in terms of the true zenith angle zt, for which he claimed a maximum error 74 (at the horizon) of 0.0037 air mass. 75 76 + Finally, Pickering (2002) developed 77 78 X = 1./sin(h+244/(165+47*h**1.1)) 79 80 where h is apparent altitude (90.-z) in degrees. Pickering claimed his 81 equation to have a tenth the error of Schaefer (1998) near the horizon. 82 83 Respectively, the keyword should be: 84 'YoungAndIrvine1967', 'Hardie1962', 'Rozenberg1966', 'KastenAndYoung1989', 85 'Young1994', 'Pickering2002' 86 87 Example usage: 88 get the different airmasses for all models between 70 and 90 degrees 89 >>> zz = np.arange(70.000, 89.999, .001) 90 >>> models = get_modelnames() 91 >>> airmasses = np.zeros((len(zz), len(models))) 92 >>> for i, model in enumerate(models): airmasses[:,i] = airmass(zz, model) 93 >>> print(airmasses[-10,:]) 94 [ -1.69573409e+08 -1.14232330e+08 3.97784409e+01 3.77562570e+01 95 3.16224873e+01 3.85395375e+01] 96 """ 97 # define the standard quantities 98 zrad = _deg2rad(zz) 99 if model.lower() == 'youngandirvine1967': 100 out = _sec(zrad)*(1.-0.0012*(_sec(zrad)**2. - 1.)) 101 elif model.lower() == 'hardie1962': 102 out = _sec(zrad) - 0.0018167*(_sec(zrad) - 1.) - 0.002875*(_sec(zrad) - 1.)**2. - 0.0008083*(_sec(zrad) - 1.)**3. 103 elif model.lower() == 'rozenberg1966': 104 out = 1./(np.cos(zrad) + 0.025*np.exp(-11.*np.cos(zrad))) 105 elif model.lower() == 'kastenandyoung1989': 106 out = 1./(np.cos(zrad) + 0.50572*(96.07995 - zz)**(-1.6364)) 107 elif model.lower() == 'young1994': 108 out = (1.002432*np.cos(zrad)**2. + 0.148386*np.cos(zrad) + 0.0096467)/(np.cos(zrad)**3. + 0.149864*np.cos(zrad)**2. + 0.0102963*np.cos(zrad) + 0.000303978) 109 elif model.lower() == 'pickering2002': 110 hh = 90.-zz 111 out = 1./np.sin(_deg2rad(hh+244/(165+47*hh**1.1))) 112 return (out)
113