From c1e6b195e843024d6147078a395fd4aa7c0f1981 Mon Sep 17 00:00:00 2001 From: Marcel Frommelt Date: Mon, 3 May 2021 17:39:57 +0900 Subject: [PATCH] fixed major bug for thermal absorption of balloon + minor restructuring --- main.py | 775 ++++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 523 insertions(+), 252 deletions(-) diff --git a/main.py b/main.py index e49633e..861f20b 100644 --- a/main.py +++ b/main.py @@ -1,245 +1,410 @@ +import sys +import warnings import numpy as np -import math import xarray as xr -import matplotlib.pyplot as plt -from scipy.spatial import cKDTree -from scipy import interpolate -from scipy.integrate import odeint, solve_ivp +import pandas as pd import cartopy.crs as ccrs +import astropy.units as unit +import matplotlib.pyplot as plt +from dask import delayed from datetime import datetime -start = datetime.now() -print(start) -begin_time = datetime.now() -from models.sun import sun_angles_analytical, tau -# from models.sun import sun_angles_astropy -from models.drag import drag, Cd_model -from models.transformation import visible_cells, transform, radii +from netCDF4 import Dataset +from astropy.time import Time +from scipy import interpolate +from scipy.spatial import cKDTree +from scipy.integrate import solve_ivp from input.user_input import * from input.natural_constants import * -from astropy.time import Time -import astropy.units as unit -# from models.thermal import -from netCDF4 import Dataset -import pandas as pd +from models.gravity import grav +from models.valving import valving +from models.ballasting import ballasting +from models.sun import sun_angles_analytical, tau +from models.drag import drag, cd_PalumboLow, cd_Palumbo, cd_PalumboHigh, cd_sphere +from models.transformation import visible_cells, transform, radii -data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Mappe1.xls', sheet_name='Tabelle3') # Tabelle3 +starttime = datetime.now() + +if not sys.warnoptions: + warnings.simplefilter("ignore") + +data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Data_PoGo2016.xls', sheet_name='Tabelle3') # Tabelle3 comp_time = pd.DataFrame(data, columns=['Time']).to_numpy().squeeze() comp_height = pd.DataFrame(data, columns=['Height']).to_numpy().squeeze() comp_lat = pd.DataFrame(data, columns=['Latitude']).to_numpy().squeeze() comp_lon = pd.DataFrame(data, columns=['Longitude']).to_numpy().squeeze() -print("Initialisiere Simulation...") -print("STARTORT: Längengrad %.4f, Breitengrad: %.4f" % (start_lon, start_lat)) -print("STARTZEIT: " + str(start_utc) + " (UTC)") -#print("SIMULATION INITIALISIERT. Simulierte Flugdauer %.2f Stunden (oder %.2f Tage)" % ( -#t_end / 3600, t_end / (3600 * 24))) -#print("GESCHÄTZE SIMULATIONSDAUER: %.2f Minuten" % (t_end / (3600 * 10))) +print("") +print("Initialising simulation...") +print("Launch location : longitude %.4f, latitude: %.4f" % (start_lon, start_lat)) +print("Launch time: " + str(start_utc) + " (UTC)") +print("") +print("Reading ERA5-datasets, please wait.") + +ascend_data = xr.open_dataset("ascend_2019_kiruna.nc") +#float_data = xr.open_dataset("float_2019.nc") + +df = xr.open_mfdataset(['float1_2019.nc', 'float2_2019.nc', 'float3_2019.nc', 'float4_2019.nc'], combine='by_coords', concat_dim="time", parallel=True) + +startNC = Dataset('float1_2019.nc') +endNC = Dataset('float4_2019.nc') +start = int(startNC.variables['time'][0]) +end = int(endNC.variables['time'][-1]) + +float_data = df.assign_coords(time=np.linspace(start, end, (end - start) + 1)) + +single_data = xr.open_dataset("single_2019.nc") + +ERAtime = float_data.variables['time'][:] # time -level_data = Dataset("level.nc", "r", format="NETCDF4") -single_data = Dataset("single.nc", "r", format="NETCDF4") +# ERA5 MULTI-LEVEL ASCENT (WIND + ATMOSPHERIC DATA DURING ASCENT) -ERAtime = level_data.variables['time'][:] # time - -ERAlat = level_data.variables['latitude'][:] # latitude (multi-level) -ERAlon = level_data.variables['longitude'][:] # longitude (multi-level) -ERAlats = single_data.variables['latitude'][:] # latitude (single level) -ERAlons = single_data.variables['longitude'][:] # longitude (single level) +ERAlat0 = ascend_data.variables['latitude'][:].values # latitude [deg] +ERAlon0 = ascend_data.variables['longitude'][:].values # longitude [deg] +ERAz_ascend = ascend_data.variables['z'][:].values / g # geopotential [m^-2/s^-2] to geopotential height [m] +ERApress_ascend = ascend_data.variables['level'][:].values # pressure level [-] +ERAtemp_ascend = ascend_data.variables['t'][:].values # air temperature in K +vw_x_ascend = ascend_data.variables['u'][:].values # v_x in [m/s] +vw_y_ascend = ascend_data.variables['v'][:].values # v_y in [m/s] +vw_z_ascend = ascend_data.variables['w'][:].values # v_z in [m/s] -ERAz = level_data.variables['z'][:] / g # geopotential to geopotential height -ERApress = level_data.variables['level'][:] # pressure level -ERAtemp = level_data.variables['t'][:] # air temperature in K +# ERA5 MULTI-LEVEL FLOAT (WIND + ATMOSPHERIC DATA DURING FLOAT) + +ERAlat1 = float_data.variables['latitude'][:].values # latitude [deg] +ERAlon1 = float_data.variables['longitude'][:].values # longitude [deg] +ERAz_float = float_data.variables['z'][:].values / g # geopotential [m^-2/s^-2] to geopotential height [m] +ERApress_float = float_data.variables['level'][:].values # pressure level [-] +ERAtemp_float = float_data.variables['t'][:].values # air temperature in [K] +vw_x_float = float_data.variables['u'][:].values # v_x in [m/s] +vw_y_float = float_data.variables['v'][:].values # v_y in [m/s] +vw_z_float = float_data.variables['w'][:].values # v_z in [m/s] +# ERA5 SINGLE-LEVEL (RADIATIVE ENVIRONMENT) -ERAtcc = single_data.variables['tcc'][:] # total cloud cover [0 to 1] -ERAal = single_data.variables['fal'][:] # forecast albedo [0 to 1] -# ERAst = single_data.variables['skt'][:] # skin (ground) temperature in K -ERAcbh = single_data.variables['cbh'][:] # cloud base height in m -ERAlcc = single_data.variables['lcc'][:] # low cloud cover -ERAmcc = single_data.variables['mcc'][:] # medium cloud cover -ERAhcc = single_data.variables['hcc'][:] # high cloud cover -ERAssr = single_data.variables['ssr'][:] # surface net solar radiation -ERAstrn = single_data.variables['str'][:] # surface net thermal radiation -ERAstrd = single_data.variables['strd'][:] # surface thermal radiation downwards -ERAssrd = single_data.variables['ssrd'][:] # surface solar radiation downwards -ERAtsr = single_data.variables['tsr'][:] # top net solar radiation -ERAttr = single_data.variables['ttr'][:] # top net thermal radiation -ERAsst = single_data.variables['sst'][:] # sea surface temperature -ERAtisr = single_data.variables['tisr'][:] # TOA incident solar radiation -ERAstrdc = single_data.variables['strdc'][:] # surface thermal radiation downward clear-sky +ERAlat2 = single_data.variables['latitude'][:].values # latitude [deg] +ERAlon2 = single_data.variables['longitude'][:].values # longitude [deg] +ERAtcc = single_data.variables['tcc'][:].values # total cloud cover [-] +ERAskt = single_data.variables['skt'][:].values # skin (ground) temperature in [K] +ERAcbh = single_data.variables['cbh'][:].values # cloud base height in [m] +ERAlcc = single_data.variables['lcc'][:].values # low cloud cover [-] +ERAmcc = single_data.variables['mcc'][:].values # medium cloud cover [-] +ERAhcc = single_data.variables['hcc'][:].values # high cloud cover [-] +ERAssr = single_data.variables['ssr'][:].values # hourly accumulated surface net solar radiation [J/m^2] +ERAstrn = single_data.variables['str'][:].values # hourly accumulated surface net thermal radiation [J/m^2] +ERAstrd = single_data.variables['strd'][:].values # hourly accumulated surface thermal radiation downwards [J/m^2] +ERAssrd = single_data.variables['ssrd'][:].values # hourly accumulated surface solar radiation downwards [J/m^2] +ERAtsr = single_data.variables['tsr'][:].values # hourly accumulated top net solar radiation [J/m^2] +ERAttr = single_data.variables['ttr'][:].values # hourly accumulated top net thermal radiation [J/m^2] +ERAtisr = single_data.variables['tisr'][:].values # hourly accumulated TOA incident solar radiation [J/m^2] +ERAstrdc = single_data.variables['strdc'][:].values # hourly accumulated surface thermal radiation downward clear-sky [J/m^2] +ERAsp = single_data.variables['sp'][:].values # surface pressure in [Pa] +print("Finished reading ERA5-datasets.") -vw_x = level_data.variables['u'][:] # v_x -vw_y = level_data.variables['v'][:] # v_y -vw_z = level_data.variables['w'][:] # v_z +lon_era2d0, lat_era2d0 = np.meshgrid(ERAlon0, ERAlat0) +lon_era2d1, lat_era2d1 = np.meshgrid(ERAlon1, ERAlat1) +lon_era2d2, lat_era2d2 = np.meshgrid(ERAlon2, ERAlat2) -lon_era2d, lat_era2d = np.meshgrid(ERAlon, ERAlat) -lon_era2ds, lat_era2ds = np.meshgrid(ERAlons, ERAlats) +xs0, ys0, zs0 = transform(lon_era2d0.flatten(), lat_era2d0.flatten()) +xs1, ys1, zs1 = transform(lon_era2d1.flatten(), lat_era2d1.flatten()) +xs2, ys2, zs2 = transform(lon_era2d2.flatten(), lat_era2d2.flatten()) -xs, ys, zs = transform(lon_era2d.flatten(), lat_era2d.flatten()) -xs2, ys2, zs2 = transform(lon_era2ds.flatten(), lat_era2ds.flatten()) - -tree = cKDTree(np.column_stack((xs, ys, zs))) +print("") +tree0 = cKDTree(np.column_stack((xs0, ys0, zs0))) +tree1 = cKDTree(np.column_stack((xs1, ys1, zs1))) tree2 = cKDTree(np.column_stack((xs2, ys2, zs2))) +print("Built kd-trees") -# INITIALISATION TIME-COUNTERS - - -# MAYBE LATER -# t_pre = int(t_epoch) -# t_post = t_pre + 1 -# t1 = np.searchsorted(ERAtime, t_pre, side='left') -# t2 = t1 + 1 - def ERA5Data(lon, lat, h, t, deltaT_ERA): - t_epoch = deltaT_ERA + t/3600 + t_epoch = deltaT_ERA + t / 3600 t_pre = int(t_epoch) - # t_post = t_pre + 1 - t_pre_ind = t_pre - ERAtime[0] + t_pre_ind = t_pre - int(ERAtime[0]) t_post_ind = t_pre_ind + 1 xt, yt, zt = transform(lon, lat) # current coordinates - d1, inds1 = tree.query(np.column_stack((xt, yt, zt)), k=4) # longitude, latitude - d2, inds2 = tree2.query(np.column_stack((xt, yt, zt)), k=visible_cells(h)) # longitude, latitude - w1 = 1.0 / d1[0] # ** 2 + d0, inds0 = tree0.query(np.column_stack((xt, yt, zt)), k=4) + d1, inds1 = tree1.query(np.column_stack((xt, yt, zt)), k=4) # longitude, latitude + d2, inds2 = tree2.query(np.column_stack((xt, yt, zt)), k=visible_cells(h)) # longitude, latitude visible_cells(h) + + w0 = 1.0 / d0[0] + w1 = 1.0 / d1[0] w2 = 1.0 / d2[0] ** 2 - lat_ind1 = np.unravel_index(inds1[0], lon_era2d.shape)[0] - lon_ind1 = np.unravel_index(inds1[0], lon_era2d.shape)[1] - lat_ind2 = np.unravel_index(inds2[0], lon_era2ds.shape)[0] - lon_ind2 = np.unravel_index(inds2[0], lon_era2ds.shape)[1] + lat_ind0 = np.unravel_index(inds0[0], lon_era2d0.shape)[0] + lon_ind0 = np.unravel_index(inds0[0], lon_era2d0.shape)[1] + lat_ind1 = np.unravel_index(inds1[0], lon_era2d1.shape)[0] + lon_ind1 = np.unravel_index(inds1[0], lon_era2d1.shape)[1] + lat_ind2 = np.unravel_index(inds2[0], lon_era2d2.shape)[0] + lon_ind2 = np.unravel_index(inds2[0], lon_era2d2.shape)[1] - interp4d_temp_pre = np.ma.dot(w1, ERAtemp[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_temp_post = np.ma.dot(w1, ERAtemp[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_temp = (interp4d_temp_post - interp4d_temp_pre) * (t_epoch - t_pre) + interp4d_temp_pre + if h >= 30000: + try: + interp4d_temp_pre = np.ma.dot(w1, ERAtemp_float[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_temp_post = np.ma.dot(w1, ERAtemp_float[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_temp = (interp4d_temp_post - interp4d_temp_pre) * (t_epoch - t_pre) + interp4d_temp_pre - interp4d_height_pre = np.ma.dot(w1, ERAz[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_height_post = np.ma.dot(w1, ERAz[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_height = (interp4d_height_post - interp4d_height_pre) * (t_epoch - t_pre) + interp4d_height_pre + interp4d_height_pre = np.ma.dot(w1, ERAz_float[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_height_post = np.ma.dot(w1, ERAz_float[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_height = (interp4d_height_post - interp4d_height_pre) * (t_epoch - t_pre) + interp4d_height_pre - interp4d_vw_x_pre = np.ma.dot(w1, vw_x[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_vw_x_post = np.ma.dot(w1, vw_x[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_vw_x = (interp4d_vw_x_post - interp4d_vw_x_pre) * (t_epoch - t_pre) + interp4d_vw_x_pre + interp4d_vw_x_pre = np.ma.dot(w1, vw_x_float[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_vw_x_post = np.ma.dot(w1, vw_x_float[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_vw_x = (interp4d_vw_x_post - interp4d_vw_x_pre) * (t_epoch - t_pre) + interp4d_vw_x_pre - interp4d_vw_y_pre = np.ma.dot(w1, vw_y[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_vw_y_post = np.ma.dot(w1, vw_y[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_vw_y = (interp4d_vw_y_post - interp4d_vw_y_pre) * (t_epoch - t_pre) + interp4d_vw_y_pre + interp4d_vw_y_pre = np.ma.dot(w1, vw_y_float[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_vw_y_post = np.ma.dot(w1, vw_y_float[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_vw_y = (interp4d_vw_y_post - interp4d_vw_y_pre) * (t_epoch - t_pre) + interp4d_vw_y_pre - interp4d_vw_z_pre = np.ma.dot(w1, vw_z[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_vw_z_post = np.ma.dot(w1, vw_z[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) - interp4d_vw_z = (interp4d_vw_z_post - interp4d_vw_z_pre) * (t_epoch - t_pre) + interp4d_vw_z_pre + interp4d_vw_z_pre = np.ma.dot(w1, vw_z_float[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_vw_z_post = np.ma.dot(w1, vw_z_float[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1) + interp4d_vw_z = (interp4d_vw_z_post - interp4d_vw_z_pre) * (t_epoch - t_pre) + interp4d_vw_z_pre - albedo_pre = np.ma.dot(w2, ERAal[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) - albedo_post = np.ma.dot(w2, ERAal[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) - al = (albedo_post - albedo_pre) * (t_epoch - t_pre) + albedo_pre + pressure_hPa = np.array([1, 2, 3, 5, 7, 10, 20]) + + pressure = 100 * pressure_hPa + + 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) + + except IndexError: + print("Error: Please check time range of ERA5 data!") + elif np.abs(lat - start_lat) <= 10.0 and np.abs(lon - start_lon) <= 10.0: + try: + interp4d_temp_pre = np.ma.dot(w0, ERAtemp_ascend[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_temp_post = np.ma.dot(w0, ERAtemp_ascend[t_post_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_temp = (interp4d_temp_post - interp4d_temp_pre) * (t_epoch - t_pre) + interp4d_temp_pre + + interp4d_height_pre = np.ma.dot(w0, ERAz_ascend[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_height_post = np.ma.dot(w0, ERAz_ascend[t_post_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_height = (interp4d_height_post - interp4d_height_pre) * (t_epoch - t_pre) + interp4d_height_pre + + interp4d_vw_x_pre = np.ma.dot(w0, vw_x_ascend[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_x_post = np.ma.dot(w0, vw_x_ascend[t_post_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_x = (interp4d_vw_x_post - interp4d_vw_x_pre) * (t_epoch - t_pre) + interp4d_vw_x_pre + + interp4d_vw_y_pre = np.ma.dot(w0, vw_y_ascend[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_y_post = np.ma.dot(w0, vw_y_ascend[t_post_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_y = (interp4d_vw_y_post - interp4d_vw_y_pre) * (t_epoch - t_pre) + interp4d_vw_y_pre + + interp4d_vw_z_pre = np.ma.dot(w0, vw_z_ascend[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_z_post = np.ma.dot(w0, vw_z_ascend[t_post_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_z = (interp4d_vw_z_post - interp4d_vw_z_pre) * (t_epoch - t_pre) + interp4d_vw_z_pre + + 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 + + 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) + + except IndexError: + print("Error: Check time range of ERA5 data!") + else: + pass tcc_pre = np.ma.dot(w2, ERAtcc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) tcc_post = np.ma.dot(w2, ERAtcc[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) tcc = (tcc_post - tcc_pre) * (t_epoch - t_pre) + tcc_pre + if isinstance(tcc, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"tcc\"!") + print("Assuming simplified value for parameter \"tcc\".") + tcc = cc + cbh_pre = np.ma.dot(w2, ERAcbh[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) cbh_post = np.ma.dot(w2, ERAcbh[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) cbh = (cbh_post - cbh_pre) * (t_epoch - t_pre) + cbh_pre + if isinstance(tcc, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"cbh\"!") + print("Assuming simplified value for parameter \"cbh\".") + cbh = 2000 + lcc_pre = np.ma.dot(w2, ERAlcc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) lcc_post = np.ma.dot(w2, ERAlcc[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) lcc = (lcc_post - lcc_pre) * (t_epoch - t_pre) + lcc_pre + if isinstance(lcc, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"lcc\"!") + print("Assuming simplified value for parameter \"lcc\".") + lcc = cc/3 + mcc_pre = np.ma.dot(w2, ERAmcc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) mcc_post = np.ma.dot(w2, ERAmcc[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) mcc = (mcc_post - mcc_pre) * (t_epoch - t_pre) + mcc_pre + if isinstance(mcc, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"mcc\"!") + print("Assuming simplified value for parameter \"mcc\".") + mcc = cc/3 + hcc_pre = np.ma.dot(w2, ERAhcc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) hcc_post = np.ma.dot(w2, ERAhcc[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) hcc = (hcc_post - hcc_pre) * (t_epoch - t_pre) + hcc_pre + if isinstance(hcc, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"hcc\"!") + print("Assuming simplified value for parameter \"hcc\".") + hcc = cc/3 + ssr_pre = np.ma.dot(w2, ERAssr[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) ssr_post = np.ma.dot(w2, ERAssr[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) - ssr = ((ssr_post - ssr_pre) * (t_epoch - t_pre) + ssr_pre)/3600 + ssr = ((ssr_post - ssr_pre) * (t_epoch - t_pre) + ssr_pre) / 3600 strn_pre = np.ma.dot(w2, ERAstrn[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) strn_post = np.ma.dot(w2, ERAstrn[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) - strn = ((strn_post - strn_pre) * (t_epoch - t_pre) + strn_pre)/3600 + strn = ((strn_post - strn_pre) * (t_epoch - t_pre) + strn_pre) / 3600 + + if isinstance(strn, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"strn\"!") + print("Assuming simplified value for parameter \"strn\".") + strn = 0 + + skt_pre = np.ma.dot(w2, ERAskt[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) + skt_post = np.ma.dot(w2, ERAskt[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) + skt = ((skt_post - skt_pre) * (t_epoch - t_pre) + skt_pre) + + if isinstance(skt, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"skt\"!") + print("Assuming simplified value for parameter \"skt\".") + skt = T_ground strd_pre = np.ma.dot(w2, ERAstrd[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) strd_post = np.ma.dot(w2, ERAstrd[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) - strd = ((strd_post - strd_pre) * (t_epoch - t_pre) + strd_pre)/3600 + strd = ((strd_post - strd_pre) * (t_epoch - t_pre) + strd_pre) / 3600 + + if isinstance(strd, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"strd\"!") + print("Assuming simplified value for parameter \"strd\".") + strd = epsilon_ground * sigma * T_ground ** 4 strdc_pre = np.ma.dot(w2, ERAstrdc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) strdc_post = np.ma.dot(w2, ERAstrdc[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) strdc = ((strdc_post - strdc_pre) * (t_epoch - t_pre) + strdc_pre) / 3600 + if isinstance(strdc, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"strdc\"!") + print("Assuming simplified value for parameter \"strdc\".") + strdc = epsilon_ground * sigma * T_ground ** 4 + ssrd_pre = np.ma.dot(w2, ERAssrd[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) ssrd_post = np.ma.dot(w2, ERAssrd[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) - ssrd = ((ssrd_post - ssrd_pre) * (t_epoch - t_pre) + ssrd_pre)/3600 + ssrd = ((ssrd_post - ssrd_pre) * (t_epoch - t_pre) + ssrd_pre) / 3600 + + if isinstance(ssrd, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"ssrd\"!") + print("Assuming simplified value for parameter \"ssrd\".") + ssrd = 1 + ssr = 1 - Albedo + + if isinstance(ssr, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"ssr\"!") + print("Assuming simplified value for parameter \"ssr\".") + ssrd = 1 + ssr = 1 - Albedo tsr_pre = np.ma.dot(w2, ERAtsr[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) tsr_post = np.ma.dot(w2, ERAtsr[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) - tsr = ((tsr_post - tsr_pre) * (t_epoch - t_pre) + tsr_pre)/3600 + tsr = ((tsr_post - tsr_pre) * (t_epoch - t_pre) + tsr_pre) / 3600 tisr_pre = np.ma.dot(w2, ERAtisr[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) tisr_post = np.ma.dot(w2, ERAtisr[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) - tisr = ((tisr_post - tisr_pre) * (t_epoch - t_pre) + tisr_pre)/3600 + tisr = ((tisr_post - tisr_pre) * (t_epoch - t_pre) + tisr_pre) / 3600 + + if isinstance(tisr, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"tisr\"!") + print("Assuming simplified value for parameter \"tisr\".") + AZ, ELV = sun_angles_analytical(lat, lon, h, utc) + 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 + tisr = I_Sun * np.sin(np.deg2rad(ELV)) + + if isinstance(tsr, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"tsr\"!") + print("Assuming simplified value for parameter \"tsr\".") + tsr = (1 - Albedo) * tisr ttr_pre = np.ma.dot(w2, ERAttr[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) ttr_post = np.ma.dot(w2, ERAttr[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) - ttr = ((ttr_post - ttr_pre) * (t_epoch - t_pre) + ttr_pre)/3600 + ttr = ((ttr_post - ttr_pre) * (t_epoch - t_pre) + ttr_pre) / 3600 - 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]) + p0_pre = np.ma.dot(w2, ERAsp[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) + p0_post = np.ma.dot(w2, ERAsp[t_post_ind, lat_ind2, lon_ind2]) / np.sum(w2) + p0 = (p0_post - p0_pre) * (t_epoch - t_pre) + p0_pre - pressure = 100 * pressure_hPa + if isinstance(p0, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"sp\"!") + print("Assuming simplified value for parameter \"sp\".") + p0 = 101325.0 - 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) + if isinstance(ttr, float) != True: + print("WARNING: Corrupt ERA5 Data for parameter \"ttr\"!") + print("Assuming simplified value for parameter \"ttr\".") + AZ, ELV = sun_angles_analytical(lat, lon, h, utc) + tau_atm, tau_atmIR = tau(ELV, h, p_air, p0) + HalfConeAngle = np.arcsin(R_E / (R_E + h)) + ViewFactor = (1 - np.cos(HalfConeAngle)) / 2 + ttr = epsilon_ground * sigma * T_ground ** 4 * tau_atmIR * ViewFactor * 2 if h > interp4d_height[0]: p_air = press_interp1d(interp4d_height[0]) T_air = temp_interp1d(interp4d_height[0]) u = vw_x_interp1d(interp4d_height[0]) v = vw_y_interp1d(interp4d_height[0]) - w = -1 / g * vw_z_interp1d(interp4d_height[0]) * R_air * T_air / p_air + w = -1 / grav(lat, h) * vw_z_interp1d(interp4d_height[0]) * R_air * T_air / p_air elif h < interp4d_height[-1]: p_air = press_interp1d(interp4d_height[-1]) T_air = temp_interp1d(interp4d_height[-1]) u = vw_x_interp1d(interp4d_height[-1]) v = vw_y_interp1d(interp4d_height[-1]) - w = -1 / g * vw_z_interp1d(interp4d_height[-1]) * R_air * T_air / p_air + w = -1 / grav(lat, h) * vw_z_interp1d(interp4d_height[-1]) * R_air * T_air / p_air else: p_air = press_interp1d(h) T_air = temp_interp1d(h) u = vw_x_interp1d(h) v = vw_y_interp1d(h) - w = -1 / g * vw_z_interp1d(h) * R_air * T_air / p_air + w = -1 / grav(lat, h) * vw_z_interp1d(h) * R_air * T_air / p_air rho_air = p_air / (R_air * T_air) - return p_air, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, ssrd, tsr, ttr, al, tisr, strdc -# p_air, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, ssrd, tsr, ttr, al, tisr + return p_air, p0, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, strdc, ssrd, tsr, ttr, tisr, skt deltaT_ERA = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 -p_air0, T_air0, rho_air0, u0, v0, w0, cbh0, tcc0, lcc0, mcc0, hcc0, ssr0, strn0, strd0, ssrd0, tsr0, ttr0, al0, tisr0, strdc0 = ERA5Data(start_lon, start_lat, start_height, 0, deltaT_ERA) +p_air0, p00, T_air0, rho_air0, u0, v0, w0, cbh0, tcc0, lcc0, mcc0, hcc0, ssr0, strn0, strd0, strdc0, ssrd0, tsr0, ttr0, tisr0, skt0 = ERA5Data( + start_lon, start_lat, start_height, 0, deltaT_ERA) + +y0 = [ + start_lon, # start longitude [deg] + start_lat, # start latitude [deg] + start_height, # start altitude [m] + 0, # initial v_x [m/s] + 0, # initial v_y [m/s] + 0, # initial v_z [m/s] + T_air0, # initial gas temperature [K] = initial air temperature [K] + T_air0, # initial film temperature [K] = initial air temperature [K] + m_gas_init, # initial lifting gas mass [kg] + 0, # initial factor c2 [-] + m_bal_init # initial ballast mass [kg] +] - - - -y0 = [start_lon, start_lat, start_height, 0, 0, 0, T_air0, T_air0, m_gas_init, 0] -## y = [longitude [deg], latitude [deg], height [m], v_x (lon) [m/s], v_y [m/s], v_z [m/s], T_gas [K], T_film [K], m_gas_init [kg], c2_init [-]] - -valving = False - -def model(t, y, m_pl, m_film, m_bal, c_virt): +def model(t, y, m_pl, m_film, c_virt): utc = Time(start_utc) + t * unit.second lon = y[0] # 1 lat = y[1] # 2 @@ -251,12 +416,36 @@ def model(t, y, m_pl, m_film, m_bal, c_virt): T_film = y[7] # 8 m_gas = y[8] # 9 c2 = y[9] # 10 + m_bal = y[10] # 11 - r_lon, r_lat = radii(lat, h) + if (lon % 360) > 180: # convert longitude to value in standard interval [-180, 180] + lon = (lon % 360) - 360 + else: + lon = (lon % 360) + + if lat > 90: # set plausible limits for latitudes + lat = 90 + elif lat < -90: + lat = -90 + else: + lat = lat + + r_lon, r_lat = radii(lat, h) # calculate radii for velocity conversion between cartesian and Earth reference frame deltaT_ERA = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 # conversion to ERA5 time format - p_air, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, ssrd, tsr, ttr, al, tisr, strdc = ERA5Data(lon, lat, h, t, deltaT_ERA) + try: + p_air, p0, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, strdc, ssrd, tsr, ttr, tisr, skt = ERA5Data(lon, lat, h, t, deltaT_ERA) + except: + # in case of solver (temporarily) exceeding interpolation area (with subsequent correction) + # or permanent drift out of interpolation area + if h >= 30000 or (np.abs(lat - start_lat) <= 10.0 and np.abs(lon - start_lon) <= 10.0): + print("solver exceeds definition area") + p_air, p0, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, strdc, ssrd, tsr, ttr, tisr, skt = 0, 101325, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 + else: + print("instable trajectory!") + p_air, p0, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, strdc, ssrd, tsr, ttr, tisr, skt = 0, 101325, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 + p_gas = p_air h_valve = 1.034 * V_design ** (1 / 3) @@ -272,36 +461,53 @@ def model(t, y, m_pl, m_film, m_bal, c_virt): rho_gas = p_gas / (R_gas * T_gas) # calculate gas density through *ideal* gas equation - dP_valve = g * (rho_air - rho_gas) * h_valve - dP_duct = g * (rho_air - rho_gas) * h_duct + dP_valve = grav(lat, h) * (rho_air - rho_gas) * h_valve + dP_duct = grav(lat, h) * (rho_air - rho_gas) * h_duct + + if m_gas > m_gas_init: # limit gas mass to plausible value + m_gas = m_gas_init # + elif m_gas < 0: + m_gas = 0 + else: + pass V_b = m_gas / rho_gas # calculate balloon volume from current gas mass and gas density if V_b > V_design: c_duct = c_ducts + elif V_b < 0: + c_duct = 0 + V_b = 0 else: c_duct = 0 - #if valving == True: # opening valve process - # if c2 == 0: - # c2 = 1.0 - # c2dot = 0 - # elif c_valve < c2 <= 1.0: - # c2dot = (c_valve - 1.0) / t_open - # else: - # c2dot = 0 - # c2 = c_valve + if ballasting(utc) == True: + if m_bal >= 0: + mdot = m_baldot + else: + mdot = 0 + else: + mdot = 0 - #if valving == False: # closing valve process - # if c2 == 0: - # c2dot = 0 - # elif c_valve <= c2 < 1.0: - # c2dot = (1.0 - c_valve) / t_close - # else: - # c2dot = 0 - # c2 = 0 + if valving(utc) == True: # opening valve process + if c2 == 0: + c2 = 1.0 + c2dot = 0 + elif c_valve < c2 <= 1.0: + c2dot = (c_valve - 1.0) / t_open + else: + c2dot = 0 + c2 = c_valve + + if valving(utc) == False: # closing valve process + if c2 == 0: + c2dot = 0 + elif c_valve <= c2 < 1.0: + c2dot = (1.0 - c_valve) / t_close + else: + c2dot = 0 + c2 = 0 - ## m_gas/dt = -(A_valv * c2 * np.sqrt(2 * dP_valve * rho_gas)) # mass loss due to valving m_gross = m_pl + m_film + m_bal m_tot = m_pl + m_film + m_gas m_virt = m_tot + c_virt * rho_air * V_b @@ -313,195 +519,260 @@ def model(t, y, m_pl, m_film, m_bal, c_virt): A_eff = 0.65 * A_surf + 0.35 * A_surf1 A_top = np.pi / 4 * d_b ** 2 - AZ, ELV = sun_angles_analytical(lat, lon, utc) + AZ, ELV = sun_angles_analytical(lat, lon, h, utc) A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV))) # projected area for sun radiation - A_drag = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * alpha)) # projected area for drag + A_drag = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * alpha)) # projected area for drag # CALCULATIONS FOR THERMAL MODEL - tau_atm, tau_atmIR = tau(ELV, h, p_air) - tau_atm0, tau_atmIR0 = tau(ELV, 0, p_0) + tau_atm, tau_atmIR = tau(ELV, h, p_air, p0) + + tau_atm0, tau_atmIR0 = tau(ELV, 0, p0, p0) 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 I_Sun0 = I_Sun * tau_atm0 - ## q_sun = I_SunZ - ## q_IRground = np.abs(str) #epsilon_ground * sigma * T_ground ** 4 - ## q_IREarth = q_IRground * tau_atmIR - tcc = 0 + HalfConeAngle = np.arcsin(R_E / (R_E + h)) + ViewFactor = (1 - np.cos(HalfConeAngle)) / 2 + + if simple == True: + + q_IREarth = epsilon_ground * sigma * T_ground ** 4 * tau_atmIR - if tcc <= 0.01: - q_IRground = strd - strn - q_IREarth = q_IRground * tau_atmIR - q_sun = I_SunZ if ELV <= 0: q_Albedo = 0 else: - q_Albedo = (1 - ssr/ssrd) * I_Sun0 * np.sin(np.deg2rad(ELV)) # ssrd - ssr + q_Albedo = Albedo * I_Sun * np.sin(np.deg2rad(ELV)) + + 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_sun = I_SunZ + else: - q_IRground_bc = (strd - strn) + (strd - strdc) - q_IREarth_bc = q_IRground_bc * tau_atmIR - q_sun_bc = I_SunZ * (1 - tcc) - q_Albedo_bc = (1 - ssr/ssrd) * (1 - tcc) * I_Sun * np.sin(np.deg2rad(ELV)) - q_IREarth_ac = np.abs(ttr) #q_IRground_bc * tau_atmIR * (1 - tcc) - q_sun_ac = I_SunZ - q_Albedo_ac = (tisr - tsr)/tisr * I_Sun * np.sin(np.deg2rad(ELV)) + if tcc <= 0.01: + + q_IREarth1 = alpha_IR * A_surf * (strd - strn) * tau_atmIR * ViewFactor * (1 + tau_VIS / (1 - r_VIS)) + q_IREarth2 = alpha_IR * A_surf * np.abs(ttr) * 0.5 * (1 + tau_IR / (1 - r_IR)) + # q_IREarth2 = alpha_IR * A_surf * (0.04321906 * np.abs(ttr) + 84.67820281) * (1 + tau_IR / (1 - r_IR)) + + if h > 40000: + Q_IREarth = q_IREarth2 + else: + Q_IREarth = (q_IREarth2 - q_IREarth1) / 40000 * h + q_IREarth1 - if h <= cbh: # "below clouds" - q_IREarth = q_IREarth_bc - q_sun = q_sun_bc if ELV <= 0: - q_Albedo = 0 + q_Albedo1 = 0 + q_Albedo2 = 0 else: - q_Albedo = q_Albedo_bc - elif h >= 12000: # "above clouds" - q_IREarth = q_IREarth_ac - q_sun = q_sun_ac - if ELV <= 0: - q_Albedo = 0 - else: - q_Albedo = q_Albedo_ac - elif h >= 6000: - if hcc >= 0.01: - q_IREarth = ((h - 6000) / 6000 * hcc + mcc + lcc) / (hcc + mcc + lcc) * (q_IREarth_ac - q_IREarth_bc) + q_IREarth_bc - q_sun = ((h - 6000) / 6000 * hcc + mcc + lcc) / (hcc + mcc + lcc) * (q_sun_ac - q_sun_bc) + q_sun_bc - if ELV <= 0: - q_Albedo = 0 + if ssrd == 0: + q_Albedo1 = 0 else: - q_Albedo = ((h - 6000) / 6000 * hcc + mcc + lcc) / (hcc + mcc + lcc) * (q_Albedo_ac - q_Albedo_bc) + q_Albedo_bc - else: - q_IREarth = q_IREarth_ac - q_sun = q_sun_ac - if ELV <= 0: - q_Albedo = 0 + q_Albedo1 = alpha_VIS * A_surf * (1 - ssr / ssrd) * I_Sun0 * tau_atmIR * np.sin( + np.deg2rad(ELV)) * ViewFactor * (1 + tau_VIS / (1 - r_VIS)) # ! + if tisr == 0: + q_Albedo2 = 0 else: - q_Albedo = q_Albedo_ac + q_Albedo2 = alpha_VIS * A_surf * (1 - tsr / tisr) * I_Sun * np.sin(np.deg2rad(ELV)) * 0.5 * ( + 1 + tau_VIS / (1 - r_VIS)) # ! + if h > 40000: + Q_Albedo = q_Albedo2 + else: + Q_Albedo = (q_Albedo2 - q_Albedo1) / 40000 * h + q_Albedo1 + + q_sun = I_SunZ - elif h >= 2000: - if mcc > 0.01 or hcc > 0.01: - q_IREarth = ((h - 2000)/4000 * mcc + lcc)/(hcc + mcc + lcc) * (q_IREarth_ac - q_IREarth_bc) + q_IREarth_bc - q_sun = ((h - 2000)/4000 * mcc + lcc)/(hcc + mcc + lcc) * (q_sun_ac - q_sun_bc) + q_sun_bc - if ELV <= 0: - q_Albedo = 0 - else: - q_Albedo = ((h - 2000)/4000 * mcc + lcc)/(hcc + mcc + lcc) * (q_Albedo_ac - q_Albedo_bc) + q_Albedo_bc - else: - q_IREarth = q_IREarth_ac - q_sun = q_sun_ac - if ELV <= 0: - q_Albedo = 0 - else: - q_Albedo = q_Albedo_ac else: - q_IREarth = (h/2000 * lcc)/(hcc + mcc + lcc) * (q_IREarth_ac - q_IREarth_bc) + q_IREarth_bc - q_sun = (h/2000 * lcc)/(hcc + mcc + lcc) * (q_sun_ac - q_sun_bc) + q_sun_bc + q_IRground_bc = (strd - strn) + (strd - strdc) + q_IREarth_bc = alpha_IR * A_surf * q_IRground_bc * tau_atmIR * ViewFactor * (1 + tau_VIS / (1 - r_VIS)) + q_sun_bc = I_SunZ * (1 - tcc) + q_Albedo_bc = alpha_VIS * A_surf * (1 - ssr / ssrd) * I_Sun0 * tau_atmIR * np.sin( + np.deg2rad(ELV)) * ViewFactor * (1 + tau_VIS / (1 - r_VIS)) - if ELV <= 0: - q_Albedo = 0 + # q_IREarth_ac = alpha_IR * A_surf * (0.04321906 * np.abs(ttr) + 84.67820281) * (1 + tau_IR / (1 - r_IR)) + q_IREarth_ac = alpha_IR * A_surf * np.abs(ttr) * 0.5 * (1 + tau_VIS / (1 - r_VIS)) + q_sun_ac = I_SunZ + q_Albedo_ac = alpha_VIS * A_surf * (1 - tsr / tisr) * I_Sun * np.sin(np.deg2rad(ELV)) * ViewFactor * ( + 1 + tau_VIS / (1 - r_VIS)) + + if h <= cbh: # "below clouds" + Q_IREarth = q_IREarth_bc + q_sun = q_sun_bc + if ELV <= 0: + Q_Albedo = 0 + else: + Q_Albedo = q_Albedo_bc + elif h >= 12000: # "above clouds" + Q_IREarth = q_IREarth_ac + q_sun = q_sun_ac + if ELV <= 0: + Q_Albedo = 0 + else: + Q_Albedo = q_Albedo_ac + elif h >= 6000: + if hcc >= 0.01: + Q_IREarth = ((h - 6000) / 6000 * hcc + mcc + lcc) / (hcc + mcc + lcc) * ( + q_IREarth_ac - q_IREarth_bc) + q_IREarth_bc + q_sun = ((h - 6000) / 6000 * hcc + mcc + lcc) / (hcc + mcc + lcc) * (q_sun_ac - q_sun_bc) + q_sun_bc + if ELV <= 0: + Q_Albedo = 0 + else: + Q_Albedo = ((h - 6000) / 6000 * hcc + mcc + lcc) / (hcc + mcc + lcc) * ( + q_Albedo_ac - q_Albedo_bc) + q_Albedo_bc + else: + Q_IREarth = q_IREarth_ac + q_sun = q_sun_ac + if ELV <= 0: + Q_Albedo = 0 + else: + Q_Albedo = q_Albedo_ac + elif h >= 2000: + if mcc > 0.01 or hcc > 0.01: + Q_IREarth = ((h - 2000) / 4000 * mcc + lcc) / (hcc + mcc + lcc) * ( + q_IREarth_ac - q_IREarth_bc) + q_IREarth_bc + q_sun = ((h - 2000) / 4000 * mcc + lcc) / (hcc + mcc + lcc) * (q_sun_ac - q_sun_bc) + q_sun_bc + if ELV <= 0: + Q_Albedo = 0 + else: + Q_Albedo = ((h - 2000) / 4000 * mcc + lcc) / (hcc + mcc + lcc) * ( + q_Albedo_ac - q_Albedo_bc) + q_Albedo_bc + else: + Q_IREarth = q_IREarth_ac + q_sun = q_sun_ac + if ELV <= 0: + Q_Albedo = 0 + else: + Q_Albedo = q_Albedo_ac else: - q_Albedo = (h/2000 * lcc)/(hcc + mcc + lcc) * (q_Albedo_ac - q_Albedo_bc) + q_Albedo_bc + Q_IREarth = (h / 2000 * lcc) / (hcc + mcc + lcc) * (q_IREarth_ac - q_IREarth_bc) + q_IREarth_bc + q_sun = (h / 2000 * lcc) / (hcc + mcc + lcc) * (q_sun_ac - q_sun_bc) + q_sun_bc + if ELV <= 0: + Q_Albedo = 0 + else: + Q_Albedo = (h / 2000 * lcc) / (hcc + mcc + lcc) * (q_Albedo_ac - q_Albedo_bc) + q_Albedo_bc my_air = (1.458 * 10 ** -6 * T_air ** 1.5) / (T_air + 110.4) 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) + Gr_air = (rho_air ** 2 * grav(lat, h) * 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_rel) * d_b * rho_air / my_air - Fr = np.abs(v_rel) / np.sqrt(g * d_b) + Fr = np.abs(v_rel) / np.sqrt(grav(lat, h) * d_b) 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)) ** ( + (rho_gas ** 2 * grav(lat, h) * 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_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) - c_d = Cd_model(Fr, Re, A_top) + if simple == True: + c_d = 0.47 # sphere + else: + c_d = cd_PalumboHigh(Fr, Re, A_top) # Ref. xx + + #c_d = 0.47 # cd_PalumboHigh(Fr, Re, A_top) #cd_sphere(Re) #0.47 # cd_Palumbo(Fr, Re, A_top) # 0.8 #cd_sphere(Re) # cd_palumbo(Fr, Re, A_top) # cd_sphere(Re) / 0.8 / 0.47 D = drag(c_d, rho_air, A_drag, v_rel) # calculate drag force if v_rel == 0: - Drag_x, Dray_y, Drag_z = 0, 0, 0 + Drag_x, Drag_y, Drag_z = 0, 0, 0 else: Drag_x, Drag_y, Drag_z = D * v_relx / v_rel, D * v_rely / v_rel, D * v_relz / v_rel - F = g * V_b * (rho_air - rho_gas) - g * m_gross + Drag_z # gross inflation - weight + drag + F = grav(lat, h) * V_b * (rho_air - rho_gas) - grav(lat, h) * m_gross + Drag_z # gross inflation - weight + drag a_x, a_y, a_z = Drag_x / m_virt, Drag_y / m_virt, F / m_virt eqn1 = np.rad2deg(y[3] / r_lon) eqn2 = np.rad2deg(y[4] / r_lat) - eqn3 = y[5] + eqn3 = v_z eqn4 = a_x eqn5 = a_y eqn6 = a_z - eqn7 = Q_ConvInt/(gamma * c_v * m_gas) - (gamma - 1)/gamma * (rho_air * g)/(rho_gas * R_gas) * y[5] + eqn7 = Q_ConvInt / (gamma * c_v * m_gas) - (gamma - 1) / gamma * (rho_air * grav(lat, h)) / (rho_gas * R_gas) * v_z eqn8 = (Q_Sun + Q_Albedo + Q_IREarth + Q_IRFilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film) - eqn9 = -(A_ducts * c_duct * np.sqrt(2 * dP_duct * rho_gas)) # - (A_valve * c2 * np.sqrt(2 * dP_valve * rho_gas)) - eqn10 = 0 #c2dot + eqn9 = -(A_ducts * c_duct * np.sqrt(np.abs(2 * dP_duct * rho_gas))) - (A_valve * c2 * np.sqrt(np.abs(2 * dP_valve * rho_gas))) + eqn10 = c2dot + if m_bal > 0: + eqn11 = -mdot + else: + eqn11 = 0 - return [eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10] + return [eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10, eqn11] -m_pl = 1728 + 172 + 500 # + 500 # mass of payload and flight chain in [kg] -m_bal = 0 # initial ballast mass in [kg] -m_film = 1897 # mass of balloon film in [kg] -FreeLift = 4 # [%] -t_max = 100000 #0 -dt = 1.0 -t = np.arange(2.0, t_max, dt) +# DEFINITION OF EVENTS FOR SOLVER -m_gas = ((m_pl + m_bal + m_film) * (FreeLift/100 + 1))/(R_gas/R_air - 1) +def at_ground(t, y, m_pl, m_film, c_virt): + return y[2] -def ending(t, y): return y[2] -sol = solve_ivp(fun=model, t_span=[2.0, t_max], y0=y0, t_eval=t, args=(m_pl, m_film, m_bal, c_virt), max_step=35.0) #, events=ending.terminal) # events= -print("hier:") +def above_float(t, y, m_pl, m_film, c_virt): + return 45000 - y[2] + + +def below_float(t, y, m_pl, m_film, c_virt): + return y[2] - 30050 + + +hit_ground = lambda t, x: at_ground(t, x, m_pl, m_film, c_virt) +hit_ground.terminal = True +hit_ground.direction = -1 +excess_ascend = lambda t, x: above_float(t, x, m_pl, m_film, c_virt) +excess_ascend.terminal = True +excess_ascend.direction = -1 +instable = lambda t, x: below_float(t, x, m_pl, m_film, c_virt) +instable.terminal = True +instable.direction = -1 + +t0 = 2.0 +tf = 30000 + +print("Beginning simulation") + +sol = solve_ivp(fun=lambda t, x: model(t, x, m_pl, m_film, c_virt), t_span=[t0, tf], y0=y0, method='BDF', events=[hit_ground, excess_ascend, instable]) + +# tnew = np.linspace(0, sol.t[-1], len(Vol_list)) + print(sol.message) -res = sol.y[2, :] +print(datetime.now() - starttime) -end = datetime.now() -print(end-start) -plt.plot(comp_time, comp_height, 'r--', label='PoGo Flight Test') -plt.plot(t, res) +plt.plot(sol.t, sol.y[2, :], 'k--', label='Simulation') +plt.plot(comp_time, comp_height, 'r-', label='PoGo Flight Test') +plt.legend() +plt.title('high factor') plt.xlabel('time in s') plt.ylabel('Balloon Altitude in m') - -plt.xlim(0, t_max) plt.show() plt.clf() -ax = plt.axes(projection=ccrs.Mollweide()) +# ax = plt.axes(projection=ccrs.Mollweide()) +ax = plt.axes(projection=ccrs.AzimuthalEquidistant(central_latitude=90)) ax.coastlines() ax.gridlines() ax.stock_img() ax.set_extent([-120, 30, 60, 80], crs=ccrs.PlateCarree()) - plt.plot(start_lon, start_lat, 'rx', transform=ccrs.Geodetic()) plt.plot(sol.y[0, :], sol.y[1, :], 'k--', transform=ccrs.Geodetic()) plt.plot(comp_lon, comp_lat, 'r-', transform=ccrs.Geodetic()) # plt.xticks() # figname = "LatLon_%.2f_%.2f.png" % (Albedo, epsilon_ground) # plt.savefig(os.path.join(rootdir, figname)) -plt.show() \ No newline at end of file +plt.show()