BASTET/main.py

import numpy as np
import math
import xarray as xr
import matplotlib.pyplot as plt
from scipy.spatial import cKDTree
from scipy import interpolate
import cartopy.crs as ccrs
from datetime import datetime
start = datetime.now()
print(start)
begin_time = datetime.now()
from models.sun import sun_angles_analytical
from models.sun import sun_angles_astropy
from models.drag import drag  #, c_d
from input.user_input import *
from input.natural_constants import *
from astropy.time import Time
import astropy.units as unit
from models.thermal import AirMass
from netCDF4 import Dataset
import pandas as pd

data = pd.read_excel(r'C:\Users\marcel\PycharmProjects\MasterThesis\Mappe1.xls', sheet_name='Tabelle3')

comp_time = pd.DataFrame(data, columns= ['Time'])
comp_height = pd.DataFrame(data, columns= ['Height'])
comp_lat = pd.DataFrame(data, columns= ['Latitude'])
comp_lon = pd.DataFrame(data, columns= ['Longitude'])

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]

cd_interp1d = interpolate.interp1d(Reynolds, drag_coeff)

Re_list = []

def cd_func(Re):
    w = np.log10(Re)

    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

    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)))  # !




lat = start_lat   # deg
lon = start_lon   # deg

t = 0  # simulation time in seconds
h = start_height
utc = Time(start_utc)
epoch_diff = (Time(start_utc).jd - Time('1900-01-01 00:00:00.0').jd) * 24.000000
t_epoch = epoch_diff



xt, yt, zt, tt = transform(lon, lat, t)  # test coordinates

d, inds = tree.query(np.column_stack((xt, yt, zt, tt)), k=8)  # longitude, latitude, time in h  # !
w = 1.0 / d[0] ** 2

lat_sel = np.unravel_index(inds[0], lon_era2d.shape)[0]
lon_sel = np.unravel_index(inds[0], lon_era2d.shape)[1]
time_sel = np.unravel_index(inds[0], lon_era2d.shape)[2]

interp4d_height = np.ma.dot(w, ERAz[time_sel, :, lat_sel, lon_sel]) / np.sum(w)
interp4d_temp = np.ma.dot(w, ERAtemp[time_sel, :, lat_sel, lon_sel]) / np.sum(w)
interp4d_vw_x = np.ma.dot(w, vw_x[time_sel, :, lat_sel, lon_sel]) / np.sum(w)
interp4d_vw_y = np.ma.dot(w, vw_y[time_sel, :, lat_sel, lon_sel]) / np.sum(w)
interp4d_vw_z = np.ma.dot(w, vw_z[time_sel, :, lat_sel, lon_sel]) / np.sum(w)

pressure_hPa = np.array([1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 125, 150, 175, 200,225, 250, 300, 350, 400,
                        450, 500, 550, 600, 650, 700, 750, 775, 800, 825, 850, 875, 900, 925, 950, 975, 1000])

pressure = 100 * pressure_hPa
# print(pressure)


# height_interp1d = interpolate.interp1d(interp4d_height, pressure)
temp_interp1d = interpolate.interp1d(interp4d_height, interp4d_temp)
press_interp1d = interpolate.interp1d(interp4d_height, pressure)
vw_x_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_x)
vw_y_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_y)
vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z)


#u = 0
#v = 0
#w = 0
T_air = temp_interp1d(h)
p_air = press_interp1d(h)
rho_air = p_air/(R_air * T_air)

u = vw_x_interp1d(h)
v = vw_y_interp1d(h)
w = -1 / g * vw_z_interp1d(h) * R_air * T_air / p_air


v_x = 0
v_y = 0
v_z = 0

v_rel = 0

T_gas = T_air
T_film = T_gas

t_list = []
h_list = []
v_list = []
lat_list = []
lon_list = []
p_list = []
Temp_list = []
rho_list = []

while t <= t_end and h >= 0:
    t_list.append(t)
    h_list.append(h)
    v_list.append(v_rel)
    lat_list.append(lat)
    lon_list.append(lon)

    p_list.append(p_air)
    Temp_list.append(T_air)
    rho_list.append(rho_air)

    t_epoch = epoch_diff + t/3600

    xt, yt, zt, tt = transform(lon, lat, t_epoch)  # current balloon coordinates in cartesian coordinates: x [m], y [m], z [m], time [h since 1900-01-01 00:00:00.0]

    d, inds = tree.query(np.column_stack((xt, yt, zt, tt)), k=8)  # longitude, latitude, time in h  # !
    w = 1.0 / d[0] ** 2

    lat_sel = np.unravel_index(inds[0], lon_era2d.shape)[0]
    lon_sel = np.unravel_index(inds[0], lon_era2d.shape)[1]
    time_sel = np.unravel_index(inds[0], lon_era2d.shape)[2]

    interp4d_height = np.ma.dot(w, ERAz[time_sel, :, lat_sel, lon_sel]) / np.sum(w)
    interp4d_temp = np.ma.dot(w, ERAtemp[time_sel, :, lat_sel, lon_sel]) / np.sum(w)
    interp4d_vw_x = np.ma.dot(w, vw_x[time_sel, :, lat_sel, lon_sel]) / np.sum(w)
    interp4d_vw_y = np.ma.dot(w, vw_y[time_sel, :, lat_sel, lon_sel]) / np.sum(w)
    interp4d_vw_z = np.ma.dot(w, vw_z[time_sel, :, lat_sel, lon_sel]) / np.sum(w)

    press_interp1d = interpolate.interp1d(interp4d_height, pressure)
    temp_interp1d = interpolate.interp1d(interp4d_height, interp4d_temp)
    vw_x_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_x)
    vw_y_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_y)
    vw_z_interp1d = interpolate.interp1d(interp4d_height, interp4d_vw_z)


    try:
        p_air = press_interp1d(h)
        T_air = temp_interp1d(h)
    except:
        if (h > interp4d_height[0]):
            h = interp4d_height[0]
        else:
            h = interp4d_height[-1]

        p_air = press_interp1d(h)
        T_air = temp_interp1d(h)

        #print("height")
        #print(h)
        #print("temperature")
        #print(temp_interp1d(h))
        #print("pressure")
        #print(press_interp1d(h))
        #print("vw_x")
        #print(vw_x_interp1d(h))

    p_gas = p_air
    rho_air = p_air / (R_air * T_air)
    u = vw_x_interp1d(h)
    v = vw_y_interp1d(h)
    w = -1/g * vw_z_interp1d(h) * R_air * T_air / p_air

    v_relx = u - v_x
    v_rely = v - v_y
    v_relz = w - v_z

    v_rel = (v_relx ** 2 + v_rely ** 2 + v_relz ** 2) ** (1/2)

    rho_gas = p_gas/(R_gas * T_gas)  # calculate gas density through ideal(!) gas equation

    V_b = m_gas/rho_gas  # calculate balloon volume from current gas mass and gas density

    if V_b > V_design:
       V_b = V_design
       m_gas = V_design * rho_gas
    else:
        pass

    m_gross = m_pl + m_film
    m_tot = m_pl + m_film + m_gas
    m_virt = m_tot + c_virt * rho_air * V_b

    d_b = 1.383 * V_b ** (1/3)  # calculate diameter of balloon from its volume
    L_goreB = 1.914 * V_b ** (1/3)
    h_b = 0.748 * d_b
    A_surf = 4.94 * V_b ** (2/3)
    A_surf1 = 4.94 * V_design ** (2/3) * (1 - np.cos(np.pi * L_goreB/L_goreDesign))
    A_eff = 0.65 * A_surf + 0.35 * A_surf1
    A_top = np.pi/4 * d_b ** 2

    AZ, ELV = sun_angles_analytical(lat, lon, utc)

    A_proj = A_top * (0.9125 + 0.0875 * np.cos(np.pi - 2 * np.deg2rad(ELV)))

    # CALCULATIONS FOR THERMAL MODEL

    if ELV >= -(180 / np.pi * np.arccos(R_E / (R_E + h))):
        tau_atm = 0.5 * (np.exp(-0.65 * AirMass(p_air, p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air, p_0, ELV, h)))
        tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air / p_0) + np.exp(-0.095 * p_air / p_0))
    else:
        tau_atm = 0
        tau_atmIR = 0

    doy = int(utc.doy)

    MA = (357.52911 + 0.98560028 * (utc.jd - 2451545)) % 360  # in degree, reference: see folder "literature"
    TA = MA + 2 * e * np.sin(np.deg2rad(MA)) + 5 / 4 * e ** 2 * np.sin(np.deg2rad(2 * MA))
    I_Sun = 1367.5 * ((1 + e * np.cos(np.deg2rad(TA))) / (1 - e ** 2)) ** 2
    I_SunZ = I_Sun * tau_atm
    q_sun = I_SunZ
    q_IRground = epsilon_ground * sigma * T_ground ** 4
    q_IREarth = q_IRground * tau_atmIR

    if ELV <= 0:
        q_Albedo = 0
    else:
        q_Albedo = Albedo * I_Sun * np.sin(np.deg2rad(ELV))

    my_air = (1.458 * 10 ** -6 * T_air ** 1.5) / (T_air + 110.4)
    my_gas = 1.895 * 10 ** -5 * (T_gas / 273.15) ** 0.647
    k_air = 0.0241 * (T_air / 273.15) ** 0.9
    k_gas = 0.144 * (T_gas / 273.15) ** 0.7
    Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air
    Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas

    Gr_air = (rho_air ** 2 * g * np.abs(T_film - T_air) * d_b ** 3) / (T_air * my_air ** 2)
    Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25
    HC_free = Nu_air * k_air / d_b
    Re = np.abs(v_relz) * d_b * rho_air / my_air

    Re_list.append(Re)

    HC_forced = k_air / d_b * (2 + 0.41 * Re ** 0.55)
    HC_internal = 0.13 * k_gas * ((rho_gas ** 2 * g * np.abs(T_film - T_gas) * Pr_gas) / (T_gas * my_air ** 2)) ** (
                1 / 3)
    HC_external = np.maximum(HC_free, HC_forced)

    HalfConeAngle = np.arcsin(R_E / (R_E + h))
    ViewFactor = (1 - np.cos(HalfConeAngle)) / 2

    Q_Sun = alpha_VIS * A_proj * q_sun * (1 + tau_VIS / (1 - r_VIS))
    Q_Albedo = alpha_VIS * A_surf * q_Albedo * ViewFactor * (1 + tau_VIS / (1 - r_VIS))
    Q_IREarth = alpha_IR * A_surf * q_IREarth * ViewFactor * (1 + tau_IR / (1 - r_IR))
    Q_IRfilm = sigma * epsilon * alpha_IR * A_surf * T_film ** 4 * 1 / (1 - r_IR)
    Q_IRout = sigma * epsilon * A_surf * T_film ** 4 * (1 * tau_IR / (1 - r_IR))
    Q_ConvExt = HC_external * A_eff * (T_air - T_film)
    Q_ConvInt = HC_internal * A_eff * (T_film - T_gas)

    try:
        c_d = cd_func(Re)
    except:
        print(Re)

    D = drag(c_d, rho_air, d_b, v_rel)  # calculate drag force

    if v_rel == 0:
        Drag_x = 0
        Drag_y = 0
        Drag_z = 0
    else:
        Drag_x = D * v_relx / v_rel
        Drag_y = D * v_rely / v_rel
        Drag_z = D * v_relz / v_rel

    I = g * V_b * (rho_air - rho_gas)  # calculate gross inflation
    W = g * m_gross  # calculate weight (force)
    F = I - W + Drag_z

    # RoC = -v_z
    # dT_gas = (Q_ConvInt / (gamma * m_gas * c_v) - (gamma - 1) / gamma * rho_air(h) * g / (rho_gas * R_gas) * RoC) * dt
    # dT_film = ((Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film)) * dt
    #  v_w = w
    #  v = v_z

    s = h

    a_x = Drag_x/m_virt
    a_y = Drag_y/m_virt
    a_z = F/m_virt

    # DIFFERENTIAL EQUATIONS

    dx = dt * v_x + 0.5 * a_x * dt ** 2   # lon
    dy = dt * v_y + 0.5 * a_y * dt ** 2   # lat

    lon = lon + np.rad2deg(np.arctan(dx/((6371229.0 + h) * np.cos(np.deg2rad(lat)))))
    lat = lat + np.rad2deg(np.arctan(dy/(6371229.0 + h)))

    v_xn = v_x + dt * a_x
    v_yn = v_y + dt * a_y

    s_n = s + dt * v_z + 0.5 * a_z * dt ** 2
    dh = s_n - s
    v_n = v_z + dt * a_z

    T_gn = T_gas + dt * (Q_ConvInt / (gamma * c_v * m_gas) - g * R_gas * T_gas * dh / (c_v * gamma * T_air * R_air))
    T_en = T_film + (Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout) / (c_f * m_film) * dt

    T_g = T_gn
    T_e = T_en
    s = s_n
    v_x = v_xn
    v_y = v_yn
    v_z = v_n

    h = s
    T_gas = T_g
    T_film = T_e
#   v_z = v

    t += dt  # time increment
    utc = Time(start_utc) + t * unit.second



#print(len(lon_list), len(lat_list))
#"""
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))





#"""
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.xlabel('time in s')
plt.ylabel('Balloon Altitude in m')
plt.legend()
plt.show()
#"""