clc
clear all
close all

datos=load('DATOS_TEMPORALES_CALENTAMIENTO.txt');

tiempo=datos(:,1);
potencia=datos(:,2);
temperatura=datos(:,3);
nMuestras=size(temperatura,1);

figure(1)
title('Datos experimentales')
subplot(2,1,1)
plot(potencia,'r','linewidth',2)
title ('potencia')
axis([0,3500,-10,60])
subplot(2,1,2)
plot(temperatura,'b','linewidth',2)
title ('temperatura')
axis([0,3500,10,100])

instante=1;
while(potencia(instante)==potencia(instante+1))
    instante=instante+1;
end
disp('El escalon se aplica en el instante:')
instanteEscalon=instante

%% Identificación del sistema con un modelo de 1er orden en el que solo hay que calcular la ganancia y la constante de tiempo
cteTiempo=0.63*(temperatura(nMuestras)-temperatura(instanteEscalon));
while temperatura(instante)<(cteTiempo+temperatura(instanteEscalon))
    instante=instante+1;
end
disp('Parámetros del sistema:')
cteTiempo=instante-instanteEscalon
amplitudEscalon=(potencia(nMuestras)-potencia(instanteEscalon));
ganancia=(temperatura(nMuestras)-temperatura(instanteEscalon))/amplitudEscalon

modelo=tf(ganancia,[cteTiempo,1])
% la funsión setpDataOptions solo esta disponible para Matlab 2013b
opciones = stepDataOptions('InputOffset',0,'StepAmplitude',amplitudEscalon);
[y,t] = step(modelo,opciones);

% con esta figura se comprueba que la respuesta no coincide con un modelo
% de primer orden
figure(2)
hold off
plot(t,y+temperatura(instanteEscalon),'linewidth',2)
hold on
plot([0:(nMuestras-instanteEscalon)],temperatura(instanteEscalon:nMuestras),'r','linewidth',2)
title('Primer orden')