From 89e9ec1db1ccd78b2090cd95dae2f2a5463fc412 Mon Sep 17 00:00:00 2001 From: Marcel Frommelt Date: Mon, 3 May 2021 17:42:26 +0900 Subject: [PATCH] included fail-safe: in case analytical equations do not yield plausible result, output from AstroPy sunangles will be used --- models/sun.py | 83 +++++++++++++++++++++++++++------------------------ 1 file changed, 44 insertions(+), 39 deletions(-) diff --git a/models/sun.py b/models/sun.py index 6f55c65..4616a23 100644 --- a/models/sun.py +++ b/models/sun.py @@ -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