function u = CalcBlock_num(x,v,d,p)
% static calculation block to solve for u = f(v,x,d,p)

qnom = 1e-2;                                     % [m3/s] nominal flow rate
Qnom = 1672;                                      % [kW] nominal heat rate

MV     = p(6);                   % switch between MVs; 1:u=q; 2=u=Q

switch MV
    case 1
        u0=qnom;
    case 2
        u0=Qnom;
    otherwise
end
opts = optimoptions('fsolve','Display','off');
u    = fsolve(@(u) rootu(x,v,p,d,u),u0, opts);          % find u numericaly

    function F = rootu(x,v,p,d,u)
        V    = p(1);                                     % [m3] tank volume
        rho  = p(2);                                     % [m3/kg] density
        cp   = p(3);                              % [kJ/kg/C] specific heat
        A    = p(4);             %[1/s] tuning param gen linear transf (vL)
        tau0 = p(5);                         % [s] tuning param IFAC. (vL0)
        
        T    = x;                                    % [C] tank temperature
        T0   = d(1);                               % inlet flow temperature
        
        
        % switch variables to change between case
        MV     = p(6);                   % switch between MVs; 1:u=q; 2=u=Q
        transf = p(7);   % switch trasform: 0=static; 1=lin;2=vL0; 3=fb lin
        
        switch MV                                  % select if q or Q is MV
            case 1                                       % q is MV, Q is DV
                q = u; % [m3/s] inflow input applied to the plant. after saturation
                Q = d(2);                             % [kW] heat flow rate
                f = q/V*(T0-T) + Q/(V*rho*cp);      % RHS of energy balance
                
                switch transf                  % select transformation type
                    case 0                          % static transformation
                        F = v - (T0+Q/(q*rho*cp));                  % u = q
                    case 1                          % linear trasnformation
                        F = v - (f-A*T);                            % u = q
                    case 2                                      % IFAC. vL0
                        F = v - (tau0*f+T);                           % u=q
                    case 3  %integrating transform. similar to feedback lin
                        F = v - f;                                  % u = q
                    otherwise
                end
                
            case 2                                       % q is DV, Q is MV
                Q = u; %[kW] heat rate input applied to the plant. after saturation
                q = d(2);                                % [m3/s] flow rate
                f = q/V*(T0-T) + Q/(V*rho*cp);      % RHS of energy balance
                switch transf
                    case 0                          % static transformation
                        F = v - (T0+Q/(q*rho*cp));                    % u=Q
                    case 1                          % linear trasnformation
                        F = v - (f-A*T);                            % u = Q
                    case 2                                      % IFAC. vL0
                        F = v - (tau0*f+T);                           % u=Q
                    case 3  %integrating transform. similar to feedback lin
                        F = v - f ;                                 % u = Q
                    otherwise
                end
            otherwise
        end
    end
end
