saturation curve plotting tool

saturation_plot.py

import CoolProp.CoolProp as CP
import matplotlib.pyplot as plt
import numpy as np

def plot_sat(Fluid, chart, ax=None, subplot=False, numpoints=200):
    if ax == None:
        fig, ax = plt.subplots(1, 1, figsize=(8,6))
    prop1, prop2 = chart[0], chart[1]
    # propriedade crítica
    p_crit = CP.PropsSI('pcrit', Fluid)
    T_crit = CP.PropsSI('Tcrit', Fluid)
    rho_crit = CP.PropsSI('rhocrit', Fluid)
    h_crit = CP.PropsSI('H', 'P', p_crit, 'T', T_crit, Fluid)
    s_crit = CP.PropsSI('S', 'P', p_crit, 'T', T_crit, Fluid)
    
    # valor mínimo
    p_min = CP.PropsSI('ptriple', Fluid) + 100
    T_min = CP.PropsSI('Ttriple', Fluid) + 0.1
    rho_min = CP.PropsSI('D', 'P', p_min, 'T', T_min, Fluid)
    h_min = CP.PropsSI('H', 'P', p_min, 'T', T_min, Fluid)
    s_min = CP.PropsSI('S', 'P', p_min, 'T', T_min, Fluid)
    
    # Dicionário de propriedades
    prop_dict = {'p': ['Pressão (kPa)', 'P', p_crit, p_min],
                 'T': ['Temperatura (K)', 'T', T_crit, T_min],
                 'v': ['Volume específico (m3/kg)', 'D', rho_crit, rho_min],
                 'h': ['Entalpia (kJ/kg)', 'H', h_crit, h_min],
                 's': ['Entropia (kJ/kg.K)', 'S', s_crit, s_min]}
    
    # cálculo
    entry = prop_dict[prop1]
    out = prop_dict[prop2]
    prop1_array = np.linspace(entry[3], entry[2], numpoints)
    prop2_array = np.zeros(2 * numpoints - 1)
    for i, var in enumerate(prop1_array[:-1]):
        prop2_array[i] = CP.PropsSI(out[1], 'Q', 0, entry[1], var, Fluid)
        prop2_array[::-1][i] = CP.PropsSI(out[1], 'Q', 1, entry[1], var, Fluid)
    prop2_array[numpoints - 1] = out[2]
    prop1_array = np.append(prop1_array, prop1_array[::-1][1:])
       
    # Plotando o gráfico
    
    if chart == 'pv':
        ax.loglog(1 / prop2_array, prop1_array / 1000)
        ax.set_xlim(min(1 / prop2_array) / 10, max(1 / prop2_array) * 10)
    elif chart == 'Ts':
        ax.plot(prop2_array / 1000, prop1_array)
    elif chart == 'ph' or chart == 'pT':
        if prop2 == 'h':
            prop2_array = prop2_array / 1000
        ax.semilogy(prop2_array, prop1_array / 1000)
    elif chart == 'Tv':
        ax.semilogx(1 / prop2_array, prop1_array)
        ax.set_xlim(min(1 / prop2_array) / 10, max(1 / prop2_array) * 10)
    ax.set_title(chart + '-diagram')
    ax.set_xlabel(out[0])
    ax.set_ylabel(entry[0])
    if not subplot:
        plt.show()