From d107aaae2c4fe35a449514d7d3f8dafcc2833f26 Mon Sep 17 00:00:00 2001 From: Marcel Frommelt Date: Mon, 15 Feb 2021 21:40:19 +0900 Subject: [PATCH] full implementation of ERA5-netCDF-file-interface --- DataRequestNew.py | 55 +++++++++ main_atmodel.py | 278 ++++++++++++++++++++++++++++++++++++++++++++++ main_netCDF.py | 277 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 610 insertions(+) create mode 100644 DataRequestNew.py create mode 100644 main_atmodel.py create mode 100644 main_netCDF.py diff --git a/DataRequestNew.py b/DataRequestNew.py new file mode 100644 index 0000000..d6447c0 --- /dev/null +++ b/DataRequestNew.py @@ -0,0 +1,55 @@ +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_atmodel.py b/main_atmodel.py new file mode 100644 index 0000000..dacd3d2 --- /dev/null +++ b/main_atmodel.py @@ -0,0 +1,278 @@ +import numpy as np +import matplotlib.pyplot as plt +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.simple_atmosphere import T_air, p_air, rho_air +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 + +#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'][:] +#vw_x = data.variables['u'][:] # v_x +#vw_y = data.variables['v'][:] # v_y +#vw_z = data.variables['w'][:] # v_z + +# t: simulation time (0 .. t_max) in [s] +# utc = Time(start_utc) + t * u.second # current time (UTC) + +lat = start_lat # deg +lon = start_lon # deg + +# date = Time('2020-12-13 12:00:00.000') + + +# INITIALISATION + +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 + +#sq_diff_lat = (ERAlat - lat) ** 2 +#sq_diff_lon = (ERAlon - lon) ** 2 +#sq_diff_time = (ERAtime - t_epoch) ** 2 + +#min_index_lat = sq_diff_lat.argmin() +#min_index_lon = sq_diff_lat.argmin() +#min_index_time = sq_diff_time.argmin() + +#sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 + +#min_index_h = sq_diff_h.argmin() + +#h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] +#p_grid = ERApress[min_index_h] * 100 +#T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] +#u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] +#v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] +#w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid + +#v = v_grid +#u = u_grid +#w = w_grid +#T_air = T_grid +#p_air = p_grid +#rho_air = p_air/(R_air * T_air) + +v_x = 0 +v_y = 0 +v_z = 0 + +T_gas = T_air(h) +T_film = T_gas + +t_list = [] +h_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) + lat_list.append(lat) + lon_list.append(lon) + + p_list.append(p_air(h)) + Temp_list.append(T_air(h)) + rho_list.append(rho_air(h)) + + + #t_epoch = epoch_diff + t/3600 + + #sq_diff_lat = (ERAlat - lat) ** 2 + #sq_diff_lon = (ERAlon - lon) ** 2 + #sq_diff_time = (ERAtime - t_epoch) ** 2 + + #min_index_lat = sq_diff_lat.argmin() + #min_index_lon = sq_diff_lat.argmin() + #min_index_time = sq_diff_time.argmin() + + #sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 + + #min_index_h = sq_diff_h.argmin() + + #h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] + #p_grid = ERApress[min_index_h] * 100 + #T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] + #u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] + #v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] + #w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid + + # p_air = p_grid + p_gas = p_air(h) + # rho_air = p_air / (R_air * T_air) + u = 0 + v = 0 + w = 0 + + 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) + + p_gas = p_air(h) + 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(h) * 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 + + D = drag(c_d, rho_air(h), 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(h) - rho_gas) # calculate gross inflation + W = g * m_gross # calculate weight (force) + F = I - W + Drag_z + + 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(h), p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air(h), p_0, ELV, h))) + tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air(h) / p_0) + np.exp(-0.095 * p_air(h) / 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(h) ** 1.5) / (T_air(h) + 110.4) + my_gas = 1.895 * 10 ** -5 * (T_gas / 273.15) ** 0.647 + k_air = 0.0241 * (T_air(h) / 273.15) ** 0.9 + k_gas = 0.144 * (T_gas / 273.15) ** 0.7 + Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air(h) + Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas + + Gr_air = (rho_air(h) ** 2 * g * np.abs(T_film - T_air(h)) * d_b ** 3) / (T_air(h) * 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(h) / my_air + + 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(h) - T_film) + Q_ConvInt = HC_internal * A_eff * (T_film - T_gas) + + # 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(s) * 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 + +end = datetime.now() +print(end) +plt.plot(t_list, h_list, 'k-') +plt.plot(t_list, p_list, 'r-') +plt.plot(t_list, Temp_list, 'b-') +plt.plot(t_list, rho_list, 'g-') +plt.show() \ No newline at end of file diff --git a/main_netCDF.py b/main_netCDF.py new file mode 100644 index 0000000..15e7def --- /dev/null +++ b/main_netCDF.py @@ -0,0 +1,277 @@ +import numpy as np +import matplotlib.pyplot as plt +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 + +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'][:] +vw_x = data.variables['u'][:] # v_x +vw_y = data.variables['v'][:] # v_y +vw_z = data.variables['w'][:] # v_z + +# t: simulation time (0 .. t_max) in [s] +# utc = Time(start_utc) + t * u.second # current time (UTC) + +lat = start_lat # deg +lon = start_lon # deg + +# date = Time('2020-12-13 12:00:00.000') + + +# INITIALISATION + +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 + +sq_diff_lat = (ERAlat - lat) ** 2 +sq_diff_lon = (ERAlon - lon) ** 2 +sq_diff_time = (ERAtime - t_epoch) ** 2 + +min_index_lat = sq_diff_lat.argmin() +min_index_lon = sq_diff_lat.argmin() +min_index_time = sq_diff_time.argmin() + +sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 + +min_index_h = sq_diff_h.argmin() + +h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] +p_grid = ERApress[min_index_h] * 100 +T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] +u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] +v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] +w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid + +v = v_grid +u = u_grid +w = w_grid +T_air = T_grid +p_air = p_grid +rho_air = p_air/(R_air * T_air) + +v_x = 0 +v_y = 0 +v_z = 0 + +T_gas = T_air +T_film = T_gas + +t_list = [] +h_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) + 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 + + sq_diff_lat = (ERAlat - lat) ** 2 + sq_diff_lon = (ERAlon - lon) ** 2 + sq_diff_time = (ERAtime - t_epoch) ** 2 + + min_index_lat = sq_diff_lat.argmin() + min_index_lon = sq_diff_lat.argmin() + min_index_time = sq_diff_time.argmin() + + sq_diff_h = (ERAz[min_index_time, :, min_index_lat, min_index_lon] - h) ** 2 + + min_index_h = sq_diff_h.argmin() + + h_grid = ERAz[min_index_time, min_index_h, min_index_lat, min_index_lon] + p_grid = ERApress[min_index_h] * 100 + T_grid = ERAtemp[min_index_time, min_index_h, min_index_lat, min_index_lon] + u_grid = vw_x[min_index_time, min_index_h, min_index_lat, min_index_lon] + v_grid = vw_y[min_index_time, min_index_h, min_index_lat, min_index_lon] + w_grid = -1/g * vw_z[min_index_time, min_index_h, min_index_lat, min_index_lon] * R_air * T_grid / p_grid + + p_air = p_grid + p_gas = p_air + T_air = T_grid + rho_air = p_air / (R_air * T_air) + u = u_grid + v = v_grid + w = w_grid + + 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) + + p_gas = p_air + 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 + + 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 + + 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 + + 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) + + # 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 + +end = datetime.now() +print(end) +plt.plot(t_list, h_list, 'k-') +plt.plot(t_list, p_list, 'r-') +plt.plot(t_list, Temp_list, 'b-') +plt.plot(t_list, rho_list, 'g-') +plt.show()