from math import sqrt, exp, log10, pi, atan, tan, asin, sin, cos, degrees as deg, radians as rad
import numpy as np

#  Используемые константы
Cpc = 1004.7  # удельная стандартная теплоемкость при постоянном давлении при Tc1
muc = 1.79 / 100000  # стандартная вязкость при Tc1
lambda_c = 23.2 / 1000  # стандартная теплопроводность при Tc2
k = 1.4  # показатель адиабаты для воздуха
Tc1 = 288.15  # стандартная температура для вязкости
Tc2 = 261  # стандартная температура для теплопроводности
fi = 0.1  # степень для определения теплоемкости
n = 0.76  # степень для определения вязкости
eps = 0.85  # степень для определения теплопроводности
r = 6356767  # условный радиус Земли
Reкрст = 5 * 10 ** 6  # критическое число Рейнольдса (при Tст/Tr = 1)

# Исходные данные
N = 11
nг = 1

if nг != 1 and nг != 2:
    while (nг != 1) and (n != 2):
        print('ошибка, введите другой номер группы')
        nг = int(input())
if nг == 1:
    H = 10 + 0.5 * N  # геометрическая высота полета в км
    Minf = 4 + 0.15 * N  # скорость полета
    Tст = 405  # температура стенки в К
if nг == 2:
    H = 10 + 0.5 * N  # геометрическая высота полета в км
    Minf = 4.5 + 0.1 * N  # скорость полета
    Tст = 373  # температура стенки в К
H = H * 1000  # геометрическая высота в м
while H > 71000:
    print('ошибка, введите другой номер варианта')
    N = int(input())
    if nг == 1:
        H = 10 + 0.5 * N  # геометрическая высота полета в км
        Minf = 4.5 + 0.15 * N  # скорость полета
        Tст = 405  # температура стенки в К
    if nг == 2:
        H = 10 + 0.4 * N  # геометрическая высота полета в км
        Minf = 5 + 0.1 * N  # скорость полета
        Tст = 373  # температура стенки в К
    H = H * 1000  # геометрическая высота в м
H = r * H / (r + H)  # переход на геопотенциальную высоту

if 1 <= N <= 6:
    c = 80  # мм
    b = 1000  # мм
if 7 <= N <= 12:
    c = 90  # мм
    b = 1100  # мм
if 13 <= N <= 18:
    c = 100  # мм
    b = 1200  # мм
if 19 <= N <= 24:
    c = 110  # мм
    b = 1300  # мм
if N > 24:
    c = 110  # мм
    b = 1300  # мм
c = c / 1000  # высота профиля в м
b = b / 1000  # длина профиля в м
L = sqrt((c / 2) ** 2 + (b / 2) ** 2)  # длина грани в м
betta = atan(c / b)  # полуугол клина
alpha = rad(2 + 0.1 * N)  # угол атаки
if alpha > betta:
    print('проверка на образование СУ не пройдена')


# Используемые функции
def Atm(H):  # ГОСТ 4401-81 для расчета параметров атмосферы до высоты 71 км
    gс = 9.80665  # ускорение свободного падения
    k = 1.4  # показатель адиабаты
    R = 287.05287  # удельная газовая постоянная, R=const для высот до 94 км
    S = 110.4  # эмпирические коэффициент Сатерлэнда
    beta_s = 1.458 * 10 ** -6  # в уравнении для определения динамической вязкости

    # значения температур и вертикальных градиентов, принятые для расчета стандартной атмосферы до высоты 120 км
    H_st = (-2000, 0, 11000, 20000, 32000, 47000, 51000, 71000)  # границы слоев атмосферы
    T_st = (301.15, 288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65)  # температуры на границах
    beta_ = (-0.0065, -0.0065, 0, 0.001, 0.0028, 0, -0.0028, -0.002)  # градиаент температуры по слоям
    P_st = (127774, 101325, 22632, 5474.87, 868.014, 110.906, 66.9384, 3.95639)  # давление на границах
    for i in range(7):
        if H_st[i] <= H < H_st[i + 1]:
            Hst = H_st[i]
            Tst = T_st[i]
            beta = beta_[i]
            Pst = P_st[i]
        elif H == 71000:
            Hst = H_st[7]
            Tst = T_st[7]
            beta = beta[7]
            Pst = P_st[7]
        else:
            i += 1
    T = Tst + beta * (H - Hst)
    if beta != 0:
        P = Pst * (1 + (beta * (H - Hst)) / Tst) ** ((-gс) / (R * beta))
    else:
        P = Pst * 10 ** ((-gс * 0.434294 * (H - Hst)) / (R * T))
    ro = P / (R * T)  # до высот 94 км
    a = sqrt(k * R * T)  # для R=const
    mu = (beta_s * T ** 1.5) / (T + S)  # динамическая вязкость
    nu = mu / ro  # кинематическая вязкость
    return P, ro, T, a, mu, nu


def thetaC(M, beta):  # Нахождение угла наклона СУ
    a = 0.0000000000001
    p = 0.001
    l = 0
    r = 1
    step = 0.000001
    theta_C = beta + step
    while abs(l - r) > a:
        while abs(l - r) > p:
            l = tan(theta_C) / tan(theta_C - beta)
            r = 0.5 * (k + 1.) * (M ** 2) * (sin(theta_C) ** 2.) \
                / (1. + 0.5 * (k - 1.) * (M ** 2.) * (sin(theta_C) ** 2.))
            theta_C += step
        theta_C -= step
        p /= 10
        step /= 10
    if theta_C >= rad(90):
        print('-------------------------')
        print('ошибка нахождения угла СУ')
    return theta_C


def parametr0(Pi, roi, Ti, Mi):  # Нахождение параметров торможения
    P0i = Pi * (1 + 0.5 * (k - 1) * Mi ** 2) ** (k / (k - 1))
    ro0i = roi * (1 + 0.5 * (k - 1) * Mi ** 2) ** (1 / (k - 1))
    T0i = Ti * (1 + 0.5 * (k - 1) * Mi ** 2)
    return P0i, ro0i, T0i


def CY(beta, Minf, H, tettac):  # Определение параметров для 1 и 2 граней
    Pinf, roinf, Tinf, ainf, muinf, nuinf = Atm(H)
    Vinf = Minf * ainf
    P = Pinf * ((2 * k * Minf ** 2 * sin(tettac) ** 2) / (k + 1) - (k - 1) / (k + 1))
    ro = roinf * (0.5 * (k + 1) * Minf ** 2 * sin(tettac) ** 2) / (1 + 0.5 * (k - 1) * Minf ** 2 * sin(tettac) ** 2)
    T = Tinf * (P * roinf) / (Pinf * ro)
    V = Vinf * (cos(tettac)) / (cos(tettac - beta))
    a = sqrt((k * P) / ro)
    M = V / a
    muu = asin(1 / M)  # Угол Маха
    P0, ro0, T0 = parametr0(P, ro, T, M)
    Cp = Cpc * (T / Tc1) ** fi  # Удельная теплоемкость при постоянном давлении
    lambda_ = lambda_c * (T / Tc2) ** eps  # Теплопроводность
    mu = muc * (T / Tc1) ** n  # Динамическая вязкость
    return M, P, T, ro, Cp, mu, lambda_, muu, V, P0, ro0, T0


def PM(Mi, P0i, ro0i, T0i):
    dM = 0.0001
    omegak = sqrt((k + 1) / (k - 1)) * atan(sqrt(((k - 1) * (Mi ** 2 - 1)) / (k + 1))) - atan(sqrt(Mi ** 2 - 1))
    omega = omegak + 2 * betta
    M = 1
    p = 0.00001
    l = 0
    r = 1
    while abs(l - r) > p:
        l = omega
        r = (sqrt((k + 1) / (k - 1)) * atan(sqrt(((k - 1) * (M ** 2 - 1)) / (k + 1))) - atan(
            sqrt(M ** 2 - 1)))
        M += dM
        if M > 20:
            print('ОШИБКА: не найдено решение уравнения задачи Прандтля-Майера')
            break
    muu = asin(1 / M)
    P = P0i * (1 + 0.5 * (k - 1) * M ** 2) ** (-k / (k - 1))
    ro = ro0i * (1 + 0.5 * (k - 1) * M ** 2) ** (-1 / (k - 1))
    T = T0i * (1 + 0.5 * (k - 1) * M ** 2) ** -1
    a = sqrt((k * P) / ro)
    V = M * a
    Cp = Cpc * (T / Tc1) ** fi  # Удельная теплоемкость при постоянном давлении
    lambda_ = lambda_c * (T / Tc2) ** eps  # Теплопроводность
    mu = muc * (T / Tc1) ** n  # Динамическая вязкость
    return M, P, T, ro, Cp, mu, lambda_, muu, V, P0i, ro0i, T0i, betta * 2


def opr_temp(T, M, r, q):
    Tr = T * (1 + r * 0.5 * (k - 1) * M ** 2)  # Температура восстановления
    Topr = 0.5 * (Tст + T) + 0.22 * (Tr - T)  # Определяющая температура
    Topr1 = 10
    while (Topr1 - Topr) >= 1:  # Пока определяющая температура изменяется больше чем на 1
        Topr1 = Topr
        Cp = Cpc * (Topr / Tc1) ** fi  # Удельная теплоемкость при постоянном давлении
        mu = muc * (Topr / Tc1) ** n  # Динамическая вязкость
        lambda_ = lambda_c * (Topr / Tc2) ** eps  # Теплопроводность
        Pr = (Cp * mu) / lambda_  # Число Прандтля
        if q == 1:
            r = sqrt(Pr)  # Коэффициент восстановления для ЛПС
        else:
            r = Pr ** (1 / 3)  # Коэффициент восстановления для ТПС
        Tr = T * (1 + r * 0.5 * (k - 1) * M ** 2)
        Topr = 0.5 * (Tст + T) + 0.22 * (Tr - T)
    return Tr, Topr


def X(Pi, Pinf):
    X = (Pi - Pinf) * (c / 2)
    return X


def Y(Pi, Pinf):
    Y = (Pi - Pinf) * (b / 2)
    return Y


def Re(Tст, Tr, M):
    Tстср = Tст / Tr
    x = (Tстср - 1) / M ** 2
    y = 1 - 16 * x - 412.5 * (x ** 2) - 35000 * (x ** 3) - 375000 * (x ** 4)
    Reкр = Reкрст * y
    return Reкр

def lam_nc(x, ro, V, mu):
    Re = (V * ro * x) / mu
    delt = (4.64 * x) / sqrt(Re)
    tau = (0.323 * ro * V ** 2) / (sqrt(Re))
    Cfx = 0.646 / sqrt(Re)
    Cf = 1.292 / sqrt(Re)
    delt_st = 1.74 * x / sqrt(Re)
    delt_st_st = 0.646 * x / sqrt(Re)
    return delt, tau, Cfx, Cf, delt_st, delt_st_st, Re


def lam(delt_nc, tau_nc, Cfx_nc, Cf_nc, delt_st_nc, delt_st_st_nc, T, Toprl):
    delt = delt_nc * (Toprl / T) ** (0.5 * (n + 1))
    tau = tau_nc * (Toprl / T) ** (0.5 * (n - 1))
    Cfx = Cfx_nc * (Toprl / T) ** (0.5 * (n - 1))
    Cf = Cf_nc * (Toprl / T) ** (0.5 * (n - 1))
    delt_st = delt_st_nc * delt / delt_nc
    delt_st_st = delt_st_st_nc * Cfx / Cfx_nc
    return delt, tau, Cfx, Cf, delt_st, delt_st_st


def tur_nc(x, ro, V, mu, x1):
    Re = (V * ro * x) / mu
    delt = (0.37 * x) / (Re ** 0.2)
    tau = (0.02899 * ro * V ** 2) / (Re ** 0.2)
    Cfx = 0.0578 / Re ** 0.2
    Cf = 0.074 / Re ** 0.2
    if x1 == 0:
        Cf = 0
    kk = 7
    delt_st = delt / (kk + 1)
    delt_st_st = (delt * kk) / ((kk + 1) * (kk + 2))

    return delt, tau, Cfx, Cf, delt_st, delt_st_st, Re


def tur(delt_nc, tau_nc, Cfx_nc, Cf_nc, delt_st_nc, delt_st_st_nc, T, Toprt):
    delt = delt_nc * (Toprt / T) ** (0.2 * (n + 1))
    tau = tau_nc * (Toprt / T) ** (0.2 * (n - 4))
    Cfx = Cfx_nc * (Toprt / T) ** (0.2 * (n - 4))
    Cf = Cf_nc * (Toprt / T) ** (0.2 * (n - 4))
    delt_st = delt_st_nc * delt / delt_nc
    delt_st_st = delt_st_st_nc * Cfx / Cfx_nc

    return delt, tau, Cfx, Cf, delt_st, delt_st_st


def Cfср12(xкр, deltaX, ro, V, mu, Reкр, Toprl, Toprt, T):
    l = L - xкр + deltaX
    Rel = ro * V * l / mu
    RedeltaX = ro * V * deltaX / mu
    Cfсрl = 1.292 * ((Toprl / T) ** (0.5 * (n - 1))) / Reкр ** 0.5
    Cfср1t = 0.074 * ((Toprt / T) ** (0.2 * (n - 4))) / Rel ** 0.2
    Cfср2t = 0.074 * ((Toprt / T) ** (0.2 * (n - 4))) / RedeltaX ** 0.2
    Cfср = Cfсрl * xкр / L + Cfср1t * l / L - Cfср2t * deltaX / L
    return Cfср


def Cfср34(deltaX, ro, V, mu, Toprt, T):
    l = L + deltaX
    Rel = ro * V * l / mu
    RedeltaX = ro3 * V3 * deltaX / mu
    Cfср1t = 0.074 * ((Toprt / T) ** (0.2 * (n - 4))) / Rel ** 0.2
    Cfср2t = 0.074 * ((Toprt / T) ** (0.2 * (n - 4))) / RedeltaX ** 0.2
    Cfср = Cfср1t * l / L - Cfср2t * deltaX / L
    return Cfср


# Основная программа

# Первая часть

# Параметры набегающего потока
Pinf, roinf, Tinf, ainf, muinf, nuinf = Atm(H)
Vinf = Minf * ainf
lambda_inf = (2.648151 * 10 ** -3 * Tinf ** 1.5) / (Tinf + 245.4 * 10 ** (-12 / Tinf))  # теплопроводность
print('\nИсходные данные')
print('-----------------')
print('номер варианта=', N, '\nномер группы=', nг, '\nTстенки=', Tст, '\nMinf=', Minf,
      '\nальфа=', '{0:.1f}'.format(deg(alpha)), '\nH(геопотенциальная)=', '{0:.5f}'.format(H),
      '\nh(геметрическая)=', (H * r) / (r - H), '\nc=', c, '\nb=', b, '\nL=', '{0:.5f}'.format(L),
      '\nbetta=', '{0:.3f}'.format(deg(betta)), '\nTinf=', Tinf, '\nPinf=', '{0:.5f}'.format(Pinf),
      '\nroinf=', '{0:.7f}'.format(roinf), '\nainf=', '{0:.5f}'.format(ainf), '\nVinf=', '{0:.5f}'.format(Vinf),
      '\nlambda_inf=', '{0:.6f}'.format(lambda_inf), '\nmuinf=', '{0:.9f}'.format(muinf))

# Расчет параметров обтекания сверхзвуковым невязким потоком

# Параметры торможения невозмущенного набегающего потока
P0, ro0, T0 = parametr0(Pinf, roinf, Tinf, Minf)
print('\nПараметры торможения набегающего потока:', '\n--------------------------------------', '\nP0=', P0, '\nro0=',
      ro0, '\nT0=', T0)

# Перед гранями 1 и 2 образуется косой СУ

# Грань №1
beta1 = betta - alpha  # Угол отклонения потока
tettac1 = thetaC(Minf, beta1)
M1, P1, T1, ro1, Cp1, mu1, lambda_1, muu1, V1, P01, ro01, T01 = CY(beta1, Minf, H, tettac1)
r1 = 0.83  # начальное приближение коэффициента восстановления для ЛПС
r2 = 0.88  # начальное приближение коэффициента восстановления для ТПС
Trl1, Toprl1 = opr_temp(T1, M1, r1, 1)
Trt1, Toprt1 = opr_temp(T1, M1, r2, 2)
print('\nГрань №1:', '\n-----------------', '\nM=', M1, '\nV=', V1, '\nP=', P1, '\nP0=', P01, '\nT=', T1, '\nro=', ro1,
      '\nCp=', Cp1, '\nmu=', mu1, '\nlambda=', lambda_1, '\nTr(л)=', Trl1, '\nT*(л)=', Toprl1, '\nTr(т)=', Trt1,
      '\nT*(т)=', Toprt1, '\nteta=', deg(tettac1), '\nmuu=', deg(muu1), '\nbeta=', deg(beta1))

# Грань №2
beta2 = betta + alpha  # Угол отклонения потока
tettac2 = thetaC(Minf, beta2)
M2, P2, T2, ro2, Cp2, mu2, lambda_2, muu2, V2, P02, ro02, T02 = CY(beta2, Minf, H, tettac2)
r1 = 0.83  # начальное приближение коэффициента восстановления для ЛПС
r2 = 0.88  # начальное приближение коэффициента восстановления для ТПС
Trl2, Toprl2 = opr_temp(T2, M2, r1, 1)
Trt2, Toprt2 = opr_temp(T2, M2, r2, 2)
print('\nГрань №2:', '\n-----------------', '\nM=', M2, '\nV=', V2, '\nP=', P2, '\nP0=', P02, '\nT=', T2, '\nro=', ro2,
      '\nCp=', Cp2, '\nmu=', mu2, '\nlambda=', lambda_2, '\nTr(л)=', Trl2, '\nT*(л)=', Toprl2, '\nTr(т)=', Trt2,
      '\nT*(т)=', Toprt2, '\nteta=', deg(tettac2), '\nmuu=', deg(muu2), '\nbeta=', deg(beta2))

# Грань №3
# невязкий поток, параметры торможенения не меняются
M3, P3, T3, ro3, Cp3, mu3, lambda_3, muu3, V3, P03, ro03, T03, beta3 = PM(M1, P01, ro01, T01)
r1 = 0.83  # начальное приближение коэффициента восстановления для ЛПС
r2 = 0.88  # начальное приближение коэффициента восстановления для ТПС
Trl3, Toprl3 = opr_temp(T3, M3, r1, 1)
Trt3, Toprt3 = opr_temp(T3, M3, r2, 2)
print('\nГрань №3:', '\n-----------------', '\nM=', M3, '\nV=', V3, '\nP=', P3, '\nP0=', P03, '\nT=', T3, '\nro=', ro3,
      '\nCp=', Cp3, '\nmu=', mu3, '\nlambda=', lambda_3, '\nTr(л)=', Trl3, '\nT*(л)=', Toprl3, '\nTr(т)=', Trt3,
      '\nT*(т)=', Toprt3, '\nmuu=', deg(muu3), '\nbeta=', deg(beta3))

# Грань №4
# невязкий поток, параметры торможенения не меняются
M4, P4, T4, ro4, Cp4, mu4, lambda_4, muu4, V4, P04, ro04, T04, beta4 = PM(M2, P02, ro02, T02)
r1 = 0.83  # начальное приближение коэффициента восстановления для ЛПС
r2 = 0.88  # начальное приближение коэффициента восстановления для ТПС
Trl4, Toprl4 = opr_temp(T4, M4, r1, 1)
Trt4, Toprt4 = opr_temp(T4, M4, r2, 2)
print('\nГрань №4:', '\n-----------------', '\nM=', M4, '\nV=', V4, '\nP=', P4, '\nP0=', P04, '\nT=', T4, '\nro=', ro4,
      '\nCp=', Cp4, '\nmu=', mu4, '\nlambda=', lambda_4, '\nTr(л)=', Trl4, '\nT*(л)=', Toprl4, '\nTr(т)=', Trt4,
      '\nT*(т)=', Toprt4, '\nmuu=', deg(muu4), '\nbeta=', deg(beta4))

# Области №5 и №6
# Нахождение из условия равенства давлений в этх облостях линии тангенсального разрыва
delta = alpha - rad(1.)
ddelta = 0.01
P5 = 1
P6 = 0
m = 10
while abs(P6 - P5) > 0.01:
    # Углы отклонения потока
    delta += ddelta
    beta5 = betta + delta  # + т.к. в общей формуле  в функии - (т.е. вычтется и betta, и detta)
    beta6 = betta - delta  # - т.к. в общей формуле  в функии - (т.е. вычтется betta, detta добавится)
    tettac5 = thetaC(M3, beta5)
    tettac6 = thetaC(M4, beta6)
    P5 = P3 * ((2 * k * M3 ** 2 * sin(tettac5) ** 2) / (k + 1) - (k - 1) / (k + 1))
    P6 = P4 * ((2 * k * M4 ** 2 * sin(tettac6) ** 2) / (k + 1) - (k - 1) / (k + 1))
    if abs(P6 - P5) / P5 * 100 < m:
        print('Прошли', m, '% разницу давлений')
        delta -= ddelta
        print('P5=', P5, '\nP6=', P6)
        ddelta /= 10
        m /= 10

# Область №5
# По фоормулам для косого СУ
beta5 = beta3 / 2 - delta
ro5 = ro3 * (0.5 * (k + 1) * M3 ** 2 * sin(tettac5) ** 2) / (1 + 0.5 * (k - 1) * M3 ** 2 * sin(tettac5) ** 2)
T5 = T3 * (P5 * ro3) / (P3 * ro5)
V5 = V3 * (cos(tettac5)) / (cos(tettac5 - beta5))
a5 = sqrt((k * P5) / ro5)
M5 = V5 / a5
muu5 = asin(1 / M5)  # Угол Маха
P05, ro05, T05 = parametr0(P5, ro5, T5, M5)
Cp5 = Cpc * (T5 / Tc1) ** fi  # Удельная теплоемкость при постоянном давлении
lambda_5 = lambda_c * (T5 / Tc2) ** eps  # Теплопроводность
mu5 = muc * (T5 / Tc1) ** n  # Динамическая вязкость
print('\nОбласть №5:', '\n-----------------', '\nM=', M5, '\nV=', V5, '\nP=', P5, '\nP0=', P05, '\nT=', T5, '\nro=', ro5,
      '\nCp=', Cp5, '\nmu=', mu5, '\nlambda=', lambda_5, '\ntettac=', deg(tettac5), '\nmuu=', deg(muu5), '\nbeta=',
      deg(beta5))

# Область №6
# По фоормулам для косого СУ
beta6 = beta3 / 2 + delta
ro6 = ro4 * (0.5 * (k + 1) * M4 ** 2 * sin(tettac6) ** 2) / (1 + 0.5 * (k - 1) * M4 ** 2 * sin(tettac6) ** 2)
T6 = T4 * (P6 * ro4) / (P4 * ro6)
V6 = V4 * (cos(tettac6)) / (cos(tettac6 - beta6))
a6 = sqrt((k * P6) / ro6)
M6 = V6 / a6
muu6 = asin(1 / M6)  # Угол Маха
P06, ro06, T06 = parametr0(P6, ro6, T6, M6)
Cp6 = Cpc * (T6 / Tc1) ** fi  # Удельная теплоемкость при постоянном давлении
lambda_6 = lambda_c * (T6 / Tc2) ** eps  # Теплопроводность
mu6 = muc * (T6 / Tc1) ** n  # Динамическая вязкость
print('\nОбласть №6:', '\n-----------------', '\nM=', M6, '\nV=', V6, '\nP=', P6, '\nP0=', P06, '\nT=', T6, '\nro=', ro6,
      '\nCp=', Cp6, '\nmu=', mu6, '\nlambda=', lambda_6, '\ntettac=', deg(tettac6), '\nmuu=', deg(muu6), '\nbeta=',
      deg(beta6))

# Расчет а.д. характеристик (идеальная среда - силовое воздействие только от действия давления на грани
qinf = (roinf * Vinf ** 2) / 2
print('\nАэродинамические характеристики')
print('--------------------------------')
print('qinf', qinf)
X1 = X(P1, Pinf)
print('X1', X1)
X2 = X(P2, Pinf)
X3 = X(P3, Pinf)
X4 = X(P4, Pinf)
Cx = (X1 + X2 - X3 - X4) / (b * qinf)
print('Cx=', Cx)
Y1 = Y(P1, Pinf)
Y2 = Y(P2, Pinf)
Y3 = Y(P3, Pinf)
Y4 = Y(P4, Pinf)
Cy = (-Y1 + Y2 - Y3 + Y4) / (b * qinf)
print('Cy=', Cy)
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) / (qinf * b ** 2)
print('mz=', mz)
Cxa = Cx * cos(alpha) + Cy * sin(alpha)
print('Cxa =', Cxa)
Cya = - Cx * sin(alpha) + Cy * cos(alpha)
print('Cya=', Cya)

# Вторая часть

# Определение типа течения и расчет xкр на 1 и 2 гранях
Reкр1 = Re(Tст, Trl1, M1)  # Возможно неверная Tr Критическое число Рейнольдса
Reкр2 = Re(Tст, Trl2, M2)  # Возможно неверная Tr Критическое число Рейнольдса
ReL1 = ro1 * V1 * L / mu1  # Среднее число Рейнольдса на грани
ReL2 = ro2 * V2 * L / mu2  # Среднее число Рейнольдса на грани

# если значения чисел Рейнолдса по граням больше критических => смешанный пограничный слой
if (ReL1 > Reкр1) and (ReL2 > Reкр2):
    print('\nСмешанный пограничный слой на 1 и 2 гранях')

# координаты точек перехода
xкр1 = Reкр1 * mu1 / (V1 * ro1)
xкр2 = Reкр2 * mu2 / (V2 * ro2)
print('\nКритические точки')
print('-------------------')
print('Reкр1=', Reкр1)
print('Reкр2=', Reкр2)
print('xкр1=', xкр1)
print('xкр2=', xкр2)


# Расчет смешанного ПС

# Грани №1 и №3

# Координаты характерных точек
x1 = xкр1 / 2  # Ламинарный ПС
x2 = xкр1  # Ламинарный ПС
x3 = xкр1  # Турбулентный ПС
x4 = xкр1 + (L - xкр1) / 3  # Турбулентный ПС
x5 = L - (L - xкр1) / 3  # Турбулентный ПС
x6 = L  # Турбулентный ПС
x7 = 0  # Турбулентный ПС
x8 = L / 2  # Турбулентный ПС
x9 = L  # Турбулентный ПС

# Грань №1

# Точки в ламинарном ПС
# Точка 1
# Параметры несжимаемого ПС
delt_nc1, tau_nc1, Cfx_nc1, Cf_nc1, delt_st_nc1, delt_st_st_nc1, Re1 = lam_nc(x1, ro1, V1, mu1)
# Параметры сжимаемого ПС
delt1, tau1, Cfx1, Cf1, delt_st1, delt_st_st1 = lam(delt_nc1, tau_nc1, Cfx_nc1, Cf_nc1, delt_st_nc1, delt_st_st_nc1, T1,
                                                    Toprl1)

# Точка 2
# Параметры несжимаемого ПС
delt_nc2, tau_nc2, Cfx_nc2, Cf_nc2, delt_st_nc2, delt_st_st_nc2, Re2 = lam_nc(x2, ro1, V1, mu1)
# Параметры сжимаемого ПС
delt2, tau2, Cfx2, Cf2, delt_st2, delt_st_st2 = lam(delt_nc2, tau_nc2, Cfx_nc2, Cf_nc2, delt_st_nc2, delt_st_st_nc2, T1,
                                                    Toprl1)

# Точки в турбулентном ПС
deltaX1 = ((delt2 ** 5 * V1 * ro1 * (T1 / Toprt1) ** (n + 1)) / (mu1 * 0.37 ** 5)) ** 0.25
x3f = deltaX1
x4f = (L - xкр1) / 3 + deltaX1
x5f = 2 * (L - xкр1) / 3 + deltaX1
x6f = L - xкр1 + deltaX1


# Точка 3
# Параметры несжимаемого ПС
delt_nc3, tau_nc3, Cfx_nc3, Cf_nc3, delt_st_nc3, delt_st_st_nc3, Re3 = tur_nc(x3f, ro1, V1, mu1, x3)
# Параметры сжимаемого ПС
delt3, tau3, Cfx3, Cf3, delt_st3, delt_st_st3 = tur(delt_nc3, tau_nc3, Cfx_nc3, Cf_nc3, delt_st_nc3, delt_st_st_nc3, T1,
                                                    Toprt1)

# Точка 4
# Параметры несжимаемого ПС
delt_nc4, tau_nc4, Cfx_nc4, Cf_nc4, delt_st_nc4, delt_st_st_nc4, Re4 = tur_nc(x4f, ro1, V1, mu1, x4)
# Параметры сжимаемого ПС
delt4, tau4, Cfx4, Cf4, delt_st4, delt_st_st4 = tur(delt_nc4, tau_nc4, Cfx_nc4, Cf_nc4, delt_st_nc4, delt_st_st_nc4, T1,
                                                    Toprt1)

# Точка 5
# Параметры несжимаемого ПС
delt_nc5, tau_nc5, Cfx_nc5, Cf_nc5, delt_st_nc5, delt_st_st_nc5, Re5 = tur_nc(x5f, ro1, V1, mu1, x5)
# Параметры сжимаемого ПС
delt5, tau5, Cfx5, Cf5, delt_st5, delt_st_st5 = tur(delt_nc5, tau_nc5, Cfx_nc5, Cf_nc5, delt_st_nc5, delt_st_st_nc5, T1,
                                                    Toprt1)

# Точка 6
# Параметры несжимаемого ПС
delt_nc6, tau_nc6, Cfx_nc6, Cf_nc6, delt_st_nc6, delt_st_st_nc6, Re6 = tur_nc(x6f, ro1, V1, mu1, x6)
# Параметры сжимаемого ПС
delt6, tau6, Cfx6, Cf6, delt_st6, delt_st_st6 = tur(delt_nc6, tau_nc6, Cfx_nc6, Cf_nc6, delt_st_nc6, delt_st_st_nc6, T1,
                                                    Toprt1)

# Грань  № 3

# Точки в турбулентном ПС
deltaX3 = ((delt6 ** 5 * V3 * ro3 * (T3 / Toprt3) ** (n + 1)) / (mu3 * 0.37 ** 5)) ** 0.25
print('deltaX3=', deltaX3)
x7f = deltaX3
x8f = x8 + deltaX3
x9f = x9 + deltaX3

# Точка 7
# Параметры несжимаемого ПС
delt_nc7, tau_nc7, Cfx_nc7, Cf_nc7, delt_st_nc7, delt_st_st_nc7, Re7 = tur_nc(x7f, ro3, V3, mu3, x7)
# Параметры сжимаемого ПС
delt7, tau7, Cfx7, Cf7, delt_st7, delt_st_st7 = tur(delt_nc7, tau_nc7, Cfx_nc7, Cf_nc7, delt_st_nc7, delt_st_st_nc7, T3,
                                                    Toprt3)

# Точка 8
# Параметры несжимаемого ПС
delt_nc8, tau_nc8, Cfx_nc8, Cf_nc8, delt_st_nc8, delt_st_st_nc8, Re8 = tur_nc(x8f, ro3, V3, mu3, x8)
# Параметры сжимаемого ПС
delt8, tau8, Cfx8, Cf8, delt_st8, delt_st_st8 = tur(delt_nc8, tau_nc8, Cfx_nc8, Cf_nc8, delt_st_nc8, delt_st_st_nc8, T3,
                                                    Toprt3)

# Точка 9
delt_nc9, tau_nc9, Cfx_nc9, Cf_nc9, delt_st_nc9, delt_st_st_nc9, Re9 = tur_nc(x9f, ro3, V3, mu3, x9)
delt9, tau9, Cfx9, Cf9, delt_st9, delt_st_st9 = tur(delt_nc9, tau_nc9, Cfx_nc9, Cf_nc9, delt_st_nc9, delt_st_st_nc9, T3,
                                                    Toprt3)

# Грани №2 и №4

# Координаты характерных точек
x10 = xкр2 / 2  # Ламинарный ПС
x11 = xкр2  # Ламинарный ПС
x12 = xкр2  # Турбулентный ПС
x13 = xкр2 + (L - xкр2) / 3  # Турбулентный ПС
x14 = L - (L - xкр2) / 3  # Турбулентный ПС
x15 = L  # Турбулентный ПС
x16 = 0  # Турбулентный ПС
x17 = L / 2  # Турбулентный ПС
x18 = L  # Турбулентный ПС

# Грань №2

# Точки в ламинарном ПС
# Точка 10
# Параметры несжимаемого ПС
delt_nc10, tau_nc10, Cfx_nc10, Cf_nc10, delt_st_nc10, delt_st_st_nc10, Re10 = lam_nc(x10, ro2, V2, mu2)
# Параметры сжимаемого ПС
delt10, tau10, Cfx10, Cf10, delt_st10, delt_st_st10 = lam(delt_nc10, tau_nc10, Cfx_nc10, Cf_nc10, delt_st_nc10,
                                                          delt_st_st_nc10, T2, Toprl2)

# Точка 11
# Параметры несжимаемого ПС
delt_nc11, tau_nc11, Cfx_nc11, Cf_nc11, delt_st_nc11, delt_st_st_nc11, Re11 = lam_nc(x11, ro2, V2, mu2)
# Параметры сжимаемого ПС
delt11, tau11, Cfx11, Cf11, delt_st11, delt_st_st11 = lam(delt_nc11, tau_nc11, Cfx_nc11, Cf_nc11, delt_st_nc11,
                                                          delt_st_st_nc11, T2, Toprl2)

# Точки в турбулентном ПС
deltaX2 = ((delt11 ** 5 * V2 * ro2 * (T2 / Toprt2) ** (n + 1)) / (mu2 * 0.37 ** 5)) ** 0.25
x12f = deltaX2
x13f = (L - xкр2) / 3 + deltaX2
x14f = 2 * (L - xкр2) / 3 + deltaX2
x15f = L - xкр2 + deltaX2

# Точка 12
# Параметры несжимаемого ПС
delt_nc12, tau_nc12, Cfx_nc12, Cf_nc12, delt_st_nc12, delt_st_st_nc12, Re12 = tur_nc(x12f, ro2, V2, mu2, x12)
# Параметры сжимаемого ПС
delt12, tau12, Cfx12, Cf12, delt_st12, delt_st_st12 = tur(delt_nc12, tau_nc12, Cfx_nc12, Cf_nc12, delt_st_nc12,
                                                          delt_st_st_nc12, T2, Toprt2)

# Точка 13
# Параметры несжимаемого ПС
delt_nc13, tau_nc13, Cfx_nc13, Cf_nc13, delt_st_nc13, delt_st_st_nc13, Re13 = tur_nc(x13f, ro2, V2, mu2, x13)
# Параметры сжимаемого ПС
delt13, tau13, Cfx13, Cf13, delt_st13, delt_st_st13 = tur(delt_nc13, tau_nc13, Cfx_nc13, Cf_nc13, delt_st_nc13,
                                                          delt_st_st_nc13, T2, Toprt2)

# Точка 14
# Параметры несжимаемого ПС
delt_nc14, tau_nc14, Cfx_nc14, Cf_nc14, delt_st_nc14, delt_st_st_nc14, Re14 = tur_nc(x14f, ro2, V2, mu2, x14)
# Параметры сжимаемого ПС
delt14, tau14, Cfx14, Cf14, delt_st14, delt_st_st14 = tur(delt_nc14, tau_nc14, Cfx_nc14, Cf_nc14, delt_st_nc14,
                                                          delt_st_st_nc14, T2, Toprt2)

# Точка 15
# Параметры несжимаемого ПС
delt_nc15, tau_nc15, Cfx_nc15, Cf_nc15, delt_st_nc15, delt_st_st_nc15, Re15 = tur_nc(x15f, ro2, V2, mu2, x15)
# Параметры сжимаемого ПС
delt15, tau15, Cfx15, Cf15, delt_st15, delt_st_st15 = tur(delt_nc15, tau_nc15, Cfx_nc15, Cf_nc15, delt_st_nc15,
                                                          delt_st_st_nc15, T2, Toprt2)

# Грань  № 4

# Точки в турбулентном ПС
deltaX4 = ((delt15 ** 5 * V4 * ro4 * (T4 / Toprt4) ** (n + 1)) / (mu4 * 0.37 ** 5)) ** 0.25
x16f = deltaX4
x17f = x17 + deltaX4
x18f = x18 + deltaX4

# Точка 16
# Параметры несжимаемого ПС
delt_nc16, tau_nc16, Cfx_nc16, Cf_nc16, delt_st_nc16, delt_st_st_nc16, Re16 = tur_nc(x16f, ro4, V4, mu4, x16)
# Параметры сжимаемого ПС
delt16, tau16, Cfx16, Cf16, delt_st16, delt_st_st16 = tur(delt_nc16, tau_nc16, Cfx_nc16, Cf_nc16, delt_st_nc16,
                                                          delt_st_st_nc16, T4, Toprt4)

# Точка 17
# Параметры несжимаемого ПС
delt_nc17, tau_nc17, Cfx_nc17, Cf_nc17, delt_st_nc17, delt_st_st_nc17, Re17 = tur_nc(x17f, ro4, V4, mu4, x17)
# Параметры сжимаемого ПС
delt17, tau17, Cfx17, Cf17, delt_st17, delt_st_st17 = tur(delt_nc17, tau_nc17, Cfx_nc17, Cf_nc17, delt_st_nc17,
                                                          delt_st_st_nc17, T4, Toprt4)

# Точка 18
delt_nc18, tau_nc18, Cfx_nc18, Cf_nc18, delt_st_nc18, delt_st_st_nc18, Re18 = tur_nc(x18f, ro4, V4, mu4, x18)
delt18, tau18, Cfx18, Cf18, delt_st18, delt_st_st18 = tur(delt_nc18, tau_nc18, Cfx_nc18, Cf_nc18, delt_st_nc18,
                                                          delt_st_st_nc18, T4, Toprt4)
# Средний коэффициент трения для граней

Cfср1 = Cfср12(xкр1, deltaX1, ro1, V1, mu1, Reкр1, Toprl1, Toprt1, T1)
Cfср2 = Cfср12(xкр2, deltaX2, ro2, V2, mu2, Reкр2, Toprl2, Toprt2, T2)

Cfср3 = Cfср34(deltaX3, ro3, V3, mu3, Toprt3, T3)
Cfср4 = Cfср34(deltaX4, ro4, V4, mu4, Toprt4, T4)

Cxsum = Cx + (Cfср1 + Cfср2 + Cfср3 + Cfср4) / 2
Cysum = Cy

Cxasum = Cxsum * cos(alpha) + Cysum * sin(alpha)
Cyasum = - Cxsum * sin(alpha) + Cysum * cos(alpha)

# Вывод по граням

print('\nГрань №1:', '\n---------')

print('\nТочка 1:', '\nx=', x1, '\nxф=', x1, '\nRe =', Re1, '\ndelta=', delt1, '\ntau=', tau1, '\nCfx=', Cfx1,
      '\nCf ср=', Cf1, '\ndelta*=', delt_st1, '\ndelta**=', delt_st_st1)

print('\nТочка 2:', '\nx=', x2, '\nxф=', x2, '\nRex=', Re2, '\ndelta=', delt2, '\ntau=', tau2, '\nCfx=', Cfx2,
      '\nCf ср=', Cf2, '\ndelta*=', delt_st2, '\ndelta**=', delt_st_st2)

print('\nТочка 3:', '\nx=', x3, '\nxф=', x3f, '\nRex=', Re3, '\ndelta=', delt3, '\ntau=', tau3, '\nCfx=', Cfx3,
      '\nCf ср=', Cf3, '\ndelta*=', delt_st3, '\ndelta**=', delt_st_st3)

print('\nТочка 4:', '\nx=', x4, '\nxф=', x4f, '\nRex=', Re4, '\ndelta=', delt4, '\ntau=', tau4, '\nCfx=', Cfx4,
      '\nCf ср=', Cf4, '\ndelta*=', delt_st4, '\ndelta**=', delt_st_st4)

print('\nТочка 5:', '\nx=', x5, '\nxф=', x5f, '\nRex=', Re5, '\ndelta=', delt5, '\ntau=', tau5, '\nCfx=', Cfx5,
      '\nCf ср=', Cf5, '\ndelta*=', delt_st5, '\ndelta**=', delt_st_st5)

print('\nТочка 6:', '\nx=', x6, '\nxф=', x6f, '\nRex=', Re6, '\ndelta=', delt6, '\ntau=', tau6, '\nCfx=', Cfx6,
      '\nCf ср=', Cf6, '\ndelta*=', delt_st6, '\ndelta**=', delt_st_st6)

print('\nГрань №3:', '\n---------')

print('\nТочка 7:', '\nx=', x7, '\nxф=', x7f, '\nRex=', Re7, '\ndelta=', delt7, '\ntau=', tau7, '\nCfx=', Cfx7,
      '\nCf ср=', Cf7, '\ndelta*=', delt_st7, '\ndelta**=', delt_st_st7)

print('\nТочка 8:', '\nx=', x8, '\nxф=', x8f, '\nRex=', Re8, '\ndelta=', delt8, '\ntau=', tau8, '\nCfx=', Cfx8,
      '\nCf ср=', Cf8, '\ndelta*=', delt_st8, '\ndelta**=', delt_st_st8)

print('\nТочка 9:', '\nx=', x9, '\nxф=', x9f, '\nRex=', Re9, '\ndelta=', delt9, '\ntau=', tau9, '\nCfx=', Cfx9,
      '\nCf ср=', Cf9, '\ndelta*=', delt_st9, '\ndelta**=', delt_st_st9)

print('\nГрань №2:', '\n---------')

print('\nТочка 10:', '\nx=', x10, '\nxф=', x10, '\nRe =', Re10, '\ndelta=', delt10, '\ntau=', tau10, '\nCfx=', Cfx10,
      '\nCf ср=', Cf10, '\ndelta*=', delt_st10, '\ndelta**=', delt_st_st10)

print('\nТочка 11:', '\nx=', x11, '\nxф=', x11, '\nRex=', Re11, '\ndelta=', delt11, '\ntau=', tau11, '\nCfx=', Cfx11,
      '\nCf ср=', Cf11, '\ndelta*=', delt_st11, '\ndelta**=', delt_st_st11)

print('\nТочка 12:', '\nx=', x12, '\nxф=', x12f, '\nRex=', Re12, '\ndelta=', delt12, '\ntau=', tau12, '\nCfx=', Cfx12,
      '\nCf ср=', Cf12, '\ndelta*=', delt_st12, '\ndelta**=', delt_st_st12)

print('\nТочка 13:', '\nx=', x13, '\nxф=', x13f, '\nRex=', Re13, '\ndelta=', delt13, '\ntau=', tau13, '\nCfx=', Cfx13,
      '\nCf ср=', Cf13, '\ndelta*=', delt_st13, '\ndelta**=', delt_st_st13)

print('\nТочка 14:', '\nx=', x14, '\nxф=', x14f, '\nRex=', Re14, '\ndelta=', delt14, '\ntau=', tau14, '\nCfx=', Cfx14,
      '\nCf ср=', Cf14, '\ndelta*=', delt_st14, '\ndelta**=', delt_st_st14)

print('\nТочка 15:', '\nx=', x15, '\nxф=', x15f, '\nRex=', Re15, '\ndelta=', delt15, '\ntau=', tau15, '\nCfx=', Cfx15,
      '\nCf ср=', Cf15, '\ndelta*=', delt_st15, '\ndelta**=', delt_st_st15)

print('\nГрань №4:', '\n---------')

print('\nТочка 16:', '\nx=', x16, '\nxф=', x16f, '\nRex=', Re16, '\ndelta=', delt16, '\ntau=', tau16, '\nCfx=', Cfx16,
      '\nCf ср=', Cf16, '\ndelta*=', delt_st16, '\ndelta**=', delt_st_st16)

print('\nТочка 17:', '\nx=', x17, '\nxф=', x17f, '\nRex=', Re17, '\ndelta=', delt17, '\ntau=', tau17, '\nCfx=', Cfx17,
      '\nCf ср=', Cf17, '\ndelta*=', delt_st17, '\ndelta**=', delt_st_st17)

print('\nТочка 18:', '\nx=', x18, '\nxф=', x18f, '\nRex=', Re18, '\ndelta=', delt18, '\ntau=', tau18, '\nCfx=', Cfx18,
      '\nCf ср=', Cf18, '\ndelta*=', delt_st18, '\ndelta**=', delt_st_st18)

print('\nСредние коэффициенты продольной силы от трения', '\n---------------------------------------------------')
print('\nCfср1', Cfср1, '\nCfср2', Cfср2, '\nCfср3', Cfср3, '\nCfср4', Cfср4)

print('\nПолные аэродинамические коэффициенты в связанной СК', '\n---------------------------------------------------')
print('\nCxsum', Cxsum, '\nCysum', Cysum)
print('\nПолные аэродинамические коэффициенты в скоростной СК', '\n---------------------------------------------------')
print('\nCxasum', Cxasum, '\nCyasum', Cyasum)