# Path: proptools/nozzle.py
# def thrust_coef(p_c, p_e, gamma, p_a=None, er=None):
# def c_star(gamma, m_molar, T_c):
# def er_from_p(p_c, p_e, gamma):
# def pressure_from_er(er, gamma):
# def throat_area(m_dot, p_c, T_c, gamma, m_molar):
# def mass_flow(A_t, p_c, T_c, gamma, m_molar):
# def thrust(A_t, p_c, p_e, gamma, p_a=None, er=None):
# def mach_from_er(er, gamma):
# def mach_from_pr(p_c, p_e, gamma):
# def is_choked(p_c, p_e, gamma):
# def mach_from_area_subsonic(area_ratio, gamma):
# def area_from_mach(M, gamma):
#     C_F = (2 * gamma**2 / (gamma - 1) \
#         * (2 / (gamma + 1))**((gamma + 1) / (gamma - 1)) \
#         * (1 - (p_e / p_c)**((gamma - 1) / gamma))
#           )**0.5
#     R = R_univ / m_molar
#     X = np.zeros((n,))
#     M = np.zeros((n,))
#     B = (y+1)/(y-1)
#     X[0] = (u*k)**(B/(1-B))
#     M[0] = X[0]
#     P = 2 / (gamma + 1)
#     Q = 1 - P
#     E = 1 / Q
#     R = area_ratio**2
#     X = fsolve(
#         lambda X: (P + Q * X)**E - R * X,
#         X_init
#         )


from scipy.misc import derivative
from scipy.integrate import ode
from proptools import nozzle
    from matplotlib import pyplot as plt
import numpy as np
    return (T_o, M, choked)


def main():
    def f_f(x):
        return 0
    def f_q(x):
        return 0
    def f_A(x):
        return 1 + x
    mdot = 100
    c_p = 2000
    gamma = 1.4
    R = c_p * (1 - 1 / gamma)
    x = np.linspace(0, 1)
    T_o_in = 300

    for M_in in [0.2, 0.8, 1.2, 2]:
        T_in = T_o_in * (1 + (gamma - 1) / 2 * M_in**2)**-1
        v_in = M_in * (gamma * R * T_in)**0.5
        rho_in = mdot / (v_in * f_A(0))
        p_in = rho_in * R * T_in
        p_o_in = p_in * (T_o_in / T_in)**(gamma / (gamma -1))
        

        A_t = nozzle.throat_area(mdot, p_o_in, T_o_in, gamma, nozzle.R_univ / R)

        if M_in < 1:
            M_area = [nozzle.mach_from_area_subsonic(f_A(xi) / A_t, gamma) for xi in x]
        else:

# Path: proptools/nozzle.py
# def thrust_coef(p_c, p_e, gamma, p_a=None, er=None):
# def c_star(gamma, m_molar, T_c):
# def er_from_p(p_c, p_e, gamma):
# def pressure_from_er(er, gamma):
# def throat_area(m_dot, p_c, T_c, gamma, m_molar):
# def mass_flow(A_t, p_c, T_c, gamma, m_molar):
# def thrust(A_t, p_c, p_e, gamma, p_a=None, er=None):
# def mach_from_er(er, gamma):
# def mach_from_pr(p_c, p_e, gamma):
# def is_choked(p_c, p_e, gamma):
# def mach_from_area_subsonic(area_ratio, gamma):
# def area_from_mach(M, gamma):
#     C_F = (2 * gamma**2 / (gamma - 1) \
#         * (2 / (gamma + 1))**((gamma + 1) / (gamma - 1)) \
#         * (1 - (p_e / p_c)**((gamma - 1) / gamma))
#           )**0.5
#     R = R_univ / m_molar
#     X = np.zeros((n,))
#     M = np.zeros((n,))
#     B = (y+1)/(y-1)
#     X[0] = (u*k)**(B/(1-B))
#     M[0] = X[0]
#     P = 2 / (gamma + 1)
#     Q = 1 - P
#     E = 1 / Q
#     R = area_ratio**2
#     X = fsolve(
#         lambda X: (P + Q * X)**E - R * X,
#         X_init
#         )


# Path: proptools/nonsimple_comp_flow.py
from scipy.misc import derivative
from scipy.integrate import ode
from proptools import nozzle
    from matplotlib import pyplot as plt
import numpy as np

    return (T_o, M, choked)


def main():
    def f_f(x):
        return 0
    def f_q(x):
        return 0
    def f_A(x):
        return 1 + x
    mdot = 100
    c_p = 2000
    gamma = 1.4
    R = c_p * (1 - 1 / gamma)
    x = np.linspace(0, 1)
    T_o_in = 300

    for M_in in [0.2, 0.8, 1.2, 2]:
        T_in = T_o_in * (1 + (gamma - 1) / 2 * M_in**2)**-1
        v_in = M_in * (gamma * R * T_in)**0.5
        rho_in = mdot / (v_in * f_A(0))
        p_in = rho_in * R * T_in
        p_o_in = p_in * (T_o_in / T_in)**(gamma / (gamma -1))
        

        A_t = nozzle.throat_area(mdot, p_o_in, T_o_in, gamma, nozzle.R_univ / R)

        if M_in < 1:
            M_area = [nozzle.mach_from_area_subsonic(f_A(xi) / A_t, gamma) for xi in x]
        else:

            M_area = [nozzle.mach_from_er(f_A(xi) / A_t, gamma) for xi in x]