function [A,C,Bu,Bv,Cm,Dv,Dvm,QQ,RR,NN] = mpc_plantupdate(...
    a,b,c,d,A,B,C,D,mvindex,mdindex,unindex,nxp,nup,ny,nmv,nv,nxid, ...
    myindex,Uscale,Yscale,Cid,Did)
% Helper function to update prediction model and nominal conditions 
% for adaptive mode.  Assumes all dimensions have been validated.

% a,b,c,d = new plant matrices in engineering units
% A,B,C,D = the MPCstruct matrices to be updated.
% mvindex, mdindex,unindex,myindex = pointers to locations in the
% original plant model.  mdindex DOES NOT include the offset.
% nv = length(mdindex) + 1 (includes offset).  NOTE:  EML doesn't like
% empty matrices so an empty mdindex is replaced by zero.  Thus, need to
% provide nv separately.  (Same for nxid).
% nxid = number of input disturbance states.
% Uscale, Yscale are plant scale factors.
% DX is updated nominal DX value, discrete time
% Cid, Did are coeffient matrices of the input disturbance model
% Bu,Bv,Cm,Dv,Dvm are subsets of the full model parameters
% QQ,RR,NN are updated covariance matrices for use in LTVKF calculations.

% Note:  only the plant model changes but the augmented [A,B,C,D] model
% includes  terms that depend on both the plant and input disturbance
% model.  If the plant includes unmeasured input disturbances, these terms
% must be updated.

% Authors:  A. Bemporad, N.L. Ricker, R. Chen
% Copyright 1994-2020 The MathWorks, Inc.

%#codegen
ONE32 = ones('int32');
ZERO32 = zeros('int32');
nym = int32(length(myindex));
nx = int32(size(B,1));
nutot = int32(size(B,2));
% Convert b,c,d to dimensionless
b = b.*Uscale(:,ones(ONE32,nxp))';
c = c./Yscale(:,ones(ONE32,nxp));
d = Uscale(:,ones(ONE32,ny))'.*(d./Yscale(:,ones(ONE32,nup)));
% Update the [A,B,C,D] model
ixp = ONE32:nxp;
A(ixp,ixp) = a;

% These are the terms involving both the plant and input disturbance
% model.
if ~isempty(unindex) && all(unindex~=ZERO32)
    nuid = int32(size(Did,2));   % Number of noise inputs to input disturbance model
    iid = ONE32:nuid;
    for ct = ONE32:nuid
        iid(ct) = iid(ct) + (nmv+nv);   % Points to noise inputs to ID model
    end
    if nxid > ZERO32
        ixid = ONE32:nxid;      % Points to input disturbance model states    
        for ct = ONE32:nxid
            ixid(ct) = ixid(ct) + nxp;   % Points to noise inputs to ID model
        end
        A(ixp,ixid)  = b(:,unindex)*Cid;
        C(:,ixid)    = d(:,unindex)*Cid;
    end
    B(ixp,iid)   = b(:,unindex)*Did;
    D(:,iid)     = d(:,unindex)*Did;
end
B(ixp,ONE32:nmv) = b(:,mvindex);
C(:,ixp)    = c;
if nv > ONE32
    % Measured disturbance terms
    imd = ONE32:(nv-ONE32);
    for ct = ONE32:(nv-ONE32)
        imd(ct) = imd(ct) + nmv;   % Points to noise inputs to ID model
    end
    B(ixp,imd) = b(:,mdindex);
    D(:,imd)   = d(:,mdindex);
end
% Return estimator parameters
Bu = B(:,ONE32:nmv);
Bv = B(:,nmv+ONE32:nmv+nv);
Cm = C(myindex,:);
Dv = D(:,nmv+ONE32:nmv+nv);
Dvm = Dv(myindex,:);
% Exclude DX Offset input from "unknown" inputs for covariance calculation
UnknownIn = zeros(ONE32,nutot-ONE32);
UnknownIn(ONE32:nmv+nv-ONE32) = ONE32:nmv+nv-ONE32;
if nutot > nmv+nv
    UnknownIn(nmv+nv:nutot-ONE32) = (nmv+nv+ONE32):nutot;
end
CovMat = [B(:,UnknownIn); D(myindex,UnknownIn)]*...
                [B(:,UnknownIn); D(myindex,UnknownIn)]';
ix = ONE32:nx;
iy = nx+ONE32:nx+nym;
QQ = CovMat(ix,ix);
RR = CovMat(iy,iy);
NN = CovMat(ix,iy);

