BASTET/main.py

1115 lines
44 KiB
Python
Raw Normal View History

import sys
2021-06-18 14:58:02 +02:00
import pickle
from pyfiglet import Figlet
import warnings
import numpy as np
import xarray as xr
import pandas as pd
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
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
from scipy.spatial import cKDTree
from scipy.integrate import solve_ivp
from input.user_input import *
from input.natural_constants import *
from models.gravity import grav
from models.valving import valving
from models.ballasting import ballasting
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
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")
2021-06-18 14:58:02 +02:00
data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Data_PoGo2016.xls', sheet_name='SuperTIGER2') # Tabelle3
2021-04-09 13:34:38 +02:00
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("")
2021-06-18 14:58:02 +02:00
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.")
2021-06-18 14:58:02 +02:00
first_file = Dataset(ERA5_float[0])
last_file = Dataset(ERA5_float[-1])
2021-06-18 14:58:02 +02:00
tstart = int(first_file.variables['time'][0])
tend = int(last_file.variables['time'][-1])
2021-06-18 14:58:02 +02:00
first_file.close()
last_file.close()
2021-06-18 14:58:02 +02:00
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)
2021-06-18 14:58:02 +02:00
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]
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]
2021-06-18 14:58:02 +02:00
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)
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]
2021-06-18 14:58:02 +02:00
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)
lon_era2d1, lat_era2d1 = np.meshgrid(ERAlon1, ERAlat1)
lon_era2d2, lat_era2d2 = np.meshgrid(ERAlon2, ERAlat2)
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())
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)))
2021-06-18 14:58:02 +02:00
print("Built kd-trees.")
print("")
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
flag_arr = np.zeros(20)
def ERA5Data(lon, lat, h, t, deltaT_ERA, flag_arr):
t_epoch = deltaT_ERA + t / 3600
2021-04-09 13:34:38 +02:00
t_pre = int(t_epoch)
t_pre_ind = t_pre - int(ERAtime[0])
2021-04-09 13:34:38 +02:00
t_post_ind = t_pre_ind + 1
2021-04-09 13:34:38 +02:00
xt, yt, zt = transform(lon, lat) # current coordinates
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_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]
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_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_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_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_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
2021-06-18 14:58:02 +02:00
pressure_hPa = np.array([1, 2, 3, 5, 7, 10, 20, 30]) # !!!
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:
2021-06-18 14:58:02 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
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(
[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:
2021-06-18 14:58:02 +02:00
if flag_arr[19] == 0:
print("Error: Check time range of ERA5 data!")
flag_arr[19] = 1
else:
flag_arr[19] = 1
else:
pass
2021-04-09 13:34:38 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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
2021-04-09 13:34:38 +02:00
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
if isinstance(strn, float) != True:
2021-06-18 14:58:02 +02:00
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)
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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
if isinstance(strd, float) != True:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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
if isinstance(ssrd, float) != True:
2021-06-18 14:58:02 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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
2021-04-09 13:34:38 +02:00
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
if isinstance(tisr, float) != True:
2021-06-18 14:58:02 +02:00
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))
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
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
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
if isinstance(p0, float) != True:
2021-06-18 14:58:02 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
else:
p_air = press_interp1d(h)
T_air = temp_interp1d(h)
2021-04-09 13:34:38 +02:00
u = vw_x_interp1d(h)
v = vw_y_interp1d(h)
w = -1 / grav(lat, h) * vw_z_interp1d(h) * R_air * T_air / p_air
2021-04-09 13:34:38 +02:00
rho_air = p_air / (R_air * T_air)
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
2021-06-18 14:58:02 +02:00
t_start = Time(start_utc)
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)
2021-04-09 13:34:38 +02:00
2021-06-18 14:58:02 +02:00
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]
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]
]
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
lon = y[0] # 1
lat = y[1] # 2
h = y[2] # 3
v_x = y[3] # 4
v_y = y[4] # 5
v_z = y[5] # 6
T_gas = y[6] # 7
T_film = y[7] # 8
m_gas = y[8] # 9
c2 = y[9] # 10
m_bal = y[10] # 11
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
2021-04-09 13:34:38 +02:00
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
try:
2021-06-18 14:58:02 +02:00
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:
2021-06-18 14:58:02 +02:00
# 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):
2021-06-18 14:58:02 +02:00
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:
2021-06-18 14:58:02 +02:00
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
2021-04-09 13:34:38 +02:00
h_valve = 1.034 * V_design ** (1 / 3)
h_duct = 0.47 * h_valve
v_relx = u - v_x # relative wind velocity x-dir (balloon frame)
v_rely = v - v_y # relative wind velocity y-dir (balloon frame)
v_relz = w - v_z # relative wind velocity z-dir (balloon frame)
v_rel = np.sqrt(v_relx ** 2 + v_rely ** 2 + v_relz ** 2) # total relative wind velocity (balloon frame)
2021-04-09 13:34:38 +02:00
alpha = np.arcsin(v_relz / v_rel) # "angle of attack": angle between longitudinal axis and rel. wind (in [rad])
2021-04-09 13:34:38 +02:00
rho_gas = p_gas / (R_gas * T_gas) # calculate gas density through *ideal* gas equation
dP_valve = grav(lat, h) * (rho_air - rho_gas) * h_valve
dP_duct = grav(lat, h) * (rho_air - rho_gas) * h_duct
2021-06-18 14:58:02 +02:00
if m_gas < 0: # limit gas mass to plausible value
m_gas = 0
2021-04-09 13:34:38 +02:00
V_b = m_gas / rho_gas # calculate balloon volume from current gas mass and gas density
2021-06-18 14:58:02 +02:00
rhog_list.append(rho_gas)
if V_b > V_design:
2021-04-09 13:34:38 +02:00
c_duct = c_ducts
elif V_b < 0:
c_duct = 0
2021-06-18 14:58:02 +02:00
V_b = 1.0
else:
2021-04-09 13:34:38 +02:00
c_duct = 0
2021-06-18 14:58:02 +02:00
V_b_list.append(V_b)
if ballasting(utc) == True:
if m_bal >= 0:
mdot = m_baldot
else:
mdot = 0
else:
mdot = 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
2021-04-09 13:34:38 +02:00
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
2021-04-09 13:34:38 +02:00
d_b = 1.383 * V_b ** (1 / 3) # calculate diameter of balloon from its volume
L_goreB = 1.914 * V_b ** (1 / 3)
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
2021-04-09 13:34:38 +02:00
A_top = np.pi / 4 * d_b ** 2
2021-06-18 14:58:02 +02:00
A_top0 = A_top0
2021-04-09 13:34:38 +02:00
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
if simple == True:
q_IREarth = epsilon_ground * sigma * T_ground ** 4 * tau_atmIR
2021-04-09 13:34:38 +02:00
if ELV <= 0:
q_Albedo = 0
else:
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:
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
2021-04-09 13:34:38 +02:00
else:
Q_IREarth = (q_IREarth2 - q_IREarth1) / 40000 * h + q_IREarth1
2021-04-09 13:34:38 +02:00
if ELV <= 0:
q_Albedo1 = 0
q_Albedo2 = 0
2021-04-09 13:34:38 +02:00
else:
if ssrd == 0:
q_Albedo1 = 0
2021-04-09 13:34:38 +02:00
else:
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
2021-04-09 13:34:38 +02:00
else:
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
2021-04-09 13:34:38 +02:00
else:
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))
# 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
2021-04-09 13:34:38 +02:00
if ELV <= 0:
Q_Albedo = 0
2021-04-09 13:34:38 +02:00
else:
Q_Albedo = q_Albedo_bc
elif h >= 12000: # "above clouds"
Q_IREarth = q_IREarth_ac
2021-04-09 13:34:38 +02:00
q_sun = q_sun_ac
if ELV <= 0:
Q_Albedo = 0
2021-04-09 13:34:38 +02:00
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
2021-04-09 13:34:38 +02:00
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
2021-04-09 13:34:38 +02:00
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 * 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
2021-04-09 13:34:38 +02:00
Re = np.abs(v_rel) * d_b * rho_air / my_air
Fr = np.abs(v_rel) / np.sqrt(grav(lat, h) * d_b)
HC_forced = k_air / d_b * (2 + 0.41 * Re ** 0.55)
2021-04-09 13:34:38 +02:00
HC_internal = 0.13 * k_gas * (
(rho_gas ** 2 * grav(lat, h) * np.abs(T_film - T_gas) * Pr_gas) / (T_gas * my_air ** 2)) ** (
2021-04-09 13:34:38 +02:00
1 / 3)
HC_external = np.maximum(HC_free, HC_forced)
Q_Sun = alpha_VIS * A_proj * q_sun * (1 + tau_VIS / (1 - r_VIS))
2021-04-09 13:34:38 +02:00
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)
2021-06-18 14:58:02 +02:00
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)
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:
2021-06-18 14:58:02 +02:00
c_d = c_d
else:
2021-06-18 14:58:02 +02:00
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)
2021-04-09 13:34:38 +02:00
D = drag(c_d, rho_air, A_drag, v_rel) # calculate drag force
if v_rel == 0:
Drag_x, Drag_y, Drag_z = 0, 0, 0
else:
2021-04-09 13:34:38 +02:00
Drag_x, Drag_y, Drag_z = D * v_relx / v_rel, D * v_rely / v_rel, D * v_relz / v_rel
F = grav(lat, h) * V_b * (rho_air - rho_gas) - grav(lat, h) * m_gross + Drag_z # gross inflation - weight + drag
2021-04-09 13:34:38 +02:00
a_x, a_y, a_z = Drag_x / m_virt, Drag_y / m_virt, F / m_virt
2021-04-09 13:34:38 +02:00
eqn1 = np.rad2deg(y[3] / r_lon)
eqn2 = np.rad2deg(y[4] / r_lat)
eqn3 = v_z
2021-04-09 13:34:38 +02:00
eqn4 = a_x
eqn5 = a_y
eqn6 = a_z
eqn7 = Q_ConvInt / (gamma * c_v * m_gas) - (gamma - 1) / gamma * (rho_air * grav(lat, h)) / (rho_gas * R_gas) * v_z
2021-04-09 13:34:38 +02:00
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)))
2021-06-18 14:58:02 +02:00
if eqn9 > 0:
eqn9 = 0
eqn10 = c2dot
2021-06-18 14:58:02 +02:00
if m_bal > 0:
eqn11 = -mdot
else:
eqn11 = 0
return [eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10, eqn11]
# DEFINITION OF EVENTS FOR SOLVER
def at_ground(t, y, m_pl, m_film, c_virt):
return y[2]
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
2021-06-18 14:58:02 +02:00
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
2021-06-18 14:58:02 +02:00
t0 = 0
tf = t_sim
2021-06-18 14:58:02 +02:00
print("")
print("BEGINNING SIMULATION")
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
2021-06-18 14:58:02 +02:00
tnew = np.linspace(0, sol.t[-1], len(V_b_list))
2021-04-09 13:34:38 +02:00
print(sol.message)
2021-06-18 14:58:02 +02:00
"""
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)
2021-06-18 14:58:02 +02:00
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')
2021-06-18 14:58:02 +02:00
plt.plot(comp_time, comp_height, 'r-', label='PoGo+ Flight Test')
plt.legend()
plt.title('high factor')
2021-04-09 13:34:38 +02:00
plt.xlabel('time in s')
plt.ylabel('Balloon Altitude in m')
plt.show()
2021-04-09 13:34:38 +02:00
plt.clf()
2021-06-18 14:58:02 +02:00
ax = plt.axes(projection=ccrs.AzimuthalEquidistant(central_latitude=-90))
2021-04-09 13:34:38 +02:00
ax.coastlines()
2021-06-18 14:58:02 +02:00
ax.gridlines(draw_labels=True, linewidth=0.25, color='black')
ax.stock_img()
2021-04-09 13:34:38 +02:00
ax.set_extent([-120, 30, 60, 80], crs=ccrs.PlateCarree())
2021-04-09 13:34:38 +02:00
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.savefig(os.path.join(rootdir, figname))
2021-06-18 14:58:02 +02:00
plt.show()