VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   3840
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   6540
   LinkTopic       =   "Form1"
   ScaleHeight     =   3840
   ScaleWidth      =   6540
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   4680
      TabIndex        =   0
      Top             =   3240
      Width           =   1095
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Dim C2, FACTR, M, N
Private Sub Command1_Click()
    'PROGRAM D15R1
    'Driver for routine SHOOT
    'Solves for eigenvalues of Spheroidal Harmonics. Both
    'Prolate and Oblate case are handled simultaneously,
    'leading to six first-order equations.  Unknown to
    'SHOOT, these are actually two independent sets of
    'three coupled equations, one set with c^2 positive
    'and the other with c^2 negative.
    NVAR = 6: N2 = 2: DELTA = 0.001: EPS = 0.000001
    Dim V(2), DELV(2), F(2), DV(2)
    DX = 0.0001
    Do
      Print "     Input M,N,C-Squared (999 to end)"
      M = 2
      N = 2
      C2 = 0.1
      If C2 = 999 Then End
    Loop While (N < M) Or (M < 0) Or (N < 0)
    Print Tab(5); M; Tab(10); N, Tab(15); C2
    FACTR = 1
    If M <> 0 Then
        Q1 = N
        For I = 1 To M
            FACTR = -0.5 * FACTR * (N + I) * (Q1 / I)
            Q1 = Q1 - 1
        Next I
    End If
    V(1) = N * (N + 1) - M * (M + 1) + C2 / 2#
    V(2) = N * (N + 1) - M * (M + 1) - C2 / 2#
    DELV(1) = DELTA * V(1)
    DELV(2) = DELV(1)
    H1 = 0.1
    HMIN = 0#
    X1 = -1# + DX
    X2 = 0#
    Print "            Prolate                   Oblate"
    Print "     Mu(M,N)      Error Est.   Mu(M,N)     Error Est."
    Do
      Call SHOOT(NVAR, V(), DELV(), N2, X1, X2, EPS, H1, HMIN, F(), DV())
      Print Tab(5); Format$(V(1), "#.#####0"); Tab(18); Format$(DV(1), "#.#####0"),
      Print Tab(31); Format$(V(2), "#.#####0"), Tab(45); Format$(DV(2), "#.#####0")
    Loop While Abs(DV(1)) > Abs(EPS * V(1)) Or Abs(DV(2)) > Abs(EPS * V(2))
End Sub
Sub Load(X1, V(), Y())
    Y(3) = V(1)
    Y(2) = -(Y(3) - C2) * FACTR / 2# / (M + 1#)
    Y(1) = FACTR + Y(2) * DX
    Y(6) = V(2)
    Y(5) = -(Y(6) + C2) * FACTR / 2# / (M + 1#)
    Y(4) = FACTR + Y(5) * DX
End Sub
Sub SCORE(X2, Y(), F())
    If ((N - M) Mod 2) = 0 Then
        F(1) = Y(2)
        F(2) = Y(5)
    Else
        F(1) = Y(1)
        F(2) = Y(4)
    End If
End Sub
Sub DERIVS(X, Y(), DYDX())
    DYDX(1) = Y(2)
    DYDX(3) = 0#
    DYDX(2) = (2# * X * (M + 1#) * Y(2) - (Y(3) - C2 * X * X) * Y(1)) / (1# - X * X)
    DYDX(4) = Y(5)
    DYDX(6) = 0#
    DYDX(5) = (2# * X * (M + 1#) * Y(5) - (Y(6) + C2 * X * X) * Y(4)) / (1# - X * X)
End Sub
Sub SHOOT(NVAR, V(), DELV(), N2, X1, X2, EPS, H1, HMIN, F(), DV())
    Dim Y(20), DFDV(20, 20), INDX(20), XP(200), YP(10, 200)
    Call Load(X1, V(), Y())
    Call ODEINT(Y(), NVAR, X1, X2, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
    Call SCORE(X2, Y(), F())
    For IV = 1 To N2
        SAV = V(IV)
        V(IV) = V(IV) + DELV(IV)
        Call Load(X1, V(), Y())
        Call ODEINT(Y(), NVAR, X1, X2, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
        Call SCORE(X2, Y(), DV())
        For I = 1 To N2
            DFDV(I, IV) = (DV(I) - F(I)) / DELV(IV)
        Next I
        V(IV) = SAV
    Next IV
    For IV = 1 To N2
        DV(IV) = -F(IV)
    Next IV
    Call LUDCMP(DFDV(), N2, INDX(), DET)
    Call LUBKSB(DFDV(), N2, INDX(), DV())
    For IV = 1 To N2
        V(IV) = V(IV) + DV(IV)
    Next IV
End Sub
Sub ODEINT(YSTART(), NVAR, X1, X2, EPS, H1, HMIN, NOK, NBAD, XP(), YP())
    MAXSTP = 10000
    TWO = 2#
    ZERO = 0#
    TINY = 1E-30
    Dim YSCAL(10), Y(10), DYDX(10)
    X = X1
    H = Abs(H1) * Sgn(X2 - X1)
    NOK = 0
    NBAD = 0
    KOUNT = 0
    For I = 1 To NVAR
        Y(I) = YSTART(I)
    Next I
    If KMAX > 0 Then XSAV = X - DXSAV * TWO
    For NSTP = 1 To MAXSTP
        Call DERIVS(X, Y(), DYDX())
        For I = 1 To NVAR
            YSCAL(I) = Abs(Y(I)) + Abs(H * DYDX(I)) + TINY
        Next I
        If KMAX > 0 Then
            If Abs(X - XSAV) > Abs(DXSAV) Then
                If KOUNT < KMAX - 1 Then
                    KOUNT = KOUNT + 1
                    XP(KOUNT) = X
                    For I = 1 To NVAR
                        YP(I, KOUNT) = Y(I)
                    Next I
                    XSAV = X
                End If
            End If
        End If
        If (X + H - X2) * (X + H - X1) > ZERO Then H = X2 - X
        Call RKQC(Y(), DYDX(), NVAR, X, H, EPS, YSCAL(), HDID, HNEXT)
        If HDID = H Then
            NOK = NOK + 1
        Else
            NBAD = NBAD + 1
        End If
        If (X - X2) * (X2 - X1) >= ZERO Then
            For I = 1 To NVAR
                YSTART(I) = Y(I)
            Next I
            If KMAX <> 0 Then
                KOUNT = KOUNT + 1
                XP(KOUNT) = X
                For I = 1 To NVAR
                    YP(I, KOUNT) = Y(I)
                Next I
            End If
            Erase DYDX, Y, YSCAL
            Exit Sub
        End If
        If Abs(HNEXT) < HMIN Then
            Print " Stepsize smaller than minimum."
            Exit Sub
        End If
        H = HNEXT
    Next NSTP
    Print " Too many steps."
End Sub
Sub RKQC(Y(), DYDX(), N, X, HTRY, EPS, YSCAL(), HDID, HNEXT)
      ONE = 1#
      SAFETY = 0.9
      ERRCON = 0.0006
      FCOR = 0.066667
      Dim YTEMP(10), YSAV(10), DYSAV(10)
      PGROW = -0.2
      PSHRNK = -0.25
      XSAV = X
      For I = 1 To N
          YSAV(I) = Y(I)
          DYSAV(I) = DYDX(I)
      Next I
      H = HTRY
      Do
        HH = 0.5 * H
        Call RK4(YSAV(), DYSAV(), N, XSAV, HH, YTEMP())
        X = XSAV + HH
        Call DERIVS(X, YTEMP(), DYDX())
        Call RK4(YTEMP(), DYDX(), N, X, HH, Y())
        X = XSAV + H
        If X = XSAV Then
            Print " Stepsize not significant in RKQC."
            Exit Sub
        End If
        Call RK4(YSAV(), DYSAV(), N, XSAV, H, YTEMP())
        ERRMAX = 0#
        For I = 1 To N
            YTEMP(I) = Y(I) - YTEMP(I)
            If ERRMAX < Abs(YTEMP(I) / YSCAL(I)) Then
                ERRMAX = Abs(YTEMP(I) / YSCAL(I))
            End If
        Next I
        ERRMAX = ERRMAX / EPS
        If ERRMAX > ONE Then
            H = SAFETY * H * (ERRMAX ^ PSHRNK)
            FLAG = 1
        Else
            HDID = H
            If ERRMAX > ERRCON Then
                HNEXT = SAFETY * H * (ERRMAX ^ PGROW)
            Else
                HNEXT = 4# * H
            End If
            FLAG = 0
        End If
      Loop While FLAG = 1
      For I = 1 To N
          Y(I) = Y(I) + YTEMP(I) * FCOR
      Next I
      Erase DYSAV, YSAV, YTEMP
End Sub
Sub RK4(Y(), DYDX(), N, X, H, YOUT())
      Dim YT(10), DYT(10), DYM(10)
      HH = H * 0.5
      H6 = H / 6#
      XH = X + HH
      For I = 1 To N
          YT(I) = Y(I) + HH * DYDX(I)
      Next I
      Call DERIVS(XH, YT(), DYT())
      For I = 1 To N
          YT(I) = Y(I) + HH * DYT(I)
      Next I
      Call DERIVS(XH, YT(), DYM())
      For I = 1 To N
          YT(I) = Y(I) + H * DYM(I)
          DYM(I) = DYT(I) + DYM(I)
      Next I
      Call DERIVS(X + H, YT(), DYT())
      For I = 1 To N
          YOUT(I) = Y(I) + H6 * (DYDX(I) + DYT(I) + 2# * DYM(I))
      Next I
      Erase DYM, DYT, YT
End Sub
Sub LUBKSB(A(), N, INDX(), B())
    II = 0
    For I = 1 To N
        LL = INDX(I)
        Sum = B(LL)
        B(LL) = B(I)
        If II <> 0 Then
            For J = II To I - 1
                Sum = Sum - A(I, J) * B(J)
            Next J
        ElseIf Sum <> 0# Then
            II = I
        End If
        B(I) = Sum
    Next I
    For I = N To 1 Step -1
        Sum = B(I)
        If I < N Then
            For J = I + 1 To N
                Sum = Sum - A(I, J) * B(J)
            Next J
        End If
        B(I) = Sum / A(I, I)
    Next I
End Sub
Sub LUDCMP(A(), N, INDX(), D)
    NMAX = 100
    TINY = 1E-20
    Dim VV(100)
    D = 1#
    For I = 1 To N
        AAMAX = 0#
        For J = 1 To N
            If Abs(A(I, J)) > AAMAX Then AAMAX = Abs(A(I, J))
        Next J
        If AAMAX = 0# Then Print "Singular matrix."
        VV(I) = 1# / AAMAX
    Next I
    For J = 1 To N
        If J > 1 Then
            For I = 1 To J - 1
                Sum = A(I, J)
                If I > 1 Then
                    For K = 1 To I - 1
                        Sum = Sum - A(I, K) * A(K, J)
                    Next K
                    A(I, J) = Sum
                End If
            Next I
        End If
        AAMAX = 0#
        For I = J To N
            Sum = A(I, J)
            If J > 1 Then
                For K = 1 To J - 1
                    Sum = Sum - A(I, K) * A(K, J)
                Next K
                A(I, J) = Sum
            End If
            DUM = VV(I) * Abs(Sum)
            If DUM >= AAMAX Then
                IMAX = I
                AAMAX = DUM
            End If
        Next I
        If J <> IMAX Then
            For K = 1 To N
                DUM = A(IMAX, K)
                A(IMAX, K) = A(J, K)
                A(J, K) = DUM
            Next K
            D = -D
            VV(IMAX) = VV(J)
        End If
        INDX(J) = IMAX
        If J <> N Then
            If A(J, J) = 0# Then A(J, J) = TINY
            DUM = 1# / A(J, J)
            For I = J + 1 To N
                A(I, J) = A(I, J) * DUM
            Next I
        End If
    Next J
    If A(N, N) = 0# Then A(N, N) = TINY
End Sub