Untitled

# -*- coding: utf-8 -*-
"""
Created on Thu Nov 11 15:23:11 2021

@author: carlo
"""

import thermo
import matplotlib.pyplot as plt
import numpy as np

#Set the properties for our substance
methane = np.array([190.6, 4.6e6, .008])

#Code for problem 1
def problem_1():
    #initial conditions
    P_0, T_0, R = .48e6, 140, 8.31446
    param = thermo.PR_param(P_0, T_0, methane)
    A, B = param[3], param[4]
    z_values=thermo.PR_Z(A,B)
    #obtain z-values, intialize our correct z
    zl, zh, z  = z_values[0] , z_values[2], 0
    #determine the smallest z
    if thermo.PR_log_phi(A,B, zl) < thermo.PR_log_phi(A,B, zh):
        z = zl
    else:
        z = zh
    print("Z of Methane = " + str(z))
    #Calculate the molar volume and enthalpy function
    mVol = (z * R * T_0) / P_0
    h_d = thermo.PR_H(P_0, T_0, param, z)
    print("Molar volume = " + str(mVol))
    print("Enthalpy departure = " + str(h_d))

#Listed below are my functions for problem 2    
#calculates difference between the fugacities    
def df(guess, t):
    #find z values
    params = thermo.PR_param(guess, t, methane)
    A, B = params[3], params[4]
    z_values = thermo.PR_Z(A,B)
    zl, zh  = z_values[0] , z_values[2]
    #calculate the fugacities
    fl = np.exp(thermo.PR_log_phi(A,B, zl)) * guess
    fh = np.exp(thermo.PR_log_phi(A,B, zh)) * guess
    #take the absolute value of the ratio of the fugacities minus one
    df = abs((fl/fh)-1)
    fugacity_ratio = fl/fh 
    #returns this quanitity and the ratio of the fugacities
    return df, fugacity_ratio

#returns the approximate vapor pressure
def vp(guess, t):
    iterations = 1
    #checks if our error is within a small value
    while abs(df(guess,t)[0]) > 0.0001:
        guess *= df(guess,t)[1]
        iterations += 1
    #returns the vapor pressure in MPa and the number of iterations
    return guess / 1e6 , iterations 

#a simple function to print vapor pressure data
def print_vp(guess, t):
    print("Temperature:", t)
    print("Approximated Vapor Pressure:",  vp(guess,t)[0])
    print("Iterations:", vp(guess,t)[1])
    print("")

#prints out the vapor pressures for the 5 given temperatures
def problem_2():
    print_vp(0.05e6,100)
    print_vp(.1e6,120)
    print_vp(.3e6,140)
    print_vp(.5e6,160)
    print_vp(2.5e6,180)
    
#makes two semilog plots: one w a t axis and another with a t/1000 axis
def problem_3():
    #initialize the x and y arrays for each plot
    temps = np.linspace(100,180, 5)
    real_pressures = np.array([0.034931,0.199173,0.686465,1.65732,3.28685])
    approx_pressures = np.array([vp(0.05e6,100)[0],
        vp(.1e6,120)[0],vp(.3e6,140)[0],vp(.5e6,160)[0],vp(2.5e6,180)[0]])
    
    #creates the semilog plot of vapor pressure with respect to t
    plt.figure(figsize=(12,4))
    plt.subplot(1,2,1)
    plt.semilogy(temps,approx_pressures, 'r' , label = "Approximation")
    plt.semilogy(temps,real_pressures, label = "Real")
    plt.title("Real vs Approximated Vapor Pressures")
    plt.xlabel("Temperature (K)")
    plt.ylabel("Vapor Pressure (MPa)")
    plt.legend()
    
    #creates the semilog plot of vapor pressure with respect to t/1000
    plt.subplot(1,2,2)
    plt.semilogy(1000 / temps,approx_pressures, 'r' , label = "Approximation")
    plt.semilogy(1000/ temps,real_pressures, label = "Real")
    plt.title("Real vs Approximated Vapor Pressures")
    plt.xlabel("Temperature (K)")
    plt.ylabel("Vapor Pressure (MPa)")
    plt.legend()

problem_1()
problem_2()
problem_3()