VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   5310
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   5310
   ScaleWidth      =   4680
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3240
      TabIndex        =   0
      Top             =   4680
      Width           =   1095
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
    'PROGRAM D10R7
    'Driver for routine ZBRENT
    N = 100
    NBMAX = 20
    X1 = 1#
    X2 = 50#
    Dim XB1(20), XB2(20)
    NB = NBMAX
    Call ZBRAK(X1, X2, N, XB1(), XB2(), NB)
    Print
    Print Tab(5); "Roots of BESSJ0:"
    Print
    Print "                 x              F(x)"
    Print
    For I = 1 To NB
        TOL = (0.000001) * (XB1(I) + XB2(I)) / 2#
        ROOT = ZBRENT(XB1(I), XB2(I), TOL)
        Print Tab(5); "Root"; Format$(I, "##");
        Print Tab(15); Format$(ROOT, "##.####00");
        Print Tab(29); Format$(BESSJ0(ROOT), ".###0E+00")
    Next I
End Sub
Function FUNC(X)
    FUNC = BESSJ0(X)
End Function
Function ZBRENT(X1, X2, TOL)
    ITMAX = 100
    EPS = 0.00000003
    A = X1
    B = X2
    FA = FUNC(A)
    FB = FUNC(B)
    If FB * FA > 0# Then
        Print "Root must be bracketed for ZBRENT"
    End If
    FC = FB
    For ITER = 1 To ITMAX
        If FB * FC > 0 Then
            C = A
            FC = FA
            D = B - A
            E = D
        End If
        If Abs(FC) < Abs(FB) Then
            A = B
            B = C
            C = A
            FA = FB
            FB = FC
            FC = FA
        End If
        TOL1 = 2# * EPS * Abs(B) + 0.5 * TOL
        XM = 0.5 * (C - B)
        If Abs(XM) <= TOL1 Or FB = 0 Then
            ZBRENT = B
            Exit Function
        End If
        If Abs(E) >= TOL1 And Abs(FA) > Abs(FB) Then
            S = FB / FA
            If A = C Then
                P = 2# * XM * S
                Q = 1 - S
            Else
                Q = FA / FC
                R = FB / FC
                P = S * (2# * XM * Q * (Q - R) - (B - A) * (R - 1#))
                Q = (Q - 1#) * (R - 1#) * (S - 1#)
            End If
            If P > 0# Then Q = -Q
            P = Abs(P)
            If 3# * XM * Q - Abs(TOL1 * Q) < Abs(E * Q) Then
                AAA = 3# * XM * Q - Abs(TOL1 * Q)
            Else
                AAA = Abs(E * Q)
            End If
            If 2# * P < AAA Then
                E = D
                D = P / Q
            Else
                D = XM
                E = D
            End If
        Else
            D = XM
            E = D
        End If
        A = B
        FA = FB
        If Abs(D) > TOL1 Then
            B = B + D
        Else
            B = B + Abs(TOL1) * Sgn(XM)
        End If
        FB = FUNC(B)
    Next ITER
    Print "ZBRENT exceeding maximum iterations."
    ZBRENT = B
End Function
Sub ZBRAK(X1, X2, N, XB1(), XB2(), NB)
      NBB = NB
      NB = 0
      X = X1
      DX = (X2 - X1) / N
      FP = FUNC(X)
      For I = 1 To N
          X = X + DX
          FC = FUNC(X)
          If FC * FP < 0# Then
              NB = NB + 1
              XB1(NB) = X - DX
              XB2(NB) = X
          End If
          FP = FC
          If NBB = NB Then Exit Sub
      Next I
End Sub
Function BESSJ0(X)
    P1# = 1#
    P2# = -0.001098628627
    P3# = 0.00002734510407
    P4# = -0.000002073370639
    P5# = 2.093887211E-07
    Q1# = -0.01562499995
    Q2# = 0.0001430488765
    Q3# = -0.000006911147651
    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#
    If Abs(X) < 8# Then
       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#))
    Else
       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 = Sqr(0.636619772 / AX#) * (Cos(XX) * AAA# - EEE#)
    End If
End Function