included fail-safe: in case analytical equations do not yield plausible result, output from AstroPy sunangles will be used
This commit is contained in:
parent
449f92da0f
commit
89e9ec1db1
|
|
@ -37,48 +37,53 @@ def sun_angles_astropy(lat, lon, h, utc): # get current sun elevation and azimu
|
|||
return az, elv
|
||||
|
||||
|
||||
def sun_angles_analytical(lat, lon, utc): # get current sun elevation and azimuth through several equations (see [xx])
|
||||
if np.abs(lat) == 90: # handling collapse of longitudes at poles by
|
||||
lat = np.sign(lat) * 89.999999 # expanding one point to a very small circle
|
||||
else:
|
||||
pass
|
||||
def sun_angles_analytical(lat, lon, h, utc): # get current sun elevation and azimuth through several equations (see [xx])
|
||||
try:
|
||||
if np.abs(lat) == 90: # handling collapse of longitudes at poles by
|
||||
lat = np.sign(lat) * 89.999999 # expanding one point to a very small circle
|
||||
else:
|
||||
pass
|
||||
|
||||
jd = utc.jd
|
||||
jc = (jd - 2451545) / 36525
|
||||
gml = (280.46646 + jc * (36000.76983 + jc * 0.0003032)) % 360
|
||||
gma = 357.52911 + jc * (35999.05029 - 0.0001537 * jc)
|
||||
eeo = 0.016708634 - jc * (0.000042037 + 0.0000001267 * jc)
|
||||
sec = np.sin(np.deg2rad(gma)) * (1.914602 - jc * (0.004817 + 0.000014 * jc)) + np.sin(np.deg2rad(2 * gma)) * (
|
||||
jd = utc.jd
|
||||
jc = (jd - 2451545) / 36525
|
||||
gml = (280.46646 + jc * (36000.76983 + jc * 0.0003032)) % 360
|
||||
gma = 357.52911 + jc * (35999.05029 - 0.0001537 * jc)
|
||||
eeo = 0.016708634 - jc * (0.000042037 + 0.0000001267 * jc)
|
||||
sec = np.sin(np.deg2rad(gma)) * (1.914602 - jc * (0.004817 + 0.000014 * jc)) + np.sin(np.deg2rad(2 * gma)) * (
|
||||
0.019993 - 0.000101 * jc) + np.sin(np.deg2rad(3 * gma)) * 0.000289
|
||||
stl = gml + sec
|
||||
sal = stl - 0.00569 - 0.00478 * np.sin(np.deg2rad(125.04 - 1934.136 * jc))
|
||||
moe = 23 + (26 + (21.448 - jc * (46.815 + jc * (0.00059 - jc * 0.001813))) / 60) / 60
|
||||
oc = moe + 0.00256 * np.cos(np.deg2rad(125.04 - 1934.136 * jc))
|
||||
sd = np.rad2deg(np.arcsin(np.sin(np.deg2rad(oc)) * np.sin(np.deg2rad(sal)))) # radian
|
||||
var_y = np.tan(np.deg2rad(oc / 2)) ** 2
|
||||
eot = 4 * np.rad2deg(
|
||||
var_y * np.sin(2 * np.deg2rad(gml)) - 2 * eeo * np.sin(np.deg2rad(gma)) + 4 * eeo * var_y * np.sin(
|
||||
np.deg2rad(gma)) * np.cos(2 * np.deg2rad(gml)) - 0.5 * var_y ** 2 * np.sin(
|
||||
4 * np.deg2rad(gml)) - 1.25 * eeo ** 2 * np.sin(2 * np.deg2rad(gma)))
|
||||
tst = (((jd - 0.5) % 1) * 1440 + eot + 4 * lon) % 1440
|
||||
stl = gml + sec
|
||||
sal = stl - 0.00569 - 0.00478 * np.sin(np.deg2rad(125.04 - 1934.136 * jc))
|
||||
moe = 23 + (26 + (21.448 - jc * (46.815 + jc * (0.00059 - jc * 0.001813))) / 60) / 60
|
||||
oc = moe + 0.00256 * np.cos(np.deg2rad(125.04 - 1934.136 * jc))
|
||||
sd = np.rad2deg(np.arcsin(np.sin(np.deg2rad(oc)) * np.sin(np.deg2rad(sal)))) # radian
|
||||
var_y = np.tan(np.deg2rad(oc / 2)) ** 2
|
||||
eot = 4 * np.rad2deg(
|
||||
var_y * np.sin(2 * np.deg2rad(gml)) - 2 * eeo * np.sin(np.deg2rad(gma)) + 4 * eeo * var_y * np.sin(
|
||||
np.deg2rad(gma)) * np.cos(2 * np.deg2rad(gml)) - 0.5 * var_y ** 2 * np.sin(
|
||||
4 * np.deg2rad(gml)) - 1.25 * eeo ** 2 * np.sin(2 * np.deg2rad(gma)))
|
||||
tst = (((jd - 0.5) % 1) * 1440 + eot + 4 * lon) % 1440
|
||||
|
||||
if tst / 4 < 0:
|
||||
ha = tst / 4 + 180
|
||||
else:
|
||||
ha = tst / 4 - 180
|
||||
if tst / 4 < 0:
|
||||
ha = tst / 4 + 180
|
||||
else:
|
||||
ha = tst / 4 - 180
|
||||
|
||||
sza = np.rad2deg(np.arccos(
|
||||
np.sin(np.deg2rad(lat)) * np.sin(np.deg2rad(sd)) + np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(sd)) * np.cos(
|
||||
np.deg2rad(ha))))
|
||||
sea = 90 - sza
|
||||
sza = np.rad2deg(np.arccos(
|
||||
np.sin(np.deg2rad(lat)) * np.sin(np.deg2rad(sd)) + np.cos(np.deg2rad(lat)) * np.cos(
|
||||
np.deg2rad(sd)) * np.cos(
|
||||
np.deg2rad(ha))))
|
||||
sea = 90 - sza
|
||||
|
||||
if ha > 0:
|
||||
saa = (np.rad2deg(np.arccos(((np.sin(np.deg2rad(lat)) * np.cos(np.deg2rad(sza))) - np.sin(np.deg2rad(sd))) / (
|
||||
np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(sza))))) + 180) % 360
|
||||
else:
|
||||
saa = (540 - np.rad2deg(np.arccos(
|
||||
((np.sin(np.deg2rad(lat)) * np.cos(np.deg2rad(sza))) - np.sin(np.deg2rad(sd))) / (
|
||||
if ha > 0:
|
||||
saa = (np.rad2deg(
|
||||
np.arccos(((np.sin(np.deg2rad(lat)) * np.cos(np.deg2rad(sza))) - np.sin(np.deg2rad(sd))) / (
|
||||
np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(sza))))) + 180) % 360
|
||||
else:
|
||||
saa = (540 - np.rad2deg(np.arccos(
|
||||
((np.sin(np.deg2rad(lat)) * np.cos(np.deg2rad(sza))) - np.sin(np.deg2rad(sd))) / (
|
||||
np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(sza)))))) % 360
|
||||
except:
|
||||
saa, sea = sun_angles_astropy(lat, lon, h, utc)
|
||||
|
||||
return saa, sea # Azimuth, Elevation
|
||||
|
||||
|
|
@ -95,11 +100,11 @@ def AirMass(p_air, p_0, ELV, h): # get atmospheric air mass over balloon
|
|||
return res
|
||||
|
||||
|
||||
def tau(ELV, h, p_air): # get atmospheric transmissivity as function of balloon altitude and sun elevation
|
||||
def tau(ELV, h, p_air, p0): # get atmospheric transmissivity as function of balloon altitude and sun elevation
|
||||
if ELV >= -(180 / np.pi * np.arccos(R_E / (R_E + h))):
|
||||
tau_atm = 0.5 * (
|
||||
np.exp(-0.65 * AirMass(p_air, p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air, p_0, ELV, h)))
|
||||
tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air / p_0) + np.exp(-0.095 * p_air / p_0))
|
||||
np.exp(-0.65 * AirMass(p_air, p0, ELV, h)) + np.exp(-0.095 * AirMass(p_air, p0, ELV, h)))
|
||||
tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air / p0) + np.exp(-0.095 * p_air / p0))
|
||||
else:
|
||||
tau_atm = 0
|
||||
tau_atmIR = 0
|
||||
|
|
|
|||
Loading…
Reference in New Issue