% Tanque mezclador con control de Temperatura y de Nivel
% En X están las variables de estado.
% En Y deben ir las variables que se requieren en las ODEs o que se quieren graficar.

clear all; close all; clc;

%=============== Modelo =================

% ODEs
function dX = ODEs(t,X)
  % En dX devuelve el vector columna de derivadas

  % Recupera variables X
  [L T Ai1] = num2cell(X'){1,:};

  % Recupera variables Y
  Y = AEs(t,X);
  [F1 F2 F3 A T1 T2 e1 Lsp Tsp x1 x3 rho Cp Wg lambda] = num2cell(Y){1,:};

  % Ecuaciones diferenciales
  dL=(F1+F2-F3)/A;
  dT=(F1*rho*Cp*(T1-T)+F2*rho*Cp*(T2-T)+Wg*lambda)/(rho*Cp*A*L);
  dAi1=e1;

  dX = [dL dT dAi1]'; % vector columna
endfunction % ODEs
%---------------------------------------

% AEs
function Y = AEs(t,X)
  % En Y devuelve el vector fila de variables requeridas por ODEs o a graficar.

  % Recupera variables X
  [L T Ai1] = num2cell(X'){1,:};

  % Parámetros
  Wg=1; %kg/s
  T1=25; T2=100; %°C
  D=1.5; %m
  Pv1=10130; %N/m2 (Pa)
  g=9.81; %m/s2

  Ab1=0.5;
  Kp1=0.05; %°C
  tau1=30; %s
  Tsp=52; %°C

  Ab3=0.5;
  Kp3=10; %m-1
  Lsp=1; %m

  Cv1=4.039E-4; Cv3=8.078E-4; %m3.5/kg0.5

  rho=1000; %kg/m3
  Cp=1; %kcal//kg*°C
  lambda=539.4; %kcal/kg


  % Ecuaciones algebraicas
   %Perturbación
   if t<100;
    Tsp=Tsp;
    else
    Tsp=Tsp+10;
    endif

    A=pi*D^2/4;
    F2=Wg/rho;

    %Control de temperatura
    e1=T-Tsp;
    Ac1=Ab1+Kp1*(e1+Ai1/tau1);
    x1=max(0,Ac1); %Acota entre 0 y 1
    x1=min(x1,1);
    F1=Cv1*x1*sqrt(Pv1);

    %Control de nivel
    e3=L-Lsp;
    Ac3=Ab3+Kp3*e3;
    x3=max(0,Ac3); %Acota entre 0 y 1
    x3=min(x3,1);
    F3=Cv3*x3*sqrt(rho*g*L);

    Y=[F1 F2 F3 A T1 T2 e1 Lsp Tsp x1 x3 rho Cp Wg lambda];

endfunction % AEs
%---------------------------------------

% Inicialización
function [tfin dt Xini LX LY] = inicializacion
  % Inicializa la simulación

  % Parámetros de simulación
  tfin = 1100; % tiempo final
  dt = 1; % paso temporal

  % Inicialización
  Lini=1;
  Tini=52;
  Ai1ini = 0;
  Xini = [Lini Tini Ai1ini]; % Inicializa la variable de estado

  % Leyendas
  LX = {"L" "T" "Ai1"}; % Leyendas de las variables X
  LY = {"F1" "F2" "F3" "A" "T1" "T2" "e1" "Lsp" "Tsp" "x1" "x3" "rho" "Cp" "Wg" "lambda"}; % Leyendas de las variables Y
endfunction % inicializar
%---------------------------------------

% Análisis
function analizar(LX,LY,tpts,X,Y)
  % Análisis de resultados. Funciones disponibles:
  % exportar('resultados.csv',[{leyendas}])
  % graficar({leyendas}, 'título', 'rótulo x', 'rótulo y', [limitesy])
  % vector(leyenda)
  exportar('TP5_Ej2.csv', {"F1" "T1" "x1" "F2" "T2" "F3" "T" "x3"  "L" "Ai1"});
  graficar({"Ai1"}, 'Efecto integral vs. tiempo', 's', '°C*s');
  graficar({"F1" "F2" "F3"}, 'Caudales vs. tiempo', 's', 'm^3/s');
  graficar({"T1" "T2" "T" "Tsp"}, 'Temperaturas vs. tiempo', 's', '°C');
  graficar({"x1" "x3"}, 'Apertura de válvulas vs. tiempo', 's', ' ');
  graficar({"L" "Lsp"}, 'Nivel vs. tiempo', 's', 'm');
endfunction % analizar
%=======================================

% Llama al resolverdor
resolvedor;