function [sys,x0,str,ts] = advancedmdl (t,x,u,flag,init,par)

switch flag,
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes(t,x,u,flag,init,par);
  case 1,
    sys=mdlDerivatives(t,x,u,flag,init,par);
  case 3,
    sys=mdlOutputs(t,x,u,flag,init,par);
  case { 2, 4, 9 }
     sys = []; % Unused flags
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);
end

function [sys,x0,str,ts]=mdlInitializeSizes(t,x,u,flag,init,par)

sizes = simsizes;
sizes.NumContStates  = 5;   
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 5;   %number of outputs
sizes.NumInputs      = 5;   %number of inputs
sizes.DirFeedthrough = 0; %Put to 1 if input is used in output
sizes.NumSampleTimes = 1; % At least one sample time is needed
sys = simsizes(sizes);

% %Initial conditions
x0  = [init];
% str is always an empty matrix
str = [];
% initialize the array of sample times
ts  = [0 0];

function sys=mdlDerivatives(t,x,u,flag,init,par)


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;

co2in=par.co2in;
co2rate=par.co2rate;

% Identify inputs
qhf = u(1);         
qhr = u(2);
qin = u(3);
To = u(4);
price=u(5);
% Identify states
Tf = x(1);                
Tr = x(2);                  
mr = x(3);                
co2=x(4);

hours=0:23;

 %data used for temperature prediction
tempmeas=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];
To=interp1(hours,tempmeas,t/3600,'linear','extrap');

% 


% 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/mr-max(0,qout)*co2/mr+max(0,-qout)*co2in/mr  + co2breath ;
dxdt(4)=qin*(co2in-co2)/mr+ co2breath ;
dxdt(5)=price*(qhf+qhr)/3600;
dxdt=dxdt(:);

%       Floor   Room    Mass
% sys = [dxdt(1) dxdt(2) dxdt(3)];
sys=dxdt;

function sys=mdlOutputs(t,x,u,flag,init,par)

% States
Tf = x(1)-273; % C                
Tr = x(2)-273; % C   
mr = x(3); 

sys = [Tf;Tr;mr;x(4);x(5)];