From 011ef6968d82eb4dfb0af2c188f4523bca6e5e43 Mon Sep 17 00:00:00 2001 From: Marcel Frommelt Date: Fri, 9 Apr 2021 20:34:38 +0900 Subject: [PATCH] Implemented ERA5 radiation properties --- DataRequest.py | 55 ---- main.py | 704 ++++++++++++++++++++++++++----------------------- 2 files changed, 379 insertions(+), 380 deletions(-) delete mode 100644 DataRequest.py diff --git a/DataRequest.py b/DataRequest.py deleted file mode 100644 index d6447c0..0000000 --- a/DataRequest.py +++ /dev/null @@ -1,55 +0,0 @@ -import cdsapi - -c = cdsapi.Client() - - - - -r = c.retrieve( - 'reanalysis-era5-pressure-levels', - { - 'product_type': 'reanalysis', - 'variable': [ - 'fraction_of_cloud_cover', 'geopotential', 'relative_humidity', - 'specific_humidity', 'temperature', 'u_component_of_wind', - 'v_component_of_wind', 'vertical_velocity', - ], - 'pressure_level': [ - '1', '2', '3', - '5', '7', '10', - '20', '30', '50', - '70', '100', '125', - '150', '175', '200', - '225', '250', '300', - '350', '400', '450', - '500', '550', '600', - '650', '700', '750', - '775', '800', '825', - '850', '875', '900', - '925', '950', '975', - '1000', - ], - 'year': '2016', - 'month': '07', - 'day': [ - '11', '12', '13', - '14', '15', '16', - '17', '18', - ], - 'time': [ - '00:00', '01:00', '02:00', - '03:00', '04:00', '05:00', - '06:00', '07:00', '08:00', - '09:00', '10:00', '11:00', - '12:00', '13:00', '14:00', - '15:00', '16:00', '17:00', - '18:00', '19:00', '20:00', - '21:00', '22:00', '23:00', - ], - 'area': [ - 72, -111, 67, - 22, - ], - 'format': 'netcdf', - }) -r.download('test2021.nc') \ No newline at end of file diff --git a/main.py b/main.py index e56564e..e49633e 100644 --- a/main.py +++ b/main.py @@ -4,325 +4,424 @@ 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 cartopy.crs as ccrs from datetime import datetime start = datetime.now() print(start) begin_time = datetime.now() -from models.sun import sun_angles_analytical -from models.sun import sun_angles_astropy -from models.drag import drag #, c_d +from 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 input.user_input import * from input.natural_constants import * from astropy.time import Time import astropy.units as unit -from models.thermal import AirMass +# from models.thermal import from netCDF4 import Dataset import pandas as pd -data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Mappe1.xls', sheet_name='Tabelle3') +data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Mappe1.xls', sheet_name='Tabelle3') # Tabelle3 -comp_time = pd.DataFrame(data, columns= ['Time']) -comp_height = pd.DataFrame(data, columns= ['Height']) -comp_lat = pd.DataFrame(data, columns= ['Latitude']) -comp_lon = pd.DataFrame(data, columns= ['Longitude']) +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() -Reynolds = [0.0, 0.002, 0.004, 0.007, 0.01] #, 0.04, 0.07, 0.1, 0.2, 0.3, 0.5, 0.7, 1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 200, -# 300, 500, 700, 1000, 2000, 3000, 5000, 7000, 10000, 20000, 30000, 50000] -drag_coeff = [20000, 12000, 6000, 3429, 2400] #, 600, 343, 240, 120, 80, 49, 36.5, 26.5, 14.4, 10.4, 6.9, 5.4, 4.1, 2.55, 2, 1.5, -# 1.27, 1.07, 0.77, 0.65, 0.55, 0.5, 0.46, 0.42, 0.4, 0.385, 0.39, 0.405, 0.45, 0.47, 0.49] +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))) -cd_interp1d = interpolate.interp1d(Reynolds, drag_coeff) -Re_list = [] +level_data = Dataset("level.nc", "r", format="NETCDF4") +single_data = Dataset("single.nc", "r", format="NETCDF4") -def cd_func(Re): - w = np.log10(Re) +ERAtime = level_data.variables['time'][:] # time - if Re <= 0.01: - c_d = 9/2 + 24/Re - elif Re <= 20: - c_d = 24/Re * (1 + 0.1315 * Re ** (0.82 - 0.05 * w)) - elif Re <= 260: - c_d = 24/Re * (1 + 0.1935 * Re ** 0.6305) - elif Re <= 1.5 * 10e3: - c_d = 10 ** (1.6435 - 1.1242 * w + 0.1558 * w ** 2) - elif Re <= 1.2 * 10e4: - c_d = 10 ** (-2.4571 + 2.5558 * w - 0.9295 * w ** 2 + 0.1049 * w ** 3) - elif Re <= 4.4 * 10e4: - c_d = 10 ** (-1.9181 + 0.6370 * w - 0.0636 * w ** 3) - elif Re <= 3.38 * 10e5: - c_d = 10 ** (-4.3390 + 1.5809 * w - 0.1546 * w ** 2) - elif Re <= 4 * 10e5: - c_d = 29.78 - 5.3 * w - elif Re <= 10 ** 6: - c_d = 0.1 * w - 0.49 - elif Re > 10e6: - c_d = 0.19 - 8 * 10e4/Re - else: - c_d = 0.47 +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) - return c_d -def cd_func_garg(Re): - if Re == 0: - c_d = 2.0 - elif Re < 2 * 10e5: - c_d = 5.4856 * 10e9 * np.tanh(4.3774 * 10e-9 / Re) + 0.0709 * np.tanh(700.6575 / Re) + 0.3894 * np.tanh(74.1539 / Re) - 0.1198 * np.tanh(7429.0843 / Re) + 1.7174 * np.tanh(9.9851 / Re + 2.3384) + 0.4744 - elif 2 * 10e5 < Re <= 10e6: - c_d = 8 * 10e-6 * ((Re/6530) ** 2 + np.tanh(Re) - 8 * np.log(Re) / np.log(10)) - 0.4119 * np.exp(-2.08 * 10 ** 43 / (Re + Re ** 2) ** 4) - 2.1344 * np.exp(-((np.log(Re ** 2 + 10.7563) / np.log(10)) ** 2 + 9.9867) / Re) + 0.1357 * np.exp(-((Re / 1620) ** 2 + 10370) / Re) - 8.5 * 10e-3 * (2 * np.log(np.tanh(np.tanh(Re))) / np.log(10) - 2825.7162) / Re + 2.4795 - elif Re > 10e6: - c_d = 0.212546 - else: - c_d = 2.0 - return 0.80 # c_d - -#cd_interp1d = interpolate.interp1d(Reynolds, drag_coeff) - -def transform(lon, lat, t): - # WGS 84 reference coordinate system parameters - A = 6378137.0 # major axis [m] - E2 = 6.69437999014e-3 # eccentricity squared - - t_s = 1000 * t - - lon_rad = np.radians(lon) - lat_rad = np.radians(lat) - # convert to cartesian coordinates - r_n = A / (np.sqrt(1 - E2 * (np.sin(lat_rad) ** 2))) - x = r_n * np.cos(lat_rad) * np.cos(lon_rad) - y = r_n * np.cos(lat_rad) * np.sin(lon_rad) - z = r_n * (1 - E2) * np.sin(lat_rad) - return x, y, z, t_s - -data = Dataset("test2021.nc", "r", format="NETCDF4") - -ERAtime = data.variables['time'][:] # time -ERAlat = data.variables['latitude'][:] # latitude -ERAlon = data.variables['longitude'][:] # longitude -ERAz = data.variables['z'][:]/g # geopotential to geopotential height -ERApress = data.variables['level'][:] # pressure level -ERAtemp = data.variables['t'][:] # temperature in K - -vw_x = data.variables['u'][:] # v_x -vw_y = data.variables['v'][:] # v_y -vw_z = data.variables['w'][:] # v_z - -lon_era2d, lat_era2d, time_era = np.meshgrid(ERAlon, ERAlat, ERAtime) - -xs, ys, zs, ts = transform(lon_era2d.flatten(), lat_era2d.flatten(), time_era.flatten()) - -tree = cKDTree(np.column_stack((xs, ys, zs, ts))) # ! +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 - -lat = start_lat # deg -lon = start_lon # deg - -t = 0 # simulation time in seconds -h = start_height -utc = Time(start_utc) -epoch_diff = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000 -t_epoch = epoch_diff +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 +vw_x = level_data.variables['u'][:] # v_x +vw_y = level_data.variables['v'][:] # v_y +vw_z = level_data.variables['w'][:] # v_z -xt, yt, zt, tt = transform(lon, lat, t) # test coordinates +lon_era2d, lat_era2d = np.meshgrid(ERAlon, ERAlat) +lon_era2ds, lat_era2ds = np.meshgrid(ERAlons, ERAlats) -d, inds = tree.query(np.column_stack((xt, yt, zt, tt)), k=8) # longitude, latitude, time in h # ! -w = 1.0 / d[0] ** 2 +xs, ys, zs = transform(lon_era2d.flatten(), lat_era2d.flatten()) +xs2, ys2, zs2 = transform(lon_era2ds.flatten(), lat_era2ds.flatten()) -lat_sel = np.unravel_index(inds[0], lon_era2d.shape)[0] -lon_sel = np.unravel_index(inds[0], lon_era2d.shape)[1] -time_sel = np.unravel_index(inds[0], lon_era2d.shape)[2] - -interp4d_height = np.ma.dot(w, ERAz[time_sel, :, lat_sel, lon_sel]) / np.sum(w) -interp4d_temp = np.ma.dot(w, ERAtemp[time_sel, :, lat_sel, lon_sel]) / np.sum(w) -interp4d_vw_x = np.ma.dot(w, vw_x[time_sel, :, lat_sel, lon_sel]) / np.sum(w) -interp4d_vw_y = np.ma.dot(w, vw_y[time_sel, :, lat_sel, lon_sel]) / np.sum(w) -interp4d_vw_z = np.ma.dot(w, vw_z[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - -pressure_hPa = np.array([1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 125, 150, 175, 200,225, 250, 300, 350, 400, - 450, 500, 550, 600, 650, 700, 750, 775, 800, 825, 850, 875, 900, 925, 950, 975, 1000]) - -pressure = 100 * pressure_hPa -# print(pressure) +tree = cKDTree(np.column_stack((xs, ys, zs))) +tree2 = cKDTree(np.column_stack((xs2, ys2, zs2))) -# height_interp1d = interpolate.interp1d(interp4d_height, pressure) -temp_interp1d = interpolate.interp1d(interp4d_height, interp4d_temp) -press_interp1d = interpolate.interp1d(interp4d_height, pressure) -vw_x_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_x) -vw_y_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_y) -vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z) +# INITIALISATION TIME-COUNTERS -#u = 0 -#v = 0 -#w = 0 -T_air = temp_interp1d(h) -p_air = press_interp1d(h) -rho_air = p_air/(R_air * T_air) +# MAYBE LATER +# t_pre = int(t_epoch) +# t_post = t_pre + 1 +# t1 = np.searchsorted(ERAtime, t_pre, side='left') +# t2 = t1 + 1 -u = vw_x_interp1d(h) -v = vw_y_interp1d(h) -w = -1 / g * vw_z_interp1d(h) * R_air * T_air / p_air +def ERA5Data(lon, lat, h, t, deltaT_ERA): + t_epoch = deltaT_ERA + t/3600 + t_pre = int(t_epoch) + # t_post = t_pre + 1 + t_pre_ind = t_pre - ERAtime[0] + t_post_ind = t_pre_ind + 1 -v_x = 0 -v_y = 0 -v_z = 0 + xt, yt, zt = transform(lon, lat) # current coordinates -v_rel = 0 + 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 + w2 = 1.0 / d2[0] ** 2 -T_gas = T_air -T_film = T_gas + 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] -t_list = [] -h_list = [] -v_list = [] -lat_list = [] -lon_list = [] -p_list = [] -Temp_list = [] -rho_list = [] + 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 -while t <= t_end and h >= 0: - t_list.append(t) - h_list.append(h) - v_list.append(v_rel) - lat_list.append(lat) - lon_list.append(lon) + 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 - p_list.append(p_air) - Temp_list.append(T_air) - rho_list.append(rho_air) + 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 - t_epoch = epoch_diff + t/3600 + 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 - xt, yt, zt, tt = transform(lon, lat, t_epoch) # current balloon coordinates in cartesian coordinates: x [m], y [m], z [m], time [h since 1900-01-01 00:00:00.0] + 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 - d, inds = tree.query(np.column_stack((xt, yt, zt, tt)), k=8) # longitude, latitude, time in h # ! - w = 1.0 / d[0] ** 2 + 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 - lat_sel = np.unravel_index(inds[0], lon_era2d.shape)[0] - lon_sel = np.unravel_index(inds[0], lon_era2d.shape)[1] - time_sel = np.unravel_index(inds[0], lon_era2d.shape)[2] + 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 - interp4d_height = np.ma.dot(w, ERAz[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - interp4d_temp = np.ma.dot(w, ERAtemp[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - interp4d_vw_x = np.ma.dot(w, vw_x[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - interp4d_vw_y = np.ma.dot(w, vw_y[time_sel, :, lat_sel, lon_sel]) / np.sum(w) - interp4d_vw_z = np.ma.dot(w, vw_z[time_sel, :, lat_sel, lon_sel]) / np.sum(w) + 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 + + 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 + + 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 + + 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 + + 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 + + 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 + + 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 + + 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 + + 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 + + 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 + + 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 + + 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 + + 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 - press_interp1d = interpolate.interp1d(interp4d_height, pressure) temp_interp1d = interpolate.interp1d(interp4d_height, interp4d_temp) + press_interp1d = interpolate.interp1d(interp4d_height, pressure) vw_x_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_x) vw_y_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_y) vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z) - - try: + 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 + 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 + else: p_air = press_interp1d(h) T_air = temp_interp1d(h) - except: - if (h > interp4d_height[0]): - h = interp4d_height[0] - else: - h = interp4d_height[-1] + u = vw_x_interp1d(h) + v = vw_y_interp1d(h) + w = -1 / g * vw_z_interp1d(h) * R_air * T_air / p_air - p_air = press_interp1d(h) - T_air = temp_interp1d(h) - - #print("height") - #print(h) - #print("temperature") - #print(temp_interp1d(h)) - #print("pressure") - #print(press_interp1d(h)) - #print("vw_x") - #print(vw_x_interp1d(h)) - - p_gas = p_air rho_air = p_air / (R_air * T_air) - u = vw_x_interp1d(h) - v = vw_y_interp1d(h) - w = -1/g * vw_z_interp1d(h) * R_air * T_air / p_air - v_relx = u - v_x - v_rely = v - v_y - v_relz = w - v_z + 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 - v_rel = (v_relx ** 2 + v_rely ** 2 + v_relz ** 2) ** (1/2) - rho_gas = p_gas/(R_gas * T_gas) # calculate gas density through ideal(!) gas equation +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) - V_b = m_gas/rho_gas # calculate balloon volume from current gas mass and gas density + + + + +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): + utc = Time(start_utc) + t * unit.second + 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 + + r_lon, r_lat = radii(lat, h) + + 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) + p_gas = p_air + + 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) + + alpha = np.arcsin(v_relz / v_rel) # "angle of attack": angle between longitudinal axis and rel. wind (in [rad]) + + 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 + + V_b = m_gas / rho_gas # calculate balloon volume from current gas mass and gas density if V_b > V_design: - V_b = V_design - m_gas = V_design * rho_gas + c_duct = c_ducts else: - pass + c_duct = 0 - m_gross = m_pl + m_film + #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 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 + + ## 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 - d_b = 1.383 * V_b ** (1/3) # calculate diameter of balloon from its volume - L_goreB = 1.914 * V_b ** (1/3) - h_b = 0.748 * d_b - A_surf = 4.94 * V_b ** (2/3) - A_surf1 = 4.94 * V_design ** (2/3) * (1 - np.cos(np.pi * L_goreB/L_goreDesign)) + 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 - A_top = np.pi/4 * d_b ** 2 + A_top = np.pi / 4 * d_b ** 2 AZ, ELV = sun_angles_analytical(lat, lon, utc) - A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV))) + 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 ELV >= -(180 / np.pi * np.arccos(R_E / (R_E + h))): - tau_atm = 0.5 * (np.exp(-0.65 * AirMass(p_air, p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air, p_0, ELV, h))) - tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air / p_0) + np.exp(-0.095 * p_air / p_0)) - else: - tau_atm = 0 - tau_atmIR = 0 - - doy = int(utc.doy) + tau_atm, tau_atmIR = tau(ELV, h, p_air) + tau_atm0, tau_atmIR0 = tau(ELV, 0, p_0) MA = (357.52911 + 0.98560028 * (utc.jd - 2451545)) % 360 # in degree, reference: see folder "literature" TA = MA + 2 * e * np.sin(np.deg2rad(MA)) + 5 / 4 * e ** 2 * np.sin(np.deg2rad(2 * MA)) I_Sun = 1367.5 * ((1 + e * np.cos(np.deg2rad(TA))) / (1 - e ** 2)) ** 2 I_SunZ = I_Sun * tau_atm - q_sun = I_SunZ - q_IRground = epsilon_ground * sigma * T_ground ** 4 - q_IREarth = q_IRground * tau_atmIR + 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 - if ELV <= 0: - q_Albedo = 0 + tcc = 0 + + 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 else: - q_Albedo = Albedo * I_Sun * np.sin(np.deg2rad(ELV)) + 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 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_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) - my_gas = 1.895 * 10 ** -5 * (T_gas / 273.15) ** 0.647 k_air = 0.0241 * (T_air / 273.15) ** 0.9 k_gas = 0.144 * (T_gas / 273.15) ** 0.7 Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas - Gr_air = (rho_air ** 2 * g * np.abs(T_film - T_air) * d_b ** 3) / (T_air * my_air ** 2) Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25 HC_free = Nu_air * k_air / d_b - Re = np.abs(v_relz) * d_b * rho_air / my_air - - Re_list.append(Re) + Re = np.abs(v_rel) * d_b * rho_air / my_air + Fr = np.abs(v_rel) / np.sqrt(g * 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)) ** ( - 1 / 3) + HC_internal = 0.13 * k_gas * ( + (rho_gas ** 2 * g * np.abs(T_film - T_gas) * Pr_gas) / (T_gas * my_air ** 2)) ** ( + 1 / 3) HC_external = np.maximum(HC_free, HC_forced) HalfConeAngle = np.arcsin(R_E / (R_E + h)) @@ -331,123 +430,78 @@ while t <= t_end and h >= 0: 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_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) - try: - c_d = cd_func(Re) - except: - print(Re) + c_d = Cd_model(Fr, Re, A_top) - D = drag(c_d, rho_air, d_b, v_rel) # calculate drag force + D = drag(c_d, rho_air, A_drag, v_rel) # calculate drag force if v_rel == 0: - Drag_x = 0 - Drag_y = 0 - Drag_z = 0 + Drag_x, Dray_y, Drag_z = 0, 0, 0 else: - Drag_x = D * v_relx / v_rel - Drag_y = D * v_rely / v_rel - Drag_z = D * v_relz / v_rel + Drag_x, Drag_y, Drag_z = D * v_relx / v_rel, D * v_rely / v_rel, D * v_relz / v_rel - I = g * V_b * (rho_air - rho_gas) # calculate gross inflation - W = g * m_gross # calculate weight (force) - F = I - W + Drag_z + F = g * V_b * (rho_air - rho_gas) - g * m_gross + Drag_z # gross inflation - weight + drag - # RoC = -v_z - # dT_gas = (Q_ConvInt / (gamma * m_gas * c_v) - (gamma - 1) / gamma * rho_air(h) * g / (rho_gas * R_gas) * RoC) * dt - # dT_film = ((Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film)) * dt - # v_w = w - # v = v_z + a_x, a_y, a_z = Drag_x / m_virt, Drag_y / m_virt, F / m_virt - s = h + eqn1 = np.rad2deg(y[3] / r_lon) + eqn2 = np.rad2deg(y[4] / r_lat) + eqn3 = y[5] + 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] + 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 - a_x = Drag_x/m_virt - a_y = Drag_y/m_virt - a_z = F/m_virt + return [eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10] - # DIFFERENTIAL EQUATIONS +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 # [%] - dx = dt * v_x + 0.5 * a_x * dt ** 2 # lon - dy = dt * v_y + 0.5 * a_y * dt ** 2 # lat +t_max = 100000 #0 +dt = 1.0 +t = np.arange(2.0, t_max, dt) - lon = lon + np.rad2deg(np.arctan(dx/((6371229.0 + h) * np.cos(np.deg2rad(lat))))) - lat = lat + np.rad2deg(np.arctan(dy/(6371229.0 + h))) +m_gas = ((m_pl + m_bal + m_film) * (FreeLift/100 + 1))/(R_gas/R_air - 1) - v_xn = v_x + dt * a_x - v_yn = v_y + dt * a_y +def ending(t, y): return y[2] - s_n = s + dt * v_z + 0.5 * a_z * dt ** 2 - dh = s_n - s - v_n = v_z + dt * a_z +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:") +print(sol.message) - T_gn = T_gas + dt * (Q_ConvInt / (gamma * c_v * m_gas) - g * R_gas * T_gas * dh / (c_v * gamma * T_air * R_air)) - T_en = T_film + (Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film) * dt +res = sol.y[2, :] - T_g = T_gn - T_e = T_en - s = s_n - v_x = v_xn - v_y = v_yn - v_z = v_n - - h = s - T_gas = T_g - T_film = T_e -# v_z = v - - t += dt # time increment - utc = Time(start_utc) + t * unit.second - - - -#print(len(lon_list), len(lat_list)) -#""" -ax = plt.axes(projection=ccrs.Mollweide()) -#ax.coastlines() -ax.stock_img() -plt.plot(start_lon, start_lat, 'rx', transform=ccrs.Geodetic()) -plt.plot(lon_list, lat_list, 'k--', transform=ccrs.Geodetic()) -plt.plot(comp_lon, comp_lat, 'r-', transform=ccrs.Geodetic()) -#plt.xticks() -plt.show() -#""" - - -### WORKS! -### -###dset = xr.open_dataset("test2021.nc") -### -###print(dset['t'][1][0]) # [time] [level] -### -#fig = plt.figure() #figsize=[120,50]) -### -#ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) -### -###dset['t'][1][0].plot(ax=ax, cmap='jet', transform=ccrs.PlateCarree()) -###ax.coastlines() -###plt.show() - -#print(len(ERAlat)) -#print(len(ERAlon)) -#print(len(ERAz)) - - - - - -#""" +end = datetime.now() +print(end-start) plt.plot(comp_time, comp_height, 'r--', label='PoGo Flight Test') -#plt.plot(t_list, Re_list, 'k-', label="Reynolds-Number") -#plt.plot(t_list, v_list, 'k--', label='Ballon v_rel') -plt.plot(t_list, h_list, 'k-', label='Balloon Altitude') -#plt.plot(t_list, p_list, 'r-', label='Air Pressure [Pa]') -#plt.plot(t_list, Temp_list, 'b-', label='Air Temperature [K]') -#plt.plot(t_list, rho_list, 'g-', label='Air Density in [kg/m$^3$]') +plt.plot(t, res) plt.xlabel('time in s') plt.ylabel('Balloon Altitude in m') -plt.legend() + +plt.xlim(0, t_max) plt.show() -#""" + +plt.clf() +ax = plt.axes(projection=ccrs.Mollweide()) +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