%% s-function sent into the simulink model.

% States: 
%   x1 = Temperature of the room
%   x2 = Temperature of the radiator
%   x3 = Mass of ice
%   x4 = Mass water

% Inputs: 
%   u1 = Outdoor temperature
%   u2 = splitfraction
%   u3 = Q_chiller

function [sys,x0,str,ts,simStateCompliance] = sfuntmpl(t,x,u,flag,a)

switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u,a);

  %%%%%%%%%%%%%%%%%%%%%%
  % Update & Terminate %
  %%%%%%%%%%%%%%%%%%%%%%
  case {2,9},
    sys=[]; %do nothing

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u);

  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag',num2str(flag));

end

% end sfuntmpl

%
%================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, 
% and sample times for the S-function.
%================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes

%
% call simsizes for a sizes structure, fill it in and convert
%it to a sizes array.
%

sizes = simsizes;

sizes.NumContStates  = 4;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 4;
sizes.NumInputs      = 3;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed

sys = simsizes(sizes);

%
% initialize the initial conditions for the four states, x1,x2,
% x3 and 4x.
%
x0  = [ 296
        288
        50
        80];

str = [];
ts  = [0 0]; %sample time
simStateCompliance = 'UnknownSimState';

% end mdlInitializeSizes

%
%================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%================================================================
%
% The equality constraints of the system:

function sys=mdlDerivatives(t,x,u,a)

sys = [a.vent*(u(1)-x(1))/a.room - a.U*a.A*(x(1)-x(2))/(a.room*a.cpair*a.densair);
        a.massys*u(2)*(a.Tw-x(2))/(a.radiator*a.denswater) + a.U*a.A*(x(1)-x(2))/(a.radiator*a.cpwater*a.denswater);
        u(3)/a.fus - a.massys*a.cpwater*u(2)*(x(2)-a.Tw)/a.fus];
        sys(4) = -sys(3);
%
% 
% end mdlDerivatives
%
%================================================================
% mdlOutputs
% Return the block outputs.
%================================================================
%
function sys=mdlOutputs(t,x,u)

sys = x;

% end mdlOutputs

