clc; clear; 
// Ejercicio 1.4

//--- Obtengo datos experimentales desde archivo de texto ---
Beta=2
Datos=read('Ej_1.4_TablaDatos.txt',22,2) // para realizar esta opreración es necesario estar trabajando en la carpeta donde se encuentra el archivo Ej_1.4_TablaDatos.txt. Verificarlo en File Browser en Scilab

T = Datos(:,1)' // Vector de temperaturas experimentales
P = Datos(:,2)' // Vector de presiones experimentales

// Otra forma es ingresar los datos manualmente:
//T = [290.1 302.4 311.2 318.7 325.1 330.5 334.9 338.9 343.0 346.2 349.9 353.5 356.2 358.9 362.3 365.2 367.9 370.5 373.2 375.8 378.5 381.3]  //K
//P = [8.6 15.4 22.5 30.5 39 47.6 55.5 63.8 73 81.3 91.3 102.0 110.9 120.1 132.8 144.5 155.8 167.6 180 193.5 207.9 223.4] //kPa

//------------------------------ Parte a -----------------------------------------
// --- Clapeyron log10(P) = A + B/T ---
//--- Método de búsqueda directa ---
// sistema lineal: variables x=1/T  y=log10(P) ; y = B0 + B1x
// minimos cuadrados (regresion lineal): B0 = Ym - B1 Xm; B1 = (n sum(XY)- sum(Y)*sum(X))/(n sum(X^2) - sum(X)^2)
X = T .^ (-1); Xm = mean(X);
Y = log10(P) ; Ym = mean(Y);
n = length(T)
B1 = (n*sum(X.*Y)-sum(Y)*sum(X))/(n*sum(X.^2)-sum(X)^2);
B0 = Ym - B1*Xm;
disp('Coeficientes de Clapeyron hallados por método de regresión lineal:')
disp(B1,'B = ',B0,'A = ')
Yest = B0 + B1 * X;
Pest_C_RL = 10.^Yest  // presión modelada (predicciones del modelo de Clapeyron hallando coeficientes por método de regresión lineal)

//--- Método de busqueda directa: fminsearch ---
function e=clapeyron(B)
    y = B(1) + B(2)*T.^(-1); Pest = 10.^y  // calculo de valores de presión teóricos
    e = sum((Pest - P).^2)  // minimos cuadrados
endfunction
Binicial = [B0,B1]; // estimación inicial de parámetros
B = fminsearch(clapeyron,Binicial);
disp('Coeficientes de Clapeyron hallados por método de búsqueda directa:')
disp(B(2),'B = ',B(1),'A = ')
Pest_C_BD = 10.^( B(1) + B(2)./ T); // presión modelada (predicciones del modelo de Clapeyron ajustando los coeficientes por método de búsqueda directa)
 
//--- Riedel log10(P) = A + B/T + C log10(T) + D T^Beta ---
// Método de busqueda directa: fminsearch
function e=riedel(B)
    y = B(1) + B(2)*T.^(-1) + B(3)*log10(T) + B(4)*T.^(Beta)
    Pest = 10.^y
    e = sum((Pest - P).^2)
endfunction
Binicial = [B0,B1,0,0];
B = fminsearch(riedel,Binicial);
disp('Coeficientes de Riedel hallados por método de búsqueda directa:')
disp(B(4),'D = ',B(3),'C = ',B(2),'B = ',B(1),'A = ')
Pest_R_BD = 10.^( B(1) + B(2)./ T + B(3)*log10(T) + B(4)*T.^(Beta));   // presión modelada (predicciones del modelo de Riedel ajustando los coeficientes por método de búsqueda directa)

//------------------------------ Parte b -----------------------------------------
clf // clear figure
plot(T,P,'x',T,Pest_C_RL,'g',T,Pest_C_BD,'r',T,Pest_R_BD,'m');
title('Datos experimentales y  predicciones de modelos')
legend('Experimental','Clapeyron RL','Clapeyron BD','Riedel BD')
xlabel('Temperatura (K)'); ylabel('Presión (kPa)')


//------------------------------ Parte c -----------------------------------------
figure(2); clf; title('Residuos')
plot(T,(P-Pest_C_RL),'g')
plot(T,(P-Pest_C_BD),'r')
plot(T,(P-Pest_R_BD),'m')
legend('Clapeyron RL','Clapeyron BD','Riedel_BD')
xlabel('Temperatura (K)'); ylabel('Residuos (kPa)')

//------------------------------ Parte d -----------------------------------------
//Mejor ajuste con Riedel, pero no varía demasiado respecto a Clapeyron cuando se hallan los valores de los parámetros con un método de búsqueda directa. No vale la pena introducir dos parámetros más para obtener una mejora tan poco significativa.