unit2.pas

用于开发税务票据管理的软件

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics,unit1, Controls, Forms, Dialogs;
Procedure AMOEBA(var P:matrx2;var Y:array of real; MP, NP:integer;
                            NDIM:integer; FTOL:real;var ITER:integer);
Function BESSJ0(X:real):real;                            

implementation
Procedure AMOEBA(var P:matrx2;var Y:array of real; MP, NP:integer;
                            NDIM:integer; FTOL:real;var ITER:integer);
const
    NMAX = 20; ALPHA = 1;  BETA = 0.5;  GAMMA = 2;  ITMAX = 500;
var
    PR, PRR, PBAR:array[0..20] of real;
    I,J,ILO,MPTS,IHI,INHI:integer;
    RTOL,YPR,YPRR:real;
begin
    MPTS:=NDIM + 1;
    ITER:=0;
    while true do
    begin
      ILO:=1;
      If Y[1] > Y[2] Then
      begin
          IHI:=1;
          INHI:=2;
      end
      Else
      begin
          IHI:=2;
          INHI:=1;
      end;
      For I:=1 To MPTS do
      begin
          If Y[I] < Y[ILO] Then ILO:=I;
          If Y[I] > Y[IHI] Then
          begin
              INHI:=IHI;
              IHI:=I;
          end
          Else If Y[I] > Y[INHI] Then
              If I <> IHI Then INHI:=I;
      end;
      RTOL:=2 * Abs(Y[IHI] - Y[ILO]) / (Abs(Y[IHI]) + Abs(Y[ILO]));
      If RTOL < FTOL Then Exit;
      If ITER = ITMAX Then
      begin
          ShowMessage('Amoeba exceeding maximum iterations.');
          Exit;
      end;
      ITER:=ITER + 1;
      For J:=1 To NDIM do
          PBAR[J]:=0;
      For I:=1 To MPTS do
      begin
          If I <> IHI Then
          begin
              For J:=1 To NDIM do
                  PBAR[J]:=PBAR[J] + P[I, J];
          end;
      end;
      For J:=1 To NDIM do
      begin
          PBAR[J]:=PBAR[J] / NDIM;
          PR[J]:=(1 + ALPHA) * PBAR[J] - ALPHA * P[IHI, J];
      end;
      YPR:=FAMOEB(PR);
      If YPR <= Y[ILO] Then
      begin
          For J:=1 To NDIM do
              PRR[J]:=GAMMA * PR[J] + (1 - GAMMA) * PBAR[J];
          YPRR:=FAMOEB(PRR);
          If YPRR < Y[ILO] Then
          begin
              For J:=1 To NDIM do
                  P[IHI, J]:=PRR[J];
              Y[IHI]:=YPRR;
          end
          Else
          begin
              For J:=1 To NDIM do
                  P[IHI, J]:=PR[J];
              Y[IHI]:=YPR;
          end;
      end
      Else If YPR >= Y[INHI] Then
      begin
          If YPR < Y[IHI] Then
          begin
              For J:=1 To NDIM do
                  P[IHI, J]:=PR[J];
              Y[IHI]:=YPR;
          end;
          For J:=1 To NDIM do
              PRR[J]:=BETA * P[IHI, J] + (1 - BETA) * PBAR[J];
          YPRR:=FAMOEB(PRR);
          If YPRR < Y[IHI] Then
          begin
              For J:=1 To NDIM do
                  P[IHI, J]:=PRR[J];
              Y[IHI]:=YPRR;
          end
          Else
          begin
              For I:=1 To MPTS do
              begin
                  If I <> ILO Then
                  begin
                      For J:=1 To NDIM do
                      begin
                          PR[J]:=0.5 * (P[I, J] + P[ILO, J]);
                          P[I, J]:=PR[J];
                      end;
                      Y[I]:=FAMOEB(PR);
                  end;
              end;
          end;
      end
      Else if (YPR > Y[ILO]) and (YPR < Y[INHI]) then
      begin
          For J:=1 To NDIM do
              P[IHI, J]:=PR[J];
          Y[IHI]:=YPR;
      end;
    end;
end;

Function BESSJ0(X:real):real;
var
   AAA,BBB,CCC,Y,AX,Z,DDD,EEE,XX:real;
const
  P1=1;                    P2=-0.001098628627;
  P3=0.2734510407e-4;      P4=-0.2073370639e-5;
  P5=2.093887211E-07;
  Q1=-0.1562499995e-1;     Q2=0.1430488765e-3;
  Q3=-0.6911147651e-5;     Q4=7.621095161E-07;
  Q5=-9.34945152E-08;
  R1=57568490574;          R2=-13362590354;
  R3=651619640.7;          R4=-11214424.18;
  R5=77392.33017;          R6=-184.9052456;
  S1=57568490411;          S2=1029532985;
  S3=9494680.718;          S4=59272.64853;
  S5=267.8532712;          S6=1;
begin
  If Abs(X) < 8 Then
    begin
      Y:=X * X;
      BBB:=Y* (R4+ Y* (R5+ Y* R6));
      AAA:=R1+ Y* (R2+ Y* (R3+ BBB));
      CCC:=Y* (S3+ Y* (S4+ Y* (S5+ Y* S6)));
      BESSJ0:= AAA / (S1+ Y* (S2+ CCC));
    end
  Else
    begin
      AX:=Abs(X);
      Z:=8/ AX;
      Y:=Z* Z;
      XX:= AX- 0.785398164;
      CCC:=Y* (P3+ Y* (P4+ Y* P5));
      AAA:=P1+ Y* (P2+ CCC);
      DDD:=Y* (Q3+ Y* (Q4+ Y* Q5));
      EEE:=Z* Sin(XX) * (Q1+ Y* (Q2+ DDD));
      BESSJ0:= Sqrt(0.636619772 / AX) * (Cos(XX) * AAA- EEE);
    End;
End;
end.