unit2.pas

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

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
Function BESSJ0(X:real):real;
Function BESSJ1(X:real):real;
Function BESSY1(X:real):real;
Function BESSY0(X:real):real;
Function BESSY(var N:integer; X:real):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;

Function BESSJ1(X:real):real;
VAR
    AAA,BBB,CCC,AX,Z,Y,XX,SGN:real;
const
    R1 = 72362614232;         R2 = -7895059235;
    R3 = 242396853.1;         R4 = -2972611.439;
    R5 = 15704.4826;          R6 = -30.16036606;
    S1 = 144725228442;        S2 = 2300535178;
    S3 = 18583304.74;         S4 = 99447.43394;
    S5 = 376.9991397;         S6 = 1;
    P1 = 1;                   P2 = 0.00183105;
    P3 = -0.3516396496e-4;    P4 = 0.000002457520174;
    P5 = -0.240337019e-6;
    Q1 = 0.04687499995;       Q2 = -0.2002690873e-3;
    Q3 = 0.8449199096e-5;     Q4 = -0.88228987e-6;
    Q5 = 0.105787412e-6;
begin
    If Abs(X) < 8 Then
      begin
        Y:=X * X;
        AAA:=R1 + Y * (R2 + Y * (R3 + Y * (R4 + Y * (R5 + Y * R6))));
        BBB:=S1 + Y * (S2 + Y * (S3 + Y * (S4 + Y * (S5 + Y * S6))));
        BESSJ1:=X * AAA / BBB;
      end
    Else
      begin
        AX:=Abs(X);
        Z:=8 / AX;
        Y:=Z * Z;
        XX:=AX - 2.356194491;
        AAA:=P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * P5)));
        BBB:=Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * Q5)));
        CCC:=Sqrt(0.636619772 / AX);
        if X > 0 THEN
          SGN:= 1
        ELSE
          SGN:= -1;
        BESSJ1:=CCC * (Cos(XX) * AAA - Z * Sin(XX) * BBB * Sgn);
      End;
End;

Function BESSY1(X:real):real;
var
    AAA,BBB,CCC,Z,Y,XX:real;
const
    P1 = 1;                    P2 = 0.183105e-2;
    P3 = -0.3516396496e-4;     P4 = 0.2457520174e-5;
    P5 = -0.240337019e-6;
    Q1 = 0.04687499995;        Q2 = -0.2002690873e-3;
    Q3 = 0.8449199096e-5;      Q4 = -0.88228987e-6;
    Q5 = 0.105787412e-6;
    R1 = -4.900604943e12;      R2 = 1.27527439e12;
    R3 = -5.153438139e10;      R4 = 7.349264551e8;
    R5 = -4237922.726;         R6 = 8511.937935;
    S1 = 2.49958057e13;        S2 = 4.244419664e11;
    S3 = 3.733650367e9;        S4 = 2.245904002e7;
    S5 = 102042.605;           S6 = 354.9632885;
    S7 = 1;
begin
    If X < 8 Then
      begin
        Y:= X * X;
        AAA:=R1 + Y * (R2 + Y * (R3 + Y * (R4 + Y * (R5 + Y * R6))));
        BBB:=S4 + Y * (S5 + Y * (S6 + Y * S7));
        BBB:=S1 + Y * (S2 + Y * (S3 + Y * BBB));
        CCC:=BESSJ1(X) * Ln(X) - 1/ X;
        BESSY1:=X * AAA / BBB + 0.636619772 * CCC;
      end
    Else
      begin
        Z:=8 / X;
        Y:=Z * Z;
        XX:=X - 2.356194491;
        AAA:=Sqrt(0.636619772 / X);
        BBB:=P1 + Y * (P2 + Y * (P3 + Y * (P4 + Y * P5)));
        CCC:=Q1 + Y * (Q2 + Y * (Q3 + Y * (Q4 + Y * Q5)));
        BESSY1:=AAA * (Sin(XX) * BBB + Z * Cos(XX) * CCC);
      end;
end;

Function BESSY0(X:real):real;
var
  AAA,BBB,Z,Y,XX,CCC:real;
const
    P1 = 1;                  P2 = -0.001098628627;
    P3 = 0.2734510407e-4;    P4 = -0.2073370639e-5;
    P5 = 2.093887211E-07;
    Q1 = -0.01562499995;     Q2 = 0.1430488765e-3;
    Q3 = -0.6911147651e-5;   Q4 = 7.621095161E-07;
    Q5 = -9.34945152E-08;
    R1 = -2.957821389E9;     R2 = 7.062834065E9;
    R3 = -512359803.6;       R4 = 10879881.29;
    R5 = -86327.92757;       R6 = 228.4622733;
    S1 = 40076544269;        S2 = 745249964.8;
    S3 = 7189466.438;        S4 = 47447.2647;
    S5 = 226.1030244;        S6 = 1;
begin
    If X < 8 Then
      begin
        Y:=X * X;
        AAA:=R1 + y * (R2 + y * (R3 + y * (R4 + y * (R5 + y * R6))));
        BBB:=S1 + y * (S2 + y * (S3 + y * (S4 + y * (S5 + y * S6))));
        BESSY0:=AAA / BBB + 0.636619772 * BESSJ0(X) * Ln(X);
      end
    Else
      begin
        Z:=8 / X;
        y:=Z * Z;
        XX:=X - 0.785398164;
        AAA:=Sqrt(0.636619772 / X);
        BBB:=P1 + y * (P2 + y * (P3 + y * (P4 + y * P5)));
        CCC:=Q1 + y * (Q2 + y * (Q3 + y * (Q4 + y * Q5)));
        BESSY0:=AAA * (Sin(XX) * BBB + Z * Cos(XX) * CCC);
      End;
End;

Function BESSY(var N:integer; X:real):real;
var
    TOX,BY,BYM,BYP:real;j:integer;
begin
    If N < 2 Then
        ShowMessage('bad argument N in BESSY');
    TOX:= 2 / X;
    By:= BESSY1(X);
    BYM:= BESSY0(X);
    For J:= 1 To N - 1 do
    begin
        BYP:= J * TOX * By - BYM;
        BYM:= By;
        By:= BYP;
    end;
    BESSY:= By;
end;

end.