unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, unit1,Forms, Dialogs;
Function GASDEV:real;
Procedure MRQMIN(var X,Y,SIG:array of real; NDATA:integer;var A:array of real;
     MA:integer;var LISTA:array of integer; MFIT:integer;var COVAR,ALPHA:matrx2;
     NCA:integer;var CHISQ,ALAMDA:real);

implementation
Function GASDEV:real;
var
   V1,V2,FAC,R:real;
begin
    If ISET^= 0 Then
    begin
      repeat
        V1:=2 * Random - 1;
        V2:=2 * Random - 1;
        R:=Sqr(V1) + Sqr(V2);
      until (R < 1);
      FAC:=Sqrt(-2 * Ln(R) / R);
      GSET^:=V1 * FAC;
      GASDEV:=V2 * FAC;
      ISET^:=1;
    end
    Else
    begin
      GASDEV:=GSET^;
      ISET^:=0;
    end;
end;

Procedure FGAUSS(X:real; A:array of real;
                  var Y:real;var DYDA:array of real; NA:integer);
var
    I,II:integer;  ARG,EX,FAC:real;
begin
    Y:=0;
    For II:=1 To (NA DIV 3) do
    begin
        I:=3*II-2;
        ARG:=(X - A[I + 1]) / A[I + 2];
        EX:=Exp(-Sqr(ARG));
        FAC:=A[I] * EX * 2 * ARG;
        Y:=Y + A[I] * EX;
        DYDA[I]:=EX;
        DYDA[I + 1]:=FAC / A[I + 2];
        DYDA[I + 2]:=FAC * ARG / A[I + 2];
    end; 
end;

Procedure MRQCOF(X,Y,SIG:array of real; NDATA:integer;var A:array of real;
     MA:integer; LISTA:array of integer; MFIT:integer;var ALPHA:matrx2;
     var BETA:array of real; NALP:integer;var CHISQ:real);
var
    DYDA:array[0..20] of real;
    J,I,K:integer;   YMOD,SIG2I,DY,WT:real;
begin
    For J:=1 To MFIT do
    begin
      For K:=1 To J do
        ALPHA[J, K]:=0;
      BETA[J]:=0;
    end; 
    CHISQ:=0; 
    For I:=1 To NDATA do
    begin
        FGAUSS(X[I], A, YMOD, DYDA, MA);
        SIG2I:=1 / (SIG[I] * SIG[I]);
        DY:=Y[I] - YMOD;
        For J:=1 To MFIT do
        begin
            WT:=DYDA[LISTA[J]] * SIG2I;
            For K:=1 To J do
                ALPHA[J, K]:=ALPHA[J, K] + WT * DYDA[LISTA[K]];
            BETA[J]:=BETA[J] + DY * WT;
        end; 
        CHISQ:=CHISQ + DY * DY * SIG2I;
    end; 
    For J:=2 To MFIT do
        For K:=1 To J - 1 do
            ALPHA[K, J]:=ALPHA[J, K];
end;

procedure GAUSSJ(VAR A:matrx2; N:integer; VAR B:array of real);
var
    IPIV,INDXR,INDXC:array[1..50] of integer;
    J,I,K,L,LL:integer;
    BIG,PIVINV,DUM:real; IROW,ICOL:integer;
begin
    For J:=1 To N do
      IPIV[J]:=0;
    For I:=1 To N do
    begin
        BIG:=0;
        For J:=1 To N do
        begin
            If IPIV[J] <> 1 Then
            begin
                For K:=1 To N do
                begin
                    If IPIV[K] = 0 Then
                    begin
                        If Abs(A[J, K]) >= BIG Then
                        begin
                            BIG:=Abs(A[J, K]);
                            IROW:=J;
                            ICOL:=K;
                        end;
                    end
                    Else if IPIV[K] > 1 Then
                        ShowMessage('Singular matrix.');
                end;
            end;
        end;
        IPIV[ICOL]:=IPIV[ICOL] + 1;
        If IROW <> ICOL Then
        begin
            For L:=1 To N do
            begin
                DUM:=A[IROW, L];
                A[IROW, L]:=A[ICOL, L];
                A[ICOL, L]:=DUM;
            end;
            DUM:=B[IROW];
            B[IROW]:=B[ICOL];
            B[ICOL]:=DUM;
        end;
        INDXR[I]:=IROW;
        INDXC[I]:=ICOL;
        If A[ICOL, ICOL] = 0 Then ShowMessage('Singular matrix.');
        PIVINV:=1 / A[ICOL, ICOL];
        A[ICOL, ICOL]:=1;
        For L:=1 To N do
          A[ICOL, L]:=A[ICOL, L] * PIVINV;
        B[ICOL]:=B[ICOL] * PIVINV;
        For LL:=1 To N do
        begin
            If LL <> ICOL Then
            begin
                DUM:=A[LL, ICOL];
                A[LL, ICOL]:=0;
                For L:=1 To N do
                    A[LL, L]:=A[LL, L] - A[ICOL, L] * DUM;
                B[LL]:=B[LL] - B[ICOL] * DUM;
            end;
        end;
    end;
    For L:=N DownTo 1 do
    begin
        If INDXR[L] <> INDXC[L] Then
        begin
            For K:=1 To N do
            begin
                DUM:=A[K, INDXR[L]];
                A[K, INDXR[L]]:=A[K, INDXC[L]];
                A[K, INDXC[L]]:=DUM;
            end;
        end;
    end;
end;

Procedure COVSRT(var COVAR:matrx2; NCVM, MA:integer;
                        LISTA:array of integer; MFIT:integer);
var
    I,J:integer;  SWAP:real;
begin
    For J:=1 To MA - 1 do
        For I:=J + 1 To MA do
            COVAR[I, J]:=0;
    For I:=1 To MFIT - 1 do
    begin
        For J:=I + 1 To MFIT do
        begin
            If LISTA[J] > LISTA[I] Then
                COVAR[LISTA[J], LISTA[I]]:=COVAR[I, J]
            Else
                COVAR[LISTA[I], LISTA[J]]:=COVAR[I, J];
        end;
    end;
    SWAP:=COVAR[1, 1];
    For J:=1 To MA do
    begin
        COVAR[1, J]:=COVAR[J, J];
        COVAR[J, J]:=0;
    end; 
    COVAR[LISTA[1], LISTA[1]]:=SWAP;
    For J:=2 To MFIT do
        COVAR[LISTA[J], LISTA[J]]:=COVAR[1, J];
    For J:=2 To MA do
        For I:=1 To J - 1 do
            COVAR[I, J]:=COVAR[J, I];
end;

Procedure MRQMIN(var X,Y,SIG:array of real; NDATA:integer;var A:array of real;
     MA:integer;var LISTA:array of integer; MFIT:integer;var COVAR,ALPHA:matrx2;
     NCA:integer;var CHISQ, ALAMDA:real);
var
    I,J,KK,K,IHIT:integer;
begin
  If ALAMDA < 0  Then
  begin
    KK:=MFIT + 1;
    For J:=1 To MA do
    begin
      IHIT:=0;
      For K:=1 To MFIT do
        If LISTA[K] = J Then IHIT:=IHIT + 1;
      If IHIT = 0 Then
      begin
        LISTA[KK]:=J;
        KK:=KK + 1;
      end
      Else If IHIT > 1 Then
      begin
        ShowMessage('Improper permutation in LISTA');
        Exit;
      end;
    end;
    If KK <> MA + 1 Then ShowMessage('Improper permutation in LISTA');
    ALAMDA:=0.001;
    MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NCA,CHISQ);
    OCHISQ:=CHISQ;
    For J:=1 To MA do
      ATRY[J]:=A[J];
  end;
  For J:=1 To MFIT do
  begin
    For K:=1 To MFIT do
      COVAR[J, K]:=ALPHA[J, K];
    COVAR[J, J]:=ALPHA[J, J] * (1 + ALAMDA);
    DA[J]:=BETA[J];
  end;
  GAUSSJ(COVAR, MFIT, DA);
  If ALAMDA = 0 Then
  begin
    COVSRT(COVAR, NCA, MA, LISTA, MFIT);
    Exit;
  end;
  For J:=1 To MFIT do
    ATRY[LISTA[J]]:=A[LISTA[J]] + DA[J];
  MRQCOF(X,Y,SIG,NDATA,ATRY,MA,LISTA,MFIT,COVAR,DA,NCA,CHISQ);
  If CHISQ < OCHISQ Then
  begin
    ALAMDA:=0.1 * ALAMDA;
    OCHISQ:=CHISQ;
    For J:=1 To MFIT do
    begin
      For K:=1 To MFIT do
        ALPHA[J, K]:=COVAR[J, K];
      BETA[J]:=DA[J];
      A[LISTA[J]]:=ATRY[LISTA[J]];
    end;
  end
  else
  begin
    ALAMDA:=10 * ALAMDA;
    CHISQ:=OCHISQ;
  end;
end;
end.