import numpy as np
from scipy.optimize import fsolve
from scipy.optimize import root
from functions import atmosphere as atm
import data as d
from openpyxl import Workbook
import datetime


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

    def ugol_SY(tettac, betta):  # Трансцендентное уравнение
        return np.tan(tettac) / np.tan(tettac - betta) - ((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 findFictitiousAngle (M: float) -> float:  # Функция поиска фиктивного угла
    FictitiousAngle = np.sqrt((k + 1) / (k - 1)) * np.arctan(np.sqrt(((k - 1) * (M ** 2 - 1)) / (k + 1))) - np.arctan(np.sqrt(M ** 2 - 1))
    return FictitiousAngle  # Фиктивный угол
def findMachThroughFictitiousAngle(w: float) -> float:  # Функция поиска Маха, через фиктивный угол

    def phiktivnyi_u(M, w) -> float:  # Трансцендентное уравнение
        return np.sqrt((k + 1) / (k - 1)) * np.arctan(np.sqrt(((k - 1) * (M ** 2 - 1)) / (k + 1))) - np.arctan(np.sqrt(M ** 2 - 1)) - w

    # Генератор листа численных решений функции, с шагом 0.1. Решение находится 50 раз. Округление до 4 знаков после запятой
    res = [round(root(phiktivnyi_u, x0=i * 1, args=(w)).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
    Mach = closestRoot
    return Mach  # Мах
def findDeflectionAngle(p): # Функция поиска углов СУ 5 и 6, и угла отклонения
    tettac5,tettac6,delta = p
    return (
        np.tan(tettac5) / np.tan(tettac5 - betta - delta) - ((k + 1) * M3 ** 2 * (np.sin(tettac5)) ** 2) / (2 + (k - 1) * M3 ** 2 * (np.sin(tettac5)) ** 2),
        np.tan(tettac6) / np.tan(tettac6 - betta + delta) - ((k + 1) * M4 ** 2 * (np.sin(tettac6)) ** 2) / (2 + (k - 1) * M4 ** 2 * (np.sin(tettac6)) ** 2),
        p3 * ((2 * k * M3**2 * np.sin(tettac5)**2 - (k-1))/(k+1)) - p4 * ((2 * k * M4**2 * np.sin(tettac6)**2 - (k-1))/(k+1))
        )
def LongitudinalAndNormalForces(pi: float) -> float:  # Функция поиска продольной и нормальной сил
    X = (pi - p) * (c / 2)
    Y = (pi - p) * (b / 2)
    return X, Y # Продольная и нормальная силы

# ПС
def CriticalReynoldsNumber(Tr, M):
    Tst_otn = Tst / Tr  # Относительная температура стенки
    x = (Tst_otn - 1) / M ** 2
    y = 1 - 16 * x - 412.5 * x**2 - 35000 * x**3 - 375000 * x**4
    return criticalReynoldsNumber * y

def SpecificHeat_Cp(T: float) -> float:  # Функция поиска удельной теплоемкости
    return Cp * (T/T0[0]) ** fi  # Удельная теплоемкость
def DynamicViscosity_mu(T: float) -> float:  # Функция поиска динамической вязкости
    return mu * (T/T0[0]) ** n  # Динамическая вязкость
def CoeffThermalConductivity_Lambda(T: float) -> float:  # Функция поиска коэффициента теплопроводности
    return Lambda * (T/T0[1]) ** z  # Коэффициент теплопроводности
def  findPressureCoefficient (pi: float) -> float: # Функция поиска коэффициента давления
    return (pi - p) / q # Коэффициент давления

def findDefiningTemperature (lam_or_turb: float,T: float, M: float) -> float: # Функция поика температуры восстановления и определяющей температуры
    T_z = 10 ** 6
    T_z_previous = 0
    rl = 0.83
    rt = 0.88
    stepen = 0.5
    r = rl
    if (lam_or_turb == 1):
        stepen = 1 / 3
        r = rt
    while abs(T_z - T_z_previous) > (T_z * 0.01):
        T_z_previous = T_z
        Tr = T * (1 + r * (k - 1) * 0.5 * M ** 2)
        T_z = 0.5 * (Tst + T) + 0.22 * (Tr - T)
        Cpi = SpecificHeat_Cp(T_z)
        mui = DynamicViscosity_mu(T_z)
        Lambdai = CoeffThermalConductivity_Lambda(T_z)
        Prl = Cpi * mui / Lambdai
        r = Prl ** stepen
    return Tr, T_z


# Параметры потока в скачках уплотнения
def findFlowParameterPressure(p: float, M: float, tettac: float) -> float:  # Функция поиска давления
    FlowParameterPressure = p * ((2 * k * M**2 * np.sin(tettac)**2 - (k-1))/(k+1))
    return FlowParameterPressure  # Давление области
def findFlowParameterDensity(ro: float, M: float,tettac: float) -> float:  # Функция поиска плотности
    FlowParameterDensity = (ro * (k+1) * M**2 * np.sin(tettac)**2)/(2 + ((k-1) * M**2 * np.sin(tettac)**2))
    return FlowParameterDensity  # Плотность области
def findFlowParameterTemperature(T: float, p: float, ro: float, pi: float, roi: float) -> float: # Функция поиска температуры
    FlowParameterTemperature = (T * pi * ro)/(p * roi)
    return FlowParameterTemperature  # Температура области
def findFlowParameterVelocity(V: float, betta: float, tettac: float) -> float: # Функция поиска скорости
    FlowParameterVelocity = (V * np.cos(tettac))/(np.cos(tettac-betta))
    return FlowParameterVelocity  # Скорость области
def findFlowParameterАcceleration(T: float) -> float: # Функция поиска скорости
    FlowParameterАcceleration = np.sqrt(k * R * T)
    return FlowParameterАcceleration  # Скорость области
def findFlowParameterMach(V: float, a: float) -> float: # Функция поиска маха
    FlowParameterMach = V / a
    return FlowParameterMach  # Число маха области

# Параметры торможения
def findBrakingParameterPressure(p: float,M: float) -> float:  # Функция поиска давления торможения
    BrakingParameterPressure = p * (((k - 1) * M ** 2) / 2 + 1) ** (k / (k - 1))
    return BrakingParameterPressure  # Давление торможения
def findBrakingParameterDensity(ro: float, M: float) -> float:  # Функция поиска плотности торможения
     BrakingParameterDensity = ro * (((k - 1) * M ** 2) / 2 + 1) ** (1 / (k - 1))
     return BrakingParameterDensity  # Плотность торможения
def findBrakingParameterTemperature(T: float, M: float) -> float:  # Функция поиска температуры торможения
    BrakingParameterTemperature = T * (((k - 1) * M ** 2) / 2 + 1)
    return BrakingParameterTemperature  # Температура торможения

# Получение данных из data.py
M = d.M
alfa = (d.alfaDegrees * np.pi) / 180
betta = d.betta
b = d.chordLength
c = d.profileThickness
L = np.sqrt((c / 2) ** 2 + (b / 2) ** 2)  # длина грани в м
h = d.h
Tst = d.TempCooledWall
criticalReynoldsNumber = d.criticalReynoldsNumber
R = d.R
fi = d.fi
n = d.n
z = d.z
T0 = d.T0
Cp = d.Cp
mu = d.mu
Lambda = d.Lambda

# Получение данных из atmosphere.py
H = atm.geopotencial_v()
g = atm.yskorenie_s_p()
T = atm.temperatura()
p = atm.davlenie()
ro = atm.plotnost()
a = atm.skorost_z()

V = a * M  # Скорость потока

# \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 1 и 2 области \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

# Углы отклонения потока
betta1 = betta - alfa
betta2 = betta + alfa
# Поиск углов СУ
k = 1.4
tettac1 = findShockWaveAngle(betta1)
tettac2 = findShockWaveAngle(betta2)
# Определяем параметры потока в скачках уплотнения
p1 = findFlowParameterPressure(p,M,tettac1)
p2 = findFlowParameterPressure(p,M,tettac2)
ro1 = findFlowParameterDensity(ro,M,tettac1)
ro2 = findFlowParameterDensity(ro,M,tettac2)
T1 = findFlowParameterTemperature(T,p,ro,p1,ro1)
T2 = findFlowParameterTemperature(T,p,ro,p2,ro2)
V1 = findFlowParameterVelocity(V,betta1,tettac1)
V2 = findFlowParameterVelocity(V,betta2,tettac2)
a1 = findFlowParameterАcceleration(T1)
a2 = findFlowParameterАcceleration(T2)
M1 = findFlowParameterMach(V1, a1)
my1 = np.arcsin(1/M1)
M2 = findFlowParameterMach(V2, a2)
my2 = np.arcsin(1/M2)
# Определяем параметры торможения для участка
p01 = findBrakingParameterPressure(p1,M1)
p02 = findBrakingParameterPressure(p2,M2)
ro01 = findBrakingParameterDensity(ro1,M1)
ro02 = findBrakingParameterDensity(ro2,M2)
T01 = findBrakingParameterTemperature(T1,M1)
T02 = findBrakingParameterTemperature(T2,M2)

# \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 3 и 4 области \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

# Прямая задача течения Прандтля – Майера
betta3 = 2 * betta
betta4 = 2 * betta
# Фиктивные углы и Мах
w1 = findFictitiousAngle(M1)
w3 = w1 + betta3
w2 = findFictitiousAngle(M2)
w4 = w2 + betta4
M3 = findMachThroughFictitiousAngle(w3)
my3 = np.arcsin(1/M3)
M4 = findMachThroughFictitiousAngle(w4)
my4 = np.arcsin(1/M4)
# Поскольку поток является невязким, давление, плотность и температура торможения в областях 1 и 3 (2 и 4) равны
p03 = p01
p04 = p02
ro03 = ro01
ro04 = ro02
T03 = T01
T04 = T02
# Найдем параметры потока
p3 = 1 / findBrakingParameterPressure(1/p03,M3) # 1/х используется, для поиска давления, но по формуле давления торможения
p4 = 1 / findBrakingParameterPressure(1/p04,M4)
ro3 = 1 / findBrakingParameterDensity(1/ro03,M3)
ro4 = 1 / findBrakingParameterDensity(1/ro04,M4)
T3 = 1 / findBrakingParameterTemperature(1/T01,M3)
T4 = 1 / findBrakingParameterTemperature(1/T02,M4)
a3 = findFlowParameterАcceleration(T3)
a4 = findFlowParameterАcceleration(T4)
V3 = findFlowParameterMach(M3, 1/a3) # используя формулу Маха, найдем скорость
V4 = findFlowParameterMach(M4, 1/a4)

# \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 5 и 6 области \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

t1, t2, u = 0.3, 0.3, 0.17 # Начальные приближенные значение
tettac5,tettac6,delta = fsolve(findDeflectionAngle, (t1, t2, u))
betta5 = betta + delta
betta6 = betta - delta
# Определяем параметры потока в скачках уплотнения
p5 = findFlowParameterPressure(p3,M3,tettac5)
p6 = findFlowParameterPressure(p4,M4,tettac6)
V5 = findFlowParameterVelocity(V3,betta5,tettac5)
V6 = findFlowParameterVelocity(V4,betta6,tettac6)
ro5 = findFlowParameterDensity(ro3,M3,tettac5)
ro6 = findFlowParameterDensity(ro4,M4,tettac6)
T5 = findFlowParameterTemperature(T3,p3,ro3,p5,ro5)
T6 = findFlowParameterTemperature(T4,p4,ro4,p6,ro6)
a5 = findFlowParameterАcceleration(T5)
a6 = findFlowParameterАcceleration(T6)
M5 = findFlowParameterMach(V5, a5)
M6 = findFlowParameterMach(V6, a6)
# Определяем параметры торможения для участка
p05 = findBrakingParameterPressure(p5,M5)
p06 = findBrakingParameterPressure(p6,M6)
ro05 = findBrakingParameterDensity(ro5,M5)
ro06 = findBrakingParameterDensity(ro6,M6)
T05 = findBrakingParameterTemperature(T5,M5)
T06 = findBrakingParameterTemperature(T6,M6)

# \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ АД коэффициенты \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

# Скоростной напор
q = (ro * V ** 2) / 2

# Коэффициенты давления
p_1 = findPressureCoefficient(p1)
p_2 = findPressureCoefficient(p2)
p_3 = findPressureCoefficient(p3)
p_4 = findPressureCoefficient(p4)

# Продольные и нормальные силы
X1, Y1 = LongitudinalAndNormalForces(p1)
X2, Y2 = LongitudinalAndNormalForces(p2)
X3, Y3 = LongitudinalAndNormalForces(p3)
X4, Y4 = LongitudinalAndNormalForces(p4)

# АД коэффициенты
Cx = (X1 + X2 - X3 - X4) / (b * q)
Cy = (-Y1 + Y2 - Y3 + Y4) / (b * q)
mz1 = Y1 * (b / 4) + X1 * (c / 4)
mz2 = - Y2 * (b / 4) - X2 * (c / 4)
mz3 = Y3 * (3 * b / 4) - X3 * (c / 4)
mz4 = - Y4 * (3 * b / 4) + X4 * (c / 4)
mz = (mz1 + mz2 + mz3 + mz4) / (q * b ** 2) # характерная площадь S = 1 * b
Cya = Cy * np.cos(alfa) - Cx * np.sin(alfa)
Cxa = Cy * np.sin(alfa) + Cx * np.cos(alfa)
K = Cya / Cxa
Cd = - (mz / Cy)

# \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Пограничный слой \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

# Удельная теплоемкость
Cp0 = SpecificHeat_Cp(T)
Cp1 = SpecificHeat_Cp(T1)
Cp2 = SpecificHeat_Cp(T2)
Cp3 = SpecificHeat_Cp(T3)
Cp4 = SpecificHeat_Cp(T4)
Cp5 = SpecificHeat_Cp(T5)
Cp6 = SpecificHeat_Cp(T6)

# Динамическая вязкость
mu0 = DynamicViscosity_mu(T)
mu1 = DynamicViscosity_mu(T1)
mu2 = DynamicViscosity_mu(T2)
mu3 = DynamicViscosity_mu(T3)
mu4 = DynamicViscosity_mu(T4)
mu5 = DynamicViscosity_mu(T5)
mu6 = DynamicViscosity_mu(T6)

# Коэффициент теплопроводности
Lambda0 = CoeffThermalConductivity_Lambda(T)
Lambda1 = CoeffThermalConductivity_Lambda(T1)
Lambda2 = CoeffThermalConductivity_Lambda(T2)
Lambda3 = CoeffThermalConductivity_Lambda(T3)
Lambda4 = CoeffThermalConductivity_Lambda(T4)
Lambda5 = CoeffThermalConductivity_Lambda(T5)
Lambda6 = CoeffThermalConductivity_Lambda(T6)

# Для граней 1-4 температура восстановления и определяющая температура
l = 0
t = 1
Trl1, Tl1_z = findDefiningTemperature(l,T1,M1)
Trl2, Tl2_z = findDefiningTemperature(l,T2,M2)
Trl3, Tl3_z = findDefiningTemperature(l,T3,M3)
Trl4, Tl4_z = findDefiningTemperature(l,T4,M4)
Trt1, Tt1_z = findDefiningTemperature(t,T1,M1)
Trt2, Tt2_z = findDefiningTemperature(t,T2,M2)
Trt3, Tt3_z = findDefiningTemperature(t,T3,M3)
Trt4, Tt4_z = findDefiningTemperature(t,T4,M4)

# Критические числа Рейнольдса для 1 и 2 грани
Re_kr1 = CriticalReynoldsNumber(Trl1, M1)
Re_kr2 = CriticalReynoldsNumber(Trl2, M2)

# Координаты точек перехода
x_kr1 = Re_kr1 * mu1 / (V1 * ro1)
x_kr2 = Re_kr2 * mu2 / (V2 * ro2)


#Вывод в эксель

wb = Workbook()
ws = wb.active
ws.append(['','№ области','betta, град','tettac, град','M, [-]','p, 10^4 Па','p0, 10^7  Па','ro, 10^-1 кг/м^3','ro0, 10^1 кг/м^3','T, 10^2 K','T0, 10^3 K', 'V, 10^3 м/с', 'a, 10^2 м/с'])
ws['A1'] = datetime.datetime.now()
ws.append(['', 1,round(betta1*180/np.pi,3),round(tettac1*180/np.pi,3),round(M1,3),round(p1/10000,3),round(p01/10000000,3),round(ro1*10,3),round(ro01/10,3),round(T1/100,3),round(T01/1000,3),round(V1/1000,3),round(a1/100,3)])
ws['A2'] = datetime.datetime.now()
ws.append(['', 2,round(betta2*180/np.pi,3),round(tettac2*180/np.pi,3),round(M2,3),round(p2/10000,3),round(p02/10000000,3),round(ro2*10,3),round(ro02/10,3),round(T2/100,3),round(T02/1000,3),round(V2/1000,3),round(a2/100,3)])
ws['A3'] = datetime.datetime.now()
ws.append(['', 3,round(betta3*180/np.pi,3),"-",round(M3,3),round(p3/10000,3),round(p03/10000000,3),round(ro3*10,3),round(ro03/10,3),round(T3/100,3),round(T03/1000,3),round(V3/1000,3),round(a3/100,3)])
ws['A4'] = datetime.datetime.now()
ws.append(['', 4,round(betta4*180/np.pi,3),"-",round(M4,3),round(p4/10000,3),round(p04/10000000,3),round(ro4*10,3),round(ro04/10,3),round(T4/100,3),round(T04/1000,3),round(V4/1000,3),round(a4/100,3)])
ws['A5'] = datetime.datetime.now()
ws.append(['', 5,round(betta5*180/np.pi,3),round(tettac5*180/np.pi,3),round(M5,3),round(p5/10000,3),round(p05/10000000,3),round(ro5*10,3),round(ro05/10,3),round(T5/100,3),round(T05/1000,3),round(V5/1000,3),round(a5/100,3)])
ws['A6'] = datetime.datetime.now()
ws.append(['', 6,round(betta6*180/np.pi,3),round(tettac6*180/np.pi,3),round(M6,3),round(p6/10000,3),round(p06/10000000,3),round(ro6*10,3),round(ro06/10,3),round(T6/100,3),round(T06/1000,3),round(V6/1000,3),round(a6/100,3)])
ws['A8'] = datetime.datetime.now()
ws.append(['', '№ области','Cp, 10^2 Дж*К/кг','mu, 10^-5 Па*с','lambda, 10^-2 Вт*К/м', '','№ области','Tr(л), 10^3 К', 'T_z(л), 10^2 К', 'Tr(т), 10^3 К', 'T_z(т), 10^2 К'])
ws['A9'] = datetime.datetime.now()
ws.append(['',1, round(Cp1/100,3), round(mu1*100000,3), round(Lambda1*100,3),'', 1,round(Trl1/1000,3), round(Tl1_z/100,3), round(Trt1/1000,3), round(Tt1_z/100,3)])
ws['A10'] = datetime.datetime.now()
ws.append(['',2, round(Cp2/100,3), round(mu2*100000,3), round(Lambda2*100,3),'', 2,round(Trl2/1000,3), round(Tl2_z/100,3), round(Trt2/1000,3), round(Tt2_z/100,3)])
ws['A11'] = datetime.datetime.now()
ws.append(['',3, round(Cp3/100,3), round(mu3*100000,3), round(Lambda3*100,3),'', 3,round(Trl3/1000,3), round(Tl3_z/100,3), round(Trt3/1000,3), round(Tt3_z/100,3)])
ws['A12'] = datetime.datetime.now()
ws.append(['',4, round(Cp4/100,3), round(mu4*100000,3), round(Lambda4*100,3),'', 4,round(Trl4/1000,3), round(Tl4_z/100,3), round(Trt4/1000,3), round(Tt4_z/100,3)])
ws['A14'] = datetime.datetime.now()
ws.append(['', 'Cy, 10^-2', 'Cx, 10^-3', 'mz, 10^-2', 'Cya, 10^2', 'Cxa, 10^3', 'Cд', 'K'])
ws['A15'] = datetime.datetime.now()
ws.append(['', round(Cy*100,3), round(Cx*1000,3), round(mz*100,3), round(Cya*100,3), round(Cxa*1000,3), round(Cd,3), round(K,3)])
ws['A17'] = datetime.datetime.now()
ws.append(['','№ области', 'Reкрит, 10^6', 'Xкрит, 10^-1'])
ws['A18'] = datetime.datetime.now()
ws.append(['', 1, round(Re_kr1/1000000,3), round(x_kr1*10,3)])
ws['A19'] = datetime.datetime.now()
ws.append(['', 2, round(Re_kr2/1000000,3), round(x_kr2*10,3)])
wb.save("sample.xlsx")


atm.vivod()
print("Угол бетта: " + str(round(betta,3))+ ' рад; '+ str(round((betta*180)/np.pi,2))+ ' град')
print("Угол альфа: " + str(round(alfa,3))+ ' рад; '+ str(round((alfa*180)/np.pi,2100))+ ' град')
print("Угол бетта1: " + str(round(betta1,3))+ ' рад; '+ str(round((betta1*180)/np.pi,2))+ ' град')
print("Угол бетта2: " + str(round(betta2,3))+ ' рад; '+ str(round((betta2*180)/np.pi,2))+ ' град')
print("Угол теттац1: " + str(round(tettac1,3))+ ' рад; '+ str(round((tettac1*180)/np.pi,2))+ ' град')
print("Угол теттац2: " + str(round(tettac2,3))+ ' рад; '+ str(round((tettac2*180)/np.pi,2))+ ' град')
print("Давление 1 области: " + str(round(p1/10000,3)) + ' * 10^4 Па')
print("Давление 2 области: " + str(round(p2/10000,3)) + ' * 10^4'+' Па')
print("Плотность 1 области: " + str(round(ro1*10,3))+' * 10^-1 кг/м^3')
print("Плотность 2 области: " + str(round(ro2*10,3))+' * 10^-1 кг/м^3')
print("Температура 1 области: " + str(round(T1/100,3)) + ' * 10^2 K')
print("Температура 2 области: " + str(round(T2/100,3)) + ' * 10^2 K')
print("Скорость 1 области: " + str(round(V1/1000,3)) + ' * 10^3 м/с')
print("Скорость 2 области: " + str(round(V2/1000,3)) + ' * 10^3 м/с')
print("Скорость звука 1 области: " + str(round(a1/100,3))+' * 10^2 м/с')
print("Скорость звука 2 области: " + str(round(a2/100,3))+' * 10^2 м/с')
print("Число Маха 1 области: " + str(round(M1,3)))
print("Число Маха 2 области: " + str(round(M2,3)))
print("Давление торможения 1 области: " + str(round(p01/10000000,3)) + ' * 10^7  Па')
print("Давление торможения 2 области: " + str(round(p02/10000000,3)) + ' * 10^7'+' Па')
print("Плотность торможения 1 области: " + str(round(ro01/10,3))+' * 10^1'+' кг/м^3')
print("Плотность торможения 2 области: " + str(round(ro02/10,3))+' * 10^1'+' кг/м^3')
print("Температура торможения 1 области: " + str(round(T01/1000,3)) + ' * 10^3 K')
print("Температура торможения 2 области: " + str(round(T02/1000,3)) + ' * 10^3 K')
print("Угол бетта3: " + str(round(betta3,3))+ ' рад; '+ str(round((betta3*180)/np.pi,2))+ ' град')
print("Угол бетта4: " + str(round(betta4,3))+ ' рад; '+ str(round((betta4*180)/np.pi,2))+ ' град')
print("Фиктивный угол 1: " + str(round(w1,3))+ ' рад; '+ str(round((w1*180)/np.pi,2))+ ' град')
print("Фиктивный угол 2: " + str(round(w2,3))+ ' рад; '+ str(round((w2*180)/np.pi,2))+ ' град')
print("Фиктивный угол 3: " + str(round(w3,3))+ ' рад; '+ str(round((w3*180)/np.pi,2))+ ' град')
print("Фиктивный угол 4: " + str(round(w4,3))+ ' рад; '+ str(round((w4*180)/np.pi,2))+ ' град')
print("Число Маха 3 области: " + str(round(M3,3)))
print("Число Маха 4 области: " + str(round(M4,3)))
print("Давление торможения 3 области: " + str(round(p03/10000000,3)) + ' * 10^7'+' Па')
print("Давление торможения 4 области: " + str(round(p04/10000000,3)) + ' * 10^7'+' Па')
print("Плотность торможения 3 области: " + str(round(ro03/10,3))+' * 10^1'+' кг/м^3')
print("Плотность торможения 4 области: " + str(round(ro04/10,3))+' * 10^1'+' кг/м^3')
print("Температура торможения 3 области: " + str(round(T03/1000,3)) + ' * 10^3 K')
print("Температура торможения 4 области: " + str(round(T04/1000,3)) + ' * 10^3 K')
print("Давление 3 области: " + str(round(p3/1000,3)) + ' * 10^3'+' Па')
print("Давление 4 области: " + str(round(p4/1000,3)) + ' * 10^3'+' Па')
print("Плотность 3 области: " + str(round(ro3*100,3))+' * 10^-2'+' кг/м^3')
print("Плотность 4 области: " + str(round(ro4*10,3))+' * 10^-1'+' кг/м^3')
print("Температура 3 области: " + str(round(T3/100,3)) + ' * 10^2 K')
print("Температура 4 области: " + str(round(T4/100,3)) + ' * 10^2 K')
print("Скорость звука 3 области: " + str(round(a3/100,3))+' * 10^2 м/с')
print("Скорость звука 4 области: " + str(round(a4/100,3))+' * 10^2 м/с')
print("Скорость 3 области: " + str(round(V3/1000,3)) + ' * 10^3 м/с')
print("Скорость 4 области: " + str(round(V4/1000,3)) + ' * 10^3 м/с')
print("Давление 5 области: " + str(round(p5 / 10000, 3)) + ' * 10^4' + ' Па')
print("Давление 6 области: " + str(round(p6 / 10000, 3)) + ' * 10^4' + ' Па')
print("Угол теттац5: " + str(round(tettac5,3))+ ' рад; '+ str(round((tettac5*180)/np.pi,2))+ ' град')
print("Угол теттац6: " + str(round(tettac6,3))+ ' рад; '+ str(round((tettac6*180)/np.pi,2))+ ' град')
print("Угол отклонения: " + str(round(delta,3))+ ' рад; '+ str(round((delta*180)/np.pi,2))+ ' град')
print("Угол бетта5: " + str(round(betta5,3))+ ' рад; '+ str(round((betta5*180)/np.pi,2))+ ' град')
print("Угол бетта6: " + str(round(betta6,3))+ ' рад; '+ str(round((betta6*180)/np.pi,2))+ ' град')
print("Плотность 5 области: " + str(round(ro5*10,3))+' * 10^-1'+' кг/м^3')
print("Плотность 6 области: " + str(round(ro6*10,3))+' * 10^-1'+' кг/м^3')
print("Температура 5 области: " + str(round(T5/100,3)) + ' * 10^2 K')
print("Температура 6 области: " + str(round(T6/100,3)) + ' * 10^2 K')
print("Скорость 5 области: " + str(round(V5/1000,3)) + ' * 10^3 м/с')
print("Скорость 6 области: " + str(round(V6/1000,3)) + ' * 10^3 м/с')
print("Скорость звука 5 области: " + str(round(a5/100,3))+' * 10^2 м/с')
print("Скорость звука 6 области: " + str(round(a6/100,3))+' * 10^2 м/с')
print("Число Маха 5 области: " + str(round(M5,3)))
print("Число Маха 6 области: " + str(round(M6,3)))
print("Давление торможения 5 области: " + str(round(p05/10000000,3)) + ' * 10^7'+' Па')
print("Давление торможения 6 области: " + str(round(p06/10000000,3)) + ' * 10^7'+' Па')
print("Плотность торможения 5 области: " + str(round(ro05/10,3))+' * 10^1'+' кг/м^3')
print("Плотность торможения 6 области: " + str(round(ro06/10,3))+' * 10^1'+' кг/м^3')
print("Температура торможения 5 области: " + str(round(T05/1000,3)) + ' * 10^3 K')
print("Температура торможения 6 области: " + str(round(T06/1000,3)) + ' * 10^3 K')
print("Угол мю1: " + str(round(my1,3))+ ' рад; '+ str(round((my1*180)/np.pi,2))+ ' град')
print("Угол мю2: " + str(round(my2,3))+ ' рад; '+ str(round((my2*180)/np.pi,2))+ ' град')
print("Угол мю3: " + str(round(my3,3))+ ' рад; '+ str(round((my3*180)/np.pi,2))+ ' град')
print("Угол мю4: " + str(round(my4,3))+ ' рад; '+ str(round((my4*180)/np.pi,2))+ ' град')

print([betta1*180/np.pi,tettac1*180/np.pi,M1,p1,p01,ro1,ro01,T1,T01,V1,a1])
print([betta2*180/np.pi,tettac2*180/np.pi,M2,p2,p02,ro2,ro02,T2,T02,V2,a2])
print([betta3*180/np.pi,"-",M3,p3,p03,ro3,ro03,T3,T03,V3,a3])
print([betta4*180/np.pi,"-",M4,p4,p04,ro4,ro04,T4,T04,V4,a4])
print([betta5*180/np.pi,tettac5*180/np.pi,M5,p5,p05,ro5,ro05,T5,T05,V5,a5])
print([betta6*180/np.pi,tettac6*180/np.pi,M6,p6,p06,ro6,ro06,T6,T06,V6,a6])

print(Cp1,Cp2,Cp3,Cp4)
print(mu1,mu2,mu3,mu4)
print(Lambda1,Lambda2,Lambda3,Lambda4)

print(Cy)
print(Cx)
print(mz)
print(Cya)
print(Cxa)
print(K)
print(Cd)

print([Trl1, Tl1_z, Trt1, Tt1_z])
print([Trl2, Tl2_z, Trt2, Tt2_z])
print([Trl3, Tl3_z, Trt3, Tt3_z])
print([Trl4, Tl4_z, Trt4, Tt4_z])