Public Class Form1

    Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
        ' Declare all the variables that we're going to use
        Dim Mass_Sun, Mass_Earth, Radius_Earth, Radius_Sun As Double
        Dim X, Y, R, Vel_X, Vel_Y As Double
        ' Dim Accel, Accel_X, Accel_Y As Double
        Dim Step_T As Double
        Dim z(4), t As Double

        ' Now define the starting values of these variables
        Mass_Sun = 100 : Radius_Sun = 40   ' sets the masses and radii
        Mass_Earth = 2 : Radius_Earth = 10 ' of the Sun and Earth
        X = 80 : Y = 0                     ' sets the Earth's coords
        Vel_X = txt_Vel_X.Text : Vel_Y = txt_Vel_Y.Text ' read in the velocity components from the text boxes
        Step_T = 1 ' This is the increment in time between each loop

        ' Now draw the Sun
        Call Plot(0, 0, Radius_Sun, Brushes.Orange)

        While True    ' start of loop

            ' calculate the accel and new position of the Earth from Newton's Laws

            R = Math.Sqrt(X ^ 2 + Y ^ 2)

            'Accel = -Mass_Sun / R ^ 2
            'Accel_X = Accel * X / R : Accel_Y = Accel * Y / R

            ' Vel_X = Vel_X + Accel_X * Step_T : Vel_Y = Vel_Y + Accel_Y * Step_T

            ' X = X + Vel_X * Step_T : Y = Y + Vel_Y * Step_T

            z(1) = X : z(2) = Y : z(3) = Vel_X : z(4) = Vel_Y
            Call RUNGE(t, z, 4, Step_T, TextBox1.Text)
            X = z(1) : Y = z(2) : Vel_X = z(3) : Vel_Y = z(4)

            ' Plots the new Earth's position
            Call Plot(X, Y, Radius_Earth, Brushes.Aquamarine)

            Call Delay(1)  ' waits for a moment

            ' Deletes the Earth from screen
            Call Plot(X, Y, Radius_Earth, Brushes.Black)

        End While          ' end of loop

    End Sub
    Private Sub Plot(ByVal x As Double, ByVal y As Double, ByVal rad As Double, ByVal col As System.Drawing.Brush)
        ' This subroutine just plots a circle at the point (x,y)

        Dim xplot, yplot, irad As Integer
        Dim xmin, xmax, ymin, ymax As Double

        xmin = -100 : xmax = +100
        ymin = -100 : ymax = +100

        xplot = (x - xmin) / (xmax - xmin) * PictureBox1.Width
        yplot = (-y - ymin) / (ymax - ymin) * PictureBox1.Height
        ' Note that VB has (0,0) at top left, so y-axis is inverted...
        irad = Int(rad)

        PictureBox1.CreateGraphics.FillEllipse _
                (col, xplot, yplot, irad, irad)
        '               (Brushes.Green, xplot, yplot, 10, 10)

    End Sub
    Private Sub Delay(ByVal t As Integer)
        Dim i As Integer
        For i = 1 To t * 10000000

        Next i
    End Sub

    Sub RUNGE(ByVal t, ByVal z, ByVal N, ByVal h, ByVal J)

        Dim z_Mid(N), dzdt(N), k1(N), k2(N), k3(N), k4(N) As Double

        Dim I As Integer

        Call ders(t, z, N, dzdt)
        For I = 1 To N
            k1(I) = h * dzdt(I)
            z_Mid(I) = z(I) + 0.5 * k1(I)
        Next I

        Call ders(t + 0.5 * h, z_Mid, N, dzdt)
        For I = 1 To N
            k2(I) = h * dzdt(I)
            z_Mid(I) = z(I) + 0.5 * k2(I)
        Next I

        Call ders(t + 0.5 * h, z_Mid, N, dzdt)
        For I = 1 To N
            k3(I) = h * dzdt(I)
            z_Mid(I) = z(I) + k3(I)
        Next I

        Call ders(t + h, z_Mid, N, dzdt)
        For I = 1 To N
            k4(I) = h * dzdt(I)
        Next I

        If J = 1 Then
            For I = 1 To N
                z(I) = z(I) + k1(I)
            Next I
        ElseIf J = 2 Then
            For I = 1 To N
                z(I) = z(I) + k2(I)
            Next I
        ElseIf J = 4 Then
            For I = 1 To N
                z(I) = z(I) + k1(I) / 6 + k2(I) / 3 + k3(I) / 3 + k4(I) / 6
            Next I
        Else
            MsgBox("Order of Runge-Kutta requested is invalid")
            End
        End If

        t = t + h

    End Sub  ' End of Runge Kutta Subroutine

    Sub ders(ByVal t, ByVal z, ByVal N, ByVal dzdt)
        Dim G, Mass_Sun As Double
        ' put all the necessary declarations here...

        Mass_Sun = 100
        G = 1

        dzdt(1) = z(3) ' put the differential equation for z(1) here
        dzdt(2) = z(4) ' ditto for z(2)
        dzdt(3) = -G * ((Mass_Sun * z(1)) / (z(1) ^ 2 + z(2) ^ 2) ^ (3 / 2)) ' ditto for z(3)
        dzdt(4) = -G * ((Mass_Sun * z(2)) / (z(1) ^ 2 + z(2) ^ 2) ^ (3 / 2)) ' ditto for z(4)
    End Sub

End Class