From 1894ced5f63ddbe7364c1b94499de49a6f74090e Mon Sep 17 00:00:00 2001 From: Marcel Frommelt Date: Fri, 9 Apr 2021 20:36:19 +0900 Subject: [PATCH] minor edits --- Super-Pressure.py | 6 +- models/drag.py | 23 ++++++- models/simple_atmosphere.py | 26 ++++--- models/sun.py | 134 +++++++++++++++++++++++++----------- models/test1.py | 5 -- models/test2.py | 3 - models/thermal.py | 106 ---------------------------- readme.txt | 6 +- 8 files changed, 136 insertions(+), 173 deletions(-) delete mode 100644 models/test1.py delete mode 100644 models/test2.py diff --git a/Super-Pressure.py b/Super-Pressure.py index 198663e..8b542ee 100644 --- a/Super-Pressure.py +++ b/Super-Pressure.py @@ -8,7 +8,7 @@ s0 = 0.00 # Initial Balloon Altitude in [m] p0 = 101325 # (Initial) Air Pressure in [Pa] rho_a = 1.2250 # (Initial) Air Density in [kg/m^3] -tmax = 50000 +tmax = 10000 T_a = 288.15 # (Initial) Air Temperature in [K] g = 9.80665 # local gravitation in [m/s^2] @@ -101,8 +101,10 @@ while t < tmax: dv_z = (I - np.sign(v_z) * D - G) / m_B * dt + # print(z, v_z) + v_z += dv_z * dt - z += v_z * dt + z += v_z * dt + 0.5 * dv_z * dt t += dt diff --git a/models/drag.py b/models/drag.py index 53a4e04..38a4b2c 100644 --- a/models/drag.py +++ b/models/drag.py @@ -1,6 +1,23 @@ import numpy as np +from input.user_input import * -c_d = 0.47 # drag coefficient balloon (spherical) [-] -def drag(c_d, rho_air, d_b, v): - return 0.125 * np.pi * c_d * rho_air * (d_b * v) ** 2 \ No newline at end of file +def Cd_model(fr, re, a_top): + k_cd = 19.8 + k1 = 1.50 * 10 ** 7 + k2 = 0.86 * 10 ** (-7) + a_top0 = np.pi / 4 * 1.383 * V_design ** (2/3) # a_top0 for zero-pressure shape + + value = 0.2 * k_cd/fr * (k1/re + k2 * re) * a_top / (a_top0 ** (3/2)) + + if value > 1.6: + value = 1.6 + + return value + + +c_d = 0.8 # 0.47 + + +def drag(c_d, rho_air, A_proj, v): + return 0.5 * c_d * rho_air * A_proj * v ** 2 diff --git a/models/simple_atmosphere.py b/models/simple_atmosphere.py index 39cf97f..672ee76 100644 --- a/models/simple_atmosphere.py +++ b/models/simple_atmosphere.py @@ -3,26 +3,30 @@ import numpy as np # ATMOSPHERE MODEL: + def T_air(h): - if h >= 0 and h <= 11000: + if 0 <= h <= 11000: res = 288.15 - 0.0065 * h - return res - elif h > 11000 and h <= 20000: + elif 11000 < h <= 20000: res = 216.65 - return res elif h >= 20000: res = 216.65 + 0.0010 * (h - 20000) - return res + return res + + def p_air(h): - if h >= 0 and h <= 11000: + if 0 <= h <= 11000: res = 101325 * ((288.15 - 0.0065 * h)/288.15) ** 5.25577 - return res - elif h > 11000 and h <= 20000: + elif 11000 < h <= 20000: res = 22632 * np.exp(-(h - 11000)/6341.62) - return res elif h > 20000: res = 5474.87 * ((216.65 + 0.0010 * (h - 20000))/216.65) ** (-34.163) - return res + return res + + def rho_air(h): res = p_air(h)/(R_air * T_air(h)) - return res \ No newline at end of file + return res + + +def \ No newline at end of file diff --git a/models/sun.py b/models/sun.py index dc5ece3..6f55c65 100644 --- a/models/sun.py +++ b/models/sun.py @@ -2,57 +2,111 @@ import astropy.units as unit import numpy as np from astropy.coordinates import EarthLocation, AltAz from astropy.coordinates import get_sun +from astropy.time import TimeISO +from input.natural_constants import * -def sun_angles_astropy(lat, lon, h, utc): + +class TimeYearDayTimeCustom(TimeISO): + """ + day-of-year as "". + The day-of-year (DOY) goes from 001 to 365 (366 in leap years). + The allowed subformat is: + - 'doy': day of year + """ + name = 'doy' # unique format name + subfmts = (('doy', + '%j', + '{yday:03d}'), + ('doy', + '%j', + '{yday:03d}'), + ('doy', + '%j', + '{yday:03d}')) + + +def sun_angles_astropy(lat, lon, h, utc): # get current sun elevation and azimuth through astropy loc = EarthLocation(lat=lat*unit.deg, lon=lon*unit.deg, height=h*unit.m) ref = AltAz(obstime=utc, location=loc) sun_pos = get_sun(utc).transform_to(ref) - AZ = sun_pos.az.degree - ELV = sun_pos.alt.degree + az = sun_pos.az.degree + elv = sun_pos.alt.degree - return AZ, ELV + return az, elv -def sun_angles_analytical(lat, lon, utc): - 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 - # STA = GMA + SEC - # SRV = (1.000001018 * (1 - EEO ** 2)) / (1 + EEO * np.cos(np.deg2rad(STA))) - 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)) - # SRA = np.rad2deg(np.arctan2(np.cos(np.deg2rad(OC)) * np.sin(np.deg2rad(SAL)), np.cos(np.deg2rad(SAL)))) # radian - 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 +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: - HA = TST / 4 - 180 + pass - 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 + 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 - 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 + if tst / 4 < 0: + ha = tst / 4 + 180 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 + 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 + + 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 + + return saa, sea # Azimuth, Elevation + + +def AirMass(p_air, p_0, ELV, h): # get atmospheric air mass over balloon + ELV_rad = np.deg2rad(ELV) # convert ELV from degree to radian + Dip = np.arccos(R_E / (R_E + h)) # geometric "dip" in radian + + if ELV_rad >= -Dip and ELV_rad < 0: + res = p_air/p_0 * (1 + ELV_rad/Dip) - 70 * ELV_rad/Dip + else: + res = (p_air/p_0) * ((1229 + (614 * np.sin(ELV_rad)) ** 2) ** (1/2) - 614 * np.sin(ELV_rad)) + + return res + + +def tau(ELV, h, p_air): # 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)) + else: + tau_atm = 0 + tau_atmIR = 0 + return tau_atm, tau_atmIR + + + + + - return SAA, SEA \ No newline at end of file diff --git a/models/test1.py b/models/test1.py deleted file mode 100644 index 75275a1..0000000 --- a/models/test1.py +++ /dev/null @@ -1,5 +0,0 @@ -import numpy as np - -def f(x): - res = np.sin(x) - return res \ No newline at end of file diff --git a/models/test2.py b/models/test2.py deleted file mode 100644 index eda5730..0000000 --- a/models/test2.py +++ /dev/null @@ -1,3 +0,0 @@ -from test1 import f - -print(f(1)) diff --git a/models/thermal.py b/models/thermal.py index b96ec77..4c82055 100644 --- a/models/thermal.py +++ b/models/thermal.py @@ -1,111 +1,5 @@ from input.natural_constants import * import numpy as np -from datetime import datetime -import time -from astropy.time import Time - -from astropy.time import TimeISO -class TimeYearDayTimeCustom(TimeISO): - """ - day-of-year as "". - The day-of-year (DOY) goes from 001 to 365 (366 in leap years). - The allowed subformat is: - - 'doy': day of year - """ - name = 'doy' # Unique format name - subfmts = (('doy', - '%j', - '{yday:03d}'), - ('doy', - '%j', - '{yday:03d}'), - ('doy', - '%j', - '{yday:03d}')) - -def AirMass(x, p_0, ELV, h): - p_air = x - ELV_rad = np.deg2rad(ELV) # convert ELV from degree to radian - Dip = np.arccos(R_E / (R_E + h)) # Dip in radian - - if ELV_rad >= -Dip and ELV_rad < 0: - res = p_air/p_0 * (1 + ELV_rad/Dip) - 70 * ELV_rad/Dip - else: - res = (p_air/p_0) * ((1229 + (614 * np.sin(ELV_rad)) ** 2) ** (1/2) - 614 * np.sin(ELV_rad)) - - return res - -#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)) -#else: -# tau_atm = 0 -# tau_atmIR = 0 -# -#doy = int(utc.doy) -# -#MA = (357.52911 + 0.98560028 * (utc.jd - 2451545)) % 360 # in degree, reference: see folder "literature" -#TA = MA + 2 * e * np.sin(np.deg2rad(MA)) + 5/4 * e ** 2 * np.sin(np.deg2rad(2 * MA)) -#I_Sun = 1367.5 * ((1 + e * np.cos(np.deg2rad(TA)))/(1 - e ** 2)) ** 2 -#I_SunZ = I_Sun * tau_atm -#q_sun = I_SunZ -# -#q_IRground = epsilon_ground * sigma * T_ground ** 4 -# -#q_IREarth = q_IRground * tau_atmIR -# -#if ELV <= 0: -# q_Albedo = 0 -#else: -# q_Albedo = Albedo * I_Sun * np.sin(np.deg2rad(ELV)) -# -### NEEDED: -# Diameter = ... -# V_relative = ... -# rho_air = ... -# rho_gas = ... -# T_film = ... -# T_air = ... -# T_gas = ... -# -#my_air = (1.458 * 10 ** -6 * T_air ** 1.5) / (T_air + 110.4) -#my_gas = 1.895 * 10 ** -5 * (T_gas/273.15) ** 0.647 -#k_air = 0.0241 * (T_air/273.15) ** 0.9 -#k_gas = 0.144 * (T_gas/273.15) ** 0.7 -#Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air -#Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas -# -#Gr_air = (rho_air ** 2 * g * np.abs(T_film - T_air) * Diameter ** 3) / (T_air * my_air ** 2) -#Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25 -#HC_free = Nu_air * k_air / Diameter -#Re = V_relative * Diameter * rho_air / my_air -#HC_forced = k_air / Diameter * (2 + 0.41 * Re ** 0.55) -#HC_internal = 0.13 * k_gas * ((rho_gas ** 2 * g * np.abs(T_film - T_gas) * Pr_gas) / (T_gas * my_air ** 2)) ** (1/3) -#HC_external = np.maximum(HC_free, HC_forced) -# -#Q_Sun = alpha_VIS * A_proj * q_sun * (1 + tau_VIS/(1 - r_VIS)) -#Q_Albedo = alpha_VIS * A_surf * q_Albedo * ViewFactor * (1 + tau_VIS/(1 - r_VIS)) -#Q_IREarth = alpha_IR * A_surf * q_IREarth * ViewFactor * (1 + tau_IR/(1 - r_IR)) -#Q_IRFilm = sigma * epsilon * alpha_IR * A_surf * T_film ** 4 * 1/(1 - r_IR) -#Q_IROut = sigma * epsilon * A_surf * T_film ** 4 * (1 * tau_IR/(1 - r_IR)) -#Q_ConvExt = HC_external * A_effective * (T_air - T_film) -#Q_ConvInt = HC_internal * A_effective * (T_film - T_gas) -# -#HalfConeAngle = np.arcsin(R_E/(R_E + h)) -#ViewFactor = (1 - np.cos(HalfConeAngle))/2 -# -# -# -#RoC = -V_z -# -#dT_gas = (Q_ConvInt/(gamma * c_v * M_gas) - (gamma - 1)/gamma * rho_air * g / (rho_gas * R_gas) * RoC) * dt -#dT_film = ((Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout)/(c_f * M_film)) * dt - - - - - - diff --git a/readme.txt b/readme.txt index aef4cc3..c33342f 100644 --- a/readme.txt +++ b/readme.txt @@ -13,10 +13,10 @@ netcdf4 1.3.3 cdsapi 0.2.7 (*) -(*) for CDS API please refer to: https://cds.climate.copernicus.eu/api-how-to - Alternatively, use ANACONDA and the environment(*.yml)-file in this repository. +(*) for proper setup of CDS API please refer to: https://cds.climate.copernicus.eu/api-how-to (in any case!) + (B) User Input: @@ -24,4 +24,4 @@ Alternatively, use ANACONDA and the environment(*.yml)-file in this repository. 1. Please define all input parameters under input/user_input.py 2. Script "DataRequest.py" needs to be run prior to generate netCDF-datafile for atmospheric data! - Please ensure sufficient disk space (1,25 GB) \ No newline at end of file + Please ensure sufficient disk space (currently: 1,25 GB) \ No newline at end of file