all_in_one_system_4_corrected.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VLE calculation for Ethanol(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",
    "    \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",
    "    \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 Ethanol(1) in range T:194.3K-248.24K\n",
    "    if T >= 273 and T <= 351.70:\n",
    "        A1 = 5.37229\n",
    "        B1 = 1670.409\n",
    "        C1 = -40.191\n",
    "    elif T > 351.70 and T <= 366.63:\n",
    "        A1 = 5.24677\n",
    "        B1 = 1598.673\n",
    "        C1 = -46.424\n",
    "    elif T > 366.63 and T <= 513.9:\n",
    "        A1 = 4.92531\n",
    "        B1 = 1432.526\n",
    "        C1 = -61.819\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 =  0.00005960611 # m3/mol from NIST at 313.15K\n",
    "#v2 = 9.1072*0.00001 # m3/mol at 313.15K\n",
    "    V1_L = 0.00005960611 #Si unit\n",
    "    V2_L = 9.1072*0.00001 # Si unit\n",
    "\n",
    "#Pure component properties\n",
    "\n",
    "    Pc_1 = 63.4*100000  # Pa from NIST\n",
    "    Vc_1 = 0.000168 # m3 per mole from NIST\n",
    "    Tc_1 = 514.7      # 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 ny 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",
    "        #Taken from final paper\n",
    "        A12 = 1.9228\n",
    "        A21 = 1.2278\n",
    "        alpha = 0.5962\n",
    "#        print(\"case is :\",case)\n",
    "    elif case == 2:\n",
    "        # Taken from final paper\n",
    "        A12 = 2.0314 \n",
    "        A21 = 1.2691\n",
    "        alpha = 0.5962\n",
    "#        print(\"case is :\",case)\n",
    "    elif case == 3:\n",
    "        #Taken from final paper\n",
    "        A12 = 1436.139\n",
    "        A21 = 51.6395\n",
    "        alpha = 0.5962\n",
    "#        print(\"case is :\",case)\n",
    "    elif case == 4:\n",
    "        #Taken from final paper\n",
    "        A12 = 501.3384\n",
    "        A21 = 1125.0699\n",
    "        alpha = 0.5962\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",
    "    X1 = [0.0195,0.0399,0.0588,0.1000,0.1494,0.1997,0.2499,0.2988,0.3491,0.4002,0.4540,0.4995,0.5493,0.5997,\n",
    "         0.6499,0.6993,0.7489,0.7993,0.8506,0.9004,0.9400,0.9595,0.9806]\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:325\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:4\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 ny Panagiotopoulos et.al 1986\n",
      "\n",
      "Enter the number corresponding to the mixing rule which is to be used:4\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:4\n",
      "Lenght of vector X1 is : 23\n",
      "P1sat is: 32156.975933143844\n",
      "P2sat is: 38778.91012495706\n",
      "Value of Phi_1_sat is : 0.9914114916920443\n",
      "Value of Phi_2_sat is : 0.9863794493493757\n",
      "phisat1 0.9914114916920443\n",
      "phisat2 0.9863794493493757\n",
      "\n",
      "Final values of P in KPa is : array('f', [40.01987838745117, 39.80637741088867, 39.610897064208984, 39.19255447387695, 38.70505905151367, 38.224464416503906, 37.76055908203125, 37.323524475097656, 36.888916015625, 36.46241760253906, 36.029136657714844, 35.67478942871094, 35.29910659790039, 34.93132400512695, 34.57699203491211, 34.239749908447266, 33.91257095336914, 33.59221267700195, 33.27955627441406, 32.990482330322266, 32.77223587036133, 32.669036865234375, 32.56087875366211])\n",
      "\n",
      "Corresponding values of X1 are : [0.0195, 0.0399, 0.0588, 0.1, 0.1494, 0.1997, 0.2499, 0.2988, 0.3491, 0.4002, 0.454, 0.4995, 0.5493, 0.5997, 0.6499, 0.6993, 0.7489, 0.7993, 0.8506, 0.9004, 0.94, 0.9595, 0.9806]\n",
      "\n",
      "Final values of Y1 are : array('f', [0.015805402770638466, 0.03250739723443985, 0.048133496195077896, 0.08270229399204254, 0.12505939602851868, 0.16919560730457306, 0.21424487233161926, 0.2590755522251129, 0.3061518669128418, 0.35496219992637634, 0.40741094946861267, 0.4526078701019287, 0.5029523372650146, 0.554832935333252, 0.6074375510215759, 0.660112738609314, 0.7139120101928711, 0.7695160508155823, 0.8270789980888367, 0.8838749527931213, 0.9296592473983765, 0.9523970484733582, 0.9771362543106079])\n",
      "\n",
      "Final values of Y2 are : array('f', [0.9841945767402649, 0.9674925804138184, 0.9518665075302124, 0.9172977209091187, 0.8749405741691589, 0.8308044075965881, 0.7857550978660583, 0.7409244775772095, 0.6938481330871582, 0.6450377702713013, 0.5925890207290649, 0.5473921298980713, 0.49704766273498535, 0.44516709446907043, 0.39256247878074646, 0.3398872911930084, 0.2860879898071289, 0.23048393428325653, 0.17292098701000214, 0.11612502485513687, 0.07034078240394592, 0.047602925449609756, 0.022863762453198433])\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",
    "    print(\"phisat1\",phi1sat)\n",
    "    print(\"phisat2\",phi2sat)\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvSElEQVR4nO3deXyU1b3H8c+PTUFRRMAqsmhFEVk1CK0b4AZiEZe2IKKtIkVF6w6t1XK1Wqq9XES0XHBrEQQXVKq41AUVtxIFUVS8QRMIuCCCgJT9d/84T+JkMhOSkMlkMt/36zWvZJ5lnt8zk8zvOec85xxzd0REJHvVSXcAIiKSXkoEIiJZTolARCTLKRGIiGQ5JQIRkSynRCAikuWUCGoRM/uVmc2rpmMdZmYLzGy9mV1RHcfcSTy/MbPx5dz2CjMbm+KQkh37QTP7UzqOXV5m9mczu7Kc244zsxEpDklSTIkgBcxsTzPLN7NzY5Y1NrNlZnZOzLL9zOwbM+sVt/8DZvZwimOca2bDduElrgfmuntjd5+Q5PU3RYlinZm9a2ajzWy3mG3GmJnHJxIzuzJaPiZ63svMdpjZhrjHT6L1DYA/AHfEvMZkM1sS7feruPAmA+eZWYtdOP9aycyaA+cD/xs972lm/zKzb81slZk9amb7x+xyB3BD9BlU5ng7TYxm9kp07HVm9r6ZnRGzrr+ZzTOztWb2pZlNMbPGca+/Je7vpm4ZxzrXzArM7Hsze9LMmsasK/E/E/1drjGzQZU595pEiSAF3H0DMBy4M/rHArgdyHX3x2K2+wq4CphiZg0BzOxEoD+Q9qvsnWgDLN7JNiPdvTGwP3ANMAiYY2YWs82nwAVx+50fLY+10t33jHu8Fa07A/jE3VfEbP8+cCnwXnxQ7r4JeDY6TilRghqzk3OrrX4FzHH3/0TP9yEkzraEz3w98EDRxu7+BfAJMCCFMf0W2N/d9yL8Xz0Uk4z2Bv4EHAAcDhxIzAVB5Pa4v5vtiQ5iZkcQEuBQYD9gI3BPkm1PAZ4ELnT3GbtycjWBEkGKuPsLwDPAhOiK/xfAZQm2mwosAW6OksH/Ale4+6pEr1t0BW1mn0WliTvMLOHnaGY/NbP5ZvZd9POn0fJbgeOAidEV0sQk+w8ws8XR1dZcMzs8Wv4y0Dtm/0N38l587+5zCV8WPyEkuiLzgUbRP2HRP2PDaHl59QNejTvm3e7+ErApyT5z4+Iol6ikd62ZLYre15lmtnsZ23czs/eiktFMYPe49aeb2cLoPX7TzDrHrGtlZrOiq+HVRZ+Tmf3YzF6Oln1jZtPMrEl5YyzrmMS9l+7+rLs/6u7r3H0jMBE4Ju4051LGexmVIr6MYnkt5rMeDgwBro/+jv6ZaH93X+Tu24qeAvWBVtG66e7+nLtvdPc1wJQE8ZXXEOCf7v5adDF3I3BWbAkjivt04BHgXHd/opLHqlGUCFLrKqAX8BhwbXT1lMgI4EJgBvBhOa4wzgRygCMJV8MXxm8QFWmfASYA+wLjgGfMbF93vwF4nXDFvqe7j0yw/6HAw8CVQHNgDvBPM2vg7n3i9o+/ek/I3ZcBuYQkFGsqP1ydXwD8ozyvF6MTIZlWxMdAlwruU+QXQF/gIKAz4Sq6lKi65EnC+TUFHgXOjll/JHA/8BvCZ/S/wGwz2y2qvngaKCBcjbck/H0AGPBnfrgKbgWMKU+MZR0z2m9n7+XxlC4J7uy9fBZoB7QglNCmAbj75Oj3oiv2nyV7ATN72sw2Ae8QEk9uBeK71ELV1rtmdnainSJHEEqSRPEtBbYAsRc6PwMeAs5x9zllvFZGUSJIoegKZTHQCJhVxnaFwE3AScAl5Xjpv7j7t9EX63hgcIJt+gP/5+5T3X2buz9MKMIn/WeL80vgGXf/l7tvBf5KuFL/aTn3T2Yl4Usx1kPAYDOrT6g+eijBfgdEV7Cxjz2idU0IVRYVsZ5QrVAZE9x9pbt/C/wT6Jpku56Eq9fx7r41qhaMLelcDPyvu7/j7tvd/e/A5mi/owlf9NdFJapN7j4PwN3zos9lc1RyHAecUM4YyzomlPFeRiWHm4Dr4latj/ZLyN3vd/f17r6ZkLC6mFmF3nt3Px1oDJwGPO/uOxLEdzLhQuKmmMUT+CEJ3Qg8aGbJSgx7At/FLfsuOm6R3oRqyzcqEn9Np0SQQmZ2HuFq7kXgLzHLn7UfGq6GRIsXA2vKKDXEWh7zewHhCyPeAdE64rZtWc7wS+wf/eMtr8D+ybQEvo1dECW0POA2QvJanmC/le7eJO7xfbRuDSX/WcujMTH/9NEV51ozWwuMBkbHJJyn4/b9Mub3jYQvkESf6wHACi85smPsZ9IGuCY2uRGu7g+IfhbEVIkUM7MWZjbDzFaY2TpC4mxWnhh3ckxI8l6a2SGEK/vfuvvrcasbA2vj94n2q2tmY81saRRrfrQqPt6dipLps8CpZlaiTcLMegLTCVfqn8bs8567r44uhuYQSiBnJTnEBmCvuGV7UTIx3khInE/GlKIynhJBili4I+V/CFdgvwF+YWbHA7h7v5iGq2mVePlWMb+3Jlxlx1tJ+KcnbtuiBtWdDTtbYn8zs+i4K5LusRNm1go4ilCtFO8fhAblilYLASyiZPG9PA6nZDXA6UUJBhgLjI1JOKeX5wUTfK5fAC2j965I65jflwO3xiW3RlHpbTnQ2szqJTjUnwmfX+eoAfU8QnVReZR1TEjwXppZG8LFzC1Rm1a8Eu9lnHMJ1ZcnEUpgbYteNvpZmeGP6wE/jomvGzCb0HD70k72dZK/V4uJqeIys4OB3Sh548L3hFLJ3sBjUSk24ykRpM5E4El3fyW6yr+ecHdQVVxFXGdm+0RfrL8FZibYZg5wqIXb4eqZ2S+BDoR6Z4CvgIPLOMYjQH8zOzH6Y7+GcCX0ZkWDNbNGZnYC8BTw7yi2eDOBU6LjVtQc4qpGzKxB1EBqQH0z291KNqqfQLjCTaW3gG3AFdFncBahyqfIFGCEmfWwYA8Lt0M2JrxPXwBjo+W7x1RpNCZcva41s5aUrqopS1nHhLj3Mnr9l4G73X1Sktcs671sTPi7WU2oIr0tbn2Zf4dm1t7M+plZQzOrH5Wyjydq0DazjsBzwOXuXqqx2czOsXA7dx0Ld/qcR0gaiUwDfmZmx0XVjjcDs9y9RFVZ9LwvoRQ13cq4HTVjuLseVfwABhKuqJvELX+JcDWWaJ9eQGE5XtsJt5Z+Rvjn+m+gbrTuV8C8mG2PBd4lVIG8Cxwbs+4nhCudNYT65ETHOhP4KNr/VeCImHVzgWFlxDmXcMfO+uixALgB2D1mmzHAQ0n2fwgYE/Pe7CB8+cU+zo7W1weWAQfEHd/jHr2idbsDhcB+SY49pujYCdblAyeV5xyi9TnRua8nJLuZwJ9i1vcltBusJXzxPwo0jta1JjQ2rwa+KfqcCI2a70bvwUJCki4sb4w7OWaz6L1pGD3/Y/TelXjvY15r/2j7BknOf0/CBcB6QrXY+dHrHRKtbxedw1rChVP8/ocTGojXR9vMB86MWf9Agr+NxTHrXyf8/a4jlFoGxb3+BuC4mOfnRn9L30dxN032N09o63qfkEDqVOd3TFU/LDohyRBm5kA7d89Ldyw1SXQrYgd3v7Ic214OtHL361MeWAYys9uAr919fDm2/W9gqbsnvN9eMoMSQYZRIhCRqqY2AhGRLKcSgYhIllOJQEQkyyW6R7lGa9asmbdt2zbdYYiIZJR33333G3dvnmhdxiWCtm3bkpubbJgRERFJxMziRxoolrKqITO738y+NrMPk6w3M5tgZnkWRkk8MlWxiIhIcqlsI3iQ0HElmX6EziTtCGOM/y2FsYiISBIpSwTu/hpxg4vFOQP4hwdvA02s5MxHIiJSDdJ511BLSo6iWciuj2wpIiIVlM5EkGgEwISdGsxsuJnlmlnuqlUJJ+4SEZFKSmciKKTkcMoHkng4Zdx9srvnuHtO8+YJ734SEam1pk2Dtm2hTp3wc1plBq8vQzoTwWzg/OjuoZ7Ad16+SVlERLLGtGkwfDgUFIB7+Dl8eNUmg1TePvowYTz2w8ys0MwuMrMRZjYi2mQOYSjlPMIY6ZemKhYRkUx1ww2wcWPJZRs3huVVJWUdytw90Ty6sesduCxVxxcRqQ2WLavY8srQWEMiIjXUm2+GdoFEWrdOvLwylAhERGqA2AbhNm1gyBA44QRo2hR2373kto0awa23Vt2xlQhERNIsvkF42TKYMQO6doVPP4V77w3JwSz8nDw5JIqqknHzEeTk5LgGnROR2qRt25AE4rVunXh5ZZjZu+6ek2idSgQiImmWrOF3+fLEy6uaEoGISBpt2AANGyZeV5UNwmVRIhARqUaxjcItW8Khh4Z+AfXrl9yuqhuEy6JEICJSTeIbhVeuhK++Cp3DHnggtQ3CZVFjsYhINUnWKNymDeTnp/bYaiwWEakBqqOXcGUoEYiIVIOpU5Ovq65G4WSUCEREUuzmm+H88+Gww0rfIVSdjcLJKBGIiKTYWWfBTTfBBx/AlCnpaxRORo3FIiIpMHUqvPMOTJyY7kgCNRaLiFSTjRvhootCVdAHH5SeS6AmUiIQEamk+Ckk77gDevSA++8PfQNeeim0AdR0KZuYRkSkNivqHFZ0xV9QAKNGwR57wHPPwamnpje+ilCJQESkEhJNIekOe++dWUkAlAhERColWSewlSurN46qoEQgIlJB7mHmsETS3TmsMpQIREQqYP16GDoUVq8uPZ9wTegcVhlKBCIiFfDeezBzZugt/OCDNa9zWGXoriERkZ1wh9xc6N49TCi/dOkPVUBDh6Y3tqqgEoGISBnWrIGzzw79A957LyzLxHaAsqhEICKSxNtvw6BBsGJF6CzWtWu6I0oNlQhERBL4n/+B444L9f/z5sE115RuHK4taulpiYjsujPOgAULQrVQbaZEICJZLXa8oB/9CK69Niy/8kp49FFo0iSNwVUTJQIRyVrxk8l/9VWoEpo2LVQJmaU7wuqhRCAiWSvReEE7doTl2USJQESyVk2dTL66KRGISNYpmpixVavE62tbP4GdUSIQkayyZAkceyx89BHcdlvpiWMydbygXZHSRGBmfc1siZnlmdnoBOv3NrN/mtn7ZrbYzH6dynhEJHu5wwMPwJFHhmTw5ZdhXKDJk2vHeEG7ImWT15tZXeBT4GSgEJgPDHb3j2K2+T2wt7uPMrPmwBLgR+6+JdnravJ6Eamo776DSy6Bhx+G3r3DxPItW6Y7quqVrsnrjwby3P2z6It9BnBG3DYONDYzA/YEvgW2pTAmEclCEybAI4+EKp9//Sv7ksDOpHKsoZbA8pjnhUB8/7yJwGxgJdAY+KW774h/ITMbDgwHaJ1trTgiUik7doTZwg48EK6/Hvr1g5yE18OSyhJBoq4Y8fVQpwILgQOArsBEM9ur1E7uk909x91zmjdvXtVxikgt8+WXYd7g446DDRtgt92UBMqSykRQCMTenHUg4co/1q+BWR7kAZ8D7VMYk4jUcs89B507wxtvwO9/D3vske6Iar5UJoL5QDszO8jMGgCDCNVAsZYBJwKY2X7AYcBnKYxJRGqR2HGC2rSB/v1DFdCPfhQmkrn44uwZJmJXpKyNwN23mdlI4HmgLnC/uy82sxHR+knALcCDZvYBoSpplLt/k6qYRKT2KBonqGiIiGXLoLAQTjoJZs+Ghg3TG18mSdnto6mi20dFBEJJoKCg9PI2bSA/v7qjqfnSdfuoiEjKaJygqqNEICIZZ9685LOF6Q7zilMiEJGM8thjcMIJ0LQp7L57yXXZOE5QVVAiEJGMUNSceeKJ8NvfwtKlcO+9GieoKqixWERqvOnT4b774NlnoUGDdEeTmdRYLCIZad06GDo0XOVv3gxr16Y7otpJiUBEaqS33oKuXcOIof/1XzB3LrRoke6oaqdUDjonIlIpO3aEYaPd4fXX4Sc/SXdEtZtKBCKSNrFDRLRtC3feGQaJq1MHHn8cFi5UEqgOSgQikhZFQ0QUFIQr/4ICuPJKGDgwrP/xj2HvvdMZYfZQIhCRtLjhhh/GCYr18cfVH0u2UyIQkbRINhTEF19UbxyiRCAiaZJsKAgNEVH9lAhEpFr93//B6NFhKIhGjUqu0xAR6aFEICLVwj0MCdGtWxgKokeP8FNDRKSf+hGISMqtWhVmC3vqqTBW0IMPhknlDzlEX/w1gRKBiKSUO/TtCx9+COPGhQHjkg0hLemhRCAiKfGf/0C9elC/PowfH/oEdO6c7qgkEeVlEalyCxbAUUfBn/4Unh93nJJATaZEICK7JHaYiDZtYNCg0BD83Xdw7LHpjk7KQ1VDIlJpRcNEFPUQXrYsPLp3D3MH7LtveuOT8lGJQEQqLdkwEV99pSSQSZQIRKTSkg0TsXx59cYhu0aJQEQq5fnnk98GqmEiMosSgYhUyMaNMHJk6Buw336w++4l12uYiMyjRCAi5bZjB5xwAtx9d5g7IC8vDBuhYSIym7l7umOokJycHM/NzU13GCJZZds2qFs3fNk/+ig0bRqGipDMYWbvuntOonUqEYhImZYsCdNFTpkSnv/850oCtY0SgYgktGMHTJwYRgv97DNo3jzdEUmqqEOZiJSyYgX8+tfwr39Bv35w332w//7pjkpSRSUCkSwXO0RE27bh+QcfwJtvwqRJ8MwzSgK1nUoEIlksfoiIgoLwfPJkyM+HZs3SGp5Uk5SWCMysr5ktMbM8MxudZJteZrbQzBab2aupjEdESko0RMTGjWG5kkD2SFmJwMzqAncDJwOFwHwzm+3uH8Vs0wS4B+jr7svMrEWq4hGR0pINEZFsudROqSwRHA3kuftn7r4FmAGcEbfNucAsd18G4O5fpzAeEYmxZUvoG5CIhojILqlMBC2B2KGnCqNlsQ4F9jGzuWb2rpmdn+iFzGy4meWaWe6qVatSFK5Idti2Lfxs0ACGDoXddiu5XkNEZJ9UJgJLsCy+G3M94CigP3AqcKOZHVpqJ/fJ7p7j7jnNdTOzSKXNnx9mCnvmmfD8/vvDraEaIiK7pfKuoUKgVczzA4GVCbb5xt2/B743s9eALsCnKYxLJOts2RKmjbzttnAraMOGP6wbMkRf/NkulSWC+UA7MzvIzBoAg4DZcds8BRxnZvXMrBHQA/g4hTGJZJ0PPghTR95yC5x3Xnjep0+6o5KaJGUlAnffZmYjgeeBusD97r7YzEZE6ye5+8dm9hywCNgB3OvuH6YqJpFsNH8+rFwJTzwBAwemOxqpiTT6qEgt9Omn8MknMGAAuMPatbDPPumOStJJo4+K1GKxQ0S0aQPnnw9du8Lll4e2ATMlASmbEoFIBisaIqKgIFz5L1sGU6dCu3bw1lvhFlGRnVEiEMlgiYaIgFAVdMAB1R6OZCglApEMlmwoiOXLEy8XSUSJQCQDuYfOYJao2yYaIkIqRolAJMMUFkL//nDRRXDooSU7h4GGiJCKUyIQySBTp0LHjvDqq3DXXbB4cZhLWENEyK5QIhCpIqtXr6Zr16507dqVH/3oR7Rs2bL4+ZYtWwDo27cvXbp04YgjjmDEiBFs374dgHHjxtGhQwc6d+7MiSeeSEFBQfHr1q1bt/h1br11AF26wKJFMHJkuGV0yBBYsmQzP//5L6lX7xAmTOhBfn4+APn5+XTs2LH4taZMmcKRRx7JmjVrqu19GTBgQIkYKmPYsGF89NFHO98wzuWXX86ee+65S8euKtu3b6dbt26cfvrp6Q6lFM1QJlJF9t13XxYuXAjAmDFj2HPPPbn22mtLbPPII4+w11574e6cc845PProowwaNIhu3bqRm5tLo0aN+Nvf/sb111/PzJkzcYf69RsyduxC+vYNI4fWqRMese677z722Wcf8vLymDFjBqNGjWLmzJkltpk6dSp33XUXL7/8MvtUU8eCWbNmVckX8b333lvhfXJzc1m7du0uH7uq3HnnnRx++OGsW7cu3aGUohKBSDXaa6+9ANi2bRtbtmzBotbe3r1706hRIwB69uxJYWEhK1fCz34GmzfDAw+E/evVK50EAJ566ikuuOACAM455xxeeuklYkcNeOSRRxg7diwvvPACzRJMPbZq1SrOPvtsunfvTvfu3XnjjTcA+Pbbbxk4cCCdO3emZ8+eLFq0CIBXX321uJTSrVs31q9fX+o1N2zYwLhx4/jDH/6Q9P0YM2YMF1xwAaeccgpt27Zl1qxZXH/99XTq1Im+ffuydetWAHr16kVubi4FBQW0a9eOb775hh07dnDcccfxwgsvlHrd7du3c91113H77bcnPfb27du59tpr6dSpE507d+auu+4C4Oabb6Z79+507NiR4cOHF7+PEyZMKC61DRo0CIDvv/+eCy+8kO7du9OtWzeeeuqphMcqLCzkmWeeYdiwYUnjSSt3z6jHUUcd5SI13R//+Ee/44473N29S5cuJdadcsop3qRJEx88eLBv27at1L6XXnqZn3XWLd6kiXvDhu516tT1o446ynv06OFPPPFEwuMdccQRvnz58uLnBx98sK9atco///xz33PPPb158+ZeWFiYNN7Bgwf766+/7u7uBQUF3r59e3d3HzlypI8ZM8bd3V966aXiczn99NN93rx57u6+fv1637p1a6lzvfLKK33WrFn++eef+xFHHJH0fTrmmGN8y5YtvnDhQm/YsKHPmTPH3d0HDhxYfL4nnHCCz58/393dp0yZ4meffbbffvvtPnz48OLX6tevn69YscLd3cePH+/jxo1zd/c99tgj4bHvueceP+uss4pjX716dYmf7u7nnXeez549293d999/f9+0aZO7u69Zs8bd3X/3u9/51KlTi5e1a9fON2zY4CtWrPB+/foVv87ZZ5/tubm5/sorr3j//v0TxpNqQK4n+V6tVInAzH5dxflIpNYqqi4q8vzzz/PFF1+wefNmXn755RJDRDRr9hBPPJHLrFnX0bEjvP8+LF++jNzcXKZPn86VV17J0qVLSx3DE4wZVlTaaN68Oa1bt+aRRx5JGuOLL77IyJEj6dq1KwMGDGDdunWsX7+eefPmMXToUAD69OnD6tWr+e677zjmmGO4+uqrmTBhAmvXrqVevXolznXhwoXk5eVx5pln7vT96devH/Xr16dTp05s376dvn37AtCpU6fito5Yw4YNY/369UyaNIm//vWvxcvnzJnDAQccwMqVK3n00Ue5/PLLyzzuiy++yIgRI4pjb9q0KQCvvPIKPXr0oFOnTrz88sssXrwYgM6dOzNkyBAeeuih4n1eeOEFxo4dS9euXenVqxebNm1i2bJlHHDAAcyZMweAp59+mhYtWnDUUUft9L1Il8q2EfwX8EBVBiKSTXbffXcGDBjAHXc8xRtvnBz1Dn6R1atvpWHDV7niit0YN65oKsnQRfjggw+mV69eLFiwgB//+MclXu/AAw9k+fLlHHjggWzbto3vvvuOpk2bsn79eho1asSzzz7LscceS4sWLRiS4JaiHTt28NZbb9Ew7l7UZAlm9OjR9O/fnzlz5tCzZ09efPFF2rdvX7zNW2+9xbvvvkvbtm3Ztm0bX3/9Nb169WLu3LmlXm+3aIq0OnXqUL9+/eIEVqdOHbYVTacWY+PGjRQWFgKh+qlx48Yl1i9YsIC8vDwOOeSQ4u0POeQQ8vLySp2bxXXE2LRpE5deeim5ubm0atWKMWPGsGnTJgCeeeYZXnvtNWbPns0tt9zC4sWLcXcef/xxDjvssFJxFnnjjTeYPXs2c+bMYdOmTaxbt47zzjuPhx56KOk+1S1picDMFiV5fADsV40xitQKGzZs4IsvvgBCG8GcOXPIzW0fJYEFwG+A2fznPy146qmQBNasWcPmzZsB+Oabb3jjjTfo0KFDqdceMGAAf//73wF47LHH6NOnT4kvuebNm/Pcc8/x+9//nueffx6AiRMnMnHiRABOOeWU4t/hhyv7448/nmnTpgEwd+5cmjVrxl577cXSpUvp1KkTo0aNIicnh08++aREPJdccgkrV64kPz+fefPmceihhxYngSeeeILf/e53lX4fR40axZAhQ7j55pu5+OKLS63v378/X375Jfn5+eTn59OoUaPiJPDvf/+b888/v/icJ02aVJxsvv322+Iv/WbNmrFhwwYee+wxICTK5cuX07t3b26//XbWrl3Lhg0bOPXUU7nrrruKE+aCBQtKxfPnP/+ZwsJC8vPzmTFjBn369KlRSQDKbizeDzgf+FmCx+rUhyZSO3Tt2hUIDYsDBgygc+fOdOnShebNW7BmzYhoq+uADcDPga4UFAwA4OOPPyYnJ4cuXbrQu3dvRo8eXZwIbrrpJmbPDnM9XXTRRaxevZpDDjmEcePGMXbs2FJxHHTQQcyePZsLL7yQd955h08++YR9990XCA2hubm5dO7cmQ4dOjBp0iQgNOYWLR89enRxshk/fjwdO3akS5cuNGzYkH79+pU417IsXbq0uNG8ol599VXmz59fnAwaNGjAA1FL+mmnncbKlfGTIJa0bNmy4lLPsGHDaN26dfHnMX36dJo0acLFF19Mp06dGDhwIN27dwdCw/J5551Hp06d6NatG1dddRVNmjThxhtvZOvWrXTu3JmOHTty4403ArBy5UpOO+20Sp1jWiRrPADuA45Nsm56sv1S/VBjsdQGW7a4n3qqexgsovSjTZvUx9C/f3/fvHlz6g8UZ8iQIf71119X+3Hd3a+99lp///3303LsdKOMxuKkE9OY2bPApe7+eXUlpfLQxDRSW1x9NXz7LTz6aMkRRBs1Uu9gqXqVnZjmAeB5M7vBzOqnJjSR7JGXByedBEXXMePGwYMPhi99DREh6ZT0riF3f8TMngFuAnLNbCphXuGi9eOqIT6RjLd9O0yYEOYOaNAgDBqXE3NdNmSIvvglvXZ2++hW4HtgN6AxMYlARHbu44/hwgvh7bfh9NNh0iRo2TLdUYmUlDQRmFlfYBwwGzjS3RPMgyQiZXn88TCR/LRpMHhw8vkDRNKprDaCG4Cfu/toJQGR5GJ7BrdtC7fdBkX9pkaNCqWCc89VEpCaq6w2guOqMxCRTFQ0eXzRXT8FBaEtoE0b+PxzqF8fWrRIb4wiO6PRR0V2QbLJ47dvVwlAMocSgcguSDZ5/IoV1RuHyK5QIhCppC++SD5JvCaPl0yiRCBSQWvWhFtCDz0Urroq9ASOpcnjJdMoEYiUkzs89hgcfjj84x9wxRXwm9+oZ7BkPs1ZLFIO27bBL34BTzwBRx4Jzz0HRQNtqmewZDqVCETKoV69UO9/++3wzjs/JAGR2kCJQCQS3zFs3Dg4+WSYPz+sHz8errsuJAWR2kR/0iIk7hh2zTXQsCEsXw7R/CQitVJKSwRm1tfMlphZnpmNLmO77ma23czOSWU8Iskk6xjWtCmcdVb1xyNSnVKWCMysLnA30A/oAAw2s1KTrUbb/QV4PlWxiOxMso5hO5n5UKRWSGWJ4Gggz90/c/ctwAzgjATbXQ48DnydwlhEknrhBWjePPE6dQyTbJDKRNASWB7zvDBaVszMWgJnApPKeiEzG25muWaWu2rVqioPVLLTqlUwdCiceirsu686hkn2SmUiSDTkVvwEyeOBUe6+vawXcvfJ7p7j7jnNk126iZSTO/z979C+PcycCTfdBO+9p45hkr1SeddQIdAq5vmBQHyNaw4ww8Iwjc2A08xsm7s/mcK4JMu98AL86ldwzDHhy75D1HKljmGSrVKZCOYD7czsIGAFMAg4N3YDdz+o6HczexB4WklAUmHrVliwAI4+Gk45BZ56KkwdWUc9aURSVzXk7tuAkYS7gT4GHnH3xWY2wsxGpOq4IvHefjsMC9G7N3z9daj6GTBASUCkSEo7lLn7HGBO3LKEDcPu/qtUxiLZZ9260D/g7rvhgAPg4Yc1W5hIIromklohfniIe++Fjh1DEhg5Ej76KJQCRKQ0DTEhGS/R8BC//W1oA7j6aujRI73xidR0KhFIxks0PMTGjWGUUCUBkZ1TIpCMl2x4iGTLRaQkJQLJaH/4Q+ggloiGhxApHyUCyUhFX/577BFuC23YsOR6DQ8hUn5KBJJRli2DM84IU0YCjB4NL78MU6ZoeAiRylIikIywdSv89a9h4vgXX4TVq8Nyi0a0GjIE8vNhx47wU0lApPx0+6jUeP/+N1x8MSxaFG4JnTgxXPWLSNVQIpAaLy8vlABmzYKBA38oBYhI1VDVkNQIsT2D27SByy4LvYMBBg+GJUvgzDOVBERSQYlA0q6oZ3BBQbgbaNkyuOceGDcuPDcLdweJSGooEUjaJZs4/vvvVQIQqQ5KBJJ2yXoAL1+eeLmIVC0lAkmbFSvggQeS9wBWz2CR6qFEINVu69ZQ/9++PVx+eegUponjRdJHiUCq1bx5Ybawa66BE06ADz6AESM0cbxIOqkfgVSb1avDfMHNm8OTT4aJYmJ7BuuLXyQ9VCKQlNq+HWbPDreB7rsvPP10mC3sjDN0R5BITaFEICmTmws9e4Yv/VdfDcv69FGfAJGaRolAqkRsz+BWreCkk+Doo6GwEKZPD+0BIlIzqY1Adln8nMGFheFx6qkwcybsvXd64xORsqlEILssWc/gTz5REhDJBEoEskvWrg1jBCWiOYNFMoMSgVTKjh1w//1w6KHJt1HPYJHMoEQglTJsGFx0EbRrB3/6k3oGi2QyNRZLua1eDfXqhXr/YcOgVy8YOjT0B2jbNrQVLFsWSgK33qoOYiKZwtw93TFUSE5Ojufm5qY7jKyyfXuYHP6GG+C88+DOO9MdkYhUlJm96+45idapakjK9NZboT/AJZdAp05h7mARqV2UCAQo2SGsbdvwfOJE+OlP4auvYMYMeOUV6Ngx3ZGKSFVTIpBSU0UWFITnW7eGIaI/+QR++UuNDSRSW6mxWBJ2CNu4MbQF5OenJSQRqUYpLRGYWV8zW2JmeWY2OsH6IWa2KHq8aWZdUhmPJJas45c6hIlkh5QlAjOrC9wN9AM6AIPNrEPcZp8DJ7h7Z+AWYHKq4pHEnnsuVAclog5hItkhlSWCo4E8d//M3bcAM4AzYjdw9zfdfU309G3gwBTGIxF3+Oab8HuPHmFk0IYNS26jDmEi2SOViaAlsDzmeWG0LJmLgGcTrTCz4WaWa2a5q1atqsIQs8+iRdC7N5x4ImzbBvvsA3Pnhn4CmipSJDulMhEkusckYSWEmfUmJIJRida7+2R3z3H3nObNm1dhiNlj9Wq49FLo1g0+/DD0C4i9C2jIkNAwvGNH+KkkIJI9UnnXUCHQKub5gcDK+I3MrDNwL9DP3VenMJ6stWhRGA5i3Tq47DIYMwaaNk13VCJSU6SyRDAfaGdmB5lZA2AQMDt2AzNrDcwChrr7pymMpdZL1CGsqB3g8MPhnHNg4UKYMEFJQERKSlmJwN23mdlI4HmgLnC/uy82sxHR+knATcC+wD0W6im2JRsLQ5KLnyGsoAAuuCDMDbx8Oey1V6jzFxFJRIPO1QJt2yaeHGbvveHLL2H33as9JBGpYTToXC2XrOPXunVKAiKyc0oEGW7lyuQdv9QhTETKQ4kgQxUUwODBcMghcNVVmiFMRCpPiSDDbNgAf/gDtG8PTz4J114bpoycPFkdwkSkcjT6aAb5/vtwK2hhIZx7Lvz5zz9U/wwZoi9+EakcJYIMsGQJHHZYuB30qqvgJz8JDxGRqqCqoRokvlPY+PHwi1+EaqC33w7bXH21koCIVC2VCGqIRJ3CrroK6tcPQ0J06pTW8ESkFlMiqCESzRIG0KIF/PGP1R+PiGQPVQ3VEMk6ha0sNUyfiEjVUiJIs8WL4bTTNEuYiKSPEkGafPUVjBgBnTvDm2+GzmHqFCYi6aBEkAbucPLJcN99cPnlsHQpTJ+uTmEikh4afbSa7NgBM2fCwIFhfuDXXoP994d27dIdmYhkA40+mmavvAI5OaE38LRpYdnxxysJiEjNoERQheI7hN1xB/zsZ9CnD3z7baj+ufDCdEcpIlKS+hFUkUQdwkaPhgYN4C9/gSuu0NwAIlIzKRFUkUQdwnbsgGbN4Prr0xOTiEh5qGqoCmzfnniqSIAVK6o3FhGRilIi2EXPPANduiRfrw5hIlLTKRHsounTYfPm0AagDmEikomUCCpo6VIYNAgWLQrPJ06Ejz6CO+9UhzARyUxKBEnE3wo6aVK46m/fHv75T/jww7DdPvuEoaIhfOnn54dG4vx8JQERyQy6ayiBRLeCXnJJuNIfPjwMC73//umNUUSkqigRJJBsboD99w8lAxGR2iTrq4biq4CmTUs+N8AXX1RnZCIi1SOrSwSJqoDOP19zA4hIdsnqEkGy3sB77KFbQUUke2R1IkhWBbRxo24FFZHskRVVQ9Omhav/ZctC9c7114d7/8uqAhoyRF/8IpIdan0iSNQOcNll4Uq/Tx946y34z39+2F5VQCKSbWp91VBZt4K+9BJMmaIqIBHJbilNBGbW18yWmFmemY1OsN7MbEK0fpGZHVnVMezsVlD1BhaRbJeyRGBmdYG7gX5AB2CwmXWI26wf0C56DAf+VtVxJLvlU7eCiogEqSwRHA3kuftn7r4FmAGcEbfNGcA/PHgbaGJmVTp4w6236lZQEZGypDIRtASWxzwvjJZVdBvMbLiZ5ZpZ7qpVqyoUxJAhuhVURKQsqbxryBIsi79hszzb4O6TgckAOTk5SW76TE63goqIJJfKEkEh0Crm+YHAykpsIyIiKZTKRDAfaGdmB5lZA2AQMDtum9nA+dHdQz2B79xdQ7uJiFSjlFUNufs2MxsJPA/UBe5398VmNiJaPwmYA5wG5AEbgV+nKh4REUkspT2L3X0O4cs+dtmkmN8duCyVMYiISNlqfc9iEREpm3mykddqKDNbBRTsZLNmwDfVEE5Nlc3nr3PPXtl8/uU59zbu3jzRioxLBOVhZrnunpPuONIlm89f556d5w7Zff67eu6qGhIRyXJKBCIiWa62JoLJ6Q4gzbL5/HXu2Subz3+Xzr1WthGIiEj51dYSgYiIlJMSgYhIlsvoRFATZkBLl3Kc+5DonBeZ2Ztm1iUdcabKzs4/ZrvuZrbdzM6pzvhSqTznbma9zGyhmS02s1erO8ZUKcff/d5m9k8zez8691ozbI2Z3W9mX5vZh0nWV/77zt0z8kEYv2gpcDDQAHgf6BC3zWnAs4ThrnsC76Q77mo8958C+0S/96st517e84/Z7mXCMCfnpDvuavzsmwAfAa2j5y3SHXc1nvvvgb9EvzcHvgUapDv2Kjr/44EjgQ+TrK/0910mlwhqxAxoabLTc3f3N919TfT0bcIQ37VFeT57gMuBx4GvqzO4FCvPuZ8LzHL3ZQDuXlvOvzzn7kBjMzNgT0Ii2Fa9YaaGu79GOJ9kKv19l8mJoMpmQMtAFT2viwhXCrXFTs/fzFoCZwKTqF3K89kfCuxjZnPN7F0zO7/aokut8pz7ROBwwrwmHwC/dfcd1RNe2lX6+y6lo4+mWJXNgJaByn1eZtabkAiOTWlE1as85z8eGOXu28PFYa1RnnOvBxwFnAg0BN4ys7fd/dNUB5di5Tn3U4GFQB/gx8C/zOx1d1+X4thqgkp/32VyIsjmGdDKdV5m1hm4F+jn7qurKbbqUJ7zzwFmREmgGXCamW1z9yerJcLUKe/f/Tfu/j3wvZm9BnQBMj0RlOfcfw2M9VBpnmdmnwPtgX9XT4hpVenvu0yuGsrmGdB2eu5m1hqYBQytBVeC8XZ6/u5+kLu3dfe2wGPApbUgCUD5/u6fAo4zs3pm1gjoAXxczXGmQnnOfRmhJISZ7QccBnxWrVGmT6W/7zK2ROBZPANaOc/9JmBf4J7oqnib15KRGct5/rVSec7d3T82s+eARcAO4F53T3jLYSYp5+d+C/CgmX1AqCoZ5e61YmhqM3sY6AU0M7NC4I9Afdj17zsNMSEikuUyuWpIRESqgBKBiEiWUyIQEclySgQiIllOiUBEJMspEYhUkJm1MrPPzaxp9Hyf6HkbM3vOzNaa2dPpjlOkvJQIRCrI3ZcDfwPGRovGApPdvQC4AxiarthEKkOJQKRy/gfoaWZXEsZx+m8Ad38JWJ/GuEQqLGN7Foukk7tvNbPrgOeAU6JhkUUykkoEIpXXD/gC6JjuQER2hRKBSCWYWVfgZMJMUFfVkgmPJEspEYhUUDT71d+AK6NZwO4A/preqEQqT4lApOIuBpa5+7+i5/cA7c3sBDN7HXgUONHMCs3s1LRFKVJOGn1URCTLqUQgIpLllAhERLKcEoGISJZTIhARyXJKBCIiWU6JQEQkyykRiIhkuf8H5/bJnBw2xdkAAAAASUVORK5CYII=\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": [
    "  \n",
    "   "
   ]
  },
  {
   "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
}