import math


def prog_1():

    # ============================================================
    # FEUILLE 1 : CHOIX D'UN DIAPHRAGME
    # ============================================================

    # ============================================================
    # 1. DONNEES A MODIFIER
    # ============================================================

    # Conditions d'essai
    p_kPa = 700.0          # pression statique amont [kPa]
    Tt = 473.0             # température totale [K]
    qref = 2.50            # débit massique imposé [kg/s]

    # Géométrie
    d0_mm = 110.0          # diamètre du diaphragme [mm]
    D0_mm = 250.0          # diamètre intérieur conduite [mm]

    # Erreurs absolues
    erreur_d_mm = 0.5
    erreur_D_mm = 0.5
    erreur_p_kPa = 0.55
    erreur_dp_kPa = 0.2
    erreur_T_K = 10.0


    # ============================================================
    # 2. CONSTANTES
    # ============================================================

    r = 287.04
    mu0 = 18.2e-6
    alpha = 17.5e-6
    pi = math.pi

    T_ref_dilatation = 288.15
    T_diaphragme = 273.0

    tol = 1.0e-10

    max_iter_C = 100
    max_iter_T = 100


    # ============================================================
    # 3. CONVERSION DES DIAMETRES
    # ============================================================

    d0 = d0_mm / 1000.0
    D0 = D0_mm / 1000.0


    # ============================================================
    # 4. VERIFICATION DE LA TEMPERATURE
    # ============================================================

    if Tt < 200.0:
        print("ERREUR : température trop basse")
        return


    # ============================================================
    # 5. PREPARATION
    # ============================================================

    p = p_kPa * 1000.0

    C_old = 0.6


    # ============================================================
    # 6. CORRECTION THERMIQUE DES DIAMETRES
    # ============================================================

    D = D0 * (
        1.0
        + alpha * (T_diaphragme - T_ref_dilatation)
    )

    d = d0 * (
        1.0
        + alpha * (T_diaphragme - T_ref_dilatation)
    )

    beta = d / D


    # ============================================================
    # 7. CAS q = 0
    # ============================================================

    if qref == 0.0:
        qref = 1.0e-5


    # ============================================================
    # 8. INITIALISATION
    # ============================================================

    delta_p = 0.1
    delta_p_old = 1.0

    T = Tt

    flag = False

    iteration_T = 0
    iteration_C = 0

    erreur_C = float("inf")
    erreur_T = float("inf")
    erreur_delta_p = float("inf")

    # Variables initialisées pour éviter tout problème
    C = C_old
    Re_D = 0.0
    V = 0.0
    Mach = 0.0
    T_new = T
    epsilon = 1.0
    rho = 0.0
    Cp = 0.0
    gamma = 1.4
    mu = 0.0
    E = 0.0
    s = 0.0


    # ============================================================
    # 9. BOUCLE EXTERIEURE SUR LA TEMPERATURE
    # ============================================================

    for iteration_T in range(1, max_iter_T + 1):

        # --------------------------------------------------------
        # 9.1 MASSE VOLUMIQUE
        # --------------------------------------------------------

        rho = p / (r * T)


        # --------------------------------------------------------
        # 9.2 SECTION DU DIAPHRAGME
        # --------------------------------------------------------

        s = pi * d**2 / 4.0


        # --------------------------------------------------------
        # 9.3 CAPACITE THERMIQUE Cp
        # --------------------------------------------------------

        exp_term = math.exp(3090.0 / T)

        Cp = r * (
            3.5
            - 2.8e-5 * T
            + 2.24e-8 * T**2
            + (3090.0 / T)**2
            * exp_term
            / (exp_term - 1.0)**2
        )


        # --------------------------------------------------------
        # 9.4 GAMMA
        # --------------------------------------------------------

        gamma = Cp / (Cp - r)


        # --------------------------------------------------------
        # 9.5 COEFFICIENT DE DETENTE epsilon
        # --------------------------------------------------------

        rapport_pression = (p - delta_p) / p

        if rapport_pression <= 0.0:
            print("ERREUR : pression aval négative ou nulle")
            return

        epsilon = 1.0 - (
            0.351
            + 0.256 * beta**4
            + 0.93 * beta**8
        ) * (
            1.0
            - rapport_pression**(1.0 / gamma)
        )


        # --------------------------------------------------------
        # 9.6 COEFFICIENT D'APPROCHE
        # --------------------------------------------------------

        E = 1.0 / math.sqrt(
            1.0 - beta**4
        )


        # --------------------------------------------------------
        # 9.7 VISCOSITE
        # --------------------------------------------------------

        mu = (
            mu0
            * (T / 293.15)**1.5
            * (113.0 + 293.15)
            / (113.0 + T)
        )


        # ========================================================
        # 10. BOUCLE INTERIEURE SUR C
        # ========================================================

        erreur_C = float("inf")

        for iteration_C in range(1, max_iter_C + 1):

            # ----------------------------------------------------
            # 10.1 CALCUL DE Delta p
            # ----------------------------------------------------

            delta_p = (
                (
                    qref
                    / (C_old * E * epsilon * s)
                )**2
                / (2.0 * rho)
            )


            # ----------------------------------------------------
            # 10.2 VERIFICATION Delta p / p
            # ----------------------------------------------------

            if 4.0 * delta_p > p:

                flag = True
                break


            # ----------------------------------------------------
            # 10.3 NOMBRE DE REYNOLDS
            # ----------------------------------------------------

            Re_D = (
                4.0 * qref
                / (pi * D * mu)
            )


            # ----------------------------------------------------
            # 10.4 TERME A1
            # ----------------------------------------------------

            A1 = (
                (19000.0 * beta) / Re_D
            )**0.8


            # ----------------------------------------------------
            # 10.5 COEFFICIENT DE DECHARGE C
            # ----------------------------------------------------

            C = (
                0.5961
                + 0.0261 * beta**2
                - 0.216 * beta**8
                + 0.000521
                * (beta * 1.0e6 / Re_D)**0.7
                + (
                    0.0188
                    + 0.0063 * A1
                )
                * beta**3.5
                * (1.0e6 / Re_D)**0.3
            )


            # ----------------------------------------------------
            # 10.6 CONVERGENCE SUR C
            # ----------------------------------------------------

            erreur_C = abs(
                1.0 - C_old / C
            )

            if erreur_C <= tol:

                C_old = C
                break


            C_old = C


        # ========================================================
        # FIN BOUCLE C
        # ========================================================

        # Il faut tester le flag AVANT la convergence de C
        # pour éviter une erreur de variable non définie.

        if flag:
            break


        if (
            iteration_C == max_iter_C
            and erreur_C > tol
        ):

            print(
                "ERREUR : "
                "la boucle sur C n'a pas convergé"
            )

            return


        # ========================================================
        # 11. VITESSE MOYENNE
        # ========================================================

        V = (
            4.0
            * qref
            * T
            * r
            / (pi * D**2 * p)
        )


        # ========================================================
        # 12. VITESSE DU SON
        # ========================================================

        a = math.sqrt(
            gamma * r * T
        )


        # ========================================================
        # 13. NOMBRE DE MACH
        # ========================================================

        Mach = V / a


        # ========================================================
        # 14. NOUVELLE TEMPERATURE STATIQUE
        # ========================================================

        T_new = Tt / (
            1.0
            + (
                (gamma - 1.0)
                / 2.0
            )
            * Mach**2
        )


        # ========================================================
        # 15. CONVERGENCE SUR LA TEMPERATURE
        # ========================================================

        erreur_T = abs(
            1.0 - T / T_new
        )

        if erreur_T <= tol:

            T = T_new
            break


        T = T_new


        # ========================================================
        # 16. CONVERGENCE SUR Delta p
        # ========================================================

        erreur_delta_p = abs(
            1.0
            - delta_p / delta_p_old
        )

        if erreur_delta_p <= tol:
            break


        delta_p_old = delta_p


    # ============================================================
    # FIN BOUCLE TEMPERATURE
    # ============================================================

    if flag:

        print()
        print("=" * 70)
        print("ERREUR")
        print("=" * 70)

        print("DeltaP / P > 0.25")

        print(
            f"Delta p = "
            f"{delta_p / 1000.0:.6f} kPa"
        )

        return


    # ============================================================
    # 17. PRESSION TOTALE
    # ============================================================

    if Mach < 0.2:

        p_total = (
            p
            + 0.5 * rho * V**2
        )

    else:

        p_total = (
            p
            * (
                1.0
                + (
                    gamma - 1.0
                )
                / 2.0
                * Mach**2
            )**(
                gamma
                / (gamma - 1.0)
            )
        )


    # ============================================================
    # 18. VERIFICATION DU NOMBRE DE REYNOLDS
    # ============================================================

    avertissement_Re = ""

    if (
        Re_D <= 5000.0
        and 0.1 <= beta <= 0.559
    ):

        avertissement_Re = (
            "Reynolds trop faible"
        )


    if (
        Re_D <= 16000.0 * beta**2
        and beta > 0.559
    ):

        avertissement_Re = (
            "Reynolds trop faible"
        )


    # ============================================================
    # 19. CHANGEMENT D'UNITES
    # ============================================================

    delta_p_kPa = (
        delta_p / 1000.0
    )

    p_bar = (
        p / 1.0e5
    )

    d_mm = (
        d * 1000.0
    )

    D_mm = (
        D * 1000.0
    )


    # ============================================================
    # 20. INCERTITUDE SUR LA GEOMETRIE
    # ============================================================

    theta_d = (
        2.0
        + 0.11 * beta**2
        + 2.0 * beta**4
        - 0.5 * beta**8
        + 7.0 * beta**16
    )

    theta_D = (
        2.0 - theta_d
    )


    contribution_d = (
        theta_d
        * erreur_d_mm
        / d_mm
    )


    contribution_D = (
        theta_D
        * erreur_D_mm
        / D_mm
    )


    # ============================================================
    # 21. INCERTITUDE SUR C
    # ============================================================

    ec = 0.0


    if 0.1 <= beta < 0.2:

        ec = (
            0.7 - beta
        ) / 100.0


    elif 0.2 <= beta <= 0.6:

        ec = (
            0.5 / 100.0
        )


    elif 0.6 < beta <= 0.75:

        ec = (
            1.667 * beta
            - 0.5
        ) / 100.0


    # Petit diamètre
    if D_mm < 71.12:

        ec = ec + (
            0.9
            * (0.75 - beta)
            * (
                2.8
                - D_mm / 25.4
            )
        ) / 100.0


    # Faible Reynolds
    if (
        beta > 0.5
        and Re_D < 10000.0
    ):

        ec = ec + 0.005


    # ============================================================
    # 22. CONTRIBUTION GLOBALE GEOMETRIE + C
    # ============================================================

    eqbeta = math.sqrt(
        contribution_d**2
        + contribution_D**2
        + ec**2
    )


    # ============================================================
    # 23. INCERTITUDE SUR Delta p
    # ============================================================

    theta_p = (
        0.5
        + (
            (
                0.44
                + 0.4 * beta**4
            )
            * delta_p_kPa
        )
        / (
            p_bar
            * 100.0
            * 1.4
        )
    )


    theta_delta_p = (
        1.0 - theta_p
    )


    if delta_p_kPa == 0.0:

        eqdp = 0.0

    else:

        eqdp = (
            theta_delta_p
            * erreur_dp_kPa
            / delta_p_kPa
        )


    # ============================================================
    # 24. INCERTITUDE SUR LA PRESSION
    # ============================================================

    eqp = (
        theta_p
        * (erreur_p_kPa / 100.0)
        / (100.0 * p_bar)
    )


    # ============================================================
    # 25. INCERTITUDE SUR LA TEMPERATURE
    # ============================================================

    eqT = (
        0.5
        * erreur_T_K
        / Tt
    )


    # ============================================================
    # 26. INCERTITUDE SUR epsilon
    # ============================================================

    eqepsilon = (
        3.5
        * delta_p_kPa
        / (
            gamma
            * p_bar
            * 100.0
        )
    ) / 100.0


    # ============================================================
    # 27. INCERTITUDE FINALE
    # ============================================================

    incertitude = math.sqrt(
        eqbeta**2
        + eqp**2
        + eqdp**2
        + eqT**2
        + eqepsilon**2
    )


    # ============================================================
    # 28. AFFICHAGE DES RESULTATS
    # ============================================================

    print()
    print("=" * 70)
    print("FEUILLE 1 - CHOIX D'UN DIAPHRAGME")
    print("=" * 70)


    print()
    print("ENTREES")
    print("-" * 70)

    print(
        f"Pression amont p           = "
        f"{p_kPa:.3f} kPa"
    )

    print(
        f"Température totale Tt      = "
        f"{Tt:.3f} K"
    )

    print(
        f"Débit imposé q             = "
        f"{qref:.6f} kg/s"
    )

    print(
        f"Diamètre diaphragme d0     = "
        f"{d0_mm:.3f} mm"
    )

    print(
        f"Diamètre conduite D0       = "
        f"{D0_mm:.3f} mm"
    )


    print()
    print("GEOMETRIE CORRIGEE")
    print("-" * 70)

    print(
        f"d                           = "
        f"{d_mm:.6f} mm"
    )

    print(
        f"D                           = "
        f"{D_mm:.6f} mm"
    )

    print(
        f"beta = d/D                  = "
        f"{beta:.8f}"
    )

    print(
        f"Section diaphragme s        = "
        f"{s:.9e} m²"
    )

    print(
        f"Coefficient approche E      = "
        f"{E:.8f}"
    )


    print()
    print("RESULTAT PRINCIPAL")
    print("-" * 70)

    print(
        f"Delta p                     = "
        f"{delta_p_kPa:.6f} kPa"
    )


    print()
    print("PROPRIETES DE L'AIR")
    print("-" * 70)

    print(
        f"rho                         = "
        f"{rho:.6f} kg/m³"
    )

    print(
        f"Cp                          = "
        f"{Cp:.3f} J/(kg.K)"
    )

    print(
        f"gamma                       = "
        f"{gamma:.8f}"
    )

    print(
        f"epsilon                     = "
        f"{epsilon:.8f}"
    )

    print(
        f"viscosité mu                = "
        f"{mu:.9e} Pa.s"
    )


    print()
    print("ECOULEMENT")
    print("-" * 70)

    print(
        f"Coefficient de décharge C   = "
        f"{C:.8f}"
    )

    print(
        f"Reynolds                    = "
        f"{Re_D:.0f}"
    )

    print(
        f"Vitesse                     = "
        f"{V:.6f} m/s"
    )

    print(
        f"Mach                        = "
        f"{Mach:.8f}"
    )

    print(
        f"Température statique        = "
        f"{T_new:.6f} K"
    )

    print(
        f"Pression totale             = "
        f"{p_total / 1000.0:.6f} kPa"
    )


    print()
    print("INCERTITUDES")
    print("-" * 70)

    print(
        f"Contribution géométrie + C  = "
        f"{eqbeta * 100:.4f} %"
    )

    print(
        f"Contribution pression       = "
        f"{eqp * 100:.4f} %"
    )

    print(
        f"Contribution Delta p        = "
        f"{eqdp * 100:.4f} %"
    )

    print(
        f"Contribution température    = "
        f"{eqT * 100:.4f} %"
    )

    print(
        f"Contribution epsilon        = "
        f"{eqepsilon * 100:.4f} %"
    )

    print()

    print(
        f"Incertitude finale          = "
        f"{incertitude * 100:.4f} %"
    )


    print()
    print("CONVERGENCE")
    print("-" * 70)

    print(
        f"Nombre d'itérations sur T   = "
        f"{iteration_T}"
    )

    print(
        f"Nombre d'itérations sur C   = "
        f"{iteration_C}"
    )


    if avertissement_Re != "":

        print()
        print(
            "ATTENTION :",
            avertissement_Re
        )







def prog_2():

    # ============================================================
    # FEUILLE 2 : CALCUL DU DEBIT
    # ============================================================
    #
    # Sens du calcul :
    #
    #       Delta_p, p, Td, Tt, d0, D0
    #                     |
    #                     v
    #                     q
    #
    # La différence de pression Delta_p est IMPOSEE.
    # Le débit massique q est CALCULE.
    #
    # ============================================================


    # ============================================================
    # 1. DONNEES A MODIFIER
    # ============================================================

    # Conditions d'essai
    delta_p_kPa = 35    # différence de pression [kPa]
    p_kPa = 1900        # pression statique amont [kPa]

    Td = 680            # température du diaphragme [K]
    Tt = 760             # température totale de l'air [K]


    # Géométrie
    d0_mm = 249.6           # diamètre du diaphragme [mm]
    D0_mm = 400.0           # diamètre intérieur conduite [mm]


    # ============================================================
    # 2. CONSTANTES
    # ============================================================

    r = 287.04

    mu0 = 18.2e-6

    alpha = 17.5e-6

    pi = math.pi

    T_ref_dilatation = 288.15

    # Dans le VBA de la feuille 2 :
    #
    # convergence sur C : 10^-6
    # convergence sur T : 10^-6

    tol = 1.0e-6

    max_iter_C = 100
    max_iter_T = 100


    # ============================================================
    # 3. CONVERSION DES DIAMETRES
    # ============================================================

    d0 = d0_mm / 1000.0

    D0 = D0_mm / 1000.0


    # ============================================================
    # 4. CONVERSION DES PRESSIONS
    # ============================================================

    delta_p = delta_p_kPa * 1000.0

    p = p_kPa * 1000.0


    # ============================================================
    # 5. VERIFICATION DES ENTREES
    # ============================================================

    if Td <= 0.0:

        print("ERREUR : température du diaphragme <= 0")
        return


    if Tt <= 0.0:

        print("ERREUR : température totale <= 0")
        return


    if p <= 0.0:

        print("ERREUR : pression amont <= 0")
        return


    # ============================================================
    # 6. CAS Delta p <= 0
    # ============================================================
    #
    # VBA :
    #
    # If dpdiaph <= 0 Then dpdiaph = 10 ^ -6
    #
    # Attention :
    # la valeur est en Pa à ce stade.

    if delta_p <= 0.0:

        delta_p = 1.0e-6


    # ============================================================
    # 7. CORRECTION THERMIQUE DES DIAMETRES
    # ============================================================

    # VBA :
    #
    # d2 = d22 * (1 + coefdil * (tmdiaph - 288.15))
    #
    # donc :
    #
    # D = D0 [1 + alpha (Td - Tref)]

    D = D0 * (
        1.0
        + alpha
        * (Td - T_ref_dilatation)
    )


    # VBA :
    #
    # d1 = d11 * (1 + coefdil * (tmdiaph - 288.15))

    d = d0 * (
        1.0
        + alpha
        * (Td - T_ref_dilatation)
    )


    # ============================================================
    # 8. RAPPORT BETA
    # ============================================================
    #
    # ATTENTION :
    #
    # le VBA écrit :
    #
    # beta = d11 / d22
    #
    # Il utilise donc les diamètres DE REFERENCE et non
    # les diamètres corrigés thermiquement.
    #
    # Comme le même coefficient de dilatation et la même
    # température sont appliqués à d et D :
    #
    # d / D = d0 / D0
    #
    # les deux méthodes donnent ici le même beta.

    beta = d0 / D0


    # ============================================================
    # 9. VERIFICATION DU DOMAINE DE BETA
    # ============================================================

    if beta <= 0.1 or beta >= 0.75:

        print()
        print("=" * 70)
        print("ERREUR")
        print("=" * 70)

        print("Beta hors tolérance")

        print(
            f"Beta = "
            f"{beta:.6f}"
        )

        return


    # ============================================================
    # 10. VERIFICATION Delta p / p
    # ============================================================
    #
    # VBA :
    #
    # If 4 * dpdiaph > ptamont Then
    #
    # équivalent à :
    #
    # Delta_p / p > 0.25

    if 4.0 * delta_p > p:

        print()
        print("=" * 70)
        print("ERREUR")
        print("=" * 70)

        print("DeltaP / P > 0.25")

        print(
            f"Delta p / p = "
            f"{delta_p / p:.6f}"
        )

        return


    # ============================================================
    # 11. INITIALISATION
    # ============================================================

    C_old = 0.6

    T = Tt

    iteration_T = 0
    iteration_C = 0

    erreur_C = float("inf")
    erreur_T = float("inf")


    # Variables initialisées
    q = 0.0

    Re_D = 0.0

    C = C_old

    V = 0.0

    Mach = 0.0

    T_new = T

    rho = 0.0

    Cp = 0.0

    gamma = 1.4

    epsilon = 1.0

    E = 0.0

    mu = 0.0

    s = 0.0


    # ============================================================
    # 12. BOUCLE EXTERIEURE SUR LA TEMPERATURE
    # ============================================================

    for iteration_T in range(
        1,
        max_iter_T + 1
    ):


        # --------------------------------------------------------
        # 12.1 MASSE VOLUMIQUE
        # --------------------------------------------------------
        #
        # VBA :
        #
        # rho = ptamont / (r * t0)
        #
        # Physique :
        #
        # rho = p / (r T)
        #
        # Le modèle considère ici Z = 1.

        rho = p / (
            r * T
        )


        # --------------------------------------------------------
        # 12.2 SECTION DU DIAPHRAGME
        # --------------------------------------------------------
        #
        # s = pi d² / 4

        s = (
            pi
            * d**2
            / 4.0
        )


        # --------------------------------------------------------
        # 12.3 CAPACITE THERMIQUE Cp
        # --------------------------------------------------------

        exp_term = math.exp(
            3090.0 / T
        )


        Cp = r * (

            3.5

            - 2.8e-5 * T

            + 2.24e-8 * T**2

            + (
                3090.0 / T
            )**2

            * exp_term

            / (
                exp_term - 1.0
            )**2
        )


        # --------------------------------------------------------
        # 12.4 GAMMA
        # --------------------------------------------------------
        #
        # gamma = Cp / Cv
        #
        # avec :
        #
        # Cv = Cp - r

        gamma = (
            Cp
            / (Cp - r)
        )


        # --------------------------------------------------------
        # 12.5 COEFFICIENT DE DETENTE epsilon
        # --------------------------------------------------------

        rapport_pression = (
            p - delta_p
        ) / p


        if rapport_pression <= 0.0:

            print(
                "ERREUR : "
                "pression aval négative ou nulle"
            )

            return


        epsilon = 1.0 - (

            0.351

            + 0.256
            * beta**4

            + 0.93
            * beta**8

        ) * (

            1.0

            - rapport_pression
            ** (1.0 / gamma)

        )


        # --------------------------------------------------------
        # 12.6 COEFFICIENT D'APPROCHE E
        # --------------------------------------------------------
        #
        # E = 1 / sqrt(1 - beta^4)

        E = (
            1.0
            / math.sqrt(
                1.0 - beta**4
            )
        )


        # --------------------------------------------------------
        # 12.7 VISCOSITE DYNAMIQUE
        # --------------------------------------------------------
        #
        # Loi de Sutherland

        mu = (

            mu0

            * (
                T / 293.15
            )**1.5

            * (
                113.0
                + 293.15
            )

            / (
                113.0
                + T
            )
        )


        # ========================================================
        # 13. BOUCLE INTERIEURE SUR C
        # ========================================================

        erreur_C = float("inf")


        for iteration_C in range(
            1,
            max_iter_C + 1
        ):


            # ----------------------------------------------------
            # 13.1 CALCUL DU DEBIT
            # ----------------------------------------------------
            #
            # VBA :
            #
            # q =
            # co * e * epsi * s
            # * sqrt(2 * dpdiaph * rho)
            #
            # Physique :
            #
            # q =
            # C E epsilon s sqrt(2 rho Delta_p)

            q = (

                C_old

                * E

                * epsilon

                * s

                * math.sqrt(
                    2.0
                    * delta_p
                    * rho
                )
            )


            # ----------------------------------------------------
            # 13.2 NOMBRE DE REYNOLDS
            # --------------------------------------------------------
            #
            # Re_D =
            #
            # 4 q / (pi D mu)

            Re_D = (

                4.0
                * q

                / (
                    pi
                    * D
                    * mu
                )
            )


            # ----------------------------------------------------
            # 13.3 TERME A1
            # --------------------------------------------------------

            A1 = (

                (
                    19000.0
                    * beta
                )

                / Re_D

            )**0.8


            # ----------------------------------------------------
            # 13.4 COEFFICIENT DE DECHARGE C
            # --------------------------------------------------------

            C = (

                0.5961

                + 0.0261
                * beta**2

                - 0.216
                * beta**8

                + 0.000521
                * (
                    beta
                    * 1.0e6
                    / Re_D
                )**0.7

                + (
                    0.0188
                    + 0.0063
                    * A1
                )

                * beta**3.5

                * (
                    1.0e6
                    / Re_D
                )**0.3
            )


            # ----------------------------------------------------
            # 13.5 CONVERGENCE SUR C
            # --------------------------------------------------------
            #
            # VBA :
            #
            # If Abs(1 - co / c) <= 10^-6 Then Exit Do

            erreur_C = abs(
                1.0
                - C_old / C
            )


            if erreur_C <= tol:

                C_old = C

                break


            C_old = C


        # ========================================================
        # FIN BOUCLE C
        # ========================================================

        if (
            iteration_C == max_iter_C
            and erreur_C > tol
        ):

            print(
                "ERREUR : "
                "la boucle sur C n'a pas convergé"
            )

            return


        # ========================================================
        # 14. VITESSE MOYENNE
        # ========================================================
        #
        # VBA :
        #
        # vitis =
        # 4*q*t0*287.04 /
        # (pi*d2²*ptamont)
        #
        # Equivalent ici à :
        #
        # V = q / (rho S)

        V = (

            4.0

            * q

            * T

            * 287.04

            / (
                pi
                * D**2
                * p
            )
        )


        # ========================================================
        # 15. VITESSE DU SON
        # ========================================================
        #
        # a = sqrt(gamma r T)

        a = math.sqrt(
            gamma
            * r
            * T
        )


        # ========================================================
        # 16. NOMBRE DE MACH
        # ========================================================

        Mach = (
            V / a
        )


        # ========================================================
        # 17. TEMPERATURE STATIQUE
        # ========================================================
        #
        # T =
        #
        # Tt /
        # [1 + (gamma - 1)/2 * M²]

        T_new = (

            Tt

            / (

                1.0

                + (
                    gamma - 1.0
                )

                / 2.0

                * Mach**2
            )
        )


        # ========================================================
        # 18. CONVERGENCE SUR LA TEMPERATURE
        # ========================================================
        #
        # VBA :
        #
        # If Abs(1 - t0/t) <= 10^-6 Then Exit Do

        erreur_T = abs(
            1.0
            - T / T_new
        )


        if erreur_T <= tol:

            T = T_new

            break


        T = T_new


    # ============================================================
    # FIN DE LA BOUCLE TEMPERATURE
    # ============================================================

    if (
        iteration_T == max_iter_T
        and erreur_T > tol
    ):

        print(
            "ERREUR : "
            "la boucle sur la température "
            "n'a pas convergé"
        )

        return


    # ============================================================
    # 19. PRESSION TOTALE
    # ============================================================

    if Mach < 0.2:

        # VBA :
        #
        # pt =
        # ptamont
        # + 1/2 rho V²

        p_total = (

            p

            + 0.5
            * rho
            * V**2
        )


    else:

        # Relation isentropique

        p_total = (

            p

            * (

                1.0

                + (
                    gamma - 1.0
                )

                / 2.0

                * Mach**2

            )**(

                gamma
                / (
                    gamma - 1.0
                )
            )
        )


    # ============================================================
    # 20. VERIFICATION DU NOMBRE DE REYNOLDS
    # ============================================================

    avertissement_Re = ""


    if (
        Re_D <= 5000.0
        and 0.1 <= beta <= 0.559
    ):

        avertissement_Re = (
            "Reynolds trop faible"
        )


    if (
        Re_D <= 16000.0 * beta**2
        and beta > 0.559
    ):

        avertissement_Re = (
            "Reynolds trop faible"
        )


    # ============================================================
    # 21. GRANDEURS POUR AFFICHAGE
    # ============================================================

    d_mm = (
        d * 1000.0
    )

    D_mm = (
        D * 1000.0
    )


    # ============================================================
    # 22. AFFICHAGE
    #     MEME ORDRE QUE LA FEUILLE EXCEL
    # ============================================================

    print()
    print("=" * 90)
    print("FEUILLE 2 - CALCUL DU DEBIT")
    print("=" * 90)


    # ------------------------------------------------------------
    # GEOMETRIE
    # ------------------------------------------------------------

    print()
    print("DEFINITION DES DIAMETRES")
    print("-" * 90)

    print(
        f"Diamètre de la veine D0         = "
        f"{D0_mm:.3f} mm"
    )

    print(
        f"Diamètre du diaphragme d0       = "
        f"{d0_mm:.3f} mm"
    )

    print(
        f"Beta                            = "
        f"{beta:.6f}"
    )


    # ------------------------------------------------------------
    # ENTREES
    # ------------------------------------------------------------

    print()
    print("CONDITIONS D'ESSAI")
    print("-" * 90)

    print(
        f"Différence de pression Delta p  = "
        f"{delta_p / 1000.0:.3f} kPa"
    )

    print(
        f"Pression statique amont p       = "
        f"{p_kPa:.3f} kPa"
    )

    print(
        f"Température diaphragme Td       = "
        f"{Td:.3f} K"
    )

    print(
        f"Température totale air Tt       = "
        f"{Tt:.3f} K"
    )


    # ------------------------------------------------------------
    # RESULTATS
    # Même ordre que la feuille Excel
    # ------------------------------------------------------------

    print()
    print("=" * 90)
    print("RESULTATS")
    print("=" * 90)

    print()

    print(
        f"{'Grandeur':<42}"
        f"{'Valeur':>18}"
        f"{'Unité':>16}"
    )

    print("-" * 90)


    # 1. Débit

    print(
        f"{'Débit air':<42}"
        f"{q:>18.6f}"
        f"{'kg/s':>16}"
    )


    # 2. Reynolds

    print(
        f"{'Nombre de Reynolds':<42}"
        f"{Re_D:>18.0f}"
        f"{'-':>16}"
    )


    # 3. Vitesse

    print(
        f"{'Vitesse':<42}"
        f"{V:>18.3f}"
        f"{'m/s':>16}"
    )


    # 4. Mach

    print(
        f"{'Mach':<42}"
        f"{Mach:>18.5f}"
        f"{'-':>16}"
    )


    # 5. Température statique

    print(
        f"{'Température statique':<42}"
        f"{T:>18.3f}"
        f"{'K':>16}"
    )


    # 6. Rho

    print(
        f"{'Rho':<42}"
        f"{rho:>18.5f}"
        f"{'kg/m3':>16}"
    )


    # 7. Viscosité

    print(
        f"{'Viscosité dynamique':<42}"
        f"{mu:>18.8e}"
        f"{'Pa.s':>16}"
    )


    # 8. Pression totale

    print(
        f"{'Pression totale':<42}"
        f"{p_total / 1000.0:>18.3f}"
        f"{'kPa':>16}"
    )


    # ------------------------------------------------------------
    # VALEURS INTERMEDIAIRES
    # ------------------------------------------------------------

    print()
    print("=" * 90)
    print("VALEURS INTERMEDIAIRES")
    print("=" * 90)

    print(
        f"d corrigé                       = "
        f"{d_mm:.6f} mm"
    )

    print(
        f"D corrigé                       = "
        f"{D_mm:.6f} mm"
    )

    print(
        f"Beta                            = "
        f"{beta:.8f}"
    )

    print(
        f"Section diaphragme s            = "
        f"{s:.9e} m²"
    )

    print(
        f"Coefficient approche E          = "
        f"{E:.9f}"
    )

    print(
        f"Coefficient détente epsilon     = "
        f"{epsilon:.9f}"
    )

    print(
        f"Coefficient décharge C          = "
        f"{C:.9f}"
    )

    print(
        f"Cp                              = "
        f"{Cp:.6f} J/(kg.K)"
    )

    print(
        f"Gamma                           = "
        f"{gamma:.9f}"
    )


    # ------------------------------------------------------------
    # CONVERGENCE
    # ------------------------------------------------------------

    print()
    print("CONVERGENCE")
    print("-" * 90)

    print(
        f"Nombre d'itérations sur T       = "
        f"{iteration_T}"
    )

    print(
        f"Nombre d'itérations sur C       = "
        f"{iteration_C}"
    )


    # ------------------------------------------------------------
    # AVERTISSEMENT
    # ------------------------------------------------------------

    if avertissement_Re != "":

        print()

        print(
            "ATTENTION :",
            avertissement_Re
        )


# ================================================================
# MAIN
# ================================================================

def main():

    print("Lancement de prog_1...")
    prog_2()


# myCompiler ne supporte pas :
#
# if __name__ == "__main__":
#     main()
#
# On appelle donc directement main()

main()

Embed on website

To embed this project on your website, copy the following code and paste it into your website's HTML: