import math


def api_psv_release_rate(P1, Pback, k, CD, T1, Z, MW, area):
    """
    PSV vapour relief rate calculated according to API 520 Part I 2014
    Eq. 5, 9, 15, 18

    Parameters
    ----------
    P1 : float
        Upstream pressure
    Pback : float 
        Downstream / backpressure
    k : float 
        Ideal gas k (Cp/Cv) 
    CD : float
        Coefficient of discharge
    T1 : float
        Upstream temperature
    Z : float 
        Compressibility
    MW : float 
        Molecular weight of the gas relieved
    area : float 
        PSV orifice area

    Returns
    ----------
        : float 
        Relief rate / mass flow
    """


    P1 = P1 / 1000
    Pback = Pback / 1000
    area = area * 1e6 
    MW = MW * 1000
    C = 0.03948 * (k * (2 / (k + 1)) ** ((k + 1) / (k - 1))) ** 0.5
    if P1 / Pback > ((k + 1) / 2) ** ((k) / (k - 1)):
        w = CD * area * C * P1 / math.sqrt(T1 * Z / MW)
    else:
        r = Pback / P1
        f2 = ((k / (k - 1)) * r ** (2 / k) * (1 - r**((k - 1) / k)) / (1 - r))**0.5
        print(f2)
        w = CD * area * f2 / (T1 * Z / (MW * P1 * (P1 - Pback)))**0.5 / 17.9
    return w/3600