unit2.pas

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

unit Unit2;

interface
uses
  unit1;
Function BESSJ0(X:real):real;
Procedure MNBRAK(var AX, BX, CX, FA, FB, FC:real);

implementation
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;

Procedure MNBRAK(var AX, BX, CX, FA, FB, FC:real);
const
    GOLD = 1.618034;  GLIMIT = 100;  TINY = 1E-20;
var 
    DUM,R,Q,U,ULIM,FU:real;   DONE:boolean;
begin
    FA:=FUNC(AX);
    FB:=FUNC(BX);
    If FB > FA Then
    begin
        DUM:=AX;
        AX:=BX;
        BX:=DUM;
        DUM:=FB;
        FB:=FA;
        FA:=DUM;
    end;
    CX:=BX + GOLD * (BX - AX);
    FC:=FUNC(CX);
    repeat
        If FB < FC Then break;
        DONE:=true;
        R:=(BX - AX) * (FB - FC);
        Q:=(BX - CX) * (FB - FA);
        DUM:=Q - R;
        If Abs(DUM) < TINY Then DUM:=TINY;
        U:=BX - ((BX - CX) * Q - (BX - AX) * R) / (2 * DUM);
        ULIM:=BX + GLIMIT * (CX - BX);
        If (BX - U) * (U - CX) > 0  Then
        begin
            FU:=FUNC(U);
            If FU < FC Then
            begin
                AX:=BX;
                FA:=FB;
                BX:=U;
                FB:=FU;
                Exit;
            end
            Else If FU > FB Then
            begin
                CX:=U;
                FC:=FU;
                Exit;
            end;
            U:=CX + GOLD * (CX - BX);
            FU:=FUNC(U);
        end
        Else If (CX - U) * (U - ULIM) > 0  Then
        begin
            FU:=FUNC(U);
            If FU < FC Then
            begin
                BX:=CX;
                CX:=U;
                U:=CX + GOLD * (CX - BX);
                FB:=FC;
                FC:=FU;
                FU:=FUNC(U);
            end;
        end
        Else If (U - ULIM) * (ULIM - CX) >= 0  Then
        begin
            U:=ULIM;
            FU:=FUNC(U);
        end
        Else
        begin
            U:=CX + GOLD * (CX - BX);
            FU:=FUNC(U);
        end;
        If DONE Then
        begin
            AX:=BX;
            BX:=CX;
            CX:=U;
            FA:=FB;
            FB:=FC;
            FC:=FU;
        end
        Else
            DONE:=false;
    until not DONE
end;

end.