function [F,G] = gradFG(Ny, Nparam, NparID, Ndata, Nu, Nx, Nzi, izhf,...
                        dt, Z, Y, RI, param, parFlag, par_del, state_eq, Uinp, ...
                        obser_eq, times, x0, iArtifStab, StabMat, integMethod)

% compute the information matrix (i.e. matrix of second gradients) F and 
% the gradient vector G
%
% 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:
%    Ny            number of output variables
%    Nparam        total number of parameters appearing in the postulated model
%    NparID        number of unknown parameters being estimated
%    Ndata         number of data points
%    Nu            number of input variables
%    Nx            number of state variables
%    Nzi           number of time segments being analyzed simultaneously
%    izhf          Cumulative index at which the concatenated time segments end
%    dt            sampling time
%    Z             measured outputs (Ndata,Ny)
%    Y             computed (model) outputs for unperturbed parameters (Ndata,Ny)
%    RI            inverse of covariance matrix
%    param         parameter vector
%    parFlag       flags for free and fixed parameters (=1, free parameter, 0: fixed)
%    par_del       parameter perturbations for numerical approximation of response gradients
%    state_eq      function coding the state equations
%    Uinp          measured inputs (Ndata,Nu) 
%    obser_eq      function coding the observation equations
%    times         time vector
%    x0            initial conditions
%    iArtifStab    flag for artificial stabilization
%    StabMat       artificial stabilization matrix
%    integMethod   integration method
%
% Outputs:
%    F             information matrix 
%    G             gradient vector 

%--------------------------------------------------------------------------
% Simulation for nominal parameter values:
% Array Y contains system responses for unperturbed parameters 

% Create cell array for perturbed system responses corresponding to Nparam parameters
YP = cell(Nparam, 1);

% Generate perturbed parameter simulations and store in cell array "YP"
iPar = 0;
for i1=1:Nparam,                            % Loop over all parameters
    
    if  parFlag(i1) > 0,                    % compute pertubed responses only for free parameters
        iPar = iPar + 1;

        paramp    = param;
        paramp(i1) = paramp(i1) + par_del(i1);

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

            while kk <= kend,               % Loop for data points of kzi-th time segment
                u1  = Uinp(kk,1:Nu)';  
                ys  = feval(obser_eq, times, xp, u1, paramp);    % outputs, Eq. (4.11), 
                ypi(kk,:) = ys';                                 % for pertubed parameters
        
                % artificial stabilization
                if  iArtifStab > 0,
                    zmy = Z(kk,:) - ys'; 
                    xp  = xp + StabMat*zmy';     % Equation (9.37)
                end
    
                if  kk < kend  &  Nx > 0,
                    u2 = Uinp(kk+1,1:Nu)';
                    % Integrate perturbed state eqs., i.e. Eq. (4.10) for pertubed parameters
                    xp = feval(integMethod, state_eq, times, xp, dt, u1, u2, paramp, 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
    
        YP{iPar} = struct('time', times, 'yps', ypi);   
        
    end;
    
end;                                        % End of i1-Loop for parameters

%--------------------------------------------------------------------------
% Initialize F, G and Q to zero
F = zeros([NparID, NparID]);
G = zeros([NparID, 1]);
Q = zeros([Ny, NparID]);

% Compute information matrix F and gradient G
for i1=1:length(times),                     % Loop over data points
    
    iPar = 0;
    for j1=1:Nparam,
        if  parFlag(j1) > 0,                % compute for free parameters only
            iPar = iPar + 1;
            q = YP{iPar}.yps(i1,:) - Y(i1,:);
            Q(:,iPar) = q'/par_del(j1);     % response gradients, Eq. (4.33)
        end;
    end;
    
    zmy = Z(i1,:) - Y(i1,:);                % response error  
    zmy = zmy(:);
    F = F + Q' * RI * Q;                    % information matrix, Eq. (4.30)
    G = G - Q' * RI * zmy;                  % gradient vector, Eq. (4.30)
    
end;

return
%end of function