unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics,unit1, Controls, Forms, Dialogs;
procedure BSSTEP(var Y, DYDX:array of real; NV:integer; var X1:real;
                HTRY, EPS:real;YSCAL:array of real; var HDID, HEND:real);
Function BESSJ0(X:real):real;
Function BESSJ1(X:real):real;
Function BESSJ(N:Integer; X:Real):Real;

implementation
procedure MMID(Y, DYDX:array of real; NVAR:integer; XS, HTOT:real;
                                    NSTEP:integer;var YOUT:array of real);
var
    YM, YN:array[0..10] of real;  I,N:integer;  H,X,H2,SWAP:real;
begin
    H:=HTOT / NSTEP;
    For I:=1 To NVAR do
    begin
        YM[I]:=Y[I];
        YN[I]:=Y[I] + H * DYDX[I];
    end; 
    X:=XS + H;
    DERIVS(X, YN, YOUT);
    H2:=2 * H;
    For N:=2 To NSTEP do
    begin
        For I:=1 To NVAR do
        begin
            SWAP:=YM[I] + H2 * YOUT[I];
            YM[I]:=YN[I];
            YN[I]:=SWAP;
        end; 
        X:=X + H;
        DERIVS(X, YN, YOUT);
    end; 
    For I:=1 To NVAR do
        YOUT[I]:=0.5 * (YM[I] + YN[I] + H * YOUT[I]);
end;

procedure RZEXTR(IEST:integer; XEST:real; YEST:array of real;
                           var YZ, DY:array of real; NV, NUSE:integer);
var
    FX:array[0..7] of real;   J,M1,K:integer;
    YY,V,C,B1,DDY,B:real;
begin
    SetLength(D, 11, 8);
    X[IEST]:=XEST;
    If IEST = 1 Then
    begin
        For J:=1 To NV do
        begin
            YZ[J]:=YEST[J];
            D[J, 1]:=YEST[J];
            DY[J]:=YEST[J];
        end;
    end 
    Else
    begin
        M1:=IEST;
        If NUSE < IEST Then M1:=NUSE;
        For K:=1 To M1 - 1 do
            FX[K + 1]:=X[IEST - K] / XEST;
        For J:=1 To NV do
        begin
            YY:=YEST[J];
            V:=D[J, 1];
            C:=YY;
            D[J, 1]:=YY;
            For K:=2 To M1 do
            begin
                B1:=FX[K] * V;
                B:=B1 - C;
                If B <> 0 Then
                begin
                    B:=(C - V) / B;
                    DDY:=C * B;
                    C:=B1 * B;
                end
                Else
                    DDY:=V;
                If K <> M1 Then V:=D[J, K];
                D[J, K]:=DDY;
                YY:=YY + DDY;
            end; 
            DY[J]:=DDY;
            YZ[J]:=YY;
        end; 
    end;
end;

procedure BSSTEP(var Y, DYDX:array of real; NV:integer; var X1:real;
                HTRY, EPS:real;YSCAL:array of real; var HDID, HEND:real);
label 1;
const
    IMAX = 11;  NUSE = 7;  ONE = 1;  SHRINK = 0.95;  GROW = 1.2;
var
    I,J:integer;   H,XSAV,XEST,ERRMAX:real;
    YERR,YSAV,DYSAV,YSEQ:array[0..10] of real;
    NSEQ:array[0..11] of integer;
begin
    NSEQ[1]:=2;   NSEQ[2]:=4;   NSEQ[3]:=6;  NSEQ[4]:=8;  NSEQ[5]:=12;
    NSEQ[6]:=16;  NSEQ[7]:=24;  NSEQ[8]:=32; NSEQ[9]:=48;
    NSEQ[10]:=64; NSEQ[11]:=96;
    H:=HTRY;
    XSAV:=X1;
    For I:=1 To NV do
    begin
        YSAV[I]:=Y[I];
        DYSAV[I]:=DYDX[I];
    end;
1:  For I:=1 To IMAX do
    begin
        MMID(YSAV, DYSAV, NV, XSAV, H, NSEQ[I], YSEQ);
        XEST:=Sqr(H / NSEQ[I]);
        RZEXTR(I, XEST, YSEQ, Y, YERR, NV, NUSE);
        ERRMAX:=0;
        For J:=1 To NV do
            If Abs(YERR[J] / YSCAL[J]) > ERRMAX Then
                ERRMAX:=Abs(YERR[J] / YSCAL[J]);
        ERRMAX:=ERRMAX / EPS;
        If ERRMAX < ONE Then
        begin
            X1:=X1 + H;
            HDID:=H;
            If I = NUSE Then
                HEND:=H * SHRINK
            Else
            begin
                If I = NUSE - 1 Then
                    HEND:=H * GROW
                Else
                    HEND:=(H * NSEQ[NUSE - 1]) / NSEQ[I];
            end;
            Exit;
        end;
    end;
    H:=0.25 * H / Exp(((IMAX - NUSE) div 2)*Ln(2));
    if X1 + H = X1  then  ShowMessage(' Step size underflow');
    goto 1;
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;

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 BESSJ(N:Integer; X:Real):Real;
var
    BJ,BJM,BJP,SUM,TOX,ANS:Real; J,JSUM,M:Integer;
const
    IACC = 40;    BIGNO = 1.0e10;    BIGNI = 1.0e-10;
begin
    If N < 2 Then
        ShowMessage('bad argument N in BASSJ');
    TOX:= 2 / X;
    if X > 1.0 * N then
    begin
        BJM:=BESSJ0(X);
        BJ:=BESSJ1(X);
        For J:=1 To N - 1 do
        begin
            BJP:=J * TOX * BJ - BJM;
            BJM:=BJ;
            BJ:=BJP;
        end;
        ANS:=BJ
    end
    Else
    begin
        M:=2 * ((N + Trunc(Sqrt(IACC * N))) div 2);
        ANS:=0.0;
        JSUM:=0;
        Sum:=0;
        BJP:=0;
        BJ:=1;
        For J:=M DownTo 1 do
          begin
            BJM:=J * TOX * BJ - BJP;
            BJP:=BJ;
            BJ:=BJM;
            If Abs(BJ) > BIGNO Then
              begin
                BJ:=BJ * BIGNI;
                BJP:=BJP * BIGNI;
                ANS:=ANS * BIGNI;
                Sum:=Sum * BIGNI;
              end;
            If JSUM <> 0 Then Sum:=Sum + BJ;
            JSUM:=1 - JSUM;
            If J = N Then ANS:=BJP;
        end;
        Sum:=2 * Sum - BJ;
        ans:=ans / Sum
    end;
    BESSJ:=ans;
end;


end.