function [currentcost, R, RI, Y] = costfun_oem (integMethod, state_eq, obser_eq, Ndata,...
                                                Ny, Nu, Nzi, Nx, dt, times, x0, Uinp, Z,...
                                                izhf, param, iArtifStab, StabMat) 

% Compute cost function for the output error method, the covariance matrix and
% its inverse, and the model outputs through simulation.
% Simulation implies computation of state variables by numerical integration.
%
% Chapter 4: Output Error Method
% "Flight Vehicle System Identification - A Time Domain Methodology"
% Author: Ravindra V. Jategaonkar
% Published by AIAA, Reston, VA 20191, USA
%
% Inputs:
%    integMethod   integration method
%    state_eq      function coding the state equations
%    obser_eq      function coding the observation equations
%    Ndata         number of data points
%    Ny            number of output variables
%    Nu            number of input variables
%    Nzi           number of time segments being analyzed simultaneously
%    Nx            number of state variables
%    dt            sampling time
%    times         time vector
%    x0            initial conditions
%    Uinp          measured inputs (Ndata,Nu) 
%    Z             measured outputs (Ndata,Ny)
%    izhf          Cumulative index at which the concatenated time segments end
%    param         parameter vector
%    iArtifStab    flag for artificial stabilization
%    StabMat       artificial stabilization matrix
%
% Outputs:
%    currentcost   Value of current cost function (determinant of R)
%    R             covariance matrix
%    RI            inverse of covariance matrix
%    Y             computed (model) outputs (Ndata,Ny)

%--------------------------------------------------------------------------
kzi = 0;
kk  = 1;
while kzi < Nzi;                       % Loop for Nzi time segments
    kzi  = kzi + 1;
    kend = izhf(kzi);
    if isempty(x0),
        x = x0;
    else
        x = x0(:,kzi);                 % get initial conditions for the time segment 
    end                                % being processed

    while kk <= kend,                  % Loop for data points of kzi-th time segment
        u1 = Uinp(kk,1:Nu)';  
        y  = feval(obser_eq, times, x, u1, param);    % compute output variables, Eq. (4.11)
        Y(kk,:) = y';
    
        % artificial stabilization
        if  iArtifStab > 0,
            zmy = Z(kk,:) - y'; 
            x   = x + StabMat*zmy';    % Equation (9.37)
        end
    
        if  kk < kend  &  Nx > 0,
            u2 = Uinp(kk+1,1:Nu)';
            % Numerical integration of state equations, Eq. (4.10)
            x  = feval(integMethod, state_eq, times, x, dt, u1, u2, param, Nx);
        end
        kk = kk + 1;
    end                                % end of KK-loop for data points of kzi-th segment
end                                    % end of Nzi-loop for time segments

%--------------------------------------------------------------------------
% estimate covariance matrix R
R = zeros(Ny,Ny);
for kk=1:Ndata,
    delta = Z(kk,:) - Y(kk,:); 
    delta = delta(:);
    R     = R + delta * delta';
end;

% determinant and inverse of R
R    = diag(diag(R))/Ndata;            % covariance matrix R, Eq. (4.15)
detr = det(R);                         % determinant of R
RI   = inv(R);                         % Inverse of R
currentcost = detr;                    % cost function, Eq. (4.17)

return