%% Model of the plant

function xend = ModelNMPC(x0,D,Qin,Betain,dt, P2, P3) 
    
addpath('C:\Users\subsea\Desktop\Amund\casadi-windows-matlabR2016a-v3.5.5')
    import casadi.*

    %% Parameters
    
    par = initHC();
    Vo = par.VRF();
    Vf = par.Vol_HC-par.VRF;
    Zu = 0.5;       % Constant underflow valve opening
   
    %% Model variables
   
    x1 = SX.sym('x1');
    x2 = SX.sym('x2');
    x = [x1;x2];
    
    %% Model inputs
    
    u = SX.sym('u'); % Overflow valve opening
     
    %% Model equations
    
    cv1 = pi*0.005^2;
    cv2 = pi*0.001^2;
    
    Qu = cv1*Zu*sqrt(2/par.Rho_w*(P3-par.Patm));
    Qo = cv2*u*sqrt(2/par.Rho_o*(P2-par.Patm));
   
    % Qsep polynomial approximation
    
    p2 =  -9.447e+07;
    p1 =        9024;
    p0 =      0.7648; 
    
    % Oil inflow
    Qino = Betain.*Qin;
    Qsep = Qino.*(p2*Qo.^2+p1*Qo+p0);
    
    %% Model equations
    
    xdot1 = (1/Vo)*(Qsep-x1*Qo); 
    xdot2 = (1/Vf)*(Qino-Qsep-x2*Qu);
    
    xdot = [xdot1;xdot2];
   
    ode = struct('x',x,'p',u,'ode',xdot);
    
    %% Integrator
    
    opts = struct('tf',dt);
    F = integrator('F', 'cvodes', ode, opts);
    sim = F('x0',x0,'p',D);

    xend = full(sim.xf);
end
