  function[a,sd_a,var,covar,c] = mrqmin(x,y,sig,ndata,ma,a,eps);

% Routine zur Durchfuehrung einer Gauss-Newton-Marquardt-Iteration 
% nach dem Skriptum, S. 116ff

% H. Sormann  November 2001

% BESCHREIBUNG DER PARAMETER:

% Eingabe-Parameter:  x  Zeilenvektor der x-Komponenten der geg. Punkte
%                     y  Zeilenvektor der y-Komponenten der geg. Punkte
%                     sig  Zeilenvektor der Standardabweichungen der geg. Punkte
%                     ndata  Zahl der geg. Punkte
%                     ma     Zahl der Fit-Parameter
%                     a      Zeilenvektor der guessed values fuer die Fit-Paramter
%                     eps    relativer Fehler der Fit-Parameter, der zu unterschreiten ist
%
% Ausgabe-Parameter:  a      Zeilenvektor der gefitteten Parameter
%                     sd_a   Zeilenvektor der Standardabweichungen der Parameter
%                     var    Varianz des Fits   Anfang und Ende des Sollbereichs
%                     covar  Kovarianz-Matrix (normiert gemaess Skriptum, S. 109)

% Achtung: zum Unterschied zur im Skriptum beschriebenen Routine MRQMIN fuehrt
%          diese Routine eine vollstaendige GNM-Iteration durch, inklusive der
%          Fehlerdiagnostik und der statistischen Auswertung.

% Dieses Programm ruft die folgende Routine  mrqcof  auf. Diese Routine arbeitet
% genau wie im Skriptum (S. 125f) beschrieben.

% Die Routine  mrqcof  ruft ihrerseits die Routine  funcs  auf:
%
%  [ymod,dyda] = funcs(x,a,ma,ndata);
%  Eingabe-Parameter:  wie oben
%  Ausgabe-Parameter:  ymod    Zeilenvektor mit den ndata Werten der Modellfunktion
%                      dyda    Matrix (1. Index=Nummer des Parameters,
%                                      2. Index=Nummer des Datenpunktes) mit den partiellen
%                              Ableitungen der Modellfunktion nach den Modellparametern

%  DIESE FUNKTION MUSS VON BENUTZER DES MRQMIN-PROGRAMMS BEIGESTELLT WERDEN.

%  global inp
  alamda=0.0003;
  [alpha,beta,chisq] = mrqcof(x,y,sig,ndata,ma,a);

  fprintf(' t     chisq            Fit-Parameter ...\n\n');
%  fprintf(inp,' t     chisq            Fit-Parameter ...\n\n');
  fprintf('%3i  %12.6e    ',0,chisq);
%  fprintf(inp,'%3i  %12.6e    ',0,chisq);
  for j=1:ma
    fprintf('   %12.6e',a(j));
%    fprintf(inp,'   %12.6e',a(j));
  end
  fprintf('\n');
%  fprintf(inp,'\n');

  beta=beta';
  ochisq=chisq;

  aalt=a;
  for iter=1:100
    normal=alpha;
    for i=1:ma
      normal(i,i)=alpha(i,i)*(1+alamda);
    end

%    da=normal\beta;

    augmented=[normal beta];
    cr=rref(augmented);
    da=cr(:,ma+1);

    atry=a+da';
    [normal,vektor,chisq] = mrqcof(x,y,sig,ndata,ma,atry);

    if(chisq <= ochisq)
      alamda=alamda/5;
      ochisq=chisq;
      alpha=normal;
      a=atry;
      beta=vektor';

      fprintf('%3i  %12.6e    ',iter,chisq);
%      fprintf(inp,'%3i  %12.6e    ',iter,chisq);
      for j=1:ma
        fprintf('   %12.6e',a(j));
%        fprintf(inp,'   %12.6e',a(j));
      end
      fprintf('\n');
%      fprintf(inp,'\n');

% Fehlerdiagnostik:

      delmax=0;
      for j=1:ma
        del=abs((a(j)-aalt(j))/a(j));
        if(del > delmax)
          delmax=del;
        end
      end
      if(delmax < eps)
         fprintf('Genauigkeit wurde erreicht!!\n');
%         fprintf(inp,'Genauigkeit wurde erreicht!!\n');
% Berechnung der Varianz:
         var=[chisq/(ndata-ma) 1-sqrt(2/(ndata-ma)) 1+sqrt(2/(ndata-ma))];
         c=(ndata-ma)/chisq;
% Berechnung der Kovarianz-Matrix:
         alamda=0.0;
         [alpha,beta,chisq] = mrqcof(x,y,sig,ndata,ma,a);

%         covar=inv(alpha);

         covar=zeros(ma,ma);
         unitmat=eye(ma,ma);
         for i=1:ma
           inhom=unitmat(:,i);
           augmented=[normal inhom];
           cr=rref(augmented);
           covar(:,i)=cr(:,ma+1);
         end

         diag_covar=diag(covar);
         sd_a=sqrt(diag_covar)';
% Normierung der Kovarianzmatrix (s. Skriptum S. 109):
         for i=1:ma
           for j=1:ma
             covar(i,j)=covar(i,j)/sqrt(diag_covar(i)*diag_covar(j));
           end
         end
         return
      else
        aalt=a;
      end  
    else
      alamda=alamda*5;
      fprintf('Konvergenzproblem -- Marquardt!\n');
%      fprintf(inp,'Konvergenzproblem -- Marquardt!\n');
    end
  end
  fprintf('\nKEINE KONVERGENZ NACH 100 ITERATIONEN !!\n');
  fprintf('Programm mit Ctrl-c abbrechen.\n');
%  fprintf(inp,'\nKEINE KONVERGENZ NACH 100 ITERATIONEN !!\n');
%  fprintf(inp,'Programm mit Ctrl-c abbrechen.\n');
