clc 
clear all

Ufl =   10e-3;  %kW/(m^2 °C)
Afl =   50;     %m^2
mfl =   600;    %kg
Cpfl=   1;      %kJ/(kg°C)
Cp=     1;      %kJ/(kg°C)
m=      180;    %kg
mv=     0.2;    %kg/s
Aw=     100;    %m^2
Uw=     2e-3;   %kW/(m^2 °C)
alpha=  0.4;     %kW/(°C)
Tin  = 21;     %C

syms Qfl QEH QAC QHW T Tfl Q

u= [QAC; QHW; QEH];
x= [T; Tfl]; 

Qfl = Ufl*Afl*(Tfl-T);
Q   = Qfl+QEH-QAC;

dTdt    = (alpha*(Tin-T)+Q)/(m*Cp);    %rate of change of room temperature
dTfldt  = (QHW-Qfl)/(mfl*Cpfl);        %rate of change of floor temperature

f   = [(alpha*(Tin-T)+Q)/(m*Cp)
        (QHW-Qfl)/(mfl*Cpfl)]; 

%dxdt = [dTdt;dTfldt];

A=jacobian(f,x);

B=jacobian(f,u);

% operating 
Qfl = 0; 
QEH= 0; 
QAC = 0; 
QHW = 0; 
T= 21;
Tfl= 21;
Q=0; 

A=double(subs(A))
B=double(subs(B))
C=eye(2);
s=tf('s');

G=minreal(C*inv(eye(2)*s-A)*[B]);

%Gd=minreal(C*inv(eye(2,2)*s-A)*E);

G2=set(zpk(G),'DisplayFormat','time constant')
%Gd2=set(zpk(Gd),'DisplayFormat','time constant')
