function [sys,x0,str,ts] = sValve(t,x,u,flag,Input2,Input3,CV,ValveR,Model,ValveType,CurvePtsPerc,CurvePtsCV,Posflow)
% Simulink S-function for representing n paralell variable speed pumps with
% nonlinear characteristics
%
% Inputs:
% u(1) = valve opening [%]
% u(2) = Input 2 (Upstream Pressure, Downstream Pressure or Flow) [bara]
% u(3) = Input 3 (Upstream Pressure, Downstream Pressure or Flow) [bara]
%
% States:
% No states. (This is a static function.)
%
% Outputs:
% sys(1) = Output depending on Input 2 and Input 3 (the one missing....) [Sm3/s]
%
% Function arguments:
% CV            = valve coefficient [m^3/(s*%*bar^k)].
%                (multiply Cv value in [gpm/sqrt(psi)] by 2.4e-4)
% ValveR        = valve rangeability, for equal percentage valves
% Model         = 1:liquid,               q=CV*f(u)*sqrt(P_us-P_ds)
%                 2:gas,                  q=CV*f(u)*sqrt(P_us^2-P_ds^2)
%                 3:pressure independent, q=CV*f(u)
% ValveType     = 1: linear
%                 2: equal percentage
%                 3: user supplied characteristic
%                    (in CurvePtsPerc and CurvePtsCV)
% CurvePtsPerc  = user supplied characteristic curve points, % open (array)
% CurvePtsCV    = user supplied characteristic curve points, CV (array)
% Posflow       = 1: for positive flow only
%
% Note:
%  gas flow model: Pu^2-Pd^2 = 4*f*(L/D)*(R*T/M)*mdot^2/A^2 (Espen S.)
%                  --> Qstd*rhoStd = K * sqrt(Pu^2-Pd^2)
%                  --> Qstd        = (K/rhoStd) * sqrt(Pu^2-Pd^2)
%                  --> Qstd        = "CV" * sqrt(Pu^2-Pd^2)
%                      where CV is in appropriate units, here: [(Sm^3/s)/bar]
%
% History
% 2006.11.17  ESR  Fixed implementation of linear valve characteristics
%                  (i.e. conversion of % to fraction).
% 2006.12.07  ESR  Added functionality for interpolating between user
%                  supplied characteristic curve points (Cv vs %).
%                  Added description of function arguments
%

switch flag,

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

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

        %%%%%%%%%%%
        % Outputs %
        %%%%%%%%%%%
    case 3,
        sys=mdlOutputs(t,x,u,Input2,Input3,CV,ValveR,Model,ValveType,...
                       CurvePtsPerc,CurvePtsCV,Posflow);

        %%%%%%%%%%%%%%%%%%%%
        % Unexpected flags %
        %%%%%%%%%%%%%%%%%%%%
    otherwise
        error(['Unhandled flag = ',num2str(flag)]);

end

% end sfuntmpl

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

%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded.  This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;

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

sys = simsizes(sizes);

%
% initialize the initial conditions
%
x0  = [];

%
% str is always an empty matrix
%
str = [];

%
% initialize the array of sample times
%
ts  = [0 0];

% end mdlInitializeSizes


%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u,Input2,Input3,CV,ValveR,Model,ValveType,...
                        CurvePtsPerc,CurvePtsCV,Posflow)
if Input2==1;
    Pu=u(2);
    if Input3==1;
        disp('Please select 2 different inputs')
        return
    elseif Input3==2
        Pd=u(3);
        solvefor='q';
    elseif Input3==3
        q=u(3);
        solvefor='Pd';
    else
        disp('Error in input assignment')
        return
    end
elseif Input2==2
    Pd=u(2);
    if Input3==2;
        disp('Please select 2 different inputs')
        return
    elseif Input3==1
        Pu=u(3);
        solvefor='q';
    elseif Input3==3
        q=u(3);
        solvefor='Pu';
    else
        disp('Error in input assignment')
        return
    end
elseif Input2==3
    q=u(2);
    if Input3==3;
        disp('Please select 2 different inputs')
        return
    elseif Input3==1
        Pu=u(3);
        solvefor='Pd';
    elseif Input3==2
        Pd=u(3);
        solvefor='Pu';
    else
        disp('Error in input assignment')
        return
    end
else
    disp('Error in input assignment')
    return
end
if ValveType==1;
    f_av_u=u(1)/100;
    CVeffective = CV*f_av_u;
elseif ValveType==2;
    if u(1)>5
        f_av_u=ValveR^(u(1)/100-1);        % 5-100% -> =%
    else
        f_av_u=(u(1)/5)*ValveR^(5/100-1);  % 0-5% -> linear 
    end
    CVeffective = CV*f_av_u;
elseif ValveType==3,
    minPerc=min(CurvePtsPerc);
    maxPerc=max(CurvePtsPerc);
    u=max(min(u(1),maxPerc),minPerc);
    CVeffective = interp1(CurvePtsPerc,CurvePtsCV,u,'linear');
else
    disp('Error with valve type')
    return
end
if isequal(solvefor,'q')
    if Model==1 %Liquid Model
        dPcontrib=sqrt(abs(Pu-Pd))*sign(Pu-Pd);
        if Posflow;dPcontrib=max(dPcontrib,0);end
    elseif Model==2 %Gas Model
        dPcontrib=sqrt(abs(Pu^2-Pd^2))*sign(Pu^2-Pd^2);
        if Posflow;dPcontrib=max(dPcontrib,0);end
    elseif Model==3 %Pressure independent
        dPcontrib=1;
    end
    sys=CVeffective*dPcontrib;
elseif isequal(solvefor,'Pd')||isequal(solvefor,'Pu')
    dP=(q/CVeffective)^2*sign(q/CVeffective);
    if Model==1 %Liquid Model
        if isequal(solvefor,'Pd')
            sys=Pu-dP;
        elseif isequal(solvefor,'Pu')
            sys=Pd+dP;
        end
    elseif Model==2 %Gas Model
        if isequal(solvefor,'Pd')
            sys=sqrt(Pu^2-dP);
        elseif isequal(solvefor,'Pu')
            sys=sqrt(Pd^2+dP);
        end 
    else %Model==3
        disp('Cannot solve pressure from pressure independent valve model....')
        return
    end
end


