From 34a5aad93edbb65f50e54db74c45f1f5d0e9cc61 Mon Sep 17 00:00:00 2001 From: Marcel Frommelt Date: Tue, 9 Mar 2021 21:20:38 +0900 Subject: [PATCH] removed obsolote/outdated files to tidy up repository --- DataRequestNew.py | 55 ------ input/Trial_2021-01-12.py | 174 ----------------- main_atmodel.py | 278 --------------------------- main_netCDF.py | 277 --------------------------- main_netCDF_new.py | 387 -------------------------------------- main_old.py | 146 -------------- placeholder.txt | 0 7 files changed, 1317 deletions(-) delete mode 100644 DataRequestNew.py delete mode 100644 input/Trial_2021-01-12.py delete mode 100644 main_atmodel.py delete mode 100644 main_netCDF.py delete mode 100644 main_netCDF_new.py delete mode 100644 main_old.py delete mode 100644 placeholder.txt diff --git a/DataRequestNew.py b/DataRequestNew.py deleted file mode 100644 index d6447c0..0000000 --- a/DataRequestNew.py +++ /dev/null @@ -1,55 +0,0 @@ -import cdsapi - -c = cdsapi.Client() - - - - -r = c.retrieve( - 'reanalysis-era5-pressure-levels', - { - 'product_type': 'reanalysis', - 'variable': [ - 'fraction_of_cloud_cover', 'geopotential', 'relative_humidity', - 'specific_humidity', 'temperature', 'u_component_of_wind', - 'v_component_of_wind', 'vertical_velocity', - ], - 'pressure_level': [ - '1', '2', '3', - '5', '7', '10', - '20', '30', '50', - '70', '100', '125', - '150', '175', '200', - '225', '250', '300', - '350', '400', '450', - '500', '550', '600', - '650', '700', '750', - '775', '800', '825', - '850', '875', '900', - '925', '950', '975', - '1000', - ], - 'year': '2016', - 'month': '07', - 'day': [ - '11', '12', '13', - '14', '15', '16', - '17', '18', - ], - 'time': [ - '00:00', '01:00', '02:00', - '03:00', '04:00', '05:00', - '06:00', '07:00', '08:00', - '09:00', '10:00', '11:00', - '12:00', '13:00', '14:00', - '15:00', '16:00', '17:00', - '18:00', '19:00', '20:00', - '21:00', '22:00', '23:00', - ], - 'area': [ - 72, -111, 67, - 22, - ], - 'format': 'netcdf', - }) -r.download('test2021.nc') \ No newline at end of file diff --git a/input/Trial_2021-01-12.py b/input/Trial_2021-01-12.py deleted file mode 100644 index d5c3307..0000000 --- a/input/Trial_2021-01-12.py +++ /dev/null @@ -1,174 +0,0 @@ -import numpy as np -from astropy.time import Time - -# SIMPLE ATMOSPHERE MODEL - -def T_air(h): - if h >= 0 and h <= 11000: - res = 288.15 - 0.0065 * h - return res - elif h > 11000 and h <= 20000: - res = 216.65 - return res - elif h >= 20000: - res = 216.65 + 0.0010 * (h - 20000) - return res -def p_air(h): - if h >= 0 and h <= 11000: - res = 101325 * ((288.15 - 0.0065 * h)/288.15) ** 5.25577 - return res - elif h > 11000 and 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 -def rho_air(h): - res = p_air(h)/(R_air * T_air(h)) - return res - - -# USER INPUT - -start_height = 0 -start_utc = '2006-01-16 08:00:00.000' -start_lat = 0 # start latitude in [deg] -start_lon = 0 # start longitude in [deg] - - -# INITIALISATION - -t = 0 -h = start_height -lat = start_lat -lon = start_lon -utc = Time(start_utc) -v_z = 0 -p_gas = p_air(h) -T_gas = T_air(h) -T_film = T_gas -p_0 = 101325 - -dT_gas = (Q_ConvInt / (gamma * m_gas * c_v) - (gamma - 1) / gamma * rho_air(h) * 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 - - - -dT_gas = (Q_g/(c_pg * m_g) - (g * M_a * m_g * T_g)/(c_pg * T_a * M_g) * dz/dt) * dt - - -V_b_new = m_gas * R_gas * T_gas_new/p_gas_new # r.h.s. with updated values! - - -v_w = 0 - -T_g1 = T_g -T_a1 = T_a -# Q_g1 = ??? -T_e1 = T_e - -# T_g2 = T_g1 + dt * (Q_g1/(c_pg * m_g1) - g * M_a * T_g1 * dh/(c_pg * T_a1 * M_g)) -# T_e2 = T_e1 + Q_e/(c_env * m_e) * dt - -T_gn = T_g + dt * (Q_g/(c_pg * m_g) - g * M_a * T_g * dh/(c_pg * T_a * M_g)) -T_en = T_e + Q_e/(c_env * m_e) * dt - -s_n = s + dt * v + 0.5 * a * dt ** 2 -v_n = v - v_w + dt * a - - - - - - - -T_g = T_gn -T_e = T_en -s = s_n -v = v_n - - -V = m_g * R_g * T_g / (p_g) - - - -dV_b = V_b_new - V_b # calculating volume increment -dT_gas = Q_ConvInt/(c_v * m_g) * dt + (gamma - 1) * T_g * (dm_g/m_g - dV_b/V_b) - - - - -V_b = V_b_new # updating current volume - - -dT_film = ((Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film)) * dt - - - - - - - - - -m_balloon = A_balloon * rho_envelope -m_gas = Vol_balloon * (M_gas * p_gas)/(R_c * T_gas) -m_gross = m_balloon + m_payload -m_tot = m_balloon + m_gas + m_payload - -GI = F_buoyant - g * (rho_a - rho_g) * Vol_balloon - -W = m_gross * g - -m_virtual = m_tot + C_virtual * rho_a * Vol_balloon - -R_z = GI + T_z - D_z - W - -dV_A = (R_z * dt)/m_virtual - -F_Drag = 0.5 * C_D * rho * np.abs(V) * V * A_cross - -Re = (rho * V * l)/my - -C_D = 24/Re - -my = my_0 * (T_0 + C)/(T + C) * (T/T_0) ** (3/2) - -C_D = 0.23175 - 0.15757 * f + 0.04744 * f ** 2 - 7.0412 * 10 ** (-3) * f ** 3 + 5.1534 * 10 ** (-4) * f ** 4 - 1.4835 * 10 ** (-5) * f ** 5 - -Vol_balloon = m_g * R_g * T_g/p_g - - -Vol_ballon_new = m_g * R_g * T_g/p_g # r.h.s. with updated values! -dVol = Vol_ballon_new - Vol_balloon # calculating volume increment -Vol_ballon = Vol_ballon_new # updating current volume - - -L_z = g * (M_atm * p_atm/(R * T_atm) - M_gas * p_gas/(R * T_gas)) * Vol_balloon - - - -dT_g = Q_ConvInt/(c_v * m_g) * dt + (gamma - 1) * T_g * (dm_g/m_g - dVol/Vol) - -dT_e = (Q_sun + Q_albedo + Q_IREarth + Q_IRsky + Q_ConvExt - Q_ConvInt + Q_IRout)/(c_e * m_e) - -Pr_g = my_g * c_pg/kappa_g -Pr_a = my_a * c_pa/kappa_a - -beta_g = 1/T_g -beta_a = 1/T_a - -Gr_g = d_b ** 3 * rho_g ** 2 * g * beta_g * (T_e - T_g)/(my_g ** 2) -Gr_a = d_b ** 3 * rho_a ** 2 * g * beta_a * (T_e - T_a)/(my_a ** 2) - - - -h = h + dt * V + 0.5 * a * dt ** 2 -V = V + dt * a - - - - - - - diff --git a/main_atmodel.py b/main_atmodel.py deleted file mode 100644 index dacd3d2..0000000 --- a/main_atmodel.py +++ /dev/null @@ -1,278 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from datetime import datetime -start = datetime.now() -print(start) -begin_time = datetime.now() -from models.sun import sun_angles_analytical -from models.sun import sun_angles_astropy -from models.simple_atmosphere import T_air, p_air, rho_air -from models.drag import drag, c_d -from input.user_input import * -from input.natural_constants import * -from astropy.time import Time -import astropy.units as unit -from models.thermal import AirMass -from netCDF4 import Dataset - -#data = Dataset("test2021.nc", "r", format="NETCDF4") - -#ERAtime = data.variables['time'][:] # time -#ERAlat = data.variables['latitude'][:] # latitude -#ERAlon = data.variables['longitude'][:] # longitude -#ERAz = data.variables['z'][:]/g # geopotential to geopotential height -#ERApress = data.variables['level'][:] # pressure level -#ERAtemp = data.variables['t'][:] -#vw_x = data.variables['u'][:] # v_x -#vw_y = data.variables['v'][:] # v_y -#vw_z = data.variables['w'][:] # v_z - -# t: simulation time (0 .. t_max) in [s] -# utc = Time(start_utc) + t * u.second # current time (UTC) - -lat = start_lat # deg -lon = start_lon # deg - -# date = Time('2020-12-13 12:00:00.000') - - -# INITIALISATION - -t = 0 # simulation time in seconds -h = start_height -utc = Time(start_utc) -#epoch_diff = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 -#t_epoch = epoch_diff - -#sq_diff_lat = (ERAlat - lat) ** 2 -#sq_diff_lon = (ERAlon - lon) ** 2 -#sq_diff_time = (ERAtime - t_epoch) ** 2 - -#min_index_lat = sq_diff_lat.argmin() -#min_index_lon = sq_diff_lat.argmin() -#min_index_time = sq_diff_time.argmin() - -#sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 - -#min_index_h = sq_diff_h.argmin() - -#h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] -#p_grid = ERApress[min_index_h] * 100 -#T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] -#u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] -#v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] -#w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid - -#v = v_grid -#u = u_grid -#w = w_grid -#T_air = T_grid -#p_air = p_grid -#rho_air = p_air/(R_air * T_air) - -v_x = 0 -v_y = 0 -v_z = 0 - -T_gas = T_air(h) -T_film = T_gas - -t_list = [] -h_list = [] -lat_list = [] -lon_list = [] -p_list = [] -Temp_list = [] -rho_list = [] - -while t <= t_end and h >= 0: - t_list.append(t) - h_list.append(h) - lat_list.append(lat) - lon_list.append(lon) - - p_list.append(p_air(h)) - Temp_list.append(T_air(h)) - rho_list.append(rho_air(h)) - - - #t_epoch = epoch_diff + t/3600 - - #sq_diff_lat = (ERAlat - lat) ** 2 - #sq_diff_lon = (ERAlon - lon) ** 2 - #sq_diff_time = (ERAtime - t_epoch) ** 2 - - #min_index_lat = sq_diff_lat.argmin() - #min_index_lon = sq_diff_lat.argmin() - #min_index_time = sq_diff_time.argmin() - - #sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 - - #min_index_h = sq_diff_h.argmin() - - #h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] - #p_grid = ERApress[min_index_h] * 100 - #T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] - #u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] - #v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] - #w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid - - # p_air = p_grid - p_gas = p_air(h) - # rho_air = p_air / (R_air * T_air) - u = 0 - v = 0 - w = 0 - - v_relx = u - v_x - v_rely = v - v_y - v_relz = w - v_z - - v_rel = (v_relx ** 2 + v_rely ** 2 + v_relz ** 2) ** (1/2) - - p_gas = p_air(h) - rho_gas = p_gas/(R_gas * T_gas) # calculate gas density through ideal(!) gas equation - - V_b = m_gas/rho_gas # calculate balloon volume from current gas mass and gas density - - if V_b > V_design: - V_b = V_design - m_gas = V_design * rho_gas - else: - pass - - m_gross = m_pl + m_film - m_tot = m_pl + m_film + m_gas - m_virt = m_tot + c_virt * rho_air(h) * V_b - - d_b = 1.383 * V_b ** (1/3) # calculate diameter of balloon from its volume - L_goreB = 1.914 * V_b ** (1/3) - h_b = 0.748 * d_b - A_surf = 4.94 * V_b ** (2/3) - A_surf1 = 4.94 * V_design ** (2/3) * (1 - np.cos(np.pi * L_goreB/L_goreDesign)) - A_eff = 0.65 * A_surf + 0.35 * A_surf1 - A_top = np.pi/4 * d_b ** 2 - - D = drag(c_d, rho_air(h), d_b, v_rel) # calculate drag force - - if v_rel == 0: - Drag_x = 0 - Drag_y = 0 - Drag_z = 0 - else: - Drag_x = D * v_relx / v_rel - Drag_y = D * v_rely / v_rel - Drag_z = D * v_relz / v_rel - - I = g * V_b * (rho_air(h) - rho_gas) # calculate gross inflation - W = g * m_gross # calculate weight (force) - F = I - W + Drag_z - - AZ, ELV = sun_angles_analytical(lat, lon, utc) - - A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV))) - - # CALCULATIONS FOR THERMAL MODEL - - if ELV >= -(180 / np.pi * np.arccos(R_E / (R_E + h))): - tau_atm = 0.5 * (np.exp(-0.65 * AirMass(p_air(h), p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air(h), p_0, ELV, h))) - tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air(h) / p_0) + np.exp(-0.095 * p_air(h) / 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)) - - my_air = (1.458 * 10 ** -6 * T_air(h) ** 1.5) / (T_air(h) + 110.4) - my_gas = 1.895 * 10 ** -5 * (T_gas / 273.15) ** 0.647 - k_air = 0.0241 * (T_air(h) / 273.15) ** 0.9 - k_gas = 0.144 * (T_gas / 273.15) ** 0.7 - Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air(h) - Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas - - Gr_air = (rho_air(h) ** 2 * g * np.abs(T_film - T_air(h)) * d_b ** 3) / (T_air(h) * my_air ** 2) - Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25 - HC_free = Nu_air * k_air / d_b - Re = np.abs(v_relz) * d_b * rho_air(h) / my_air - - HC_forced = k_air / d_b * (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) - - HalfConeAngle = np.arcsin(R_E / (R_E + h)) - ViewFactor = (1 - np.cos(HalfConeAngle)) / 2 - - 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_eff * (T_air(h) - T_film) - Q_ConvInt = HC_internal * A_eff * (T_film - T_gas) - - # RoC = -v_z - # dT_gas = (Q_ConvInt / (gamma * m_gas * c_v) - (gamma - 1) / gamma * rho_air(h) * 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 - # v_w = w - # v = v_z - - s = h - - a_x = Drag_x/m_virt - a_y = Drag_y/m_virt - a_z = F/m_virt - - # DIFFERENTIAL EQUATIONS - - dx = dt * v_x + 0.5 * a_x * dt ** 2 # lon - dy = dt * v_y + 0.5 * a_y * dt ** 2 # lat - - lon = lon + np.rad2deg(np.arctan(dx/((6371229.0 + h) * np.cos(np.deg2rad(lat))))) - lat = lat + np.rad2deg(np.arctan(dy/(6371229.0 + h))) - - v_xn = v_x + dt * a_x - v_yn = v_y + dt * a_y - - s_n = s + dt * v_z + 0.5 * a_z * dt ** 2 - dh = s_n - s - v_n = v_z + dt * a_z - - T_gn = T_gas + dt * (Q_ConvInt / (gamma * c_v * m_gas) - g * R_gas * T_gas * dh / (c_v * gamma * T_air(s) * R_air)) - T_en = T_film + (Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film) * dt - - T_g = T_gn - T_e = T_en - s = s_n - v_x = v_xn - v_y = v_yn - v_z = v_n - - h = s - T_gas = T_g - T_film = T_e -# v_z = v - - t += dt # time increment - utc = Time(start_utc) + t * unit.second - -end = datetime.now() -print(end) -plt.plot(t_list, h_list, 'k-') -plt.plot(t_list, p_list, 'r-') -plt.plot(t_list, Temp_list, 'b-') -plt.plot(t_list, rho_list, 'g-') -plt.show() \ No newline at end of file diff --git a/main_netCDF.py b/main_netCDF.py deleted file mode 100644 index 15e7def..0000000 --- a/main_netCDF.py +++ /dev/null @@ -1,277 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from datetime import datetime -start = datetime.now() -print(start) -begin_time = datetime.now() -from models.sun import sun_angles_analytical -from models.sun import sun_angles_astropy -from models.drag import drag, c_d -from input.user_input import * -from input.natural_constants import * -from astropy.time import Time -import astropy.units as unit -from models.thermal import AirMass -from netCDF4 import Dataset - -data = Dataset("test2021.nc", "r", format="NETCDF4") - -ERAtime = data.variables['time'][:] # time -ERAlat = data.variables['latitude'][:] # latitude -ERAlon = data.variables['longitude'][:] # longitude -ERAz = data.variables['z'][:]/g # geopotential to geopotential height -ERApress = data.variables['level'][:] # pressure level -ERAtemp = data.variables['t'][:] -vw_x = data.variables['u'][:] # v_x -vw_y = data.variables['v'][:] # v_y -vw_z = data.variables['w'][:] # v_z - -# t: simulation time (0 .. t_max) in [s] -# utc = Time(start_utc) + t * u.second # current time (UTC) - -lat = start_lat # deg -lon = start_lon # deg - -# date = Time('2020-12-13 12:00:00.000') - - -# INITIALISATION - -t = 0 # simulation time in seconds -h = start_height -utc = Time(start_utc) -epoch_diff = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 -t_epoch = epoch_diff - -sq_diff_lat = (ERAlat - lat) ** 2 -sq_diff_lon = (ERAlon - lon) ** 2 -sq_diff_time = (ERAtime - t_epoch) ** 2 - -min_index_lat = sq_diff_lat.argmin() -min_index_lon = sq_diff_lat.argmin() -min_index_time = sq_diff_time.argmin() - -sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 - -min_index_h = sq_diff_h.argmin() - -h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] -p_grid = ERApress[min_index_h] * 100 -T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] -u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] -v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] -w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid - -v = v_grid -u = u_grid -w = w_grid -T_air = T_grid -p_air = p_grid -rho_air = p_air/(R_air * T_air) - -v_x = 0 -v_y = 0 -v_z = 0 - -T_gas = T_air -T_film = T_gas - -t_list = [] -h_list = [] -lat_list = [] -lon_list = [] -p_list = [] -Temp_list = [] -rho_list = [] - -while t <= t_end and h >= 0: - t_list.append(t) - h_list.append(h) - lat_list.append(lat) - lon_list.append(lon) - - p_list.append(p_air) - Temp_list.append(T_air) - rho_list.append(rho_air) - - t_epoch = epoch_diff + t/3600 - - sq_diff_lat = (ERAlat - lat) ** 2 - sq_diff_lon = (ERAlon - lon) ** 2 - sq_diff_time = (ERAtime - t_epoch) ** 2 - - min_index_lat = sq_diff_lat.argmin() - min_index_lon = sq_diff_lat.argmin() - min_index_time = sq_diff_time.argmin() - - sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 - - min_index_h = sq_diff_h.argmin() - - h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] - p_grid = ERApress[min_index_h] * 100 - T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] - u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] - v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] - w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid - - p_air = p_grid - p_gas = p_air - T_air = T_grid - rho_air = p_air / (R_air * T_air) - u = u_grid - v = v_grid - w = w_grid - - v_relx = u - v_x - v_rely = v - v_y - v_relz = w - v_z - - v_rel = (v_relx ** 2 + v_rely ** 2 + v_relz ** 2) ** (1/2) - - p_gas = p_air - rho_gas = p_gas/(R_gas * T_gas) # calculate gas density through ideal(!) gas equation - - V_b = m_gas/rho_gas # calculate balloon volume from current gas mass and gas density - - if V_b > V_design: - V_b = V_design - m_gas = V_design * rho_gas - else: - pass - - m_gross = m_pl + m_film - m_tot = m_pl + m_film + m_gas - m_virt = m_tot + c_virt * rho_air * V_b - - d_b = 1.383 * V_b ** (1/3) # calculate diameter of balloon from its volume - L_goreB = 1.914 * V_b ** (1/3) - h_b = 0.748 * d_b - A_surf = 4.94 * V_b ** (2/3) - A_surf1 = 4.94 * V_design ** (2/3) * (1 - np.cos(np.pi * L_goreB/L_goreDesign)) - A_eff = 0.65 * A_surf + 0.35 * A_surf1 - A_top = np.pi/4 * d_b ** 2 - - D = drag(c_d, rho_air, d_b, v_rel) # calculate drag force - - if v_rel == 0: - Drag_x = 0 - Drag_y = 0 - Drag_z = 0 - else: - Drag_x = D * v_relx/v_rel - Drag_y = D * v_rely/v_rel - Drag_z = D * v_relz/v_rel - - I = g * V_b * (rho_air - rho_gas) # calculate gross inflation - W = g * m_gross # calculate weight (force) - F = I - W + Drag_z - - AZ, ELV = sun_angles_analytical(lat, lon, utc) - - A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV))) - - # CALCULATIONS FOR THERMAL MODEL - - 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)) - - 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) * d_b ** 3) / (T_air * my_air ** 2) - Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25 - HC_free = Nu_air * k_air / d_b - Re = np.abs(v_relz) * d_b * rho_air / my_air - - HC_forced = k_air / d_b * (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) - - HalfConeAngle = np.arcsin(R_E / (R_E + h)) - ViewFactor = (1 - np.cos(HalfConeAngle)) / 2 - - 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_eff * (T_air - T_film) - Q_ConvInt = HC_internal * A_eff * (T_film - T_gas) - - # RoC = -v_z - # dT_gas = (Q_ConvInt / (gamma * m_gas * c_v) - (gamma - 1) / gamma * rho_air(h) * 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 - # v_w = w - # v = v_z - - s = h - - a_x = Drag_x/m_virt - a_y = Drag_y/m_virt - a_z = F/m_virt - - # DIFFERENTIAL EQUATIONS - - dx = dt * v_x + 0.5 * a_x * dt ** 2 # lon - dy = dt * v_y + 0.5 * a_y * dt ** 2 # lat - - lon = lon + np.rad2deg(np.arctan(dx/((6371229.0 + h) * np.cos(np.deg2rad(lat))))) - lat = lat + np.rad2deg(np.arctan(dy/(6371229.0 + h))) - - v_xn = v_x + dt * a_x - v_yn = v_y + dt * a_y - - s_n = s + dt * v_z + 0.5 * a_z * dt ** 2 - dh = s_n - s - v_n = v_z + dt * a_z - - T_gn = T_gas + dt * (Q_ConvInt / (gamma * c_v * m_gas) - g * R_gas * T_gas * dh / (c_v * gamma * T_air * R_air)) - T_en = T_film + (Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film) * dt - - T_g = T_gn - T_e = T_en - s = s_n - v_x = v_xn - v_y = v_yn - v_z = v_n - - h = s - T_gas = T_g - T_film = T_e -# v_z = v - - t += dt # time increment - utc = Time(start_utc) + t * unit.second - -end = datetime.now() -print(end) -plt.plot(t_list, h_list, 'k-') -plt.plot(t_list, p_list, 'r-') -plt.plot(t_list, Temp_list, 'b-') -plt.plot(t_list, rho_list, 'g-') -plt.show() diff --git a/main_netCDF_new.py b/main_netCDF_new.py deleted file mode 100644 index 792abd2..0000000 --- a/main_netCDF_new.py +++ /dev/null @@ -1,387 +0,0 @@ -import numpy as np -import math -import xarray as xr -import matplotlib.pyplot as plt -from scipy.spatial import cKDTree -from scipy import interpolate -import cartopy.crs as ccrs -from datetime import datetime -start = datetime.now() -print(start) -begin_time = datetime.now() -from models.sun import sun_angles_analytical -from models.sun import sun_angles_astropy -from models.drag import drag, c_d -from input.user_input import * -from input.natural_constants import * -from astropy.time import Time -import astropy.units as unit -from models.thermal import AirMass -from netCDF4 import Dataset -import pandas as pd - -data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Mappe1.xls', sheet_name='Tabelle3') - -comp_time = pd.DataFrame(data, columns= ['Time']) -comp_height = pd.DataFrame(data, columns= ['Height']) - - -def transform(lon, lat, t): - # WGS 84 reference coordinate system parameters - A = 6378137.0 # major axis [m] - E2 = 6.69437999014e-3 # eccentricity squared - - t_s = 1000 * t - - lon_rad = np.radians(lon) - lat_rad = np.radians(lat) - # convert to cartesian coordinates - r_n = A / (np.sqrt(1 - E2 * (np.sin(lat_rad) ** 2))) - x = r_n * np.cos(lat_rad) * np.cos(lon_rad) - y = r_n * np.cos(lat_rad) * np.sin(lon_rad) - z = r_n * (1 - E2) * np.sin(lat_rad) - return x, y, z, t_s - -data = Dataset("test2021.nc", "r", format="NETCDF4") - -ERAtime = data.variables['time'][:] # time -ERAlat = data.variables['latitude'][:] # latitude -ERAlon = data.variables['longitude'][:] # longitude -ERAz = data.variables['z'][:]/g # geopotential to geopotential height -ERApress = data.variables['level'][:] # pressure level -ERAtemp = data.variables['t'][:] # temperature in K - -vw_x = data.variables['u'][:] # v_x -vw_y = data.variables['v'][:] # v_y -vw_z = data.variables['w'][:] # v_z - -lon_era2d, lat_era2d, time_era = np.meshgrid(ERAlon, ERAlat, ERAtime) - -xs, ys, zs, ts = transform(lon_era2d.flatten(), lat_era2d.flatten(), time_era.flatten()) - -tree = cKDTree(np.column_stack((xs, ys, zs, ts))) # ! - - - - -lat = start_lat # deg -lon = start_lon # deg - -t = 0 # simulation time in seconds -h = start_height -utc = Time(start_utc) -epoch_diff = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 -t_epoch = epoch_diff - - - -xt, yt, zt, tt = transform(lon, lat, t) # test coordinates - -d, inds = tree.query(np.column_stack((xt, yt, zt, tt)), k=8) # longitude, latitude, time in h # ! -w = 1.0 / d[0] ** 2 - -lat_sel = np.unravel_index(inds[0], lon_era2d.shape)[0] -lon_sel = np.unravel_index(inds[0], lon_era2d.shape)[1] -time_sel = np.unravel_index(inds[0], lon_era2d.shape)[2] - -interp4d_height = np.ma.dot(w, ERAz[time_sel, :, lat_sel, lon_sel]) / np.sum(w) -interp4d_temp = np.ma.dot(w, ERAtemp[time_sel, :, lat_sel, lon_sel]) / np.sum(w) -interp4d_vw_x = np.ma.dot(w, vw_x[time_sel, :, lat_sel, lon_sel]) / np.sum(w) -interp4d_vw_y = np.ma.dot(w, vw_y[time_sel, :, lat_sel, lon_sel]) / np.sum(w) -interp4d_vw_z = np.ma.dot(w, vw_z[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - -pressure_hPa = np.array([1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 125, 150, 175, 200,225, 250, 300, 350, 400, - 450, 500, 550, 600, 650, 700, 750, 775, 800, 825, 850, 875, 900, 925, 950, 975, 1000]) - -pressure = 100 * pressure_hPa -print(pressure) - - -# height_interp1d = interpolate.interp1d(interp4d_height, pressure) -temp_interp1d = interpolate.interp1d(interp4d_height, interp4d_temp) -press_interp1d = interpolate.interp1d(interp4d_height, pressure) -vw_x_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_x) -vw_y_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_y) -vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z) - - -#u = 0 -#v = 0 -#w = 0 -T_air = temp_interp1d(h) -p_air = press_interp1d(h) -rho_air = p_air/(R_air * T_air) - -u = vw_x_interp1d(h) -v = vw_y_interp1d(h) -w = -1 / g * vw_z_interp1d(h) * R_air * T_air / p_air - - -v_x = 0 -v_y = 0 -v_z = 0 - -v_rel = 0 - -T_gas = T_air -T_film = T_gas - -t_list = [] -h_list = [] -v_list = [] -lat_list = [] -lon_list = [] -p_list = [] -Temp_list = [] -rho_list = [] - -while t <= t_end and h >= 0: - t_list.append(t) - h_list.append(h) - v_list.append(v_rel) - lat_list.append(lat) - lon_list.append(lon) - - p_list.append(p_air) - Temp_list.append(T_air) - rho_list.append(rho_air) - - t_epoch = epoch_diff + t/3600 - - xt, yt, zt, tt = transform(lon, lat, t_epoch) # current balloon coordinates in cartesian coordinates: x [m], y [m], z [m], time [h since 1900-01-01 00:00:00.0] - - d, inds = tree.query(np.column_stack((xt, yt, zt, tt)), k=8) # longitude, latitude, time in h # ! - w = 1.0 / d[0] ** 2 - - lat_sel = np.unravel_index(inds[0], lon_era2d.shape)[0] - lon_sel = np.unravel_index(inds[0], lon_era2d.shape)[1] - time_sel = np.unravel_index(inds[0], lon_era2d.shape)[2] - - interp4d_height = np.ma.dot(w, ERAz[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - interp4d_temp = np.ma.dot(w, ERAtemp[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - interp4d_vw_x = np.ma.dot(w, vw_x[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - interp4d_vw_y = np.ma.dot(w, vw_y[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - interp4d_vw_z = np.ma.dot(w, vw_z[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - - press_interp1d = interpolate.interp1d(interp4d_height, pressure) - temp_interp1d = interpolate.interp1d(interp4d_height, interp4d_temp) - vw_x_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_x) - vw_y_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_y) - vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z) - - - try: - p_air = press_interp1d(h) - T_air = temp_interp1d(h) - except: - if (h > interp4d_height[0]): - h = interp4d_height[0] - else: - h = interp4d_height[-1] - - p_air = press_interp1d(h) - T_air = temp_interp1d(h) - - #print("height") - #print(h) - #print("temperature") - #print(temp_interp1d(h)) - #print("pressure") - #print(press_interp1d(h)) - #print("vw_x") - #print(vw_x_interp1d(h)) - - p_gas = p_air - rho_air = p_air / (R_air * T_air) - u = vw_x_interp1d(h) - v = vw_y_interp1d(h) - w = -1/g * vw_z_interp1d(h) * R_air * T_air / p_air - - v_relx = u - v_x - v_rely = v - v_y - v_relz = w - v_z - - v_rel = (v_relx ** 2 + v_rely ** 2 + v_relz ** 2) ** (1/2) - - rho_gas = p_gas/(R_gas * T_gas) # calculate gas density through ideal(!) gas equation - - V_b = m_gas/rho_gas # calculate balloon volume from current gas mass and gas density - - if V_b > V_design: - V_b = V_design - #m_gas = V_design * rho_gas - else: - pass - - m_gross = m_pl + m_film - m_tot = m_pl + m_film + m_gas - m_virt = m_tot + c_virt * rho_air * V_b - - d_b = 1.383 * V_b ** (1/3) # calculate diameter of balloon from its volume - L_goreB = 1.914 * V_b ** (1/3) - h_b = 0.748 * d_b - A_surf = 4.94 * V_b ** (2/3) - A_surf1 = 4.94 * V_design ** (2/3) * (1 - np.cos(np.pi * L_goreB/L_goreDesign)) - A_eff = 0.65 * A_surf + 0.35 * A_surf1 - A_top = np.pi/4 * d_b ** 2 - - D = drag(c_d, rho_air, d_b, v_rel) # calculate drag force - - if v_rel == 0: - Drag_x = 0 - Drag_y = 0 - Drag_z = 0 - else: - Drag_x = D * v_relx/v_rel - Drag_y = D * v_rely/v_rel - Drag_z = D * v_relz/v_rel - - I = g * V_b * (rho_air - rho_gas) # calculate gross inflation - W = g * m_gross # calculate weight (force) - F = I - W + Drag_z - - AZ, ELV = sun_angles_analytical(lat, lon, utc) - - A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV))) - - # CALCULATIONS FOR THERMAL MODEL - - 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)) - - 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) * d_b ** 3) / (T_air * my_air ** 2) - Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25 - HC_free = Nu_air * k_air / d_b - Re = np.abs(v_relz) * d_b * rho_air / my_air - - HC_forced = k_air / d_b * (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) - - HalfConeAngle = np.arcsin(R_E / (R_E + h)) - ViewFactor = (1 - np.cos(HalfConeAngle)) / 2 - - 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_eff * (T_air - T_film) - Q_ConvInt = HC_internal * A_eff * (T_film - T_gas) - - # RoC = -v_z - # dT_gas = (Q_ConvInt / (gamma * m_gas * c_v) - (gamma - 1) / gamma * rho_air(h) * 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 - # v_w = w - # v = v_z - - s = h - - a_x = Drag_x/m_virt - a_y = Drag_y/m_virt - a_z = F/m_virt - - # DIFFERENTIAL EQUATIONS - - dx = dt * v_x + 0.5 * a_x * dt ** 2 # lon - dy = dt * v_y + 0.5 * a_y * dt ** 2 # lat - - lon = lon + np.rad2deg(np.arctan(dx/((6371229.0 + h) * np.cos(np.deg2rad(lat))))) - lat = lat + np.rad2deg(np.arctan(dy/(6371229.0 + h))) - - v_xn = v_x + dt * a_x - v_yn = v_y + dt * a_y - - s_n = s + dt * v_z + 0.5 * a_z * dt ** 2 - dh = s_n - s - v_n = v_z + dt * a_z - - T_gn = T_gas + dt * (Q_ConvInt / (gamma * c_v * m_gas) - g * R_gas * T_gas * dh / (c_v * gamma * T_air * R_air)) - T_en = T_film + (Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film) * dt - - T_g = T_gn - T_e = T_en - s = s_n - v_x = v_xn - v_y = v_yn - v_z = v_n - - h = s - T_gas = T_g - T_film = T_e -# v_z = v - - t += dt # time increment - utc = Time(start_utc) + t * unit.second - - - -#print(len(lon_list), len(lat_list)) -#""" -#plt.plot(start_lon, start_lat, 'rx') -#plt.plot(ax=ax, lon_list, lat_list, transform=ccrs.PlateCarree()) -#plt.show() -#""" - - -### WORKS! -### -###dset = xr.open_dataset("test2021.nc") -### -###print(dset['t'][1][0]) # [time] [level] -### -#fig = plt.figure() #figsize=[120,50]) -### -#ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) -### -###dset['t'][1][0].plot(ax=ax, cmap='jet', transform=ccrs.PlateCarree()) -###ax.coastlines() -###plt.show() - -#print(len(ERAlat)) -#print(len(ERAlon)) -#print(len(ERAz)) - - - - - -#""" -plt.plot(comp_time, comp_height, 'r--', label='PoGo Flight Test') -#plt.plot(t_list, v_list, 'k--', label='Ballon v_rel') -plt.plot(t_list, h_list, 'k-', label='Balloon Altitude') -#plt.plot(t_list, p_list, 'r-', label='Air Pressure [Pa]') -#plt.plot(t_list, Temp_list, 'b-', label='Air Temperature [K]') -#plt.plot(t_list, rho_list, 'g-', label='Air Density in [kg/m$^3$]') -plt.xlabel('time in s') -plt.ylabel('Balloon Altitude in m') -plt.legend() -plt.show() -#""" diff --git a/main_old.py b/main_old.py deleted file mode 100644 index 79542e9..0000000 --- a/main_old.py +++ /dev/null @@ -1,146 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from datetime import datetime -begin_time = datetime.now() -from models.sun import sun_angles_analytical -from models.sun import sun_angles_astropy -from models.drag import drag, c_d -from models.simple_atmosphere import T_air, p_air, rho_air -from input.user_input import * -from input.natural_constants import * -from astropy.time import Time -import astropy.units as u -from models.thermal import AirMass - -# t: simulation time (0 .. t_max) in [s] -# utc = Time(start_utc) + t * u.second # current time (UTC) - -lat = start_lat # deg -lon = start_lon # deg - -date = Time('2020-12-13 12:00:00.000') - - -# INITIALISATION - -t = 0 -h = start_height -utc = Time(start_utc) -v_z = 0 -p_gas = p_air(h) -T_gas = T_air(h) -T_film = T_gas - -t_list = [] -h_list = [] - -while t <= t_end and h >= 0: - t_list.append(t) - h_list.append(h) - - rho_gas = p_gas/(R_gas * T_gas) # calculate gas density through ideal(!) gas equation - - V_b = m_gas/rho_gas # calculate balloon volume from current gas mass and gas density - - if V_b > V_design: - V_b = V_design - m_gas = V_design * rho_gas - else: - pass - - m_gross = m_pl + m_film - m_tot = m_pl + m_film + m_gas - m_virt = m_tot + c_virt * rho_air(h) * V_b - - d_b = 1.383 * V_b ** (1/3) # calculate diameter of balloon from its volume - L_goreB = 1.914 * V_b ** (1/3) - h_b = 0.748 * d_b - A_surf = 4.94 * V_b ** (2/3) - A_surf1 = 4.94 * V_design ** (2/3) * (1 - np.cos(np.pi * L_goreB/L_goreDesign)) - A_eff = 0.65 * A_surf + 0.35 * A_surf1 - A_top = np.pi/4 * d_b ** 2 - - D = drag(c_d, rho_air(h), d_b, v_z) # calculate drag force - I = g * V_b * (rho_air(h) - rho_gas) # calculate gross inflation - W = g * m_gross # calculate weight (force) - F = I - W - D * np.sign(v_z) - - v_zrel = v_z - - AZ, ELV = sun_angles_analytical(lat, lon, utc) - - A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV))) - - # CALCULATIONS FOR THERMAL MODEL - - if ELV >= -(180 / np.pi * np.arccos(R_E / (R_E + h))): - tau_atm = 0.5 * (np.exp(-0.65 * AirMass(p_air(h), p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air(h), p_0, ELV, h))) - tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air(h) / p_0) + np.exp(-0.095 * p_air(h) / 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)) - - my_air = (1.458 * 10 ** -6 * T_air(h) ** 1.5) / (T_air(h) + 110.4) - my_gas = 1.895 * 10 ** -5 * (T_gas / 273.15) ** 0.647 - k_air = 0.0241 * (T_air(h) / 273.15) ** 0.9 - k_gas = 0.144 * (T_gas / 273.15) ** 0.7 - Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air(h) - Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas - - Gr_air = (rho_air(h) ** 2 * g * np.abs(T_film - T_air(h)) * d_b ** 3) / (T_air(h) * my_air ** 2) - Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25 - HC_free = Nu_air * k_air / d_b - Re = np.abs(v_zrel) * d_b * rho_air(h) / my_air - - HC_forced = k_air / d_b * (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) - - HalfConeAngle = np.arcsin(R_E / (R_E + h)) - ViewFactor = (1 - np.cos(HalfConeAngle)) / 2 - - 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_eff * (T_air(h) - T_film) - Q_ConvInt = HC_internal * A_eff * (T_film - T_gas) - - RoC = -v_z - - dT_gas = (Q_ConvInt / (gamma * m_gas * c_v) - (gamma - 1) / gamma * rho_air(h) * 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 - - ### - - dv_z = F/m_virt * dt # velocity increment - - v_z += dv_z * dt # velocity differential - h += v_z * dt # altitude differential - - p_gas = p_air(h) - T_gas = T_gas + dT_gas * dt - T_film = T_film + dT_film * dt - - t += dt # time increment - utc = Time(start_utc) + t * u.second - -plt.plot(t_list, h_list) -plt.show() diff --git a/placeholder.txt b/placeholder.txt deleted file mode 100644 index e69de29..0000000