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


switch flag,

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

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

  %%%%%%%%%%%%%%%%%%%%%
  % Update &terminate %
  %%%%%%%%%%%%%%%%%%%%%
  case {2,4,9}
    sys=[]; %Unused flags

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

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

end

% end sfunctionmodel


%=====================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=====================================================

function [sys,x0,str,ts]=mdlInitializeSizes


sizes = simsizes;

sizes.NumContStates  = 3;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 3;   %Number of outputs
sizes.NumInputs      = 6;   %Number of inputs
sizes.DirFeedthrough = 0;   %Put to one if input is used as output
sizes.NumSampleTimes = 1;   %at least one sample time is needed

sys = simsizes(sizes);
x0 = [200 273+90 0];
str = [];       %str is always an empty matrix
ts  = [0 0];    %initializes the array of sample times



% end mdlInitializeSizes


%=====================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=====================================================

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

%Identify inputs
Q = u(1);
qin = u(2);
qhw = u(3);
Thws = u(4);
Tin = u(5);
p = u(6);

%Identify states
V = x(1);
T = x(2);
J = x(3);

%Equations

Q_demand = par.ro.water*par.cp.water*qhw*(Thws-Tin);

p1 = 0.005;
p2 = 0.05;

if T >= Thws
pp = 0;
else
    pp = p1*(Thws-T)^2 + p2*(Thws-T);
end


%qout = qhw*((Thws-par.Tcw)/(T-par.Tin))
qins = qin/60;
qhws = qhw/60;

%Volume
dxdt(1) = qins-qhws*((Thws-Tin)/(T-Tin));
%Temperature
dxdt(2) = 1/V*((qins*(Tin-T))+Q/(par.ro.water*par.cp.water));
%Price
dxdt(3) = p*Q/3600 + pp*Q_demand/3600; %Divides with 3600 to get in hours

sys = [dxdt(1) dxdt(2) dxdt(3)];
    

%=====================================================
% mdlOutputs
% Return the block outputs.
%=====================================================
%
function sys=mdlOutputs(t,x,u,flag,par)

%States
%V = x(1);
%T = x(2);
%J = x(3)

sys = x; %[V; T];

%end mdlOutputs