# 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
#         )


import numpy as np
import skaero.atmosphere.coesa as atmo
from matplotlib import pyplot as plt
from proptools import nozzle
"""Plot C_F vs altitude."""

p_c = 10e6   # Chamber pressure [units: pascal]
gamma = 1.2    # Exhaust heat capacity ratio [units: dimensionless]
p_e_1 = 100e3    # Nozzle exit pressure, 1st stage [units: pascal]
exp_ratio_1 = nozzle.er_from_p(p_c, p_e_1, gamma)    # Nozzle expansion ratio [units: dimensionless]
p_e_2 = 15e3    # Nozzle exit pressure, 2nd stage [units: pascal]
exp_ratio_2 = nozzle.er_from_p(p_c, p_e_2, gamma)    # Nozzle expansion ratio [units: dimensionless]

alt = np.linspace(0, 84e3)    # Altitude [units: meter]
p_a = atmo.pressure(alt)    # Ambient pressure [units: pascal]

# Compute the thrust coeffieicient of the fixed-area nozzle, 1st stage [units: dimensionless]

# 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: docs/source/examples/plots/cf_alt.py
import numpy as np
import skaero.atmosphere.coesa as atmo
from matplotlib import pyplot as plt
from proptools import nozzle

"""Plot C_F vs altitude."""

p_c = 10e6   # Chamber pressure [units: pascal]
gamma = 1.2    # Exhaust heat capacity ratio [units: dimensionless]
p_e_1 = 100e3    # Nozzle exit pressure, 1st stage [units: pascal]
exp_ratio_1 = nozzle.er_from_p(p_c, p_e_1, gamma)    # Nozzle expansion ratio [units: dimensionless]
p_e_2 = 15e3    # Nozzle exit pressure, 2nd stage [units: pascal]
exp_ratio_2 = nozzle.er_from_p(p_c, p_e_2, gamma)    # Nozzle expansion ratio [units: dimensionless]

alt = np.linspace(0, 84e3)    # Altitude [units: meter]
p_a = atmo.pressure(alt)    # Ambient pressure [units: pascal]

# Compute the thrust coeffieicient of the fixed-area nozzle, 1st stage [units: dimensionless]

C_F_fixed_1 = nozzle.thrust_coef(p_c, p_e_1, gamma, p_a=p_a, er=exp_ratio_1)