unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Memo1: TMemo;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;
type
  matrx2=array of array of real;
var
  Form1: TForm1;
  ISET:^integer; GSET:^real;
  ATRY, BETA, DA:array[0..20] of real;
  OCHISQ:real;
implementation
//PROGRAM D7R8A
//Driver for routine MRQMIN
uses
  unit2;
{$R *.DFM}
procedure TForm1.Button1Click(Sender: TObject);
Label 1;
const
  s1='%10.4f'; s2='#'; s3='.##E+00';
  NPT = 100; MA = 6; SPREAD = 0.001;
var
  COVAR,ALPHA:matrx2;
  X,Y,SIG:array[0..100] of real;   A,GUES:array[0..6] of real;
  LISTA:array[0..6] of integer;
  I,K,ITST,MFIT:integer;    CHISQ,ALAMDA,OCHISQ:real;
  F:TextFile;
begin
  SetLength(COVAR,7,7);  SetLength(ALPHA,7,7);
  New(ISET); New(GSET);
  Iset^:=0; GSET^:=0;
  Randomize;
  A[1]:=5;  A[2]:=2;  A[3]:=3;
  A[4]:=2;  A[5]:=5;  A[6]:=3;
  GUES[1]:=4.5;  GUES[2]:=2.2;  GUES[3]:=2.8;
  GUES[4]:=2.5;  GUES[5]:=4.9;  GUES[6]:=2.8;
  //First try a sum of two Gaussians
  For I:= 1 To NPT do
  begin
    X[I]:= 0.1 * I;
    Y[I]:= 0;
    Y[I]:= Y[I] + A[1] * Exp(-Sqr((X[I] - A[1 + 1]) / A[1 + 2]));
    Y[I]:= Y[I] + A[4] * Exp(-Sqr((X[I] - A[4 + 1]) / A[4 + 2]));
    Y[I]:= Y[I] * (1 + SPREAD * GASDEV);
    SIG[I]:= SPREAD * Y[I];
  end;
  MFIT:=MA;
  For I:=1 To MFIT do
    LISTA[I]:=I;
  ALAMDA:=-1;
  For I:=1 To MA do
    A[I]:=GUES[I];
  MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,MA,CHISQ,ALAMDA);
  K:=1;
  ITST:=0;
  //輸出計算結果到文件
  AssignFile(F, 'd:\delphi_shu\p7\d7r8A.dat');
  Rewrite(F);
1:  Writeln(F,'Iteration # ',FormatFloat(s2,K),' Chi-squared: ',
              format(s1,[CHISQ]),'  ALAMDA: ',FormatFloat(s3,ALAMDA));
  Writeln(F,'     A(1)      A(2)      A(3)      A(4)      A(5)     A(6)');
  Writeln(F,Format(s1,[A[1]]),Format(s1,[A[2]]),Format(s1,[A[3]]),
            Format(s1,[A[4]]),Format(s1,[A[5]]),Format(s1,[A[6]]));
  Writeln(F);
  K:=K + 1;
  OCHISQ:=CHISQ;
  MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,MA,CHISQ,ALAMDA);
  If CHISQ > OCHISQ Then
    ITST:=0
  Else If Abs(OCHISQ - CHISQ) < 0.1 Then
    ITST:=ITST + 1;
  If ITST < 2 Then  GoTo 1;
  ALAMDA:=0;
  MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA,MA,CHISQ,ALAMDA);
  Writeln(F,'Uncertainties:');
  Writeln(F,Format(s1,[Sqrt(COVAR[1,1])]),Format(s1,[Sqrt(COVAR[2,2])]),
            Format(s1,[Sqrt(COVAR[3,3])]),Format(s1,[Sqrt(COVAR[4,4])]),
            Format(s1,[Sqrt(COVAR[5,5])]),Format(s1,[Sqrt(COVAR[6,6])]));
  Writeln(F);
  Writeln(F,'Expected results:');
  Writeln(F,'     5.0       2.0       3.0       2.0       5.0       3.0');
  CloseFile(F);
  //屏幕顯示計算結果
  memo1.Lines.LoadFromFile('d:\delphi_shu\p7\d7r8A.dat');
end;

end.