function [dParam, lamda] = LMmethod(integMethod, state_eq, obser_eq, Ndata, Ny, Nu, Nzi, Nx,...
                                    Nparam, NparID, lamda, param, parFlag, dt, times, ...
                                    x0, Uinp, Z, izhf, F, G, eps_zf, cCost0,...
                                    iArtifStab, StabMat)

% Optimization of cost function applying Levenberg-Marquardt method
%
% 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
%    Nparam        total number of parameters appearing in the postulated model
%    NparID        number of unknown parameters being estimated
%    lamda         Levenberg-Marquardt parameter 
%    param         parameter vector
%    parFlag       flags for free and fixed parameters (=1, free parameter, 0: fixed)
%    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
%    F             information matrix 
%    G             gradient vector 
%    eps_zf        convergence limit in terms of relative change in the cost function
%    cCost0        previous cost function value
%    iArtifStab    flag for artificial stabilization
%    StabMat       artificial stabilization matrix
%
% Outputs:
%    dParam        parameter update (delta-Theta) 
%    lamda         Levenberg-Marquardt parameter 

%--------------------------------------------------------------------------
% compute the scale factors: square root of diagonal elements
sfac = 1./sqrt(diag(F));

% Scale G and F: Eq. (4.69) and (4.70)
G = G.*sfac;
for i1=1:NparID, for j1=1:NparID, Fscal(i1,j1) = F(i1,j1)*sfac(i1)*sfac(j1); end, end

% Add lamda to the diagonal elements of scaled F, Eq. (4.68)
F1 = Fscal + lamda*eye(NparID);


%--------------------------------------------------------------------------
idum = -1;                                            % Switch for jump back condition
% jump back address
while idum <= 0;

    % Compute parameter update vector 
    uptF1 = chol(F1);                                 % Cholesky factorization: F=uptF'*uptF
    dPar1 = -uptF1 \ ( uptF1' \ G );

    % Unscale the parameter update corresponding to lamda and compute new parameters
    dPar1 = dPar1.*sfac;
    iPar  = 0;
    for i1=1:Nparam,
        if  parFlag(i1) > 0,
            iPar = iPar + 1;
            paramStp1(i1) = param(i1) + dPar1(iPar);  % for free parameters only
        else
            paramStp1(i1) = param(i1);                % fixed parameters
        end
    end

    % Compute the cost function at lamda
    [cCost1, R, RI, Y] = costfun_oem (integMethod, state_eq, obser_eq, Ndata, Ny, Nu, Nzi,...
                                      Nx, dt, times, x0, Uinp, Z, izhf, paramStp1,...
                                      iArtifStab, StabMat);

                                  
    % Add lamda/10 to the diagonal of Fscal
    F2 = Fscal + lamda*0.1*eye(NparID);
    
    % Compute parameter update vector for reduced LM-parameter
    uptF2  = chol(F2);                               % Cholesky factorization: F2=uptF'*uptF
    dParam = -uptF2 \ ( uptF2' \ G );
    
    % Unscale the parameter update corresponding to lamda and compute new parameters
    dParam = dParam.*sfac;
    iPar   = 0;
    for i1=1:Nparam,
        if  parFlag(i1) > 0,
            iPar = iPar + 1;
            paramStp2(i1) = param(i1) + dParam(iPar); % for free parameters only
        else
            paramStp2(i1) = param(i1);                % fixed parameters
        end
    end

    % Compute the cost function at lamda/10
    [cCost2, R, RI, Y] = costfun_oem (integMethod, state_eq, obser_eq, Ndata, Ny, Nu, Nzi,...
                                      Nx, dt, times, x0, Uinp, Z, izhf, paramStp2,...
                                      iArtifStab, StabMat);

   % Evaluate the results
    costDiff2 = cCost2 - cCost0;
    costDiff1 = cCost1 - cCost0;

    if (costDiff2 <= 0),                              % condition 1
    
        % reduce lamda-parameter and accept the new point  
        currentCost = cCost2;
        %cCost0 = cCost2;
        lamda = lamda*0.1;
        i_lm  = 1;
        par_prnt = sprintf('%3i %-12s  %7.2e   %10.4e', i_lm, 'LM1: Lamda=',lamda, cCost2);
        disp(par_prnt)
        % Overwrite F with unscaled F2 prior to returning 
        for i1=1:NparID, for j1=1:NparID, F(i1,j1) = F2(i1,j1)/sfac(i1)/sfac(j1); end, end
        idum = 1;                                     % set switch to 1 for return
        
    else
        dum4 = costDiff1 - eps_zf*cCost0;
        
        if (costDiff2 > 0  & costDiff1 <= 0) | dum4 <= 0,
            
            % Keep current lamda and accept new point
            currentCost = cCost1;
            %cCost0 = cCost1;
            i_lm = 2;
            par_prnt = sprintf('%3i %-12s  %7.2e   %10.4e', i_lm, 'LM2: Lamda=',lamda, cCost1);
            disp(par_prnt)
            % Prior to returning, get dParam for the current lamda and unscale F
            dParam = dPar1;
            for i1=1:NparID, for j1=1:NparID, F(i1,j1) = F1(i1,j1)/sfac(i1)/sfac(j1); end, end
            idum = 1;                                 % set switch to 1 for return

        else
        
            % Increase lamda and reject new point and loop back
            lamda = lamda*10;
            i_lm  = 3;
            par_prnt = sprintf('%3i %-12s  %7.2e   %10.4e', i_lm, 'LM3: Lamda=',lamda, cCost2);
            disp(par_prnt)
            % Add new lamda to the diagonal of scaled F and jump back  
            F1 = Fscal + lamda*eye(NparID);

        end
    
    end
    
end

return
% end of function