function [dxdt, y] = houseModel(t,x,u,flag,varargin)


% UAfr = par.UAfr;
% UAro = par.UAro;
% Cpr = par.Cpr;
% Cpf = par.Cpf;
% Po = par.Po;
% mf = par.mf;
% Mr = par.Mr;
% R = par.R;
% k = par.k;
% Vr=par.Vr;

UAfr =    0.1801;
UAro =    0.0216;
Cpr =   1.0050 ;
Cpf =  0.6300;
Po =   1.0100;
mf =   3000 ;
Mr =  0.0290 ;
R =  8.3140e-05;
k =  100;
Vr = 60;
co2in = 6.1601e-04;
co2rate = 9.0264e-06;

% Identify inputs
qhf = u(1);         
qhr = u(2);
qin = u(3);


% To=283+dist+10*sin(2*pi*(t-6*3600)/(24*3600)-pi/2); 


hours=linspace(0,24,24);
 %data used for temperature prediction
temppred=273+[-1.2 -1.1 -2.2 -1.7 -2.9 -3.6 -3.3 -4.2 -3.4 -1.9 -2.3 -3.5 -2.2 0.9 1.1 1.1 1.0 1.4 1.2 1.9 2.9 3.5 3.1 3.3];

%data simulating measured temperature: min temp
% tempreal=273+[-1.3 -1.2 -2.2 -2.8 -2.9 -3.6 -4.2 -4.5 -4.4 -4.1 -3.5 -3.5 -3.5 -2.0 0.8 0.9 0.6 1.0 0.9  -0.2 1.9 2.9 3.1 2.9];
%data simulating measured temperature: max temp
tempreal=273+5+[-0.4 -0.5 -1.1 -1.6 -1.5 -2.3 -2.7 -3.4 -2.7 -1.9 -1.6 -1.4 -0.6 0.9 1.1 1.3 1.1 1.5 1.4 2.3 3.0 3.6 3.7 3.5];

%if prediction
if flag
    
   temperature=temppred;
   dist=1;
else
    temperature=tempreal;
    dist=1.2;
   
end



To=interp1(hours,temperature,t/3600,'linear','extrap');

price=2+1*sign(sin(2*pi*(t-8*3600)/(12*3600)+dist*pi/2));

% Identify states
Tf = x(1);                
Tr = x(2);                  
mr = x(3);                
co2=x(4);



% ODE Equations
%Floor
dxdt(1) = (qhf/(mf*Cpf))-(UAfr/(mf*Cpf))*(Tf-Tr);
%Room
dxdt(2) = (qin/mr)*(To-Tr)+(qhr/(mr*Cpr))+(UAfr/(Cpr*mr))*(Tf-Tr)-...
    (UAro/(Cpr*mr))*(Tr-To);
%Mass
%dxdt(3) = qin-k*sqrt((((mr/Mr)^2)*(R^2)*(Tr^2))-(Po^2));  
qout=k*((mr/Mr/Vr)*R*Tr-Po);
dxdt(3) = qin -qout;
co2breath=0.5*co2rate/mr;
dxdt(4)=qin*(co2in-co2)/mr+ co2breath ;
dxdt(5)=price*(qhf+qhr)/3600;
dxdt=dxdt(:);

Energy=x(5);
%Output equations (full state measurement)
y=[Tf ; Tr; mr;co2;Energy;To;price];

