
function  [sys,x0,str,ts] = system(t,x,u,flag)
% This S-function is based on a genral example written by Vinicius de
% Oliveira.
% S-functions store all the information behind a simulink block and 
% therefore behaves differently appearing on the flag sent from 
% Simulink. At first the function will be called to initialize the number 
% of outputs/inputs/states that the block should have.(flag = 0)
% Then Simulink will call the function with different flags to find the
% derivatives and the updataed states. (flag = 1 and 3 respectively)
switch flag
    
    % Initialization
    
    case 0
        global par   % Load par for updated initial conditions
        sys = [3,    % Number of continuous states
            0,       % Number of discrete states
            3,       % Number of outputs--NUMBER OF CVS
            6,       % Number of inputs (to the block)
            0,       % reserved must be zero
            0,       % direct feedthrough flag
            1];      % number of sample times
        
        %INITIAL CONDITIONS FOR THE STATES      
        x0  = [par.V, par.T, 0];
        str = [];
        ts  = [0 0];                      % sample time: [period, offset]
 
    case 1
        global par
        % Model parameters
        c_p    = par.c_p;                         %heat capacity of water
        rho    = par.rho;                               %density of water
        Q_loss = par.Q_loss;                  %loss of heat from the tank

        
        % Extract information from Simulink model
        % inputs
        Q     = u(1);    % kW
        q_in  = u(2)/60; % l/min to l/s
        q_hw  = u(3)/60; % l/min to l/s
        T_hws = u(4);    % degrees celcius
        T_cw  = u(5);    % degrees celcius
        p     = u(6)/1000;%EUR/MWh to EUR/kWh
        % States
        V = x(1);
        T = x(2);
        J = x(3); % not used in any of the formulas

        % Differential equations
        dxdt(1) = q_in - q_hw/(1+( (T-T_hws)/(T_hws-T_cw) ) );      %[l/s]
        dxdt(2) = (1/V)*(q_in*(T_cw-T) + (Q-Q_loss)/(rho*c_p));     %[C/s]
        dxdt(3) = (p*Q)/3.6e3;                                    %[EUR/s]
        % divided by 3600 to transform from [EUR/h] to [EUR/s]
        sys = dxdt;
        
    case 2                                         % Discrete state update
        
        sys = []; % do nothing
        
    case 3
        
        % --------OUTPUTS------------
          
          %all the states are outputs in this case - it could be subset 
          sys=x;
        
    case 9                                                % Terminate
        
        sys = []; % do nothing
        
    otherwise
        error(['unhandled flag = ',num2str(flag)]);
end

