diff --git a/main.py b/main.py index 861f20b..c66ed43 100644 --- a/main.py +++ b/main.py @@ -1,13 +1,25 @@ import sys +import pickle +from pyfiglet import Figlet import warnings import numpy as np import xarray as xr import pandas as pd +import pickle as pkl import cartopy.crs as ccrs import astropy.units as unit import matplotlib.pyplot as plt from dask import delayed from datetime import datetime +starttime = datetime.now() +print('----------------------------------------') +ascii_banner = Figlet(font="slant") +print(ascii_banner.renderText("BASTET")) +print("Ver. 1.0, 2021 by Marcel Frommelt") +print('----------------------------------------') +print("") +print("") + from netCDF4 import Dataset from astropy.time import Time from scipy import interpolate @@ -18,16 +30,19 @@ from input.natural_constants import * from models.gravity import grav from models.valving import valving from models.ballasting import ballasting +from dask.diagnostics import ProgressBar +from models.simple_atmosphere import T_air_simple, p_air_simple, rho_air_simple 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 +from models.drag import drag, cd_PalumboLow, cd_Palumbo, cd_PalumboHigh, cd_PalumboMC, cd_sphere +from models.transformation import visible_cells, transform, radii, transform2 +from multiprocessing import Process 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 +data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Data_PoGo2016.xls', sheet_name='SuperTIGER2') # Tabelle3 comp_time = pd.DataFrame(data, columns=['Time']).to_numpy().squeeze() comp_height = pd.DataFrame(data, columns=['Height']).to_numpy().squeeze() @@ -35,43 +50,31 @@ comp_lat = pd.DataFrame(data, columns=['Latitude']).to_numpy().squeeze() comp_lon = pd.DataFrame(data, columns=['Longitude']).to_numpy().squeeze() print("") -print("Initialising simulation...") -print("Launch location : longitude %.4f, latitude: %.4f" % (start_lon, start_lat)) +print("INITIALISING SIMULATION...") +print("") +print("Launch location:") +print("longitude: %.4f deg" % (start_lon)) +print("latitude: %.4f deg" % (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") +first_file = Dataset(ERA5_float[0]) +last_file = Dataset(ERA5_float[-1]) -df = xr.open_mfdataset(['float1_2019.nc', 'float2_2019.nc', 'float3_2019.nc', 'float4_2019.nc'], combine='by_coords', concat_dim="time", parallel=True) +tstart = int(first_file.variables['time'][0]) +tend = int(last_file.variables['time'][-1]) -startNC = Dataset('float1_2019.nc') -endNC = Dataset('float4_2019.nc') -start = int(startNC.variables['time'][0]) -end = int(endNC.variables['time'][-1]) +first_file.close() +last_file.close() -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 - - -# ERA5 MULTI-LEVEL ASCENT (WIND + ATMOSPHERIC DATA DURING ASCENT) - -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] +df1 = xr.open_mfdataset(ERA5_float, combine='by_coords', engine='netcdf4', concat_dim="time", parallel=True) +float_data = df1.assign_coords(time=np.linspace(tstart, tend, (tend - tstart) + 1)) # ERA5 MULTI-LEVEL FLOAT (WIND + ATMOSPHERIC DATA DURING FLOAT) +ERAtime = float_data.variables['time'][:] # time 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] @@ -81,6 +84,18 @@ 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] +first_file = Dataset(ERA5_single[0]) +last_file = Dataset(ERA5_single[-1]) + +tstart = int(first_file.variables['time'][0]) +tend = int(last_file.variables['time'][-1]) + +first_file.close() +last_file.close() + +df2 = xr.open_mfdataset(ERA5_single, combine='by_coords', engine='netcdf4', concat_dim="time", parallel=True) +single_data = df2.assign_coords(time=np.linspace(tstart, tend, (tend - tstart) + 1)) + # ERA5 SINGLE-LEVEL (RADIATIVE ENVIRONMENT) @@ -102,6 +117,32 @@ ERAtisr = single_data.variables['tisr'][:].values # hourly accumulated TOA 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] +first_file = Dataset(ERA5_ascent[0]) +last_file = Dataset(ERA5_ascent[-1]) + +tstart = int(first_file.variables['time'][0]) +tend = int(last_file.variables['time'][-1]) + +first_file.close() +last_file.close() + +df3 = xr.open_mfdataset(ERA5_ascent, combine='by_coords', engine='netcdf4', concat_dim="time", parallel=True) +ascent_data = df3.assign_coords(time=np.linspace(tstart, tend, (tend - tstart) + 1)) + + +# ERA5 MULTI-LEVEL ASCENT (WIND + ATMOSPHERIC DATA DURING ASCENT) + +ERAlat0 = ascent_data.variables['latitude'][:].values # latitude [deg] +ERAlon0 = ascent_data.variables['longitude'][:].values # longitude [deg] +ERAz_ascent = ascent_data.variables['z'][:].values / g # geopotential [m^-2/s^-2] to geopotential height [m] +ERApress_ascent = ascent_data.variables['level'][:].values # pressure level [-] +ERAtemp_ascent = ascent_data.variables['t'][:].values # air temperature in K +vw_x_ascent = ascent_data.variables['u'][:].values # v_x in [m/s] +vw_y_ascent = ascent_data.variables['v'][:].values # v_y in [m/s] +vw_z_ascent = ascent_data.variables['w'][:].values # v_z in [m/s] + +ascent_data.close() + print("Finished reading ERA5-datasets.") lon_era2d0, lat_era2d0 = np.meshgrid(ERAlon0, ERAlat0) @@ -116,10 +157,14 @@ 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") +print("Built kd-trees.") +print("") +wflag1, wflag2, wflag3, wflag4, wflag5, wflag6, wflag7, wflag8, wflag9, wflag10, wflag11, wflag12, wflag13, wflag14, wflag15, wflag16, wflag17, wflag18, wflag19 = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 -def ERA5Data(lon, lat, h, t, deltaT_ERA): +flag_arr = np.zeros(20) + +def ERA5Data(lon, lat, h, t, deltaT_ERA, flag_arr): t_epoch = deltaT_ERA + t / 3600 t_pre = int(t_epoch) @@ -165,7 +210,7 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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 - pressure_hPa = np.array([1, 2, 3, 5, 7, 10, 20]) + pressure_hPa = np.array([1, 2, 3, 5, 7, 10, 20, 30]) # !!! pressure = 100 * pressure_hPa @@ -176,27 +221,32 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z) except IndexError: - print("Error: Please check time range of ERA5 data!") + if flag_arr[18] == 0: + print("Error: Please check time range of ERA5 data!") + flag_arr[18] = 1 + else: + flag_arr[18] = 1 + 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_pre = np.ma.dot(w0, ERAtemp_ascent[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_temp_post = np.ma.dot(w0, ERAtemp_ascent[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_pre = np.ma.dot(w0, ERAz_ascent[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_height_post = np.ma.dot(w0, ERAz_ascent[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_pre = np.ma.dot(w0, vw_x_ascent[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_x_post = np.ma.dot(w0, vw_x_ascent[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_pre = np.ma.dot(w0, vw_y_ascent[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_y_post = np.ma.dot(w0, vw_y_ascent[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_pre = np.ma.dot(w0, vw_z_ascent[t_pre_ind, :, lat_ind0, lon_ind0]) / np.sum(w0) + interp4d_vw_z_post = np.ma.dot(w0, vw_z_ascent[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( @@ -212,7 +262,12 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z) except IndexError: - print("Error: Check time range of ERA5 data!") + if flag_arr[19] == 0: + print("Error: Check time range of ERA5 data!") + flag_arr[19] = 1 + else: + flag_arr[19] = 1 + else: pass @@ -221,8 +276,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[1] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"tcc\"!") + print("Assuming simplified value for parameter \"tcc\".") + flag_arr[1] = 1 + else: + flag_arr[1] = 1 + tcc = cc cbh_pre = np.ma.dot(w2, ERAcbh[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -230,8 +290,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[2] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"cbh\"!") + print("Assuming simplified value for parameter \"cbh\".") + flag_arr[2] = 1 + else: + flag_arr[2] = 1 + cbh = 2000 lcc_pre = np.ma.dot(w2, ERAlcc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -239,8 +304,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[3] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"lcc\"!") + print("Assuming simplified value for parameter \"lcc\".") + flag_arr[3] = 1 + else: + flag_arr[3] = 1 + lcc = cc/3 mcc_pre = np.ma.dot(w2, ERAmcc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -248,8 +318,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[4] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"mcc\"!") + print("Assuming simplified value for parameter \"mcc\".") + flag_arr[4] = 1 + else: + flag_arr[4] = 1 + mcc = cc/3 hcc_pre = np.ma.dot(w2, ERAhcc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -257,8 +332,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[5] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"hcc\"!") + print("Assuming simplified value for parameter \"hcc\".") + flag_arr[5] = 1 + else: + flag_arr[5] = 1 + hcc = cc/3 ssr_pre = np.ma.dot(w2, ERAssr[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -270,8 +350,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[6] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"strn\"!") + print("Assuming simplified value for parameter \"strn\".") + flag_arr[6] = 1 + else: + flag_arr[6] = 1 + strn = 0 skt_pre = np.ma.dot(w2, ERAskt[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -279,8 +364,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[7] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"skt\"!") + print("Assuming simplified value for parameter \"skt\".") + flag_arr[7] = 1 + else: + flag_arr[7] = 1 + skt = T_ground strd_pre = np.ma.dot(w2, ERAstrd[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -288,8 +378,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[8] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"strd\"!") + print("Assuming simplified value for parameter \"strd\".") + flag_arr[8] = 1 + else: + flag_arr[8] = 1 + strd = epsilon_ground * sigma * T_ground ** 4 strdc_pre = np.ma.dot(w2, ERAstrdc[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -297,8 +392,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[9] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"strdc\"!") + print("Assuming simplified value for parameter \"strdc\".") + flag_arr[9] = 1 + else: + flag_arr[9] = 1 + strdc = epsilon_ground * sigma * T_ground ** 4 ssrd_pre = np.ma.dot(w2, ERAssrd[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -306,14 +406,24 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[10] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"ssrd\"!") + print("Assuming simplified value for parameter \"ssrd\".") + flag_arr[10] = 1 + else: + flag_arr[10] = 1 + 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\".") + if flag_arr[11] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"ssr\"!") + print("Assuming simplified value for parameter \"ssr\".") + flag_arr[11] = 1 + else: + flag_arr[11] = 1 + ssrd = 1 ssr = 1 - Albedo @@ -326,8 +436,14 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[12] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"tisr\"!") + print("Assuming simplified value for parameter \"tisr\".") + flag_arr[12] = 1 + else: + flag_arr[12] = 1 + + utc = deltaT_ERA * unit.second * 3600 + Time('1900-01-01 00:00:00.0') 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)) @@ -335,8 +451,13 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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\".") + if flag_arr[13] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"tsr\"!") + print("Assuming simplified value for parameter \"tsr\".") + flag_arr[13] = 1 + else: + flag_arr[13] = 1 + tsr = (1 - Albedo) * tisr ttr_pre = np.ma.dot(w2, ERAttr[t_pre_ind, lat_ind2, lon_ind2]) / np.sum(w2) @@ -348,31 +469,55 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): p0 = (p0_post - p0_pre) * (t_epoch - t_pre) + p0_pre if isinstance(p0, float) != True: - print("WARNING: Corrupt ERA5 Data for parameter \"sp\"!") - print("Assuming simplified value for parameter \"sp\".") + if flag_arr[14] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"sp\"!") + print("Assuming simplified value for parameter \"sp\".") + flag_arr[14] = 1 + else: + flag_arr[14] = 1 p0 = 101325.0 if isinstance(ttr, float) != True: - print("WARNING: Corrupt ERA5 Data for parameter \"ttr\"!") - print("Assuming simplified value for parameter \"ttr\".") + if flag_arr[15] == 0: + print("WARNING: Corrupt or missing ERA5 Data for parameter \"ttr\"!") + print("Assuming simplified value for parameter \"ttr\".") + flag_arr[15] = 1 + else: + flag_arr[15] = 1 + + utc = deltaT_ERA * unit.second * 3600 + Time('1900-01-01 00:00:00.0') 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 / 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 / grav(lat, h) * vw_z_interp1d(interp4d_height[-1]) * R_air * T_air / p_air + if h > np.amax(interp4d_height): + if flag_arr[16] == 0: + print("WARNING: Balloon altitude above interpolation area!") + flag_arr[16] = 1 + else: + flag_arr[16] = 1 + + p_air = press_interp1d(np.amax(interp4d_height)) + T_air = temp_interp1d(np.amax(interp4d_height)) + u = vw_x_interp1d(np.amax(interp4d_height)) + v = vw_y_interp1d(np.amax(interp4d_height)) + w = -1 / grav(lat, h) * vw_z_interp1d(np.amax(interp4d_height)) * R_air * T_air / p_air + + elif h < np.amin(interp4d_height): + if flag_arr[17] == 0: + print("WARNING: Balloon altitude below interpolation area!") + flag_arr[17] = 1 + else: + flag_arr[17] = 1 + + p_air = press_interp1d(np.amin(interp4d_height)) + T_air = temp_interp1d(np.amin(interp4d_height)) + u = vw_x_interp1d(np.amin(interp4d_height)) + v = vw_y_interp1d(np.amin(interp4d_height)) + w = -1 / grav(lat, h) * vw_z_interp1d(np.amin(interp4d_height)) * R_air * T_air / p_air + else: p_air = press_interp1d(h) T_air = temp_interp1d(h) @@ -384,10 +529,15 @@ def ERA5Data(lon, lat, h, t, deltaT_ERA): 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 +t_start = Time(start_utc) -deltaT_ERA = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 -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) +m_gas_init = ((m_pl + m_film + m_bal_init) * (FreeLift / 100 + 1)) / (R_gas / R_air - 1) + +deltaT_ERA = (t_start.jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 +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, flag_arr) + + +A_top0 = np.pi/4 * 1.383 ** 2 * (m_gas_init * R_gas * T_air0 / p_air0) ** (2/3) y0 = [ start_lon, # start longitude [deg] @@ -404,8 +554,32 @@ y0 = [ ] -def model(t, y, m_pl, m_film, c_virt): - utc = Time(start_utc) + t * unit.second +t_list, h_list, v_list = [], [], [] +lat_list, lon_list = [], [] +p_list, rho_list = [], [] +Temp_list, Tgas_list, T_film_list = [], [], [] +rhog_list = [] +V_b_list = [] +Q_Albedo_list = [] +Q_IREarth_list = [] +Q_Sun_list = [] +Q_IRFilm_list = [] +Q_IRout_list = [] +Q_ConvExt_list = [] +Q_ConvInt_list = [] +utc_list = [] +ssr_list = [] +ssrd_list = [] +ttr_list = [] +strd_list = [] +strn_list = [] +tisr_list = [] +tsr_list = [] + + + +def model(t, y, m_pl, m_film, c_virt, A_top0, t_start): + utc = t_start + t * unit.second lon = y[0] # 1 lat = y[1] # 2 h = y[2] # 3 @@ -430,21 +604,61 @@ def model(t, y, m_pl, m_film, c_virt): else: lat = lat + if h > 53700: + h = 53700 + elif h < 0: + h = 0 + else: + h = h + + h_list.append(h) + utc_list.append(utc) + lat_list.append(lat) + lon_list.append(lon) + Tgas_list.append(T_gas) + T_film_list.append(T_film) + + 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 + deltaT_ERA = (t_start.jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 # conversion to ERA5 time format + + 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 + + HalfConeAngle = np.arcsin(R_E / (R_E + h)) + ViewFactor = (1 - np.cos(HalfConeAngle)) / 2 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) + 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, flag_arr) + tau_atm, tau_atmIR = tau(ELV, h, p_air, p0) + tau_atm0, tau_atmIR0 = tau(ELV, 0, p0, p0) + I_SunZ = I_Sun * tau_atm + I_Sun0 = I_Sun * tau_atm0 except: - # in case of solver (temporarily) exceeding interpolation area (with subsequent correction) + # in case of solver (temporarily) exceeding interpolation area (with subsequent correction by the solver itself) # 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 + p0 = 101325 + p_air = p_air_simple(h) + tau_atm, tau_atmIR = tau(ELV, h, p_air, p0) + tau_atm0, tau_atmIR0 = tau(ELV, 0, p0, p0) + I_SunZ = I_Sun * tau_atm + I_Sun0 = I_Sun * tau_atm0 + + p_air, p0, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, strdc, ssrd, tsr, ttr, tisr, skt = p_air_simple(h), 101325, T_air_simple(h), rho_air_simple(h), 0, 0, 0, 2000, cc, cc/3, cc/3, cc/3, (1 - Albedo), 0, (epsilon_ground * sigma * T_ground ** 4), (epsilon_ground * sigma * T_ground ** 4), 1, (1 - Albedo) * (I_Sun * np.sin(np.deg2rad(ELV))), (epsilon_ground * sigma * T_ground ** 4 * tau_atmIR * ViewFactor * 2), (I_Sun * np.sin(np.deg2rad(ELV))), T_ground 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 + p0 = 101325 + p_air = p_air_simple(h) + tau_atm, tau_atmIR = tau(ELV, h, p_air, p0) + tau_atm0, tau_atmIR0 = tau(ELV, 0, p0, p0) + I_SunZ = I_Sun * tau_atm + I_Sun0 = I_Sun * tau_atm0 + p_air, p0, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, strdc, ssrd, tsr, ttr, tisr, skt = p_air_simple(h), 101325, T_air_simple(h), rho_air_simple(h), 0, 0, 0, 2000, cc, cc/3, cc/3, cc/3, (1 - Albedo), 0, (epsilon_ground * sigma * T_ground ** 4), (epsilon_ground * sigma * T_ground ** 4), 1, (1 - Albedo) * (I_Sun * np.sin(np.deg2rad(ELV))), (epsilon_ground * sigma * T_ground ** 4 * tau_atmIR * ViewFactor * 2), (I_Sun * np.sin(np.deg2rad(ELV))), T_ground p_gas = p_air @@ -464,23 +678,22 @@ def model(t, y, m_pl, m_film, c_virt): 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: + if m_gas < 0: # limit gas mass to plausible value m_gas = 0 - else: - pass V_b = m_gas / rho_gas # calculate balloon volume from current gas mass and gas density + rhog_list.append(rho_gas) if V_b > V_design: c_duct = c_ducts elif V_b < 0: c_duct = 0 - V_b = 0 + V_b = 1.0 else: c_duct = 0 + V_b_list.append(V_b) + if ballasting(utc) == True: if m_bal >= 0: mdot = m_baldot @@ -518,27 +731,13 @@ def model(t, y, m_pl, m_film, c_virt): 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 - - AZ, ELV = sun_angles_analytical(lat, lon, h, utc) + A_top0 = A_top0 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 # CALCULATIONS FOR THERMAL MODEL - 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 - - 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 @@ -680,12 +879,33 @@ def model(t, y, m_pl, m_film, c_virt): Q_ConvExt = HC_external * A_eff * (T_air - T_film) Q_ConvInt = HC_internal * A_eff * (T_film - T_gas) - if simple == True: - c_d = 0.47 # sphere - else: - c_d = cd_PalumboHigh(Fr, Re, A_top) # Ref. xx + Q_Albedo_list.append(Q_Albedo) + Q_IREarth_list.append(Q_IREarth) + Q_Sun_list.append(Q_Sun) + Q_IRFilm_list.append(Q_IRFilm) + Q_IRout_list.append(Q_IRout) + Q_ConvExt_list.append(Q_ConvExt) + Q_ConvInt_list.append(Q_ConvInt) - #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 + ssr_list.append(ssr) + ssrd_list.append(ssrd) + ttr_list.append(ttr) + strd_list.append(strd) + strn_list.append(strn) + tisr_list.append(tisr) + tsr_list.append(tsr) + + if simple == True: + c_d = c_d + else: + if drag_model == 'PalumboHigh': + c_d = cd_PalumboHigh(Fr, Re, A_top, A_top0) + elif drag_model == 'Palumbo': + c_d = cd_Palumbo(Fr, Re, A_top, A_top0) + elif drag_model == 'PalumboLow': + c_d = cd_PalumboLow(Fr, Re, A_top, A_top0) + else: + c_d = cd_sphere(Re) D = drag(c_d, rho_air, A_drag, v_rel) # calculate drag force @@ -707,7 +927,12 @@ def model(t, y, m_pl, m_film, c_virt): 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(np.abs(2 * dP_duct * rho_gas))) - (A_valve * c2 * np.sqrt(np.abs(2 * dP_valve * rho_gas))) + + if eqn9 > 0: + eqn9 = 0 + eqn10 = c2dot + if m_bal > 0: eqn11 = -mdot else: @@ -733,28 +958,143 @@ def below_float(t, y, m_pl, m_film, c_virt): 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 +excess_ascent = lambda t, x: above_float(t, x, m_pl, m_film, c_virt) +excess_ascent.terminal = True +excess_ascent.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") +t0 = 0 +tf = t_sim -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]) +print("") +print("BEGINNING SIMULATION") -# tnew = np.linspace(0, sol.t[-1], len(Vol_list)) +sol = solve_ivp(fun=lambda t, x: model(t, x, m_pl, m_film, c_virt, A_top0, t_start), t_span=[t0, tf], y0=y0, method='RK45', events=[hit_ground, excess_ascent, instable]) #, t_eval=comp_time + +tnew = np.linspace(0, sol.t[-1], len(V_b_list)) print(sol.message) + +""" +lonsol = sol.y[0, :] +latsol = sol.y[1, :] +hsol = sol.y[2, :] + +x_sol, y_sol, z_sol = transform2(lonsol, latsol, hsol) +x_test, y_test, z_test = transform2(comp_lon, comp_lat, comp_height) + +delta = ((x_sol - x_test)**2 + (y_sol - y_test)**2 + (z_sol - z_test)**2)**(0.5) + +val = 0 +i = 0 + +for x in delta: + print(latsol[i]) + print(comp_lat[i]) + print(latsol[i] - comp_lat[i]) + print(x) + val += x ** 2 + i += 1 + +RMS = np.sqrt(val/i) + +print('RMS') +print(RMS) +""" + + print(datetime.now() - starttime) + +arr0 = np.linspace(0, sol.t[-1], len(V_b_list)) +arr1 = np.asarray(utc_list) +arr2 = np.asarray(h_list) +arr3 = np.asarray(lat_list) +arr4 = np.asarray(lon_list) +arr5 = np.asarray(Tgas_list) +arr6 = np.asarray(T_film_list) +arr7 = np.asarray(rhog_list) +arr8 = np.asarray(V_b_list) +arr9 = np.asarray(Q_Albedo_list) +arr10 = np.asarray(Q_IREarth_list) +arr11 = np.asarray(Q_Sun_list) +arr12 = np.asarray(Q_IRFilm_list) +arr13 = np.asarray(Q_IRout_list) +arr14 = np.asarray(Q_ConvExt_list) +arr15 = np.asarray(Q_ConvInt_list) +arr16 = np.asarray(ssr_list) +arr17 = np.asarray(ssrd_list) +arr18 = np.asarray(ttr_list) +arr19 = np.asarray(strd_list) +arr20 = np.asarray(strn_list) +arr21 = np.asarray(tisr_list) +arr22 = np.asarray(tsr_list) + +ind_list = [] +for i in range(len(arr0)): + if arr0[i - 1] == arr0[i]: + ind_list.append(i) + +arr0 = np.delete(arr0, ind_list) +arr1 = np.delete(arr1, ind_list) +arr2 = np.delete(arr2, ind_list) +arr3 = np.delete(arr3, ind_list) +arr4 = np.delete(arr4, ind_list) +arr5 = np.delete(arr5, ind_list) +arr6 = np.delete(arr6, ind_list) +arr7 = np.delete(arr7, ind_list) +arr8 = np.delete(arr8, ind_list) +arr9 = np.delete(arr9, ind_list) +arr10 = np.delete(arr10, ind_list) +arr11 = np.delete(arr11, ind_list) +arr12 = np.delete(arr12, ind_list) +arr13 = np.delete(arr13, ind_list) +arr14 = np.delete(arr14, ind_list) +arr15 = np.delete(arr15, ind_list) +arr16 = np.delete(arr16, ind_list) +arr17 = np.delete(arr17, ind_list) +arr18 = np.delete(arr18, ind_list) +arr19 = np.delete(arr19, ind_list) +arr20 = np.delete(arr20, ind_list) +arr21 = np.delete(arr21, ind_list) +arr22 = np.delete(arr22, ind_list) + + + +df1 = pd.DataFrame(data={ + 'time [s]': arr0, + 'UTC': arr1, + 'Altitude [m]': arr2, + 'Latitude [deg]': arr3, + 'Longitude [deg]': arr4, + 'T_gas [K]': arr5, + 'T_film [K]': arr6, + 'rho_gas [kg/m^3]': arr7, + 'V_balloon [m^3]': arr8, + 'Q_Albedo [W/m^2]': arr9, + 'Q_IR_Earth [W/m^2]': arr10, + 'Q_Sun [W/m^2]': arr11, + 'Q_IRFilm [W/m^2]': arr12, + 'Q_IRout [W/m^2]': arr13, + 'Q_ConvExt [W/m^2]': arr14, + 'Q_ConvInt [W/m^2]': arr15, + 'SSR [W/m^2]': arr16, + 'SSRD [W/m^2]': arr17, + 'TTR [W/m^2]': arr18, + 'STRD [W/m^2]': arr19, + 'STRN [W/m^2]': arr20, + 'TISR [W/m^2]': arr21, + 'TSR [W/m^2]': arr22 +}) + +df1.to_excel("output.xlsx") + plt.plot(sol.t, sol.y[2, :], 'k--', label='Simulation') -plt.plot(comp_time, comp_height, 'r-', label='PoGo Flight Test') +plt.plot(comp_time, comp_height, 'r-', label='PoGo+ Flight Test') plt.legend() plt.title('high factor') plt.xlabel('time in s') @@ -762,17 +1102,14 @@ plt.ylabel('Balloon Altitude in m') plt.show() plt.clf() -# ax = plt.axes(projection=ccrs.Mollweide()) -ax = plt.axes(projection=ccrs.AzimuthalEquidistant(central_latitude=90)) +ax = plt.axes(projection=ccrs.AzimuthalEquidistant(central_latitude=-90)) ax.coastlines() -ax.gridlines() +ax.gridlines(draw_labels=True, linewidth=0.25, color='black') 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() +plt.show() \ No newline at end of file