112 lines
3.4 KiB
Python
112 lines
3.4 KiB
Python
from input.natural_constants import *
|
|
import numpy as np
|
|
from datetime import datetime
|
|
import time
|
|
from astropy.time import Time
|
|
|
|
from astropy.time import TimeISO
|
|
class TimeYearDayTimeCustom(TimeISO):
|
|
"""
|
|
day-of-year as "<DOY>".
|
|
The day-of-year (DOY) goes from 001 to 365 (366 in leap years).
|
|
The allowed subformat is:
|
|
- 'doy': day of year
|
|
"""
|
|
name = 'doy' # Unique format name
|
|
subfmts = (('doy',
|
|
'%j',
|
|
'{yday:03d}'),
|
|
('doy',
|
|
'%j',
|
|
'{yday:03d}'),
|
|
('doy',
|
|
'%j',
|
|
'{yday:03d}'))
|
|
|
|
def AirMass(x, p_0, ELV, h):
|
|
p_air = x
|
|
ELV_rad = np.deg2rad(ELV) # convert ELV from degree to radian
|
|
Dip = np.arccos(R_E / (R_E + h)) # Dip in radian
|
|
|
|
if ELV_rad >= -Dip and ELV_rad < 0:
|
|
res = p_air/p_0 * (1 + ELV_rad/Dip) - 70 * ELV_rad/Dip
|
|
else:
|
|
res = (p_air/p_0) * ((1229 + (614 * np.sin(ELV_rad)) ** 2) ** (1/2) - 614 * np.sin(ELV_rad))
|
|
|
|
return res
|
|
|
|
#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))
|
|
#
|
|
### NEEDED:
|
|
# Diameter = ...
|
|
# V_relative = ...
|
|
# rho_air = ...
|
|
# rho_gas = ...
|
|
# T_film = ...
|
|
# T_air = ...
|
|
# T_gas = ...
|
|
#
|
|
#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) * Diameter ** 3) / (T_air * my_air ** 2)
|
|
#Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25
|
|
#HC_free = Nu_air * k_air / Diameter
|
|
#Re = V_relative * Diameter * rho_air / my_air
|
|
#HC_forced = k_air / Diameter * (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)
|
|
#
|
|
#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_effective * (T_air - T_film)
|
|
#Q_ConvInt = HC_internal * A_effective * (T_film - T_gas)
|
|
#
|
|
#HalfConeAngle = np.arcsin(R_E/(R_E + h))
|
|
#ViewFactor = (1 - np.cos(HalfConeAngle))/2
|
|
#
|
|
#
|
|
#
|
|
#RoC = -V_z
|
|
#
|
|
#dT_gas = (Q_ConvInt/(gamma * c_v * M_gas) - (gamma - 1)/gamma * rho_air * 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
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|