// Autor: Carlos F. Pfeiffer
// Carlos.Pfeiffer@gmail.com

// Para usar este archivo, cambia la extensión de .txt a .sci,
// y ejecútalo desde Scilab.



//=============  Ecuación de Estado Cúbica Genérica   ==============
R = 8.31451; // Constante de los gases: Pa m3 / (mol K)

//======================= Constantes de la substancia =================
//--------------------------- n-butano --------------------------------
Pc = 37.96; // bars
Tc = 425.1; // K
w = 0.2;    // Factor de correlación de Pitzer
//==================== Condiciones de operación =======================
P = 9.4573; // bar
T = 350.0;    // K
//=====================================================================
// --------------- Cálculo de propiedades reducidas -------------------
Pr = P/Pc;
Tr = T/Tc;

// ========== Función Genérica para el Factor de Compresibilidad Z ================
function Z = CubicaGenerica(Pr,Tr,w,opcion)
// Pr = Presion reducida, Tr = Temperatura reducida, w = factor acéntrico de Pitzer
// opcion = "vdW" : vand der Waals
// opcion = "RK"  : Redlich-Kwong
// opcion = "SRK" : Soave-Redlich-Kwong
// opcion = "PR"  : Peng-Robinson

  valido = %T;   // Se supone que la opción es válida. %T representa la constante booleana "verdadero".
  if opcion == "vdW"  then // van der Waals
    sigma = 0.0;
    epsilon = 0.0;
    Omega = 1.0/8.0;
    Psi = 27.0/64.0;
    alfa = 1;
  elseif opcion == "RK" then // Redlich-Kwong
    sigma = 1.0;
    epsilon = 0.0;
    Omega = 0.08664;
    Psi = 0.42748;
    alfa = 1/sqrt(Tr)
  elseif opcion == "SRK" then // Soave-Redlich-Kwong
    sigma = 1.0;
    epsilon = 0.0;
    Omega = 0.08664;
    Psi = 0.42748;
    alfa = (1 + (0.480 + 1.57*w - 0.176*w*w)*(1 - sqrt(Tr)))^2;
  elseif opcion == "PR" then // Peng-Robinson
    sigma = 1.0 + sqrt(2);
    epsilon = 1.0 - sqrt(2);
    Omega = 0.07780;
    Psi = 0.45724;
    alfa = (1 + (0.37464 + 1.5422*w - 0.26992*w*w)*(1 - sqrt(Tr)))^2;
  else
  valido = %F;  // La opción no es válida! 
  end
       
   if valido == %T then
    // ---------- Cálculo de los parámetros para la ecuación cúbica ----------
    Beta = Omega * Pr/Tr;
    q = (Psi * alfa)/(Omega*Tr);

    // ---------- Coeficientes del polinomio cúbico en Z ---------------------
    a0 = -sigma*epsilon*(Beta^3 + Beta^2) - q*Beta^2;
    a1 = (Beta^2*epsilon-Beta^2-Beta)*sigma + Beta*q - (Beta^2 + Beta)*epsilon;
    a2 = (epsilon + sigma -1)*Beta - 1;
    a3 = 1; 
    RK = poly([a0,a1,a2,a3],'Z','coeff')  // Polinomio cúbico en Z

    //------------- Obtención de las raíces del polinomio -----------------
    Z = roots(RK) // Raíces del polinomio cúbico en Z
  else
    Z = 0;
    printf("OPCION DE ECUACIÓN NO VÁLIDA \n")
  end

endfunction

//=============================================================================================
// Prueba de la función con el ejemplo 3.9 del libro de Smith, Van Ness y Abbot.
printf(" Volumen específico del gas @ 350 K y 9.4573 Pa en m3/mol \n")
printf("van der Waals \n");
Z = CubicaGenerica(Pr,Tr,w,"vdW");
V = Z(3)*R*T/(P*100000)


printf("Redlich-Kwong \n");
Z = CubicaGenerica(Pr,Tr,w,"RK");
V = Z(3)*R*T/(P*100000)

printf("Soave-Redlich-Kwong \n");
Z = CubicaGenerica(Pr,Tr,w,"SRK");
V = Z(3)*R*T/(P*100000)

printf("Peng-Robinson \n");
Z = CubicaGenerica(Pr,Tr,w,"PR");
V = Z(3)*R*T/(P*100000)


printf(" Volumen específico del líquido @ 350 K y 9.4573 Pa en m3/mol \n")
printf("van der Waals \n");
Z = CubicaGenerica(Pr,Tr,w,"vdW");
V = Z(1)*R*T/(P*100000)


printf("Redlich-Kwong \n");
Z = CubicaGenerica(Pr,Tr,w,"RK");
V = Z(1)*R*T/(P*100000)

printf("Soave-Redlich-Kwong \n");
Z = CubicaGenerica(Pr,Tr,w,"SRK");
V = Z(1)*R*T/(P*100000)

printf("Peng-Robinson \n");
Z = CubicaGenerica(Pr,Tr,w,"PR");
V = Z(1)*R*T/(P*100000)