all_in_one_system_2_corrected.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VLE calculation for DME(1)/Benzene(2) binary mixture using Bubble-P algorithm\n",
    "\n",
    "import numpy as np\n",
    "import math\n",
    "import array as arr\n",
    "import matplotlib.pyplot as plt\n",
    "import csv\n",
    "import pandas as pd\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def volumecal (a,b,P,T,eos,aplha_eos):\n",
    "    \n",
    "    Delv = 1\n",
    "    Vold = (8.314*T)/(P) #initial guess\n",
    "    if eos == 1:    \n",
    "        while Delv > 0.001:\n",
    "            Vnew = ((8.314*T)/(P+(a/(Vold*Vold))))+b\n",
    "            Delv = abs(Vnew - Vold)\n",
    "            Vold = Vnew\n",
    "        V = Vold\n",
    "#    print(\"Value of Volume in function is\", V)\n",
    "        return V\n",
    "    elif eos == 2:\n",
    "        while Delv > 0.001:\n",
    "            Vnew = b + ((8.314*T)/(P+(a/(Vold*(Vold+b)*(T**0.5)))))\n",
    "            Delv = abs(Vnew - Vold)\n",
    "            Vold = Vnew\n",
    "        V = Vold\n",
    "#    print(\"Value of Volume in function is\", V)\n",
    "        return V\n",
    "    elif eos == 3:\n",
    "        while Delv > 0.001:\n",
    "            Vnew = b + ((8.314*T)/(P+(a/(Vold*(Vold+b)))))\n",
    "            Delv = abs(Vnew - Vold)\n",
    "            Vold = Vnew\n",
    "        V = Vold\n",
    "#    print(\"Value of Volume in function is\", V)\n",
    "        return V\n",
    "    elif eos == 4:\n",
    "        while Delv > 0.001:\n",
    "            Vnew = b + ((8.314*T)/(P+(a/((Vold*Vold)+(2*b*Vold)-(b*b)))))\n",
    "            Delv = abs(Vnew - Vold)\n",
    "            Vold = Vnew\n",
    "        V = Vold\n",
    "#    print(\"Value of Volume in function is\", V)\n",
    "        return V\n",
    "#V = volumecal(2,3,10,325.15)\n",
    "#print(\"Volume is :\",V)\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mixingrule (a1,b1,a2,b2,Y1,Y2,eos,mix,alpha_eos_1,alpha_eos_2,T,case,A12,A21,V1_L,V2_L,alpha):\n",
    "    if mix == 1:\n",
    "        b12 = 0.5*(b1+b2)\n",
    "        a12 = (a1*a2)**0.5\n",
    "        bmix = (Y1*Y1*b1) + (Y2*Y2*b2) + (2*Y1*Y2*b12)\n",
    "        amix = (Y1*Y1*a1) + (Y2*Y2*a2) + (2*Y1*Y2*a12)\n",
    "    #   print(\"Value of amix in function is : \",amix)\n",
    "    #   print(\"Value of bmix in function is : \", bmix)\n",
    "        return (amix,bmix,a12,b12)\n",
    "    elif mix == 2:\n",
    "        a12 = (a1*alpha_eos_1*a2*alpha_eos_2)**0.5\n",
    "        b12 = 0.5*(b1+b2)\n",
    "        amix = (Y1*Y1*a1*alpha_eos_1) + (Y2*Y2*a2*alpha_eos_2) + (2*Y1*Y2*a12)\n",
    "        bmix = (Y1*b1) + (Y2*b2)\n",
    "        return(amix,bmix,a12,b12)\n",
    "    \n",
    "    elif mix == 3:\n",
    "        a12 = (a1*alpha_eos_1*a2*alpha_eos_2)**0.5\n",
    "        b12 = 0.5*(b1+b2) \n",
    "        A1 = ((Y1 * a1) / b1) + ((Y2 * a2) / b1)\n",
    "        A2 = ((Y1 * Y1 *a1) + (2 * Y1 * Y2 * a12) + (Y2 * Y2 *a2)) / ((Y1 * Y1 * b1) + (2 * Y1 * Y2 * b12) + (Y2 * Y2 *b2))\n",
    "        A_E_Inf = A1 - A2\n",
    "        kij = 0.21 #guess\n",
    "        p11 = b1 - ((a1)/(8.314*T))\n",
    "        p22 = b2 - ((a2)/(8.314*T))\n",
    "        p12 = (p11 + p22) * 0.5 * (1-0.21)\n",
    "        bmix = ((Y1 * Y1 * p11) + (2*Y1*Y2*p12) + (Y2*Y2*p22)) / ((1) + ((A_E_Inf)/(8.314*T)) - ((Y1*a1)/(b1*8.314*T)) - ((Y2*a2)/(b2*8.314*T)))\n",
    "        amix = bmix * (((Y1*a1)/b1) + ((Y2*a2)/b2) - A_E_Inf)\n",
    "        return(amix,bmix,a12,b12)\n",
    "    \n",
    "    elif mix == 4:    \n",
    "        k12 = 0.21\n",
    "        k21 = 0.50\n",
    "        a12 = ((a1*a2)**0.5) * ((1-k12) + ((k12 - k21) * Y1))\n",
    "        a21 = ((a1*a2)**0.5) * ((1-k21) + ((k21 - k12) * Y2))\n",
    "        amix = (Y1*Y1*a1) + (Y2*Y2*a2) + ((Y1*Y2) * (a12 + a21))\n",
    "        b12 = (b1 + b2) * 0.5\n",
    "        bmix = (Y1 * b1) + (Y2 * b2)\n",
    "        return(amix,bmix,a12,b12)\n",
    "    \n",
    "#    elif mix == 4 and eos == 1 and case == 1:\n",
    "#        kij = 0.21 #guess\n",
    "#        a12 = (a1 + a2) *0.5 *(1-kij)\n",
    "#        b12 = (b1 + b2)* 0.5 * (1-kij)\n",
    "#        print(\"a12 is :\",a12)\n",
    "#        print(\"b12 is :\",b12)\n",
    "#        A_E_Inf = (Y1*Y2) *((A21*Y1) + (A12*Y2)) * 8.314 * T\n",
    "#        p11 = b1 - ((a1)/(8.314*T))\n",
    "#        p22 = b2 - ((a2) / (8.314*T))\n",
    "#        p12 = (p11 + p22)*0.5*(1-kij)\n",
    "#        p21 = (p22 + p11)*0.5*(1-kij)\n",
    "#        bmix = ((Y1*Y1*p11) + (Y1*Y2*p12) +(Y1*Y2*p21) + (Y2*Y2*p22)) / ((1) + ((A_E_Inf) / (8.314*T)) - ((Y1*a1) / (b1*8314*T)) - ((Y2*a2) / (b2*8.314*T)))\n",
    "#        amix = bmix * ( ( (Y1*a1) / (b1) ) + ( (Y2*a2) / (b2) ) - A_E_Inf )\n",
    "#        mix = mix +1\n",
    "#        print(\"mix is :\",mix)\n",
    "#        return(amix,bmix,a12,b12)\n",
    "        \n",
    "#    elif mix == 4 and eos == 1 and case == 2: \n",
    "#        kij = 0.21 #guess\n",
    "#        a12 = (a1 + a2) *0.5 *(1-kij)\n",
    "#        b12 = (b1 + b2)* 0.5 * (1-kij)\n",
    "#        print(\"asd\",a12)\n",
    "#        A_E_Inf = (8.314 * T * Y1 * Y2) * ((A12 * A21) / ((A21 * Y1) + (A12 * Y2)))\n",
    "#        p11 = b1 - ((a1)/(8.314*T))\n",
    "#        p22 = b2 - ((a2) / (8.314*T))\n",
    "#        p12 = (p11 + p22)*0.5*(1-kij)\n",
    "#        p21 = (p22 + p11)*0.5*(1-kij)\n",
    "#        bmix = ((Y1*Y1*p11) + (Y1*Y2*p12) +(Y1*Y2*p21) + (Y2*Y2*p22)) / ((1) + ((A_E_Inf) / (8.314*T)) - ((Y1*a1) / (b1*8314*T)) - ((Y2*a2) / (b2*8.314*T)))\n",
    "#        amix = bmix * ( ( (Y1*a1) / (b1) ) + ( (Y2*a2) / (b2) ) - A_E_Inf )\n",
    "#        return(amix,bmix,a12,b12)\n",
    "        \n",
    "#    elif mix == 4 and eos == 1 and case == 3:\n",
    "#        kij = 0.21 #guess\n",
    "#        a12 = (a1 + a2) *0.5 *(1-kij)\n",
    "#        b12 = (b1 + b2)* 0.5 * (1-kij)\n",
    "#        G12 = (V2_L/V1_L)*(np.exp(-A12/(8.314*T)))\n",
    "#        G21 = (V1_L/V2_L)*(np.exp(-A21/(8.314*T)))\n",
    "#        D1 = (np.log(Y1 + (G12 * Y2))) * (-Y1)\n",
    "#        D2 = (np.log(Y2 + (G12 * Y1))) * (-Y2)\n",
    "#        A_E_Inf = (8.314 * T) * (D1 + D2)\n",
    "#        p11 = b1 - ((a1)/(8.314*T))\n",
    "#        p22 = b2 - ((a2) / (8.314*T))\n",
    "#        p12 = (p11 + p22)*0.5*(1-kij)\n",
    "#        p21 = (p22 + p11)*0.5*(1-kij)\n",
    "#        bmix = ((Y1*Y1*p11) + (Y1*Y2*p12) +(Y1*Y2*p21) + (Y2*Y2*p22)) / ((1) + ((A_E_Inf) / (8.314*T)) - ((Y1*a1) / (b1*8314*T)) - ((Y2*a2) / (b2*8.314*T)))\n",
    "#        amix = bmix * ( ( (Y1*a1) / (b1) ) + ( (Y2*a2) / (b2) ) - A_E_Inf )\n",
    "#        return(amix,bmix,a12,b12)\n",
    "        \n",
    "#    elif mix == 4 and eos == 1 and case == 4: \n",
    "#        kij = 0.21 #guess\n",
    "#        a12 = (a1 + a2) *0.5 *(1-kij)\n",
    "#        b12 = (b1 + b2)* 0.5 * (1-kij)\n",
    "#        G12 = np.exp((-1)*alpha*A12)\n",
    "#        G21 = np.exp((-1)*alpha*A21)\n",
    "#        D1 = ((A21 * G21) / (Y1 + (Y2 * G21)))\n",
    "#        D2 = ((A12 * G12) / (Y2 + (Y1 * G12)))\n",
    "#        A_E_Inf = (8.314 * T * Y1 *Y2) * (D1 + D2)\n",
    "#        p11 = b1 - ((a1)/(8.314*T))\n",
    "#        p22 = b2 - ((a2) / (8.314*T))\n",
    "#        p12 = (p11 + p22)*0.5*(1-kij)\n",
    "#        p21 = (p22 + p11)*0.5*(1-kij)\n",
    "#        bmix = ((Y1*Y1*p11) + (Y1*Y2*p12) +(Y1*Y2*p21) + (Y2*Y2*p22)) / ((1) + ((A_E_Inf) / (8.314*T)) - ((Y1*a1) / (b1*8314*T)) - ((Y2*a2) / (b2*8.314*T)))\n",
    "#        amix = bmix * ( ( (Y1*a1) / (b1) ) + ( (Y2*a2) / (b2) ) - A_E_Inf )\n",
    "#        return(amix,bmix,a12,b12)\n",
    "        \n",
    "# a,b = mixingruke (1,1,1,1,1,1)\n",
    "# print (a,b)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialvalue ():\n",
    "\n",
    "# STEP1: Input Temperature\n",
    "    T = float(input(\"Enter the Temperature of system in K:\"))\n",
    "\n",
    "# Antoine constant input for DME(1) in range T:194.3K-248.24K\n",
    "\n",
    "    A1 = 4.11475\n",
    "    B1 = 894.669\n",
    "    C1 = -30.604\n",
    "\n",
    "# Antoine constant input for Benzene (2) Values from NIST Webbook\n",
    "    if T >= 287 and T <= 354.07:\n",
    "        A2 = 4.01814\n",
    "        B2 = 1203.835\n",
    "        C2 = -53.226\n",
    "    elif T > 354.07 and T <= 373.5:\n",
    "        A2 = 4.72583\n",
    "        B2 = 1660.652\n",
    "        C2 = -1.461\n",
    "    elif T >= 421.56 and T <= 554.8:\n",
    "        A2 = 4.60362\n",
    "        B2 = 1701.073\n",
    "        C2 = 20.806\n",
    "    else:\n",
    "        A2 = 4.72583\n",
    "        B2 = 1660.652\n",
    "        C2 = -1.461\n",
    "\n",
    "\n",
    "# Pure component saturated volume at T = T input\n",
    "#v1 = 62.68 # ml/mol from NIST at 248K\n",
    "#v2 = 201.119 # ml/mol at 323.15K\n",
    "    V1_L = 0.00006268 #Si unit\n",
    "    V2_L = 0.000093320 # Si unit\n",
    "\n",
    "#Pure component properties\n",
    "\n",
    "    Pc_1 = 5400000  # Pa from NIST\n",
    "    Vc_1 = 0.000164 # m3 per mole from NIST\n",
    "    Tc_1 = 401      # K from NIST\n",
    "\n",
    "    Pc_2 = 4890000  # Pa from NIST\n",
    "    Vc_2 = 0.00025 # m3 per mole from NIST\n",
    "    Tc_2 = 562    # K from NIST\n",
    "\n",
    "#Input constants of 3 suffix Margulus equation\n",
    "    print(\"\\nSelect the Equation of State to be used from below:\")\n",
    "    print(\"\\nEnter 1 for Van der Waals Equation of State\")\n",
    "    print(\"Enter 2 for Redlich - Kwong Equation of State\")\n",
    "    print(\"Enter 3 for Redlich - Kwong- Soave Equation of State\")\n",
    "    print(\"Enter 4 for Peng - Robinson Equation of State\")\n",
    "    eos = int(input(\"\\nEnter the number corresponding to the Equation of State which is to be used:\"))\n",
    "    \n",
    "    # Finding EOS constants\n",
    "    if eos == 1:\n",
    "        a1 = (27*8.314*8.314*Tc_1*Tc_1)/(64*Pc_1)\n",
    "        b1 = (8.314*Tc_1)/(8*Pc_1)\n",
    "        a2 = (27*8.314*8.314*Tc_2*Tc_2)/(64*Pc_2)\n",
    "        b2 = (8.314*Tc_2)/(8*Pc_2)\n",
    "        omega_eos_1 = 1\n",
    "        alpha_eos_1 = 1\n",
    "        omega_eos_2 = 1\n",
    "        alpha_eos_2 = 1\n",
    "    elif eos == 2:\n",
    "        sigmaa = 1/(9*((2**0.3333)-1))\n",
    "        sigmab = ((2**0.3333)-1) / 3\n",
    "        a1 = (sigmaa*8.314*8.314*(Tc_1**2.5))/(Pc_1)\n",
    "        b1 = (sigmab*8.314*Tc_1)/(Pc_1)\n",
    "        a2 = (sigmaa*8.314*8.314*(Tc_2**2.5))/(Pc_2)\n",
    "        b2 = (sigmab*8.314*Tc_2)/(Pc_2)\n",
    "        omega_eos_1 = 1\n",
    "        alpha_eos_1 = 1\n",
    "        omega_eos_2 = 1\n",
    "        alpha_eos_2 = 1\n",
    "    elif eos == 3:\n",
    "        a1 = (0.42747*8.314*8.314*Tc_1*Tc_1)/Pc_1\n",
    "        b1 = (0.08664*8.314*Tc_1)/Pc_1\n",
    "        a2 = (0.42747*8.314*8.314*Tc_2*Tc_2)/Pc_2\n",
    "        b2 = (0.08664*8.314*Tc_2)/Pc_2\n",
    "        \n",
    "        temperature_1 = 0.7*Tc_1\n",
    "        Pr_sat_1 = ((10**(A1-(B1/(temperature_1+C1))))*100000)/(Pc_1)\n",
    "        omega_eos_1 = -1-math.log10(Pr_sat_1)\n",
    "        m_1 = 0.48508+(1.55171*omega_eos_1)-(0.15613*omega_eos_1*omega_eos_1)\n",
    "        n_1 = (1-((T/Tc_1)**0.5))\n",
    "        alpha_eos_1 = (1+(m_1*n_1))**2\n",
    "        \n",
    "        temperature_2 = 0.7*Tc_2\n",
    "        Pr_sat_2 = ((10**(A2-(B2/(temperature_2+C2))))*100000)/(Pc_2)\n",
    "        omega_eos_2 = -1-math.log10(Pr_sat_2)\n",
    "        m_2 = 0.48508+(1.55171*omega_eos_2)-(0.15613*omega_eos_2*omega_eos_2)\n",
    "        n_2 = (1-((T/Tc_2)**0.5))\n",
    "        alpha_eos_2 = (1+(m_2*n_2))**2\n",
    "    elif eos == 4:\n",
    "        a1 = (0.45724*8.314*8.314*Tc_1*Tc_1)/(Pc_1)\n",
    "        b1 = (0.07780*8.314*Tc_1)/(Pc_1)\n",
    "        a2 = (0.45724*8.314*8.314*Tc_2*Tc_2)/(Pc_2)\n",
    "        b2 = (0.07780*8.314*Tc_2)/(Pc_2)\n",
    "        \n",
    "        temperature_1 = 0.7*Tc_1\n",
    "        Pr_sat_1 = ((10**(A1-(B1/(temperature_1+C1))))*100000)/Pc_1\n",
    "        omega_eos_1 = -1-math.log10(Pr_sat_1)\n",
    "        m_1 = 0.37464+(1.54226*omega_eos_1)-(0.26992*omega_eos_1*omega_eos_1)\n",
    "        n_1 = (1-((T/Tc_1)**0.5))\n",
    "        alpha_eos_1 = (1+(m_1*n_1))**2\n",
    "        \n",
    "        temperature_2 = 0.7*Tc_2\n",
    "        Pr_sat_2 = ((10**(A2-(B2/(temperature_2+C2))))*100000)/Pc_2\n",
    "        omega_eos_2 = -1-math.log10(Pr_sat_2)\n",
    "        m_2 = 0.37464+(1.54226*omega_eos_2)-(0.26992*omega_eos_2*omega_eos_2)\n",
    "        n_2 = (1-((T/Tc_2)**0.5))\n",
    "        alpha_eos_2 = (1+(m_2*n_2))**2\n",
    "    #else: \n",
    "        #print(\"Using dfault Vander waals Equation of State\")\n",
    "        #a1 = (27*8.314*8.314*Tc_1*Tc_1)/(64*Pc_1)\n",
    "        #b1 = (8.314*Tc_1)/(8*Pc_1)\n",
    "        #a2 = (27*8.314*8.314*Tc_2*Tc_2)/(64*Pc_2)\n",
    "        #b2 = (8.314*Tc_2)/(8*Pc_2) \n",
    "        #omega_eos_1 = 0\n",
    "        #alpha_eos_1 = 0\n",
    "        #omega_eos_2 = 0\n",
    "        #alpha_eos_2 = 0\n",
    "    \n",
    "    print(\"\\nSelect mixing rule from below:\")\n",
    "    print(\"\\nEnter 1 for arithmatic, Geometric and Berthelot mixing rule for Van der Waals or Redlich - Kwong Equation of State\")\n",
    "    print(\"Enter 2 for normal mixing rule for Redlich - Kwong- Soave or Peng - Robinson Equation of State  \")\n",
    "    print(\"Enter 3 for mixing rule proposed by Huron and Vidal 1979\")\n",
    "    print(\"Enter 4 for mixing rule proposed by Panagiotopoulos et.al 1986\")\n",
    "    mix = int(input(\"\\nEnter the number corresponding to the mixing rule which is to be used:\"))\n",
    "    \n",
    "    print(\"\\nSelect activity coefficient model from below:\")\n",
    "    print(\"\\nEnter 1 for activity coefficients calculation by 3 Suffix Margulus equation\")\n",
    "    print(\"Enter 2 for activity coefficients calculation by Van Laar equation\")\n",
    "    print(\"Enter 3 for activity coefficients calculation by Wilson equation\")\n",
    "    print(\"Enter 4 for activity coefficients calculation by NRTL equation\")\n",
    "    #print(\"Enter constants of activity coefficient model\")\n",
    "    case = int(input(\"\\nEnter the number corresponding to the activity coefficient calculation model:\"))\n",
    "    if case == 1:\n",
    "        A12 = 0.0183\n",
    "        A21 = -0.2788\n",
    "        alpha = 3.53284\n",
    "        print(\"case is :\",case)\n",
    "    elif case == 2:\n",
    "        A12 = -0.0665 \n",
    "        A21 = -0.8466\n",
    "        alpha = 3.53284\n",
    "        print(\"case is :\",case)\n",
    "    elif case == 3:\n",
    "        A12 = -246.1952\n",
    "        A21 = 245.9898\n",
    "        alpha = 3.53284\n",
    "        print(\"case is :\",case)\n",
    "    elif case == 4:\n",
    "        A12 = 46.6072\n",
    "        A21 = 46.2807\n",
    "        alpha = 3.53284\n",
    "        print(\"case is :\",case)\n",
    "    #else:\n",
    "        #print(\"Using default of 3- Suffix Margulus equation\")\n",
    "        #A12 = -0.0487\n",
    "        #A21 = -0.0492\n",
    "        #alpha = 0.3048\n",
    "        #print(\"case is :\",case)\n",
    "        \n",
    "\n",
    "#Input x values\n",
    "    #X1 = []\n",
    "    # number of elements as input\n",
    "    #n = int(input(\"Enter number of elements : \"))\n",
    "    # iterating till the range\n",
    "    #for i in range(0, n):\n",
    "        #x1 = float(input())\n",
    "        #X1.append(x1)\n",
    "    #X1 = arr.array('f', [0,0.0361,0.0923,0.1674,0.2328,0.2655,0.2900,0.3215,0.4133,0.4820,0.5219,0.5573,0.5831,0.6065,0.6285,0.6535,0.6787,0.6970,0.7043,0.7296,0.7614,0.7919,0.8206,0.8476,0.8725,0.8952,0.9331,0.9476,0.9601,0.9702,0.97830.9848,0.9892,0.9601,0.9923,0.9944,0.9966,0.9987,1])\n",
    "    X1 =[0,0.0528,0.0907,0.1182,0.1551,0.1962,0.2180,0.2531,0.3157,0.3735,0.4248,0.4708,0.5122,0.5152,0.55,0.5556,0.5975,\n",
    "         0.6408,0.6842,0.7268,0.7680,0.8067,0.8422,0.8745,0.9006,0.9232,0.9423,0.9578,0.9701,0.9784,0.9842,0.9894,0.9929,\n",
    "         0.9964,0.9981,1]\n",
    "    print(\"Lenght of vector X1 is :\",len(X1))\n",
    "    X2 = arr.array('f',[])\n",
    "   # Finding X2\n",
    "\n",
    "    for i in range(len(X1)):\n",
    "        x2 = 1 - X1[i]\n",
    "        X2.append(x2)\n",
    "\n",
    "#Fugacity coefficients equal to 1\n",
    "\n",
    "    Phi_1 = 1\n",
    "    Phi_2 = 1\n",
    "\n",
    "#Poynting correction factor equal to 1\n",
    "\n",
    "    Pcf_1 = 1\n",
    "    Pcf_2 = 1\n",
    "\n",
    "# STEP 2\n",
    "\n",
    "    P1sat = 100000*(10**(A1-(B1/(T+C1)))) #Pressure in si unit\n",
    "    P2sat = 100000*(10**(A2-(B2/(T+C2)))) #Pressure in si unit\n",
    "    print(\"P1sat is:\",P1sat)\n",
    "    print(\"P2sat is:\",P2sat)\n",
    "\n",
    "\n",
    "    #print(a1,b1,a2,b2,X2)\n",
    "   # print(\"X1 isequal to :\",X1)\n",
    "   # print(\"X2 isequal to :\",X2)\n",
    "    return (T,P1sat,P2sat,X1,X2,a1,b1,a2,b2,A12,A21,V1_L,V2_L,Phi_1,Phi_2,Pcf_1,Pcf_2,case,alpha,eos,omega_eos_1,alpha_eos_1,omega_eos_2,alpha_eos_2,mix)\n",
    "#T,P1sat,P2sat,X1,X2,a1,b1,a2,b2,A12,A21,V1_L,V2_L,Phi_1,Phi_2,Pcf_1,Pcf_2,n,case,alpha,eos = initialvalue()\n",
    "#print(a1,b1,a2,b2,X1,X2,P1sat,P2sat,T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Zcal (Psat,Vsat,T):\n",
    "    Z = (Psat*Vsat)/(8.314*T)\n",
    "    #print(\"Value of Z in funtion is :\",Z)\n",
    "    return Z\n",
    "#Z = Zcal(10,10,320)\n",
    "#print(Z)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phisatcal (Z1,Z2,a1,b1,a2,b2,V1sat,V2sat,T,eos,alpha_eos_1,alpha_eos_2):\n",
    "    if eos == 1:\n",
    "        Ln_phi_1_sat = (-1)*(math.log(1-(b1/V1sat))+(a1/(V1sat*8.314*T))+(math.log(Z1))+(1-Z1))\n",
    "       # print(\"Ln_phi_1_sat\",Ln_phi_1_sat)\n",
    "        Ln_phi_2_sat = (-1)*(math.log(1-(b2/V1sat))+(a2/(V2sat*8.314*T))+(math.log(Z2))+(1-Z2))\n",
    "        #print(\"Ln_phi_2_sat\",Ln_phi_2_sat)\n",
    "        Phi_1_sat = math.exp(Ln_phi_1_sat)\n",
    "        #print(\"Phi_1_sat\",Phi_1_sat)\n",
    "        Phi_2_sat = math.exp(Ln_phi_2_sat)\n",
    "        #print(\"Phi_2_sat\",Phi_2_sat)\n",
    "        print(\"Value of Phi_1_sat is :\",Phi_1_sat)\n",
    "        print(\"Value of Phi_2_sat is :\",Phi_2_sat)\n",
    "        return (Phi_1_sat,Phi_2_sat)\n",
    "    elif eos == 2:\n",
    "        Ln_phi_1_sat = (-1)*(1+(a1/(b1*8.314*(T**1.5))))*(math.log(1-(b1/V1sat))) - math.log(Z1) + (Z1-1)\n",
    "        Ln_phi_2_sat = (-1)*(1+(a2/(b2*8.314*(T**1.5))))*(math.log(1-(b2/V2sat))) - math.log(Z2) + (Z2-1)\n",
    "        Phi_1_sat = math.exp(Ln_phi_1_sat)\n",
    "        Phi_2_sat = math.exp(Ln_phi_2_sat)\n",
    "        print(\"Value of Phi_1_sat is :\",Phi_1_sat)\n",
    "        print(\"Value of Phi_2_sat is :\",Phi_2_sat)\n",
    "        return (Phi_1_sat,Phi_2_sat)\n",
    "    elif eos == 3:\n",
    "        Ln_phi_1_sat = (-1)*(1+((a1*alpha_eos_1)/(b1*8.314*T)))*(math.log(1-(b1/V1sat)))-math.log(Z1)+(Z1-1)\n",
    "        Ln_phi_2_sat = (-1)*(1+((a2*alpha_eos_2)/(b2*8.314*T)))*(math.log(1-(b2/V2sat)))-math.log(Z2)+(Z2-1)\n",
    "        Phi_1_sat = math.exp(Ln_phi_1_sat)\n",
    "        Phi_2_sat = math.exp(Ln_phi_2_sat)\n",
    "        print(\"Value of Phi_1_sat is :\",Phi_1_sat)\n",
    "        print(\"Value of Phi_2_sat is :\",Phi_2_sat)\n",
    "        return (Phi_1_sat,Phi_2_sat)\n",
    "    elif eos ==4:\n",
    "        Ln_phi_1_sat = (-1)*(math.log(1-(b1/V1sat))) - ((a1*alpha_eos_1)/(2*1.414*8.314*b1*T))*math.log(abs(((2.414*b1)+V1sat)/((0.414*b1)-V1sat))) - math.log(Z1) + (Z1-1)\n",
    "        Ln_phi_2_sat = (-1)*(math.log(1-(b2/V2sat))) - ((a2*alpha_eos_2)/(2*1.414*8.314*b2*T))*math.log(abs(((2.414*b2)+V2sat)/((0.414*b2)-V2sat))) - math.log(Z2) + (Z2-1)\n",
    "        Phi_1_sat = math.exp(Ln_phi_1_sat)\n",
    "        Phi_2_sat = math.exp(Ln_phi_2_sat)\n",
    "        print(\"Value of Phi_1_sat is :\",Phi_1_sat)\n",
    "        print(\"Value of Phi_2_sat is :\",Phi_2_sat)\n",
    "        return(Phi_1_sat,Phi_2_sat)\n",
    "    \n",
    "                                                                                              \n",
    "                                                                    \n",
    "#phi1sat,phi2sat = phisatcal(1,1,0.5,0.4,0.2,0.54,14,15,320)\n",
    "#print(phi1sat,phi2sat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gammacal (x1,x2,A12,A21,case,alpha,V1_L,V2_L):\n",
    "    \n",
    "    if case == 1:\n",
    "        gamma1 = math.exp(((x2**2)*(2*(A21-A12)*x1)) + A12)\n",
    "        gamma2 = math.exp(((x1**2)*(2*(A12-A21)*x2)) + A21)\n",
    "        #print(\"Value of gamma1 in function is :\",gamma1)\n",
    "        #print(\"Value of gamma2 in function is :\",gamma2)\n",
    "        #print(\"Valus of case inside the function is :\",case)\n",
    "        return (gamma1,gamma2)\n",
    "    elif case == 2:\n",
    "        A = ((A21*x2)/((A12*x1)+(A21*x2)))**2\n",
    "        B = ((A12*x1)/((A12*x1)+(A21*x2)))**2\n",
    "        gamma1 = math.exp(A12*A)\n",
    "        gamma2 = math.exp(A21*B)\n",
    "        #print(\"Valus of case inside the function is :\",case)\n",
    "        return (gamma1,gamma2)\n",
    "    elif case == 3:\n",
    "        G12 = (V2_L/V1_L)*(np.exp(-A12/(8.314*T)))\n",
    "        G21 = (V1_L/V2_L)*(np.exp(-A21/(8.314*T)))\n",
    "        D1 = (-1)*(np.log(x1+(G12*x2)))\n",
    "        D2 = (G12)/(x1+(G12*x2))\n",
    "        D3 = (G21)/(x2+(G21*x1))\n",
    "        gamma1 = np.exp(D1+(x2*(D2-D3)))\n",
    "        #print(\"gamma1 is clear with value of :\",gamma1)\n",
    "        E1 = (-1)*(np.log(x2+(G21*x1)))\n",
    "        E2 = (G12)/(x1+(G12*x2))\n",
    "        E3 = (G21)/(x2+(G21*x1))\n",
    "        gamma2 = np.exp(E1+(x1*(E2-E3)))\n",
    "        #gamma1 = math.exp((-1)*(math.log(x1+(x2*A12))) + (x2)*(((A12)/(x1+(x2*A12))) - ((A21)/(x2+(x1*A21)))))\n",
    "        #gama2 =  math.exp((-1)*(math.log(x2+(x1*A21))) + (x1)*(((A12)/(x1+(x2*A12))) - ((A21)/(x2+(x1*A21)))))\n",
    "        #print(\"Valus of case inside the function is :\",case)\n",
    "        return (gamma1,gamma2)\n",
    "    elif case == 4:\n",
    "        G12 = np.exp((-1)*alpha*A12)\n",
    "        G21 = np.exp((-1)*alpha*A21)\n",
    "        A1 = ((G21)/(x1+(x2*G21)))**2\n",
    "        A2 = (((G12)/(x2+(x1*G12)))**2)/(G12)\n",
    "        gamma1 = np.exp(x2*x2*((A21*A1) + (A12*A2)))\n",
    "        B1 = ((G12)/(x2+(x1*G12)))**2\n",
    "        B2 = (((G21)/(x1+(x2*G21)))**2)/(G21)\n",
    "        gamma2 = np.exp(x1*x1*((A12*B1) + (A21*B2)))\n",
    "        #print(\"Valus of case inside the function is :\",case)\n",
    "        return(gamma1,gamma2)\n",
    "    #else:\n",
    "        #gamma1 = math.exp(((x2**2)*(2*(A21-A12)*x1)) + A12)\n",
    "        #gamma2 = math.exp(((x1**2)*(2*(A12-A21)*x2)) + A21)\n",
    "        #print(\"Valus of case inside the function is :\",case)\n",
    "        #print(\"Value of gamma1 in function is :\",gamma1)\n",
    "        #print(\"Value of gamma2 in function is :\",gamma2)\n",
    "        return (gamma1,gamma2)\n",
    "\n",
    "#gamma1,gamma2 = gammacal(0.5,0.5,-0.0487,-0.0492)\n",
    "#print(gamma1,gamma2)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pcal (x1,x2,gamma1,gamma2,P1sat,P2sat,Phi_1_sat,Phi_2_sat,Phi_1,Phi_2,Pcf_1,Pcf_2):\n",
    "    P = ((x1*gamma1*P1sat*Phi_1_sat*Pcf_1)/(Phi_1)) + ((x2*gamma2*P2sat*Phi_2_sat*Pcf_2)/(Phi_2))\n",
    "#    print(\"Value of Pcal in function is :\",P)\n",
    "    return P\n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ycal (x1,x2,gamma1,gamma2,P1sat,P2sat,Phi_1_sat,Phi_2_sat,Pcf_1,Pcf_2,Phi_1,Phi_2,P):\n",
    "    import numpy as np\n",
    "    import math\n",
    "    import array as arr\n",
    "    y1 = (x1*gamma1*P1sat*Phi_1_sat*Pcf_1)/(Phi_1*P)\n",
    "    y2 = (x2*gamma2*P2sat*Phi_2_sat*Pcf_2)/(Phi_2*P)     \n",
    "#    print(\"Value of Y1 in function is :\",y1)\n",
    "#    print(\"Value of Y2 in function is :\",y2)\n",
    "    return (y1,y2)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PCF (P,P1sat,P2sat,V1_L,V2_L,T):\n",
    "    import numpy as np\n",
    "    import math\n",
    "    import array as arr\n",
    "    Pcf_1 = math.exp(((P-P1sat)*V1_L)/(8.314*T))\n",
    "    Pcf_2 = math.exp(((P-P2sat)*V2_L)/(8.314*T))\n",
    "#    print(\"Value of Pcf_1 in function is :\",Pcf_1)\n",
    "#    print(\"Value of Pcf_2 in function is :\",Pcf_2)\n",
    "    return (Pcf_1,Pcf_2)\n",
    "#Pcf_1,Pcf_2 = PCF(10,12,14,0.25,0.254,325)\n",
    "#print (Pcf_1,Pcf_2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Phical (Vmix,bmix,a12,y1,y2,T,Pold,eos,a1,b1,a2,b2,a11,a22,amix):\n",
    "  \n",
    "    if eos == 1:\n",
    "        ln_phi_1 = math.log(Vmix/(Vmix-bmix)) + (b1/(Vmix-bmix)) - math.log((Pold*Vmix)/(8.314*T)) - ((2*(y1*a1+y2*a12))/(Vmix*8.314*T))\n",
    "        ln_phi_2 = math.log(Vmix/(Vmix-bmix)) + (b2/(Vmix-bmix)) - math.log((Pold*Vmix)/(8.314*T)) - ((2*(y1*a12+y2*a2))/(Vmix*8.314*T))\n",
    "        phi_1 = math.exp(ln_phi_1)\n",
    "        phi_2 = math.exp(ln_phi_2)\n",
    "    #    print(\"Value of phi_1 in the function is :\",phi_1)\n",
    "    #    print(\"Value of phi_2 in the function is :\",phi_2)\n",
    "        return (phi_1,phi_2)\n",
    "    \n",
    "    elif eos == 2:\n",
    "        zmix = (Pold*Vmix)/(8.314*T)\n",
    "        ln_phi_1 = (b1/bmix)*(zmix-1) - math.log(zmix*(1-(bmix/Vmix))) + ((bmix*8.314*(T**1.5))**(-1))*(math.log(1+(bmix/Vmix)))*((amix*b1)/(bmix)-2*((y1*a1)+(y2*a12)))\n",
    "        ln_phi_2 = (b2/bmix)*(zmix-1) - math.log(zmix*(1-(bmix/Vmix))) + ((bmix*8.314*(T**1.5))**(-1))*(math.log(1+(bmix/Vmix)))*((amix*b2)/(bmix)-2*((y1*a12)+(y2*a2)))\n",
    "        phi_1 = math.exp(ln_phi_1)\n",
    "        phi_2 = math.exp(ln_phi_2)\n",
    "        return(phi_1,phi_2)\n",
    "    \n",
    "    elif eos == 3:\n",
    "        zmix = (Pold*Vmix)/(8.314*T)\n",
    "        ln_phi_1 = (b1/bmix)*(zmix-1) - math.log(zmix*(1-(bmix/Vmix))) + ((bmix*8.314*(T**1.5))**(-1))*(math.log(1+(bmix/Vmix)))*((amix*b1)/(bmix)-2*((y1*a11)+(y2*a12)))\n",
    "        ln_phi_2 = (b2/bmix)*(zmix-1) - math.log(zmix*(1-(bmix/Vmix))) + ((bmix*8.314*(T**1.5))**(-1))*(math.log(1+(bmix/Vmix)))*((amix*b2)/(bmix)-2*((y1*a12)+(y2*a22)))\n",
    "        phi_1 = math.exp(ln_phi_1)\n",
    "        phi_2 = math.exp(ln_phi_2)\n",
    "        return(phi_1,phi_2)\n",
    "    \n",
    "    elif eos == 4:\n",
    "        zmix = (Pold*Vmix)/(8.314*T)\n",
    "#        ln_phi_1 = 24\n",
    "#        ln_phi_2 = 14\n",
    "        ln_phi_1 = (b1/bmix)*(zmix-1) - math.log(zmix*(1-(bmix/Vmix))) + ((amix)/(2*1.414*bmix*8.314*T))*(math.log((Vmix+(2.414*bmix))/(Vmix-(0.414*bmix))))*((b1/bmix)-2*((y1*a11)+(y2*a12)))\n",
    "        ln_phi_2 = (b2/bmix)*(zmix-1) - math.log(zmix*(1-(bmix/Vmix))) + ((amix)/(2*1.414*bmix*8.314*T))*(math.log((Vmix+(2.414*bmix))/(Vmix-(0.414*bmix))))*((b2/bmix)-2*((y1*a12)+(y2*a22)))\n",
    "        phi_1 = math.exp(ln_phi_1)\n",
    "        phi_2 = math.exp(ln_phi_2)\n",
    "        return(phi_1,phi_2)\n",
    "    \n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter the Temperature of system in K:324\n",
      "\n",
      "Select the Equation of State to be used from below:\n",
      "\n",
      "Enter 1 for Van der Waals Equation of State\n",
      "Enter 2 for Redlich - Kwong Equation of State\n",
      "Enter 3 for Redlich - Kwong- Soave Equation of State\n",
      "Enter 4 for Peng - Robinson Equation of State\n",
      "\n",
      "Enter the number corresponding to the Equation of State which is to be used:1\n",
      "\n",
      "Select mixing rule from below:\n",
      "\n",
      "Enter 1 for arithmatic, Geometric and Berthelot mixing rule for Van der Waals or Redlich - Kwong Equation of State\n",
      "Enter 2 for normal mixing rule for Redlich - Kwong- Soave or Peng - Robinson Equation of State  \n",
      "Enter 3 for mixing rule proposed by Huron and Vidal 1979\n",
      "Enter 4 for mixing rule proposed by Panagiotopoulos et.al 1986\n",
      "\n",
      "Enter the number corresponding to the mixing rule which is to be used:1\n",
      "\n",
      "Select activity coefficient model from below:\n",
      "\n",
      "Enter 1 for activity coefficients calculation by 3 Suffix Margulus equation\n",
      "Enter 2 for activity coefficients calculation by Van Laar equation\n",
      "Enter 3 for activity coefficients calculation by Wilson equation\n",
      "Enter 4 for activity coefficients calculation by NRTL equation\n",
      "\n",
      "Enter the number corresponding to the activity coefficient calculation model:1\n",
      "case is : 1\n",
      "Lenght of vector X1 is : 36\n",
      "P1sat is: 1162501.438293401\n",
      "P2sat is: 37345.37067541462\n",
      "Value of Phi_1_sat is : 0.8947171856559626\n",
      "Value of Phi_2_sat is : 1.0496879833835415\n",
      "\n",
      "Final values of P in KPa is : array('f', [29.85002899169922, 82.10055541992188, 118.5301513671875, 144.62403869628906, 179.4171905517578, 218.1605987548828, 238.8129119873047, 272.3533935546875, 333.52410888671875, 392.1610412597656, 446.4122314453125, 497.11041259765625, 544.5655517578125, 548.0748901367188, 589.5025024414062, 596.29443359375, 648.237548828125, 704.0198974609375, 762.0611572265625, 821.0311279296875, 879.811767578125, 936.4125366210938, 989.315673828125, 1038.0694580078125, 1077.7474365234375, 1112.199951171875, 1141.3109130859375, 1164.8839111328125, 1183.5308837890625, 1196.0733642578125, 1204.814697265625, 1212.63330078125, 1217.8851318359375, 1223.1278076171875, 1225.6707763671875, 1228.5101318359375])\n",
      "\n",
      "Corresponding values of X1 are : [0, 0.0528, 0.0907, 0.1182, 0.1551, 0.1962, 0.218, 0.2531, 0.3157, 0.3735, 0.4248, 0.4708, 0.5122, 0.5152, 0.55, 0.5556, 0.5975, 0.6408, 0.6842, 0.7268, 0.768, 0.8067, 0.8422, 0.8745, 0.9006, 0.9232, 0.9423, 0.9578, 0.9701, 0.9784, 0.9842, 0.9894, 0.9929, 0.9964, 0.9981, 1]\n",
      "\n",
      "Final values of Y1 are : array('f', [0.0, 0.6507737040519714, 0.7651349306106567, 0.8116505742073059, 0.8526362180709839, 0.882902979850769, 0.8950293660163879, 0.9108129143714905, 0.9314519762992859, 0.945228099822998, 0.9547836780548096, 0.9618562459945679, 0.9673063158988953, 0.9676727056503296, 0.9716772437095642, 0.9722821712493896, 0.9765012860298157, 0.9803602695465088, 0.9837961196899414, 0.986807644367218, 0.9894235730171204, 0.9916452169418335, 0.9935021996498108, 0.9950541257858276, 0.9962189197540283, 0.9971668124198914, 0.9979259371757507, 0.9985147714614868, 0.9989652037620544, 0.9992608428001404, 0.9994636178016663, 0.9996427297592163, 0.9997618198394775, 0.9998798370361328, 0.9999366998672485, 1.0])\n",
      "\n",
      "Final values of Y2 are : array('f', [1.0, 0.34922632575035095, 0.23486509919166565, 0.1883493959903717, 0.14736376702785492, 0.11709700524806976, 0.10497063398361206, 0.08918709307909012, 0.06854801625013351, 0.054771896451711655, 0.04521632194519043, 0.03814375400543213, 0.032693710178136826, 0.03232727199792862, 0.02832273580133915, 0.027717838063836098, 0.023498687893152237, 0.019639713689684868, 0.016203902661800385, 0.01319234911352396, 0.010576397180557251, 0.008354774676263332, 0.006497781723737717, 0.004945902619510889, 0.003781079314649105, 0.0028331896755844355, 0.0020740581676363945, 0.0014852298190817237, 0.0010348194045946002, 0.0007391284452751279, 0.000536383013240993, 0.00035729631781578064, 0.00023817465989850461, 0.00012018559937132522, 6.328328890958801e-05, 0.0])\n",
      "omega1 1\n",
      "omega2 1\n"
     ]
    }
   ],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    T,P1sat,P2sat,X1,X2,a1,b1,a2,b2,A12,A21,V1_L,V2_L,Phi_1,Phi_2,Pcf_1,Pcf_2,case,alpha,eos,omega_eos_1,alpha_eos_1,omega_eos_2,alpha_eos_2,mix = initialvalue() #reading and evaluating all the constants\n",
    "    a11 = a1*alpha_eos_1\n",
    "    a22 = a2*alpha_eos_2\n",
    "    \n",
    "    V1sat = volumecal(a11,b1,P1sat,T,eos,alpha_eos_1)\n",
    "    V2sat = volumecal(a22,b2,P2sat,T,eos,alpha_eos_2)\n",
    "        \n",
    "    Z1sat = Zcal(P1sat,V1sat,T)\n",
    "    Z2sat = Zcal(P2sat,V2sat,T)\n",
    "#phi1sat = 0.89\n",
    "#phi2sat = 0.94\n",
    "    phi1sat,phi2sat = phisatcal(Z1sat,Z2sat,a1,b1,a2,b2,V1sat,V2sat,T,eos,alpha_eos_1,alpha_eos_2)\n",
    "#gamma1,gamma2 = gammacal(X1[0],X2[0],A12,A21)\n",
    "    gamma1 = arr.array('f',[])\n",
    "    gamma2 = arr.array('f',[])\n",
    "    P = arr.array('f',[])\n",
    "    Y1 = arr.array('f',[])\n",
    "    Y2 = arr.array('f',[])\n",
    "    for j in range(0, len(X1)):\n",
    "        deltaP = 1\n",
    "        gam1,gam2 = gammacal(X1[j],X2[j],A12,A21,case,alpha,V1_L,V2_L)\n",
    "        gamma1.append(gam1)\n",
    "        gamma2.append(gam2)\n",
    "        Pold = Pcal(X1[j],X2[j],gam1,gam2,P1sat,P2sat,phi1sat,phi2sat,Phi_1,Phi_2,Pcf_1,Pcf_2)\n",
    "        y1,y2 = Ycal(X1[j],X2[j],gam1,gam2,P1sat,P2sat,phi1sat,phi2sat,Pcf_1,Pcf_2,Phi_1,Phi_2,Pold)\n",
    "        #print(\"y1 :\",y1)\n",
    "        #print(\"y2 :\",y2)\n",
    "        while deltaP>0.001:\n",
    "            amix,bmix,a12,b12 = mixingrule (a1,b1,a2,b2,y1,y2,eos,mix,alpha_eos_1,alpha_eos_2,T,case,A12,A21,V1_L,V2_L,alpha)\n",
    "#        print(\"amix :\",amix)\n",
    "            Vmix = volumecal(amix,bmix,Pold,T,eos,alpha_eos_1)\n",
    "#        print(\"Vmix :\",Vmix)\n",
    "        #a12 = (a1*a2)**0.5\n",
    "#       print(\"a12 :\",a12)\n",
    "            Phi_1,Phi_2 = Phical (Vmix,bmix,a12,y1,y2,T,Pold,eos,a1,b1,a2,b2,a11,a22,amix)\n",
    "#        print(\"Value of phi_1  :\",Phi_1)\n",
    "#        print(\"Value of phi_2  :\",Phi_2)\n",
    "            Pcf_1,Pcf_2 = PCF(Pold,P1sat,P2sat,V1_L,V2_L,T)\n",
    "#        print(\"Value of pcf_1 :\",Pcf_1)\n",
    "#        print(\"Value of pcf_2is :\",Pcf_2)\n",
    "            Pnew = Pcal(X1[j],X2[j],gam1,gam2,P1sat,P2sat,phi1sat,phi2sat,Phi_1,Phi_2,Pcf_1,Pcf_2)\n",
    "            y1,y2 = Ycal(X1[j],X2[j],gam1,gam2,P1sat,P2sat,phi1sat,phi2sat,Pcf_1,Pcf_2,Phi_1,Phi_2,Pnew)\n",
    "            deltaP = abs(Pnew - Pold)\n",
    "            Pold = Pnew #P in Pa\n",
    "        p = Pnew / 1000\n",
    "    \n",
    "        P.append(p)\n",
    "        Y1.append(y1)\n",
    "        Y2.append(y2)\n",
    "    \n",
    "    print(\"\\nFinal values of P in KPa is :\",P)\n",
    "    print(\"\\nCorresponding values of X1 are :\",X1)\n",
    "    print(\"\\nFinal values of Y1 are :\",Y1)\n",
    "    print(\"\\nFinal values of Y2 are :\",Y2)\n",
    "    print(\"omega1\",omega_eos_1)\n",
    "    print(\"omega2\",omega_eos_2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAudklEQVR4nO3deZwU1bn/8c8jO4KCAioMiwpRQTYFxX2LBFDRqySKo0YJIkaMMXHBS1TUa64RY/y5BXGPYIi7KAiKRqNxY4yIopALyjKikUVgANmf3x+nZujp6Z7pGbqnZ6a/79erXjNVdarqqe6ZeuqcqlNl7o6IiOSuXbIdgIiIZJcSgYhIjlMiEBHJcUoEIiI5TolARCTHKRGIiOQ4JYI6xMwuNLN3qmlbB5jZx2ZWZGa/qo5tVhDPJWZ2V4plf2Vmt2U4pGTbfszM/icb206Vmf2vmf06xbJ3mtnIDIckGaZEkAFm1szMFpnZuTHTmpvZEjMbEjNtLzNbYWbHxy3/qJn9NcMxvmlmw3diFdcAb7p7c3e/O8n6N0aJYq2ZfWRmo82sUUyZsWbm8YnEzH4dTR8bjR9vZtvNbF3ccEQ0vyHwO2BczDommNn8aLkL48KbAJxnZm12Yv/rJDNrDVwAPBCN9zOz18xslZktN7OnzWyfmEXGAWOi76Aq26swMZrZ36NtrzWzT8zs9Jh5p5jZO2a22sy+NbMHzax5gnXsEa2j3BMlMzvXzBab2Xoze8HM9oiZV+p/Jvq7/N7MzqncXtc8SgQZ4O7rgBHA/4v+sQBuBwrc/ZmYcv8BrgQeNLMmAGZ2EnAKkPWz7Ap0BOZWUGaUuzcH9gF+C5wDTDMziynzb+DncctdEE2Ptczdm8UN70XzTgfmufvXMeU/AX4J/Cs+KHffCLwSbaeMKEGNrWDf6qoLgWnu/kM03pKQODsRvvMi4NHiwu7+DTAPGJzBmK4A9nH33Qj/VxNjktHuwP8AbYGDgDxiTghi/AH4oryNmFk3QgI8H9gL2ADcn6Rsf+AFYJi7T67k/tQ4SgQZ4u6vAlOBu6Mz/p8BlyUo9wQwH7g5SgYPAL9y9+WJ1lt8Bm1mX0a1iXFmlvB7NLMjzWyWma2Jfh4ZTb8VOAa4NzqzvjfJ8oPNbG50tvWmmR0UTX8DOCFm+R9V8Fmsd/c3CQeLIwiJrtgsoGn0T1j8z9gkmp6qgcBbcdu8z91fBzYmWebNuDhSEtX0rjKzOdHn+jcza1xO+d5m9q+oZvQ3oHHc/FPNbHb0Gb9rZj1i5rU3s+eiM9mVxd+Tme1vZm9E01aY2SQza5FqjOVtk7jP0t1fcfen3X2tu28A7gWOitvNNynns4xqEd9Gsfwj5rseAeQD10R/Ry8lWt7d57j71uJRoAHQPpr3pLtPd/cN7v498GB8fBZqjgcTk8CSyAdecvd/RCdz1wNnxtcwzOxU4CngXHd/voJ11gpKBJl1JXA88AxwVXT2lMhIYBgwGfgshTOM/wL6AIcQzoaHxReIqrRTgbuBPYE7galmtqe7jwHeJpyxN3P3UQmW/xHwV+DXQGtgGvCSmTV09xPjlo8/e0/I3ZcABYQkFOsJdpyd/xz4Syrri9GdkEwr4wugZyWXKfYzYACwL9CDcBZdRtRc8gJh//YAngbOipl/CPAIcAnhO3oAmGJmjcysHvAysJhwNt6O8PcBYMD/suMsuD0wNpUYy9tmtFxFn+WxlK0JVvRZvgJ0AdoQamiTANx9QvT77dHf0WnJVmBmL5vZRuADQuIpSCW+6HO8DxhFSCLl6UaoSRLFtxDYDMSe6JwGTASGuPu0CtZXaygRZFB0hjIXaAo8V065QuAG4MfApSms+g/uvio6sN4FDE1Q5hTg/9z9CXff6u5/JVThk/6zxTkbmOrur7n7FuAOwpn6kSkun8wywkEx1kRgqJk1IDQfTUywXNvoDDZ22DWa14LQZFEZRYRmhaq4292Xufsq4CWgV5Jy/Qhnr3e5+5aoWTC2pnMx8IC7f+Du29z9cWBTtNxhhAP91VGNaqO7vwPg7gui72VTVHO8EzguxRjL2yaU81lGNYcbgKvjZhVFyyXk7o+4e5G7byIkrJ5mVqnP3t1PBZoDg4AZ7r49QXwnE04kboiZ/CvgA3f/KIXNNAPWxE1bE2232AmEZst/ph59zadEkEFmdh7hbG4moY2yePortuOCZ340eS7wfTm1hlhLY35fTDhgxGsbzSOubLsUwy+1fPSPt7QSyyfTDlgVOyFKaAuA3xOS19IEyy1z9xZxw/po3veU/mdNRXNi/umjM87VZrYaGA2Mjkk4L8ct+23M7xsIB5BE32tb4Gsv/WTH2O+kI/Db2ORGOLtvG/1cHNMkUsLM2pjZZDP72szWEhJnq1RirGCbkOSzNLPOhDP7K9z97bjZzYHV8ctEy9Uzs9vMbGEU66JoVny8FYqS6SvAT8ys1DUJM+sHPEk4U/93NK0tIRGMSXET64Dd4qbtRunEeD0hcb4QU4uq9epnO4C6ysIdKX8iVNHnAXPN7Mmo/XHgTq6+PTuqvx0IZ9nxlhH+6WN1AKZHv1dUTV5GaCYAwMws2u7XSZeogJm1Bw4lJinG+AuhyeKiKqx6DqWr76k4iNLNAKcW/27RhWJ3H1uZFcZ/r2Z2HNDOzCwmGXQAFka/LwVudfdb49cVtWt3MLP6CZLB/xK+vx7uvtLMziC03aci6TYjxZ9lSc3FzDoSTmZuia5pxSv1WcY5l9B8+WNCEtidkGyKbxioyuOP6wP7x8TXG5hCuHD7eky5wwg3Knwe/nxpAjQxs2+Bdu6+LW69c4lp4jKz/YBGlL5xYT2hVvIa8IyZnRnVmGs11Qgy517gBXf/e3SWfw3h7qB0nEVcbWYtowPrFcDfEpSZBvzIwu1w9c3sbKArod0Z4D/AfuVs4yngFDM7KWqy+S3hTOjdygZrZk2jg+KLwIdRbPH+BvSPtltZ04hrGjGzhtEFUgMamFljK31R/TjCGW4mvQdsBX4VfQdnEg5OxR4ERprZ4RbsauF2yOaEz+kb4LZoemMzK74I2pxw9rrazNpRtqmmPOVtE+I+y2j9bwD3ufv4JOss77NsTvi7WUloIv193Pxy/w7N7EAzG2hmTcysQVTLPpbograZHUw4ubnc3eMvNr9CqJH3ioYbgI+BXgmSAITrFaeZ2TFRs+PNwHPuXqqpLBofQKhFPRldh6jVlAgyIDpDO5qYf1B3fwgovhaws14EPgJmEy4IPxxfwN1XAqcSDuArCYnoVHdfERX5f8AQC/dBl+kH4O7zgfOAe4AVhGsLp7n75krEea+ZFRH+2e8CngUGJGrfdfcf3H1mzG2L8dpa2X4ExRdeXwIOjJoCir0K/EC4pjEh+v1YgChBDAIer8S+VFr0WZ1JuFD7PeG6y3Mx8wsIbfb3RvMXRGWJDlSnAZ2BJYS/nbOjRW8i3CiwhvD9J73+lCCmpNuM/AUYZNHtzMBwwoH6xtjPvriwhds4uxIuiifyF0Jz2NfA58D7cfMfBrpGzVSJ1mGE6wrfAcsJJz5nu3vxbcG/JdzM8HBMfHOjfd3k7t8WD4TPa0v0e3H868zsmKj8XMKNG5Oi7TUn3IJchruvBk4m1J7+Yknu3KstzPVimlrFzBzo4u4Lsh1LTWLhVsSu7v7rFMpeDrR392syHlgtZGa/B75z97tSKPtHYKG7J7zfXmoHJYJaRolARNKtVldnRERk56lGICKS41QjEBHJcbWuH0GrVq28U6dO2Q5DRKRW+eijj1a4e+tE82pdIujUqRMFBckeMyIiIomYWfyTBkqoaUhEJMcpEYiI5DglAhGRHKdEICKS45QIRERyXMYSgZk9YmbfmdlnSeabmd1tZgssvFLvkEzFIiJSU02aBJ06gRnUr1/65y67hJ+xQ/PmYZl0ymSN4DHCo1qTGUh4fV0Xwgup/5zBWEREKlR8UN5ll/Az/oCbbH5lpsdOa9UKhg2DxdGNndu2lf6Z6MEP69bBhRemORm4e8YGwrPAP0sy7wFgaMz4fGCfitZ56KGHuojUPRMnunfs6G4Wfk6cmL7yqZSdONG9aVP3cPgNQ9OmO8omm3/ppYmnjxjh3qRJ6en167s3bFh6WlWHjh0r9/kCBZ7sWJ1sRjqGChLBy8DRMeOvA32SlB1BeFl1QYcOHSq39yKy0yp7kK7sMhUdhKtSfsMG98WL3e+4w71x47Jlx41zf+wx9/Hj3e+6y71Fi8QH3KZN3c88s+xBvaLBLD0H/PLWXxk1NRFMTZAIDq1onaoRiFRdVQ/olTlIl7fMhAnhAO3uXlTk/s9/ur/6qnvr1okPdrvt5v7LX7pfdJH72We7Dx7sPmtWiD3ZAfKdd8L6H3us/ANpy5apH3S7dcvsQT3bNYJsPmKikPAO3GJ5JH73rogkMGkSjBkDS5ZAhw5w662Qn19++REjYMOGML54cRiHxMtt3w4LF8JVV+1YptiGDXDddfDZZ6HNunhYvz60X48Zk3iZESPCYWzECJg/H446inKtXQt/+xs0aQJNm4af69eHfU5mn33CzyOPhAcfhIsvTlxu9Wr46ito1AgaN4aePWHp0rLlOnYM+9mp0462/Fj16u1o009lejrUrx++77RJliHSMVB+jeAUwjtFDegHfJjKOlUjkLqm+Cwd3OvV23G2V5Uz7gcfdF+1ake5115znzzZ/YEHkjd9NGzoftBB7nl57rvv7v7rX4dlf/ih/DNSs7Bsy5bu7duHdfTt6/744+U3i8yZE9a/Zo379Onub7/tvvfelTvrTVYjSFQ+1bLpvkaQaHqDBmWvETRo4L7nnqW//+KfiT7HZs1Sq8nFIxtNQ8BfCS/f3kI4+/8F4X2gI6P5BtwHLAQ+Jcn1gfhBiUBquokTd/xjQ/i9Mm3dsQeT3/3O/Ykn3O+5x/3mm91/85vQtl1e08ipp+5Y/z77lH8wLx6GDAnNL7/6lfuzz4Zlt28P207WbFNe00RlDtTJPoedvUZQ1bLlNZ0lm1+Z6VVpnkuHrCSCTA1KBFLdKnvRs0GDsgfA+vXdL7nE/c473f/4xx3ld9ut/AN0o0alx3fd1X3QoPLPuF96acf6P/7Yfe5c98LCcNZe2QN68T6l6xrBzt7ZU9Xy2Tr41iRKBCKevguljRvvWHbatHCWft557v37J04C8cOee+5Yf0Vlzdznz3f/7jv3zZt3LFfZM+5k+1LRwXlnP7tcP/jWJEoEkvPKOwiuWeP+xRfub7wRxseNc1+7Niy3xx6JD7h5eWH+6NHhLH3ffd0PP7zig/qqVe7btu2Iq7wmnnQ2pcQup4NzblIikDqhKgex9evd//Wv5O3cyQ70n34alk/WBFN8D/fWraW3V96BPdFBvaJrBGrukHQpLxHUujeUSW5Kduvj5s1w4IFhPHa47jo4+mh4/XUYPDj5er//Hm6/Hdq2Dbcdtm0bhubNw/wOHRLfMtihQ/hZr17p6bfeChddBFu2lJ7esGHi2/2Kb9scMyZsp/iWw44dK74dND+//PkiKUuWIWrqoBpBbtm2zX3BguRn9PF3xbRs6d6rV2i7d3f/z3/cn3mm8rcnFqvqRc9U7xoSqS6oRiDZVlHnp6Ki0MFo3rww9OwJP/0prFoFnTsnX++338LLL4cz6I4dd5zJF2vTBs46CzZuLF2jgNBBqaJOObFn7Kl23NKZutQ6yTJETR1UI6h9Ep1VN2oUpm/f7t65c+l59eq5X3HFjuUfe8x9r72qdkYfH4fa1CVXUU6NwML82qNPnz5eUFCQ7TCkElq2DN3543XsCIsWwfXXw667hrb+Aw6A/fcPbeqx4q8RQDijnzBBZ98iqTCzj9y9T6J5ahqStFi7FmbNgg8+CMOKFfDPf4Z5iZIA7HhezC23VLz+qjTRiEhq9KpKSai8F3Rs3QoffxweSgZwww3QogX8+MfhQD1vXmjX37o1zC++wyZesunJ5OeHGsT27eGnkoBIeqhGIGUkulVz2LAwvagIPvoIfvgBPv8cDjoIjjkGbroJDj8c+vYNTUGxfv/7ql2oFZHqoWsEUkb79lBYmHjeEUeEA/7hh8PAgbD77qmts7KPTBaR9NI1AimR6IB87rkwdy5MnRqGZEnADN59t2rb1S2VIjWXEkEOKa937qWXwqZN0KsX7LZbuPgbr7Jt+iJSO+hicR2U7ELvddclfmvUTTfB88+HmsDHH8P994c2/Fhq0xepu1QjqGOSnfVPnZr4NXwQmokGDtwxrls1RXKLLhbXMcneq9qqVbjTZ/36svOKO3aJSN1V3sViNQ3VMcle6r1yJTzwgJp8RKQsJYI6YvNmOPXU8ASeRDp0CE07EyaEGoBZ+KlHNIiIEkEt5g6ffhp+b9gwNP8MGQJNmpQuF3vWr965IhJPiaAWiL8LaOLEcJfPYYeF2z0XLgzlHnsMnn4aHnxQZ/0ikjpdLK7hEj110yzUBvbbD0aPhgsugEaNshejiNR86llci40ZU/bef/fQDDR/PtTXNygiO0lNQzVceXcBKQmISDooEdRgU6aEZqBE9LgHEUkXJYIa6IcfYNQoOP10yMuDxo1Lz9e9/yKSTkoENYx76A9w333wm9/Av/8NDz2ku4BEJHN011ANUfw69l12genTw7QBA7Ibk4jUHXrERA0T3y9g/PjQDHT77WH+gAFKAiJSfXTfSTVL9HTQSy+FevXg5JOzG5uI5CbVCKpZon4BAG3awOWXV388IiIZTQRmNsDM5pvZAjMbnWD+7mb2kpl9YmZzzeyiTMZTEyTrF/Dtt9Ubh4hIsYwlAjOrB9wHDAS6AkPNrGtcscuAz929J3A88Ecza5ipmGqCZPf/q1+AiGRLJmsEhwEL3P1Ld98MTAZOjyvjQHMzM6AZsArYmsGYsu7II8OTQmOpX4CIZFMmE0E7IPbliIXRtFj3AgcBy4BPgSvcfXv8isxshJkVmFnB8uXLMxVvxj36KPz1r9Cnj/oFiEjNkcm7hhI9HCG+08JPgNnAicD+wGtm9ra7ry21kPsEYAKEfgTpDzXzXnoJLr4Y+vcPv8fXCkREsiWTNYJCoH3MeB7hzD/WRcBzHiwAvgIOzGBMWfHOO/Czn8Ehh8CzzyoJiEjNkslEMAvoYmb7RheAzwGmxJVZApwEYGZ7AQcAX2YwpqyYOTM0AU2dCs2aZTsaEZHSMpYI3H0rMAqYAXwBPOXuc81spJmNjIrdAhxpZp8CrwPXuvuKTMVUXeJ7DnfpAh9+CK1bZzsyEZGy9KyhNEv0RrGmTXVBWESyS88aqkaJeg5v2BCmi4jUREoEaZas53Cy6SIi2aZEkGbt2yeerp7DIlJTKRGk2fnnl52mnsMiUpMpEaTZTTeF10x26KCewyJSO+iuIRGRHKC7hqqBe3ir2EMPZTsSEZHKUSJIk5deghkzwpvGRERqEyWCNNi+HW64ATp3TnyxWESkJtM7i9Pguefgk0/giSegvj5REallVCPYSdu3w9ixcOCBMHRotqMREak8nb/upF12gfHjYcsWXR8QkdpJiSANjj462xGIiFSdmoZ2wpNPwqWXln3InIhIbaIaQRVt2QLXXw+77w5NmmQ7GhGRqlMiqKLHH4cvvwz9ByzR25lFRGoJNQ1VwebNcMstcNhhcMop2Y5GRGTnqEZQBQ8/HN4vMGGCagMiUvupRlAF/fuHGkH//tmORERk56lGUAX77w+/+122oxARSQ/VCFIwaRJ06hQ6jzVrBn/4Q7YjEhFJHyWCCkyaBCNGwOLF4VHT69fDjTeG6SIidYESQQXGjCnbYWzTpjBdRKQuUCKowJIllZsuIlLbKBFUoEOHyk0XEaltlAgqcOut0LRp6WlNm4bpIiJ1gRJBBfLzQ8exjh1D57GOHcN4fn62IxMRSQ/1I0hBfj6ce656EYtI3aQaQYratw/vJRYRqWuUCFJQVARffw277prtSERE0i+jicDMBpjZfDNbYGajk5Q53sxmm9lcM3srk/FU1eLF4WenTlkNQ0QkIzJ2jcDM6gH3AScDhcAsM5vi7p/HlGkB3A8McPclZtYmU/HsjEWLwk8lAhGpizJZIzgMWODuX7r7ZmAycHpcmXOB59x9CYC7f5fBeKpMiUBE6rJMJoJ2wNKY8cJoWqwfAS3N7E0z+8jMLki0IjMbYWYFZlawfPnyDIWb3I9+BMOGQZsaWV8REdk5mbx9NNHNlp5g+4cCJwFNgPfM7H13/3ephdwnABMA+vTpE7+OjOvfX+8eEJG6K5M1gkKgfcx4HrAsQZnp7r7e3VcA/wB6ZjCmKlm9Ojx5VESkLspkIpgFdDGzfc2sIXAOMCWuzIvAMWZW38yaAocDX2Qwpirp0gUuvzzbUYiIZEbGmobcfauZjQJmAPWAR9x9rpmNjOaPd/cvzGw6MAfYDjzk7p9lKqaqWLcOVqyAvLxsRyIikhkZfcSEu08DpsVNGx83Pg4Yl8k4dob6EIhIXaeexRXQraMiUtcpEVRAiUBE6jolggr07RseNrfXXtmOREQkM/QY6gocdlgYRETqKtUIKvDFF6EfgYhIXaVEUIFjjoHrrst2FCIimaNEUI6iIli5UheKRaRuUyIoh/oQSGWtXLmSXr160atXL/bee2/atWtXMr5582YABgwYQM+ePenWrRsjR45k27ZtANx555107dqVHj16cNJJJ7G4+A8wsnbtWtq1a8eoUaMSbnvTpk2cffbZdO7cmcMPP5xF0S1vixYt4uCDDy4p9+CDD3LIIYfw/fffZ+ATKG3MmDG0b9+eZs2a7fS6hg8fzueff15xwci9995L586dMTNWrFix09vfWcOGDaNNmzalvouaQomgHLp1VCprzz33ZPbs2cyePZuRI0dy5ZVXlow3bNgQgKeeeopPPvmEzz77jOXLl/P0008D0Lt3bwoKCpgzZw5DhgzhmmuuKbXu66+/nuOOOy7pth9++GFatmzJggULuPLKK7n22mvLlHniiSe45557ePXVV2nZsmUa9zyx0047jQ8//DAt63rooYfo2rVryuWPOuooZs6cSceOHdOy/Z114YUXMn369GyHkZASQTmUCCQTdtttNwC2bt3K5s2bMQsP6j3hhBNo2rQpAP369aOwsLBkmY8++oj//Oc/9C/nMbgvvvgiP//5zwEYMmQIr7/+Oh7ztMSnnnqK2267jVdffZVWrVqVWX758uWcddZZ9O3bl759+/LPf/4TgFWrVnHGGWfQo0cP+vXrx5w5cwB46623Smo7vXv3pqioqMw6+/Xrxz777FPu5zF27Fh+/vOf079/fzp16sRzzz3HNddcQ/fu3RkwYABbtmwB4Pjjj6egoIDFixfTpUsXVqxYwfbt2znmmGN49dVXy6y3d+/edKrgn3fbtm1cddVVdO/enR49enDPPfcAcPPNN9O3b18OPvhgRowYUfI53n333SW1tnPOOQeA9evXM2zYMPr27Uvv3r158cUXE27r2GOPZY899ig3nqxx90oPwEVVWS4dw6GHHurVZd489wkT3Ldvr7ZNSh1y4403+rhx49zdvWfPnqXm9e/f31u0aOFDhw71rVu3lln2sssu81tuucXd3bdt2+bHHXecL1myxB999FG/7LLLEm6vW7duvnTp0pLx/fbbz5cvX+5fffWVN2vWzFu3bu2FhYVJ4x06dKi//fbb7u6+ePFiP/DAA93dfdSoUT527Fh3d3/99ddL9uXUU0/1d955x93di4qKfMuWLQn31d191113TbrdG2+80Y866ijfvHmzz54925s0aeLTpk1zd/czzjjDn3/+eXd3P+6443zWrFnu7v7ggw/6WWed5bfffruPGDGiZF0DBw70r7/+utT6O3bs6MuXL0+47fvvv9/PPPPMkthXrlxZ6qe7+3nnnedTpkxxd/d99tnHN27c6O7u33//vbu7X3fddf7EE0+UTOvSpYuvW7fOv/76ax84cGCp7X311VferVu3pJ9FJgEFnuS4WtUawU3pS0U11wEHwMUXgyV6s4JIJcyePbvU+IwZM/jmm2/YtGkTb7zxRql5EydOpKCggKuvvhqA+++/n0GDBtG+fXvK4wmelV5c22jdujUdOnTgqaeeSrr8zJkzGTVqFL169WLw4MGsXbuWoqIi3nnnHc4//3wATjzxRFauXMmaNWs46qij+M1vfsPdd9/N6tWrqV+/fsJ9TcXAgQNp0KAB3bt3Z9u2bQwYMACA7t27l1zriDV8+HCKiooYP348d9xxR8n0adOm0bZt25S3O3PmTEaOHFkSe/EZ+9///ncOP/xwunfvzhtvvMHcuXMB6NGjB/n5+UycOLFkmVdffZXbbruNXr16cfzxx7Nx40aWLFlC27ZtmTZtWuIN1zBJO5SZ2Zxks4Cc6Gf75pvQvj3sv3+2I5G6qHHjxgwePJgXX3yRk08+GQgHpltvvZW33nqLRo0aAfDee+/x9ttvc//997Nu3To2b95Ms2bNuO2220qtLy8vj6VLl5KXl8fWrVtZs2YNe+yxB0VFRTRt2pRXXnmFo48+mjZt2pCfn18mnu3bt/Pee+/RpEmTUtOTJZjRo0dzyimnMG3aNPr168fMmTM58MADq/RZFO/rLrvsQoMGDUoS2C677MLWrVvLlN+wYUNJ09m6deto3rx5lbbr7iXbKrZx40Z++ctfUlBQQPv27Rk7diwbN24EYOrUqfzjH/9gypQp3HLLLcydOxd359lnn+WAAw6oUgw1QXk1gr2AC4DTEgwrMx9a9g0ZAjEnGyI7bd26dXzzzTdAuEYwbdq0koPnxx9/zCWXXMKUKVNoE/Ne1EmTJrFkyRIWLVrEHXfcwQUXXFAmCQAMHjyYxx9/HIBnnnmGE088sdRBrnXr1kyfPp3//u//ZsaMGUC4s+bee+8FoH///iW/w44z+2OPPZZJkyYB8Oabb9KqVSt22203Fi5cSPfu3bn22mvp06cP8+bNS/lzeP7557luJzroXHvtteTn53PzzTdz8cUXV2rZDz/8kAsuCG/F7d+/P+PHjy9JNqtWrSo56Ldq1Yp169bxzDPPACFRLl26lBNOOIHbb7+d1atXs27dOn7yk59wzz33lCTMjz/+uMr7lS3lJYKXgWbuvjhuWAS8WS3RZZH6EEg69erVCwgXFgcPHkyPHj3o2bMnbdq0YeTIkQBcffXVrFu3jp/+9KclzTMVueGGG5gyJbzv6Re/+AUrV66kc+fO3HnnnQmTxb777suUKVMYNmwYH3zwAfPmzWPPPfcEwoXQgoICevToQdeuXRk/PjwxfuzYsSXTR48eXZJs7rrrLg4++GB69uxJkyZNGDhwYKl9BbjmmmvIy8tjw4YN5OXlMXbsWAAWLlxYctG8st566y1mzZpVkgwaNmzIo48+CsCgQYNYtmxZyf7k5eVRWFhIjx49GD58OABLliwpqfUMHz6cDh06lHwfTz75JC1atODiiy+me/funHHGGfTt2xcIF5bPO+88unfvTu/evbnyyitp0aIF119/PVu2bKFHjx4cfPDBXH/99QAsW7aMQYMGlcQ9dOhQjjjiCObPn09eXh4PP/xwlfY/I5JdPABeAfZNNj9bQ3VdLP70U3dwnzy5WjYnkhWnnHKKb9q0qdq3m5+f79999121b9fd/aqrrvJPPvkkK9vOJsq5WFzeQ+ceBWaY2ePA7e6+pRryUo2hzmSSC15++eWsbHfixIlZ2S7AuHE19j1YWZM0Ebj7U2Y2FbgBKDCzJwivkyyef2c1xJc16kMgIrmiosdQbwHWA42A5sQkgrrurLNgv/0g5pqdiEidVN7towOAO4EpwCHuvqHaoqoB9t4bomtfIiJ1Wnk1gjHAT919bnUFU5NMnhyahfr1y3YkIiKZlfT2UXc/JleTAMDll8Njj2U7ChGRzNND5xJYtw5WrNCFYhHJDUoECejWURHJJUoECejWURHJJUoECSgRiEguUSJI4KKL4LPP1IdARHJDRR3KclLTptCtW7ajEBGpHqoRJHD33TB1arajEBGpHkoECdxyC7z0UrajEBGpHkoEcdSHQERyTUYTgZkNMLP5ZrbAzEaXU66vmW0zsyGZjCcV6kMgIrkmY4nAzOoB9wEDga7AUDPrmqTcH4AZmYqlMnTrqIjkmkzWCA4DFrj7l+6+GZgMnJ6g3OXAs8B3GYwlZdH7sJUIRCRnZDIRtAOWxowXRtNKmFk74L+A8eWtyMxGmFmBmRUsX7487YHGuuQSWL0a9toro5sREakxMpkILME0jxu/C7jW3beVtyJ3n+Dufdy9T+vWrdMVX1K77w6WKHoRkTookx3KCoH2MeN5wLK4Mn2AyRaOuq2AQWa21d1fyGBc5RozBrp0gQsvzFYEIiLVK5M1gllAFzPb18waAucQ3nZWwt33dfdO7t4JeAb4ZTaTAMCECfD++9mMQESkemWsRuDuW81sFOFuoHrAI+4+18xGRvPLvS6QDepDICK5KKPPGnL3acC0uGkJE4C7X5jJWFKhPgQikovUsziG+hCISC5SIoixdi00b65EICK5RYkgxtChsGaN+hCISG7R+wjiqP+AiOQa1QhiXHwxjBuX7ShERKqXEkGMF16AhQuzHYWISPVSIoioD4GI5Colgoj6EIhIrlIiiKgPgYjkKiWCiDt066ZEICK5R4kgcuqp8NlnsPfe2Y5ERKR6KRGIiOQ4JYLIGWfAVVdlOwoRkeqnRBB5910oKsp2FCIi1U+JAFi/HpYv14ViEclNSgSoD4GI5DYlAtSHQERymxIBsOuuMGAA7LdftiMREal+egw1cNxxYRARyUWqERB6FYuI5ColAuCoo+DCC7MdhYhIdigRAP/3f9C4cbajEBHJjpxPBHoPgYjkupxPBOpDICK5LucTgfoQiEiuy/lEsNdeMHw4dO6c7UhERLIj5/sR9OkTBhGRXJXzNYLVq2H79mxHISKSPTmfCE4+GU47LdtRiIhkT84ngkWLoH37bEchIpI9GU0EZjbAzOab2QIzG51gfr6ZzYmGd82sZybjibd+fehD0LFjdW5VRKRmyVgiMLN6wH3AQKArMNTMusYV+wo4zt17ALcAEzIVTyLqQyAiktkawWHAAnf/0t03A5OB02MLuPu77v59NPo+kJfBeMpQHwIRkcwmgnbA0pjxwmhaMr8AXkk0w8xGmFmBmRUsX748bQHuvz/cdBMccEDaVikiUutksh+BJZiW8IHPZnYCIREcnWi+u08gajbq06dP2h4afcABcMMN6VqbiEjtlMkaQSEQez9OHrAsvpCZ9QAeAk5395UZjKeMBQvCS+tFRHJZJhPBLKCLme1rZg2Bc4ApsQXMrAPwHHC+u/87g7EkdO65cN551b1VEZGaJWNNQ+6+1cxGATOAesAj7j7XzEZG88cDNwB7AvebGcBWd6+2Bz4sWgS9e1fX1kREaqaMPmvI3acB0+KmjY/5fTgwPJMxJLN+fWgW0h1DIpLrcrZnsfoQiIgEOZsI1IdARCTI2UTQowc88ggcdFC2IxERya6cfR9BXh5cdFG2oxARyb6crRG89x58+mm2oxARyb6crRFccQW0bAkzZmQ7EhGR7MrZGsGiRbpQLCICOZoI1IdARGSHnEwE6kMgIrJDTiYC9SEQEdkhJxPBEUfAa69B9+7ZjkREJPty8q6hli3hxz/OdhQiIjVDTtYIXn4ZZs7MdhQiIjVDTtYIbr5ZtQIRkWI5WSNQHwIRkR1yLhGoD4GISGk5lwjUh0BEpLScSwTqQyAiUlrOXSw+6SSYNw86dMh2JCIiNUPOJYJGjeCAA7IdhYhIzZFzTUOPPw4TJ2Y7ChGRmiPnEsF998ETT2Q7ChGRmiPnEoH6EIiIlJZTiUB9CEREysqpRKA+BCIiZeVUIigsDD+VCEREdsip20f794d166Bhw2xHIiJSc+RUIgDYdddsRyAiUrPkVNPQn/4E48ZlOwoRkZolpxLB5MnhFZUiIrJDTiUC9SEQESkro4nAzAaY2XwzW2BmoxPMNzO7O5o/x8wOyUQckyaFh8x99x08/XQYFxGRIGMXi82sHnAfcDJQCMwysynu/nlMsYFAl2g4HPhz9DNtJk2CESNgw4Ywvnp1GAfIz0/nlkREaqdM1ggOAxa4+5fuvhmYDJweV+Z04C8evA+0MLN90hnEmDE7kkCxDRvCdBERyWwiaAcsjRkvjKZVtgxmNsLMCsysYPny5ZUKYsmSyk0XEck1mUwElmCaV6EM7j7B3fu4e5/WrVtXKohkL6DRi2lERIJMJoJCoH3MeB6wrApldsqtt0LTpqWnNW0apouISGYTwSygi5nta2YNgXOAKXFlpgAXRHcP9QPWuPs36QwiPx8mTICOHcEs/JwwQReKRUSKZeyuIXffamajgBlAPeARd59rZiOj+eOBacAgYAGwAbgoE7Hk5+vALyKSTEafNeTu0wgH+9hp42N+d+CyTMYgIiLly6mexSIiUpYSgYhIjlMiEBHJcUoEIiI5zsL12trDzJYDi6u4eCtgRRrDqQ20z7lB+5wbdmafO7p7wh65tS4R7AwzK3D3PtmOozppn3OD9jk3ZGqf1TQkIpLjlAhERHJcriWCCdkOIAu0z7lB+5wbMrLPOXWNQEREysq1GoGIiMRRIhARyXF1MhGY2QAzm29mC8xsdIL5ZmZ3R/PnmNkh2YgznVLY5/xoX+eY2btm1jMbcaZTRfscU66vmW0zsyHVGV8mpLLPZna8mc02s7lm9lZ1x5huKfxt725mL5nZJ9E+Z+QpxtXFzB4xs+/M7LMk89N//HL3OjUQHnm9ENgPaAh8AnSNKzMIeIXwhrR+wAfZjrsa9vlIoGX0+8Bc2OeYcm8QnoI7JNtxV8P33AL4HOgQjbfJdtzVsM//Dfwh+r01sApomO3Yd2KfjwUOAT5LMj/tx6+6WCM4DFjg7l+6+2ZgMnB6XJnTgb948D7Qwsz2qe5A06jCfXb3d939+2j0fcLb4GqzVL5ngMuBZ4HvqjO4DElln88FnnP3JQDuXtv3O5V9dqC5mRnQjJAItlZvmOnj7v8g7EMyaT9+1cVE0A5YGjNeGE2rbJnapLL78wvCGUVtVuE+m1k74L+A8dQNqXzPPwJamtmbZvaRmV1QbdFlRir7fC9wEOE1t58CV7j79uoJLyvSfvzK6ItpssQSTIu/RzaVMrVJyvtjZicQEsHRGY0o81LZ57uAa919WzhZrPVS2ef6wKHASUAT4D0ze9/d/53p4DIklX3+CTAbOBHYH3jNzN5297UZji1b0n78qouJoBBoHzOeRzhTqGyZ2iSl/TGzHsBDwEB3X1lNsWVKKvvcB5gcJYFWwCAz2+ruL1RLhOmX6t/2CndfD6w3s38APYHamghS2eeLgNs8NKAvMLOvgAOBD6snxGqX9uNXXWwamgV0MbN9zawhcA4wJa7MFOCC6Op7P2CNu39T3YGmUYX7bGYdgOeA82vx2WGsCvfZ3fd1907u3gl4BvhlLU4CkNrf9ovAMWZW38yaAocDX1RznOmUyj4vIdSAMLO9gAOAL6s1yuqV9uNXnasRuPtWMxsFzCDccfCIu881s5HR/PGEO0gGAQuADYQzilorxX2+AdgTuD86Q97qtfjJjSnuc52Syj67+xdmNh2YA2wHHnL3hLch1gYpfs+3AI+Z2aeEZpNr3b3WPp7azP4KHA+0MrNC4EagAWTu+KVHTIiI5Li62DQkIiKVoEQgIpLjlAhERHKcEoGISI5TIhARyXFKBCKVZGbtzewrM9sjGm8ZjXc0s+lmttrMXs52nCKpUiIQqSR3Xwr8GbgtmnQbMMHdFwPjgPOzFZtIVSgRiFTNn4B+ZvZrwnOb/gjg7q8DRVmMS6TS6lzPYpHq4O5bzOxqYDrQP3pEskitpBqBSNUNBL4BDs52ICI7Q4lApArMrBdwMuENUVfW8hcbSY5TIhCppOhNWH8Gfh29CWwccEd2oxKpOiUCkcq7GFji7q9F4/cDB5rZcWb2NvA0cJKZFZrZT7IWpUiK9PRREZEcpxqBiEiOUyIQEclxSgQiIjlOiUBEJMcpEYiI5DglAhGRHKdEICKS4/4/UkWCG0T7yE8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "    plt.plot(X1,Y1,linestyle='--',marker ='o',color = 'b')\n",
    "    plt.title('X-Y plot of DME(1)+n-decane(2) at {}K'.format(T))\n",
    "    plt.text(0.5,0.5,\"T:{}K,eos:{},mix:{},case:{}\".format(T,eos,mix,case))\n",
    "    plt.xlabel('X1')\n",
    "    plt.ylabel('Y1')\n",
    "    plt.savefig(\"X1-Y1-1st-Mix,{},{},{},{},.PNG\".format(T,eos,mix,case))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}