Untitled

import math
import numpy as np
from scipy.optimize import root
# Исходные данные для КР
N = 11
c = 90
b = 1100
h = (10 + 0.5 * N) * 1000
M = 4+0.15 * N
alfa = 2 + 0.1 * N
betta = math.tan(c/b)
Tct = 405
Rekr = 5 * 10 ** 6
# Исходные данные для атмосферы
r = 6356767
g0 = 9.80665
R = 287.05287
i = 0
# Функции атмосферы
def geopotencial_v():
    H = (r * h) / (r + h)  # Геопотенциальная высота
    return H
def yskorenie_s_p():
    g = g0 * (r / (r + h)) ** 2  # Ускорение свободного падения
    return g
def parametry_t():
    # Параметры, для определения температуры
    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 [Hx, batm, Tx, px]
def temperatura():
    T = Tx + batm * (H - Hx)       # Температура
    return T
def davlenie():
    if batm != 0:  # Давление
        p = 10 ** (math.log10(px) - (g0 / (batm * R)) * math.log10(T / Tx))
    else:
        p = 10 ** (math.log10(px) - ((0.434294 * g0 * (H - Hx)) / (T * R)))
    return p
def plotnost():
    ro = p / (R * T)  # Плотность
    return ro
def skorost_z():
    a = 20.046796 * math.sqrt(T)  # Скорость звука
    return a
def vivod():
    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(' ')
    return
H = geopotencial_v()
g = yskorenie_s_p()
list = parametry_t()
Hx = list [0]
batm = list [1]
Tx = list [2]
px = list [3]
T = temperatura()
p = davlenie()
ro = plotnost()
a = skorost_z()
vivod = vivod()
#начало КР
Vinf = a * M # Скорость потока

# 1 и 2 области
# Углы отклонения потока
betta1 = betta - alfa
betta2 = betta + alfa
# найдем угол СУ

'''def ugol_SY():

    return tettaci'''
tettac1 = 0
tettac2 = 0
k = 1.4

def ugol_SY(tettac1):
    return np.tan(tettac1)/np.tan(tettac1-betta1)-((k+1)*M**2*(math.sin(tettac1))**2)/(2+(k-1)*M**2*(math.sin(tettac1))**2)
crange = range(1, 10)
res = [root(ugol_SY, 0.3, args=(ci, )).x[0] for ci in crange]

print(res)