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
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, 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
from netCDF4 import Dataset
import pandas as pd

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

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

print("Initialisiere Simulation...")
print("STARTORT: Längengrad %.4f, Breitengrad: %.4f" % (start_lon, start_lat))
print("STARTZEIT: " + str(start_utc) + " (UTC)")
#print("SIMULATION INITIALISIERT. Simulierte Flugdauer %.2f Stunden (oder %.2f Tage)" % (
#t_end / 3600, t_end / (3600 * 24)))
#print("GESCHÄTZE SIMULATIONSDAUER: %.2f Minuten" % (t_end / (3600 * 10)))


level_data = Dataset("level.nc", "r", format="NETCDF4")
single_data = Dataset("single.nc", "r", format="NETCDF4")

ERAtime = level_data.variables['time'][:]      # time

ERAlat = level_data.variables['latitude'][:]   # latitude (multi-level)
ERAlon = level_data.variables['longitude'][:]  # longitude (multi-level)
ERAlats = single_data.variables['latitude'][:]   # latitude (single level)
ERAlons = single_data.variables['longitude'][:]  # longitude (single level)


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



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

lon_era2d, lat_era2d = np.meshgrid(ERAlon, ERAlat)
lon_era2ds, lat_era2ds = np.meshgrid(ERAlons, ERAlats)

xs, ys, zs = transform(lon_era2d.flatten(), lat_era2d.flatten())
xs2, ys2, zs2 = transform(lon_era2ds.flatten(), lat_era2ds.flatten())

tree = cKDTree(np.column_stack((xs, ys, zs)))
tree2 = cKDTree(np.column_stack((xs2, ys2, zs2)))


# INITIALISATION TIME-COUNTERS


# MAYBE LATER
# t_pre = int(t_epoch)
# t_post = t_pre + 1
# t1 = np.searchsorted(ERAtime, t_pre, side='left')
# t2 = t1 + 1

def ERA5Data(lon, lat, h, t, deltaT_ERA):
    t_epoch = deltaT_ERA + t/3600

    t_pre = int(t_epoch)
    # t_post = t_pre + 1
    t_pre_ind = t_pre - ERAtime[0]
    t_post_ind = t_pre_ind + 1

    xt, yt, zt = transform(lon, lat)  # current coordinates

    d1, inds1 = tree.query(np.column_stack((xt, yt, zt)), k=4)  # longitude, latitude
    d2, inds2 = tree2.query(np.column_stack((xt, yt, zt)), k=visible_cells(h))  # longitude, latitude
    w1 = 1.0 / d1[0]  # ** 2
    w2 = 1.0 / d2[0] ** 2

    lat_ind1 = np.unravel_index(inds1[0], lon_era2d.shape)[0]
    lon_ind1 = np.unravel_index(inds1[0], lon_era2d.shape)[1]
    lat_ind2 = np.unravel_index(inds2[0], lon_era2ds.shape)[0]
    lon_ind2 = np.unravel_index(inds2[0], lon_era2ds.shape)[1]

    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

    interp4d_height_pre = np.ma.dot(w1, ERAz[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1)
    interp4d_height_post = np.ma.dot(w1, ERAz[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1)
    interp4d_height = (interp4d_height_post - interp4d_height_pre) * (t_epoch - t_pre) + interp4d_height_pre

    interp4d_vw_x_pre = np.ma.dot(w1, vw_x[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1)
    interp4d_vw_x_post = np.ma.dot(w1, vw_x[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1)
    interp4d_vw_x = (interp4d_vw_x_post - interp4d_vw_x_pre) * (t_epoch - t_pre) + interp4d_vw_x_pre

    interp4d_vw_y_pre = np.ma.dot(w1, vw_y[t_pre_ind, :, lat_ind1, lon_ind1]) / np.sum(w1)
    interp4d_vw_y_post = np.ma.dot(w1, vw_y[t_post_ind, :, lat_ind1, lon_ind1]) / np.sum(w1)
    interp4d_vw_y = (interp4d_vw_y_post - interp4d_vw_y_pre) * (t_epoch - t_pre) + interp4d_vw_y_pre

    interp4d_vw_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

    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

    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

    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

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

    if 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)
        u = vw_x_interp1d(h)
        v = vw_y_interp1d(h)
        w = -1 / g * vw_z_interp1d(h) * R_air * T_air / p_air

    rho_air = p_air / (R_air * T_air)

    return p_air, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, ssrd, tsr, ttr, al, tisr, strdc
# p_air, T_air, rho_air, u, v, w, cbh, tcc, lcc, mcc, hcc, ssr, strn, strd, ssrd, tsr, ttr, al, tisr


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)





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:
        c_duct = c_ducts
    else:
        c_duct = 0

    #if valving == True:  # opening valve process
    #    if c2 == 0:
    #        c2 = 1.0
    #        c2dot = 0
    #    elif c_valve < c2 <= 1.0:
    #        c2dot = (c_valve - 1.0) / t_open
    #    else:
    #        c2dot = 0
    #        c2 = c_valve

    #if 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)
    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)))  # 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)
    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
    I_Sun0 = I_Sun * tau_atm0
    ## q_sun = I_SunZ
    ## q_IRground = np.abs(str)  #epsilon_ground * sigma * T_ground ** 4
    ## q_IREarth = q_IRground * tau_atmIR

    tcc = 0

    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_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)
    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_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_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)

    c_d = Cd_model(Fr, Re, A_top)

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

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

    F = g * V_b * (rho_air - rho_gas) - g * m_gross + Drag_z  # gross inflation - weight + drag

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

    eqn1 = np.rad2deg(y[3] / r_lon)
    eqn2 = np.rad2deg(y[4] / r_lat)
    eqn3 = y[5]
    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

    return [eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10]

m_pl = 1728 + 172 + 500 # + 500  # mass of payload and flight chain in [kg]
m_bal = 0  # initial ballast mass in [kg]
m_film = 1897  # mass of balloon film in [kg]
FreeLift = 4  # [%]

t_max = 100000  #0
dt = 1.0
t = np.arange(2.0, t_max, dt)

m_gas = ((m_pl + m_bal + m_film) * (FreeLift/100 + 1))/(R_gas/R_air - 1)

def ending(t, y): return y[2]

sol = solve_ivp(fun=model, t_span=[2.0, t_max], y0=y0, t_eval=t, args=(m_pl, m_film, m_bal, c_virt), max_step=35.0)  #, events=ending.terminal)  # events=
print("hier:")
print(sol.message)

res = sol.y[2, :]

end = datetime.now()
print(end-start)
plt.plot(comp_time, comp_height, 'r--', label='PoGo Flight Test')
plt.plot(t, res)
plt.xlabel('time in s')
plt.ylabel('Balloon Altitude in m')

plt.xlim(0, t_max)
plt.show()

plt.clf()
ax = plt.axes(projection=ccrs.Mollweide())
ax.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()