function botijo  % botijo en una funcion en vez de un script
% Este fichero m de MATLAB: botijo.m resuelve un modelo dado por dos
% ecuaciones diferenciales (balance de materia y energia de un botijo) 
% y estima los valores optimos de los parametros en las ecuaciones 
% diferenciales para que se haga minima la suma de las desviaciones
% al cuadrado entre los valores experimantales y los calculados,
% resolviendo reiteradamente y numericamente el modelo.

% Ultima edicion 5/8/2005

%   Copyright (C) 1994, 2005 Jose Ignacio Zubizarreta
%   Prohibida la reproduccion total o parcial sin permiso del autor

disp '        '
disp 'Copyright (C) 1994,2005 Jose Ignacio Zubizarreta'
disp '        '
%
global k u lamb y0 tiempo vol temp resid
% 
% Los datos experimentales que me dio Pinto estan en la matriz data: tiempo
% en minutos -delta peso reduccion en kg desde 3,161 -delta theta en grados C
% disminucion de temperatura desde 39 ºC.

data =[0	0	0.1
15	20	2
27	35	3.8
36	48	5.2
60	75	8
73	90	9.4
85	103	10.2
100	118	10.8
115	134	11.7
132	150	12.2
170	184	13.2
191	203	13.5
211	222	13.8
276	281	14.3
355	357	14.6
432	429	14.8
1596	1496	13.7		
1630	1519	13.7];

% Cambio los datos para que aparezca el tiempo en horas, el volumen en
% litros en el botijo multiplicado por 10 para representar las dos curvas
% en un solo grafico y la temperatura real del agua en grados C.
tiempo = data(:,1)./60;
vol = 31.61*ones(18,1)-data(:,2)./100;
temp = 39*ones(18,1)-data(:,3);

% Valores iniciales de h (en dm de altura del segmento esferico seco;  
% 0,64112 dm de altura de segmento corresponden a un volumen de 3,161
% litros) y theta (temperatura inicial del agua del botijo)
y0=[0.64112 38.9]';

%Estimacion de los valores de los parametros
k=0.8;
u=0.22;
lamb=1.02;
estima = [k
         u
         lamb];

options = [];
opt = fminsearch(@f2, [estima]);
k = opt(1)
u = opt(2)
lamb = opt(3)
residuals = [tiempo resid];

% Solucion del modelo una vez optimizados los parametros
t0 = 0; tf=72;
[t, Y]=ode45(@f1, [t0, tf], y0);
%
[n m]=size(Y(:,1));

% Determino el volumen remanente en el botijo en funcion de h=Y(:,1)
Z=10*(4/3*pi*ones(n,1)-pi.*Y(:,1).*Y(:,1)+pi/3.*Y(:,1).*Y(:,1).*Y(:,1));

str1 = ['k = ' num2str(100*k,4) ' kg/h m2'];
str2 = ['U = ' num2str(100*u,4) ' kcal/h m2'];
str3 = ['lamb = ' num2str(100*lamb,5) ' kcal/h m2'];

% Se construye la grafica de los resultados comparando la solucion del
% modelo y los datos experimentales
plot(t,Z,t,Y(:,2),tiempo,vol,'bo',tiempo,temp,'gx')
ylabel('Volumen en litros x 10 y temperatura en ºC del agua del botijo')
xlabel('Tiempo en horas');
title('Evolucion del volumen y temperatura frente al tiempo para el valor optimo de los parametros')
legend('Volumen','temperatura', 'datos experimentales Vol','datos experimentales ºC', ...
    'Location','north')
text(40,22,str1)
text(40,20,str2)
text(40,18,str3)

dummy = input('Pulse <enter> para ver una grafica con las desviaciones. ');

figure, plot(tiempo,resid,'d')
ylabel('Residuos')
xlabel('Tiempo en horas')
title('Desviaciones con los parametros optimos')
grid

% Esta funcion presenta las ecuaciones diferenciales
function dydt = f1(t,y)
global k u lamb
dydt(1)=-0.007*k*(4*pi-pi*y(1)*y(1))/(pi*y(1)*y(1)-2*pi*y(1));
dydt(2)=(-0.233*k*(4*pi-pi*y(1)*y(1)) + lamb*2*pi*y(1)-u*(4*pi-pi*y(1)*y(1))*(y(2)-24.2))/(4/3*pi-pi*y(1)*y(1)+pi/3*y(1)*y(1)*y(1));
dydt=dydt';

%Esta funcion calcula la suma de cuadrados de las desviaciones de los datos
%experimentales y los calculados, asi como calcula los residuos o
%diferencias de datos experimentales y valores del modelo
function sqerror = f2(x)
global k u lamb y0 tiempo vol temp resid
k = x(1);
u = x(2);
lamb = x(3);
[t Y] = ode45(@f1, tiempo, y0);
Z=10*(4/3*pi*ones(18,1)-pi.*Y(:,1).*Y(:,1)+pi/3.*Y(:,1).*Y(:,1).*Y(:,1));
diff1 = vol-Z;
diff2 = temp-Y(:,2);
sqerror = diff1'*diff1 + diff2'*diff2;
resid = [diff1 diff2];
