дз

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Используемые библиотеки \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

import numpy as np
from scipy.optimize import root

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Используемые функции \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

# Функции атмосферы
def geopotencial_v():
    return (r * h) / (r + h)  # Геопотенциальная высота
def yskorenie_s_p():
    return g0 * (r / (r + h)) ** 2  # Ускорение свободного падения
def parametry_t():
    H = geopotencial_v()
    # Параметры, для определения температуры и давления
    if (H >= -2001) and (H < 0):     # 3-4 стр. ГОСТа
        Hx = -2000
        batm = -0.0065
        Tx = 301.15
        px = 127774
    elif (H >= 0) and (H <= 11000):  # 5-15 стр. ГОСТа
        Hx = 0
        batm = -0.0065
        Tx = 288.15
        px = 101325
    elif (H > 11000) and (H < 20000):  # 15-24 стр. ГОСТа
        Hx = 11000
        Tx = 216.65
        batm = 0
        px = 22632
    elif (H >= 20000) and (H < 32000):  # 25-36 стр. ГОСТа
        Hx = 20000
        batm = 0.001
        Tx = 216.65
        px = 5474.87
    elif (H >= 32000) and (H < 47000):  # 37-43 стр. ГОСТа
        Hx = 32000
        batm = 0.0028
        Tx = 228.65
        px = 868.014
    elif (H >= 47000) and (H < 51000):  # 43-46 стр. ГОСТа
        Hx = 47000
        Tx = 270.65
        batm = 0
        px = 110.906
    elif (H >= 51000) and (H < 71000):  # 45-52 стр. ГОСТа
        Hx = 51000
        batm = -0.0028
        Tx = 270.650
        px = 66.9384
    elif (H >= 71000) and (H < 85000):  # 51-53 и 159 стр. ГОСТа
        Hx = 71000
        batm = -0.002
        Tx = 214.65
        px = 3.95639
    return dict(Hx = Hx, batm = batm, Tx = Tx, px = px)
def temperatura():
    p = parametry_t()
    H = geopotencial_v()
    T = p["Tx"] + p["batm"] * (H - p["Hx"])       # Температура
    return T
def davlenie():
    p = parametry_t()
    H = geopotencial_v()
    T = temperatura()
    if p["batm"] != 0:  # Давление
        p = 10 ** (np.log10(p["px"]) - (g0 / (p["batm"] * R)) * np.log10(T / p["Tx"]))
    else:
        p = 10 ** (np.log10(p["px"]) - ((0.434294 * g0 * (H - p["Hx"])) / (T * R)))
    return p
def plotnost():
    p = davlenie()
    T = temperatura()
    ro = p / (R * T)  # Плотность
    return ro
def skorost_z():
    T = temperatura()
    a = 20.046796 * np.sqrt(T)  # Скорость звука
    return a
def vivod():
    H = geopotencial_v()
    g = yskorenie_s_p()
    T = temperatura()
    p = davlenie()
    ro = plotnost()
    a = skorost_z()
    print("Геометрическая высота: " + "{:.0f}".format(h) + ' м')
    print("Геопотенциальная высота: " + "{:.0f}".format(H) + ' м')
    print("Ускорение свободного падения: "+str(round(g,3))+' м/с^2')
    print("Температура: "+str(round(T/100,3))+' * 10^2 K')
    print("Давление: "+str(round(p/10000,3))+' * 10^4'+' Па')
    print('Плотность: '+str(round(ro*10,3))+' * 10^-1'+' кг/м^3')
    print("Скорость звука: "+str(round(a/100,3))+' * 10^2 м/с')
    print(' ')
# Функции для расчета:
# плоского СУ

def findShockWaveAngle(betta: float) -> float:  # Функция поиска угла скачка уплотнения

    def ugol_SY(tettac, bet):  # Трансцендентное уравнение
        return np.tan(tettac) / np.tan(tettac - bet) - ((k + 1) * M ** 2 * (np.sin(tettac)) ** 2) / (
                    2 + (k - 1) * M ** 2 * (np.sin(tettac)) ** 2)

    # Генератор листа численных решений функции, с шагом 0.1. Решение находится 50 раз. Округление до 4 знаков после запятой
    res = [round(root(ugol_SY, x0=i * 0.1, args=(betta)).x[0], 4) for i in range(50)]

    # Удалены повторяющиеся, отсортировано по возрастанию
    listRoots = sorted(list(set(res)))

    #print(res) # Раскомментировать, чтобы посмотреть вывод всех корней
    # Поиск наименьшего положительного корня
    closestRoot = None
    for listRoot in listRoots:
        if listRoot > 0 and (closestRoot is None or listRoot < closestRoot):
            closestRoot = listRoot
    shockWaweAngle = closestRoot
    return shockWaweAngle  # Угол СУ
# Касательной скорости
def findVr(θ_c):
    return M * skorost_z() * np.cos(θ_c)
# Нормальной скорости
def findVt(θ_c):
    return - findVr(θ_c) * np.tan(θ_c) * (
            (2 + (k - 1) * M ** 2 * np.sin(θ_c) ** 2) / ((k + 1) * M ** 2 * np.sin(θ_c) ** 2))
# Метода последовательных приближений
def V_r_next(V_r_previous, V_θ_previous, shag):
    return V_r_previous + V_θ_previous * shag
def V_tetta_next(θ_c, V_r_previous, V_θ_previous, shag):
    T_0 = T * (1 + ((k - 1) / 2) * M ** 2)
    a_kr = (k * R * T_0) - ((k - 1) / 2) * (V_r_previous ** 2 - V_θ_previous ** 2)
    dV0_dtetta = (-V_θ_previous / np.tan(θ_c) + V_r_previous * ((V_θ_previous ** 2) / a_kr - 2)) / (
            1 - (V_θ_previous ** 2) / a_kr)
    return V_θ_previous + dV0_dtetta * shag
def calc_block(θ_c, V_r, V_θ, shag):
    V_r = V_r_next(V_r, V_θ, shag)
    V_θ = V_tetta_next(θ_c, V_r, V_θ, shag)
    θ_c += shag
    return V_r, V_θ, θ_c
def ConeParametrs():
    coef = 0.85
    shag = -0.00001
    θ_c_plosk = findShockWaveAngle(β_k)
    θ_c = coef * θ_c_plosk
    while round(np.degrees(θ_c), 2) != np.degrees(β_k):
        θ_c = coef * findShockWaveAngle(β_k)
        V_r = findVr(θ_c)
        V_θ = findVt(θ_c)
        #print(θ_c, ' ',V_r, ' ', V_θ)
        while V_θ < 0:
            V_r = V_r_next(V_r, V_θ, shag)
            V_θ = V_tetta_next(θ_c, V_r, V_θ, shag)
            θ_c += shag
        #print('coef', coef, '|| β_n = ', round(np.degrees(θ_c), 2), '| V_r = ', V_r, '| V_θ = ', V_θ)
        if round(θ_c, 2) > β_k:
            coef -= 0.0001
        elif round(θ_c, 2) < β_k:
            coef += 0.0001
    θ_k = coef * θ_c_plosk
    β_n = θ_c
    V_k = V_r
    return θ_k, β_n, V_k # Угол θ корнуса, Расчетный угол конуса, # Нормальная скорость

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Исходные данные \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

N = 11
dm = 1.4
P = 0
if (N > 0) and (N < 6):
    P += 1
elif (N >= 6) and (N < 11):
    P += 2
elif (N >= 11) and (N < 16):
    P += 3
elif (N >= 16) and (N < 20):
    P += 4
M0 = 2
H0 = 3.5
βk0 = 7.5
x_cm = 0.65
β_k = np.radians(βk0 + 2.5 * (N - 2 * P))  # Угол конуса
lk = dm / (2 * np.tan(np.radians(β_k)))             # Длина конуса
Xcm = x_cm * lk                                     # Координата центра масс
l_s_chertoy_k = (2 * lk) / 3                        # Длина затупленного конуса

# Данные для атмосферы
h = (H0 + 0.2 * N) * 1000                           # Геометрическая высота
M = M0 + (P - 1)                                    # Мах
r = 6356767
g0 = 9.80665
R = 287.05287
k = 1.4

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Определение параметров стандартной атмосферы \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

H = geopotencial_v()
g = yskorenie_s_p()
T = temperatura()
p = davlenie()
ro = plotnost()
a = skorost_z()
V = a * M  # Скорость потока
vivod()
print("Мах: " + str(round(M ,3)) )
print("Скорость V: " + str(round(V,3)) + ' м/с')
print("Угол полураствора конуса: " + str(round(β_k, 3)) + ' рад; ' + str(round(np.degrees(β_k), 2)) + ' град')
print("Угол θ плоск: " + str(round(findShockWaveAngle(β_k), 3)) + ' рад; ' + str(round(np.degrees(findShockWaveAngle(β_k)), 2)) + ' град')

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Метод последовательных приближений \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

list_parametrov = ConeParametrs()
θ_c = list_parametrov[0] # Угол СУ θ корнуса
β_n = list_parametrov[1] # Расчетный угол конуса
V_k = list_parametrov[2] # Нормальная скорость
print("Угол θ корнуса: " + str(round(θ_c, 3)) + ' рад; ' + str(round(np.degrees(θ_c), 3)) + ' град')
print("Расчетный угол конуса: " + str(round(β_n, 3)) + ' рад; ' + str(round(np.degrees(β_n), 3)) + ' град')
print("Нормальная скорость V_r: " + str(round(V_k,3)) + ' м/с')

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Расчёт p̅к, Mк, Cx \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

T_0 = T * (1 + ((k - 1) / 2) * M ** 2)
T_k = T_0 * (1 - ((k - 1) * V_k ** 2 / (2 * k * R * T_0)))
M_k = V_k / np.sqrt(k * R * T_k)
p0_shtrih = p * ((k+1)/(2*k*M**2*(np.sin(θ_c))**2-(k-1)))**(1/(k-1)) * ((k+1)*M**2*(np.sin(θ_c))**2*((k-1)*M**2+2)/(2 * ((k-1)*M**2*(np.sin(θ_c))**2 + 2)))**(k/(k-1)) 
P_k = p0_shtrih * (1 + ((k - 1) * M_k ** 2)/ 2) ** ((- k)/(k - 1))
p̅k = (P_k - p) * 2 / (ro * V ** 2)
#print("T_0",T_0)
#print("T",T)
#print("k",k)
#print("M",M)
#print("T_k",T_k)
#print("M_k",M_k)
#print("p0_shtrih",p0_shtrih)
#print("P_k",P_k)
print("Коэффициент давления p̅k: " + str(round(p̅k ,3)))
print("Мах за СУ расчетный: " + str(round(M_k ,3)) )

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Расчёт по приближенным зависимостям \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

# Погрешность по 1 способу
θ_c_emp = np.arcsin(np.sqrt((1 / M ** 2) + (k + 1) * (np.sin(β_k) ** 2) / 2))
p̅k_1sposob = 0.002 * (0.8 + (1 / M ** 2)) * np.degrees(β_k) ** 1.7
# Погрешность по 2 способу
a_shtrih = np.sqrt(M-1)
p̅k_2sposob = 4 * (np.sin(β_k) ** 2) * (2.5 + 8*a_shtrih*np.sin(β_k))/(1+16*a_shtrih*np.sin(β_k))
# Погрешность по 3 способу
a_shtrih = np.sqrt(M-1)
y = 0.1*np.log10(a_shtrih*np.sin(β_k))
x = 0.18145 - 2.0923 * y + 9.092 * y**2 + 6.876 * y**3 - 52.225 * y**4 - 97.1 * y**5
p̅k_3sposob = 2 * (np.sin(β_k) ** 2) * np.exp(x)
V_c = np.sqrt(findVr(θ_c_emp)**2 + findVt(θ_c_emp)**2)
a_c = np.sqrt((k * R * T_0) - ((k-1)/2)*V_c**2)
M_c = V_c / a_c
M_k_emp = M_c * 0.98
#print("Vcr",findVr(θ_c_emp))
#print("Vcθ",findVt(θ_c_emp))
#print("Vc",V_c)
#print("a со звездой",np.sqrt(a_z_2))
#print("a_c",a_c)
#print("M_c",M_c)

print("Угол θ приближенный: " + str(round(θ_c_emp, 3)) + ' рад; ' + str(round(np.degrees(θ_c_emp), 3)) + ' град')
print("Коэффициент давления p̅k (1 способ): " + str(round(p̅k_1sposob ,10)) )
print("Коэффициент давления p̅k (2 способ): " + str(round(p̅k_2sposob ,10)) )
print("Коэффициент давления p̅k (3 способ): " + str(round(p̅k_3sposob ,10)) )
print("Мах за СУ по эмпирический: " + str(round(M_k_emp ,3)) )