Epicycles

[bobalooie]

A few years ago a friend of mine and I were talking about epicycles, which had been used in an attempt to explain planetary motion (yes, we are both nerds.) I had decided to experiment with animating epicycle orbits. The original version of this program was written in FreeBASIC, this is my QB64PE translation of that program as an exercise to learn about the graphics capability of QB64PE. (interestingly, the two BASICs have fairly similar graphics facilities.)



Hopefully the comments in the code provide enough explanation of how I approached the problem.

Code: (Select All)

'Program:       Epicycles.bas
'Purpose:       A QB64PE version of Epicycles
'Version:       0.1
'Create Date:   09/23/2022
'Rev Date:      10/28/2022

OPTION _EXPLICIT
CONST PI2 = 6.2831853

TYPE ScreenPoint
    x AS LONG
    y AS LONG
END TYPE

DIM AS INTEGER ix 'general purpose use
DIM AS STRING sx 'general purpose use
DIM SHARED AS LONG lw, lh 'desktop width and height
DIM AS LONG MinX, MinY, MaxX, MaxY 'Cartesian limits of the images
DIM AS LONG r1, r2, r3 'radii of the epicycle circles
DIM AS LONG rot1, rot2, rot3 'rotation direction
DIM AS LONG step1, step2, step3 'rotation speed
DIM AS ScreenPoint sp1, sp2, sp3 'center points of the epicycle circles
DIM AS DOUBLE Angle1, Angle2, Angle3, AngleStep(1 TO 3)
DIM AS LONG lWin1, lWin2, lWin3 'handles for the three images
' The first image is the visible one. The second image plots each successive
' endpoint of the epicycle to build the pattern. The third image is where the
' epicycles are plotted. Put the second image on the third image, draw the
' epicycle on the third image, put the third image on the first image.

'set up the images and coords
lw = _DESKTOPWIDTH
lh = _DESKTOPHEIGHT
lWin1 = _NEWIMAGE(lw, lh, 32)
lWin2 = _NEWIMAGE(lw, lh, 32)
lWin3 = _NEWIMAGE(lw, lh, 32)
MaxX = lw \ 2: MinX = -MaxX
MaxY = lh \ 2: MinY = -MaxY
r1 = r2 = r3 = 0

_DEST lWin1: WINDOW (MinX, MinY)-(MaxX, MaxY) 'set to Cartesian coords
_DEST lWin2: WINDOW (MinX, MinY)-(MaxX, MaxY) 'set to Cartesian coords
_DEST lWin3: WINDOW (MinX, MinY)-(MaxX, MaxY) 'set to Cartesian coords
SCREEN lWin1: _DEST lWin1: _FULLSCREEN

'Get the user input
CLS
COLOR _RGB(255, 0, 255): PRINT "EPICYCLES DEMONSTRATION"
COLOR _RGB(0, 255, 255)
PRINT "Screen resolution is "; lw; " wide x "; lh; " high."

DO
    LOCATE 3, 1: PRINT "                                        "
    LOCATE 3, 1: COLOR _RGB(0, 255, 255)
    PRINT "Number of epicycles (1 or 2)";: INPUT ix
    IF ix = 1 OR ix = 2 THEN EXIT DO
    COLOR _RGB(255, 0, 0): PRINT "Try again please."
LOOP

' show how this works
DrawExamples ix

DO
    LOCATE 5, 1: COLOR _RGB(0, 255, 255)
    PRINT "Enter a value for the main circle radius (1 to "; STR$(MaxY * 0.5); " )"
    PRINT "or 0 to quit:"
    INPUT r1
    IF r1 = 0 THEN END
    IF r1 > 0 AND r1 <= MaxY * 0.5 THEN EXIT DO
    LOCATE 7, 1: PRINT "                 " 'clear the input
    COLOR _RGB(255, 0, 0): PRINT "Try again please.";
LOOP

DO
    LOCATE 8, 1: COLOR _RGB(0, 255, 255)
    PRINT "Rotate clockwise (CW) or counterclockwise (CCW)"
    INPUT sx
    SELECT CASE sx
        CASE "CW", "cw"
            rot1 = -1
            EXIT DO
        CASE "CCW", "ccw"
            rot1 = 1
            EXIT DO
        CASE ELSE
            LOCATE 9, 1: PRINT "                              "
            COLOR _RGB(255, 0, 0): PRINT "Try again please."
    END SELECT
LOOP

DO
    LOCATE 10, 1: COLOR _RGB(0, 255, 255): PRINT "Rotational speed (1, 2 or 3)?"
    INPUT sx
    SELECT CASE sx
        CASE "1"
            step1 = 1
            EXIT DO
        CASE "2"
            step1 = 2
            EXIT DO
        CASE "3"
            step1 = 3
            EXIT DO
        CASE ELSE
            LOCATE 11, 1: PRINT "             " 'clear the input
            COLOR _RGB(255, 0, 0): PRINT "Try again please."
    END SELECT
LOOP

DO
    LOCATE 12, 1: COLOR _RGB(0, 255, 255): PRINT "Enter a value for the orbiting circle radius r2 (0 to "; STR$(MaxY * 0.50); ")"
    INPUT r2
    IF r2 > 0 AND r2 <= (MaxY * 0.50) THEN EXIT DO
    LOCATE 13, 1: PRINT "            "
    COLOR _RGB(255, 0, 0): PRINT "Try again please."
LOOP

DO
    LOCATE 14, 1: COLOR _RGB(0, 255, 255): PRINT "Rotate clockwise (CW) or counterclockwise (CCW)"
    INPUT sx
    SELECT CASE sx
        CASE "CW", "cw":
            rot2 = -1
            EXIT DO
        CASE "CCW", "ccw"
            rot2 = 1
            EXIT DO
        CASE ELSE
            LOCATE 15, 1: PRINT "             " 'clear the input
            COLOR _RGB(255, 0, 0): PRINT "Try again please."
    END SELECT
LOOP

DO
    LOCATE 16, 1: COLOR _RGB(0, 255, 255): PRINT "Rotational speed (1, 2 or 3)?"
    INPUT sx
    SELECT CASE sx
        CASE "1"
            step2 = 1
            EXIT DO
        CASE "2"
            step2 = 2
            EXIT DO
        CASE "3"
            step2 = 3
            EXIT DO
        CASE ELSE
            LOCATE 17, 1: PRINT "             " 'clear the input
            COLOR _RGB(255, 0, 0): PRINT "Try again please."
    END SELECT
LOOP

step3 = step1 'set a default value

IF (ix = 2) THEN
    DO
        LOCATE 18, 1: COLOR _RGB(0, 255, 255)
        PRINT "Enter a value for the orbiting circle radius r2 (0 to "; STR$(MaxY * 0.25); ")"
        INPUT r3
        IF r3 > 0 AND r3 <= (MaxY * 0.25) THEN EXIT DO
        LOCATE 19, 1: PRINT "            "
        COLOR _RGB(255, 0, 0): PRINT "Try again please."
    LOOP

    DO
        LOCATE 20, 1: COLOR _RGB(0, 255, 255): PRINT "Rotate clockwise (CW) or counterclockwise (CCW)"
        INPUT sx
        SELECT CASE sx
            CASE "CW", "cw":
                rot3 = -1
                EXIT DO
            CASE "CCW", "ccw"
                rot3 = 1
                EXIT DO
            CASE ELSE
                LOCATE 21, 1: PRINT "             " 'clear the input
                COLOR _RGB(255, 0, 0): PRINT "Try again please."
        END SELECT
    LOOP
    DO
        LOCATE 22, 1: COLOR _RGB(0, 255, 255): PRINT "Rotational speed (1, 2 or 3)?"
        INPUT sx
        SELECT CASE sx
            CASE "1"
                step3 = 1
                EXIT DO
            CASE "2"
                step3 = 2
                EXIT DO
            CASE "3"
                step3 = 3
                EXIT DO
            CASE ELSE
                LOCATE 23, 1: PRINT "             " 'clear the input
                COLOR _RGB(255, 0, 0): PRINT "Try again please."
        END SELECT
    LOOP

END IF

PRINT "Press any key to begin."
SLEEP

'-- now the fun stuff
'Use the horizontal screen size as the step size to orbit the satellite.
'use the vertical screen size to orbit the epicycle
AngleStep(1) = PI2 / lw
AngleStep(2) = PI2 / lh
AngleStep(3) = AngleStep(1) * 3
Angle1 = Angle2 = Angle3 = 0
_LIMIT 100

'Screen lWin2 tracks the epicycle points, make sure it is cleared
_DEST lWin2: CLS
'Draw a couple axes
LINE (MinX, 0)-(MaxX, 0), _RGB(64, 64, 64)
LINE (0, MinY)-(0, MaxY), _RGB(64, 64, 64)
COLOR _RGB(0, 255, 255): LOCATE 1, 1: PRINT "Press any key to exit."


DO
    WHILE INKEY$ <> "": WEND 'clear the key buffer

    Angle1 = Angle1 + AngleStep(step1) * rot1
    IF Angle1 > PI2 THEN Angle1 = 0 'gone around one full revolution
    FindCirclePoint r1, Angle1, sp1
    Angle2 = Angle2 + AngleStep(step2) * rot2
    IF Angle2 > PI2 THEN Angle2 = 0 'gone around one full revolution
    FindCirclePoint r2, Angle2, sp2
    sp2.x = sp2.x + sp1.x: sp2.y = sp2.y + sp1.y
    Angle3 = Angle3 + AngleStep(step3) * rot3
    IF Angle3 > PI2 THEN Angle3 = 0 'gone around one full revolution
    FindCirclePoint r3, Angle3, sp3
    sp3.x = sp3.x + sp2.x: sp3.y = sp3.y + sp2.y

    'track the epicycle
    _DEST lWin2
    PSET (sp3.x, sp3.y), _RGB(0, 0, 255)
    _PUTIMAGE , lWin2, lWin3

    'draw the epicycles
    _DEST lWin3

    'Circles
    CIRCLE (0, 0), 2, _RGB(255, 255, 255)
    CIRCLE (sp1.x, sp1.y), 2, _RGB(255, 255, 255)
    CIRCLE (sp2.x, sp2.y), 2, _RGB(255, 255, 255)
    CIRCLE (sp3.x, sp3.y), 2, _RGB(255, 255, 255)

    'Radius lines
    LINE (0, 0)-(sp1.x, sp1.y), _RGB(255, 0, 255)
    LINE (sp1.x, sp1.y)-(sp2.x, sp2.y), _RGB(255, 255, 0)
    LINE (sp2.x, sp2.y)-(sp3.x, sp3.y), _RGB(0, 255, 0)

    _PUTIMAGE , lWin3, lWin1

LOOP WHILE INKEY$ = ""

' clean up
_FULLSCREEN _OFF
SCREEN 0: _DEST 0
_FREEIMAGE lWin1: _FREEIMAGE lWin2: _FREEIMAGE lWin3

END
'-------------------------- end of program -----------------------------------

SUB DrawExamples (num AS INTEGER)

    DIM AS ScreenPoint sp1, sp2, sp3
    DIM AS INTEGER ix, iy

    '-- draw the example circles
    CIRCLE (0, 0), 2, _RGB(255, 255, 255) 'center of main circle
    CIRCLE (0, 0), lh \ 4, _RGB(255, 0, 0) 'main circle
    ix = (lh \ 4) * COS(PI2 \ 8)
    sp1.x = ix: sp1.y = ix
    CIRCLE (sp1.x, sp1.y), 2, _RGB(255, 255, 255) 'center of orbiting circle
    CIRCLE (sp1.x, sp1.y), lh \ 6, _RGB(255, 0, 0) 'orbiting circle
    sp2.x = sp1.x + lh \ 6
    sp2.y = sp1.y
    CIRCLE (sp2.x, sp2.y), 2, _RGB(255, 255, 255)
    LINE (0, 0)-(sp1.x, sp1.y), _RGB(255, 0, 255) 'main circle radius
    LINE (sp1.x, sp1.y)-(sp2.x, sp2.y), _RGB(255, 255, 0) 'orbiting circle radius

    IF (num = 2) THEN
        ix = (lh \ 10) * COS(PI2 \ 8)
        sp3.x = sp2.x + ix
        sp3.y = sp2.y - ix
        CIRCLE (sp2.x, sp2.y), _HYPOT(sp3.x - sp2.x, sp3.y - sp2.y), _RGB(255, 0, 0) 'orbiting circle
        CIRCLE (sp3.x, sp3.y), 2, _RGB(255, 255, 255) 'center of orbiting circle
        LINE (sp2.x, sp2.y)-(sp3.x, sp3.y), _RGB(255, 255, 0) 'orbiting circle radius
        ix = (sp3.x \ 8): iy = sp3.y \ 8
        COLOR _RGB(255, 255, 255)
    END IF

END SUB
'-----------------------------------------------------------------------------
SUB FindCirclePoint (r AS INTEGER, a AS DOUBLE, st AS ScreenPoint)
    ' calculate the offset X and Y of a point given a radius and angle
    ' Assume the offset is from 0,0. Add the returned offsets to the previous point.
    ' Angle must be in radians

    st.x = INT(r * COS(a)): st.y = INT(r * SIN(a))
END SUB

[triggered]

Say, this seems pretty cool. Never heard of epicycles. Is there a way we can see a screenshot of this program for those of us who can't install QB64PE on work computers?

[bplus]

This, epicycles, is same as this:
https://qb64phoenix.com/forum/showthread...99#pid7499
Replys #35 #36

Ha, triggered! ;-)) you were drawing really fancy curves with multiple epicycles or pendulum or whatever you want to call it.
Not to mention this: https://qb64.boards.net/thread/24/epicycles
Very fine, BTW!

[Kernelpanic]

Say, this seems pretty cool. Never heard of epicycles. Is there a way we can see a screenshot of this program for those of us who can't install QB64PE on work computers?

Who has ever been interested in astronomy knows the term. The epicycle theory, intended to explain the otherwise inexplicable movements of the planets, was taught up to the times of Copernicus and Kepler. By then, however, this theory had become so complicated that hardly any astronomer really understood it, because inexplicable movements of the planets had to be explained again and again with new epicycles of the epicycle theory.

After Copernicus, Johannes Kepler finally cleared up the epicycle theory -> Kepler's laws. Deferent and epicycle

Curves of planetary motion in geocentric perspective: Epitrochoids

Geozentrisches Weltsystem

[bobalooie]

Say, this seems pretty cool. Never heard of epicycles. Is there a way we can see a screenshot of this program for those of us who can't install QB64PE on work computers?

My program simulates one or two epicycles. The user specifies the radius of the circles, the rotation direction (CW or CCW) and the rotational speed of each segment. The second image shows the result of the settings depicted in the first image. In theory, the number of possible patterns is huge.

I couldn't take screen shots because any keystroke terminates the program. So, photo images is it.

[bobalooie]

This, epicycles, is same as this:
https://qb64phoenix.com/forum/showthread...99#pid7499
Replys #35 #36

Ha, triggered! ;-)) you were drawing really fancy curves with multiple epicycles or pendulum or whatever you want to call it.
Not to mention this: https://qb64.boards.net/thread/24/epicycles
Very fine, BTW!

Same concept, different approach. My program gives the user more control of what sort of pattern is produced. It was actually quite a fun exercise.

[triggered]

Ah cool, I think I get it. Here is the version we had floating around. It lets you click & drag to draw whatever you want. You got some catching up to do ;-)

Code: (Select All)

Option _Explicit

Do Until _ScreenExists: Loop
_Title "Epicycles"

Screen _NewImage(800, 600, 32)

Type Vector
    x As Double
    y As Double
End Type

Dim RunningMode As Integer ' 0=mouse, 1=lissajoux, 2=etc
RunningMode = 0

Dim As Long N, j, kh, i, m
Dim As Double x0, y0, xp, yp, t

Dim As Vector RawDataPoint(0 To 100000)

Do
    Cls
    Line (0, 0)-(_Width, _Height), _RGB(255, 255, 255, 255), BF
    Color _RGB(0, 0, 0), _RGB(255, 255, 255)
    Locate 1, 1
    Print "Click & Drag to draw a curve."
    Line (0, _Height / 2)-(_Width, _Height / 2), _RGB32(0, 0, 0, 255)
    Line (_Width / 2, 0)-(_Width / 2, _Height), _RGB32(0, 0, 0, 255)
    _Display

    N = 0
    If (Command$ = "") Then
        If (RunningMode = 0) Then
            If (N = 0) Then
                N = GatherMousePoints&(RawDataPoint(), 5)
            End If
        End If
        If (RunningMode <> 0) Then
            N = 250
            For j = 0 To N - 1
                RawDataPoint(j).x = 0
                RawDataPoint(j).y = 0
            Next
        End If
    Else
        RunningMode = 0
        Open Command$ For Input As #1
        Do While Not EOF(1)
            Input #1, x0, y0
            RawDataPoint(N).x = x0
            RawDataPoint(N).y = y0
            N = N + 1
        Loop
        Close #1
        Sleep
    End If

    ReDim GivenPoint(0 To N) As Vector
    For j = 0 To N - 1
        GivenPoint(j) = RawDataPoint(j)
    Next

    ReDim As Vector Q(N - 1)
    ReDim As Double rad(N - 1)
    ReDim As Double phase(N - 1)
    ReDim As Long omega(N - 1)
    ReDim As Double urad(N - 1)
    ReDim As Double uphase(N - 1)
    ReDim As Long uomega(N - 1)
    ReDim As Vector CalculatedPath(N)
    ReDim As Vector ProtoPath(N * 60)
    ReDim As Vector PathSegmentsA(N, N)
    ReDim As Vector PathSegmentsB(N, N)

    For j = 0 To N - 1
        omega(j) = j
        If (j > N / 2) Then omega(j) = j - N
        DFT Q(j).x, Q(j).y, GivenPoint(), j
        rad(j) = Sqr(Q(j).x * Q(j).x + Q(j).y * Q(j).y)
        phase(j) = _Atan2(Q(j).y, Q(j).x)
    Next

    If (RunningMode = 1) Then
        For j = 0 To (N - 1)
            rad(j) = 0
            phase(j) = 0
        Next
        Dim aa
        Dim bb
        '''
        aa = 1
        bb = 3
        '''
        rad(aa) = 30
        rad(bb) = 30
        phase(aa) = -_Pi
        phase(bb) = 3 * _Pi / 4
        rad(N - aa) = 30
        rad(N - bb) = 30
        phase(N - aa) = 0
        phase(N - bb) = -3 * _Pi / 4
    End If

    If (RunningMode = 2) Then
        For j = 0 To (N - 1)
            rad(j) = 0
            phase(j) = 0
        Next
        '''
        aa = 1
        bb = 1
        '''
        rad(aa) = 50
        phase(aa) = -_Pi / 2
        rad(N - bb) = 25
        phase(N - bb) = _Pi / 2
    End If

    For j = 0 To N - 1
        uomega(j) = omega(j)
        urad(j) = rad(j)
        uphase(j) = phase(j)
    Next

    Call QuickSort(0, N - 1, rad(), phase(), omega())

    ''''
    'Open "epi-out.txt" For Output As #1
    'For j = 0 To N - 1
    '    Print #1, omega(j), Chr$(9), rad(j), Chr$(9), phase(j)
    'Next
    'Close #1
    ''''
    'Open "epi-uout.txt" For Output As #1
    'For j = 0 To N - 1
    '    Print #1, uomega(j), Chr$(9), urad(j), Chr$(9), uphase(j)
    'Next
    'Close #1
    ''''

    m = 0
    CalculatedPath(0).x = GivenPoint(0).x
    CalculatedPath(0).y = GivenPoint(0).y
    For t = 0 To 2 * _Pi Step 2 * _Pi / N
        x0 = 0
        y0 = 0
        For j = 0 To N - 1
            PathSegmentsA(m, j).x = x0
            PathSegmentsA(m, j).y = y0
            xp = rad(j) * Cos(phase(j) + t * omega(j))
            yp = -rad(j) * Sin(phase(j) + t * omega(j))
            x0 = x0 + xp
            y0 = y0 + yp
            PathSegmentsB(m, j).x = x0
            PathSegmentsB(m, j).y = y0
        Next
        CalculatedPath(m).x = x0
        CalculatedPath(m).y = y0
        m = m + 1
    Next

    i = 0

    _KeyClear
    Do

        m = 0
        For t = 0 To (2 * _Pi) Step (2 * _Pi / (N * 60))
            x0 = 0
            y0 = 0
            For j = 0 To N - 1
                xp = rad(j) * Cos(phase(j) + t * omega(j))
                yp = -rad(j) * Sin(phase(j) + t * omega(j))
                x0 = x0 + xp
                y0 = y0 + yp
                If (j = i) Then
                    ProtoPath(m).x = x0
                    ProtoPath(m).y = y0
                    Exit For
                End If
            Next
            m = m + 1
        Next

        For j = 0 To N - 1
            kh = _KeyHit
            Cls
            Locate 1, 1
            Print "Approximation "; _Trim$(Str$(i)); " of "; _Trim$(Str$(N - 1)); ". Press any key to restart."
            Line (_Width / 2, _Height - 100)-(_Width / 2, _Height - 40), _RGB32(0, 0, 0, 55)
            For m = 0 To N - 1
                Call CCircleF(GivenPoint(m).x, GivenPoint(m).y, 3, _RGB(155, 155, 155, 255))
            Next
            For m = 0 To N - 2
                Call LineSmooth(CalculatedPath(m).x, CalculatedPath(m).y, CalculatedPath(m + 1).x, CalculatedPath(m + 1).y, _RGB32(0, 0, 0, 75))
            Next
            For m = 0 To i
                Call CCircle(PathSegmentsA(j, m).x, PathSegmentsA(j, m).y, rad(m), _RGB32(0, 127, 255, 155))
                Call LineSmooth(PathSegmentsA(j, m).x, PathSegmentsA(j, m).y, PathSegmentsB(j, m).x, PathSegmentsB(j, m).y, _RGB32(28, 28, 255, 155))
            Next
            For m = 0 To (j - 0) * 60
                Call CCircleF(ProtoPath(m).x, ProtoPath(m).y, 1, _RGB32(255, 0, 255 * m / j, 255))
                'CALL LineSmooth(ProtoPath(m).x, ProtoPath(m).y, ProtoPath(m + 60).x, ProtoPath(m + 60).y, _RGB32(255, 0, 255 * m / j, 255))
            Next
            Dim nn
            nn = N
            'IF nn > 100 THEN nn = 100
            For m = 0 To N - 1
                y0 = .9 * _Width / nn
                x0 = uomega(m) * y0 - y0 / 2
                Call CLineBF(x0, -(_Height) / 2 + 40, x0 + y0, -(_Height) / 2 + 40 + 20 * Log(1 + urad(m)), _RGB32(0, 0, 0, 55))
                If (urad(m) > .001) Then
                    Call CLineBF(x0, -(_Height) / 2 + 40, x0 + y0, -(_Height) / 2 + 40 + 10 * (uphase(m)), _RGB32(255, 0, 0, 55))
                End If
            Next
            For m = 0 To i
                y0 = .9 * _Width / nn
                x0 = omega(m) * y0 - y0 / 2
                Call CLineBF(x0, -(_Height) / 2 + 40, x0 + y0, -(_Height) / 2 + 40 + 20 * Log(1 + rad(m)), _RGB32(0, 0, 0, 155))
                If (rad(m) > .001) Then
                    Call CLineBF(x0, -(_Height) / 2 + 40, x0 + y0, -(_Height) / 2 + 40 + 10 * (phase(m)), _RGB32(255, 0, 0, 105))
                End If
            Next

            '_DELAY .5
            _Delay .1 / N
            _Display
            _Limit 60

            If (kh <> 0) Then Exit Do
        Next

        If (i = N - 1) Then
            'i = 0
        Else
            i = i + 1 + Int(Sqr(i))
            If (i >= N) Then i = N - 1
        End If

        If (kh <> 0) Then
            kh = 0
            Exit Do
        End If

        'SLEEP 15
        _Delay 1

    Loop
Loop

Sleep
System

Sub DFT (re As Double, im As Double, arr() As Vector, j0 As Long)
    Dim As Long n, k
    Dim As Double u, v, arg
    n = UBound(arr)
    re = 0
    im = 0
    For k = 0 To n
        arg = 2 * _Pi * k * j0 / n
        cmul u, v, Cos(arg), Sin(arg), arr(k).x, arr(k).y
        re = re + u
        im = im - v
    Next
    re = re / n
    im = im / n
End Sub

Sub cmul (u As Double, v As Double, xx As Double, yy As Double, aa As Double, bb As Double)
    u = xx * aa - yy * bb
    v = xx * bb + yy * aa
End Sub

Sub CCircle (x0 As Double, y0 As Double, rad As Double, shade As _Unsigned Long)
    Circle (_Width / 2 + x0, -y0 + _Height / 2), rad, shade
End Sub

Sub CPset (x0 As Double, y0 As Double, shade As _Unsigned Long)
    PSet (_Width / 2 + x0, -y0 + _Height / 2), shade
End Sub

Sub CLine (x0 As Double, y0 As Double, x1 As Double, y1 As Double, shade As _Unsigned Long)
    Line (_Width / 2 + x0, -y0 + _Height / 2)-(_Width / 2 + x1, -y1 + _Height / 2), shade
End Sub

Sub CLineBF (x0 As Double, y0 As Double, x1 As Double, y1 As Double, shade As _Unsigned Long)
    Line (_Width / 2 + x0, -y0 + _Height / 2)-(_Width / 2 + x1, -y1 + _Height / 2), shade, BF
End Sub

Sub CCircleF (x As Long, y As Long, r As Long, c As Long)
    Dim As Long xx, yy, e
    xx = r
    yy = 0
    e = -r
    Do While (yy < xx)
        If (e <= 0) Then
            yy = yy + 1
            Call CLineBF(x - xx, y + yy, x + xx, y + yy, c)
            Call CLineBF(x - xx, y - yy, x + xx, y - yy, c)
            e = e + 2 * yy
        Else
            Call CLineBF(x - yy, y - xx, x + yy, y - xx, c)
            Call CLineBF(x - yy, y + xx, x + yy, y + xx, c)
            xx = xx - 1
            e = e - 2 * xx
        End If
    Loop
    Call CLineBF(x - r, y, x + r, y, c)
End Sub

Sub LineSmooth (x0 As Single, y0 As Single, x1 As Single, y1 As Single, c As _Unsigned Long)
    ' source: https://en.wikipedia.org/w/index.php?title=Xiaolin_Wu%27s_line_algorithm&oldid=852445548
    ' translated: FellippeHeitor @ qb64.org
    ' bugfixed for alpha channel

    Dim plX As Integer, plY As Integer, plI

    Dim steep As _Byte
    steep = Abs(y1 - y0) > Abs(x1 - x0)

    If steep Then
        Swap x0, y0
        Swap x1, y1
    End If

    If x0 > x1 Then
        Swap x0, x1
        Swap y0, y1
    End If

    Dim dx, dy, gradient
    dx = x1 - x0
    dy = y1 - y0
    gradient = dy / dx

    If dx = 0 Then
        gradient = 1
    End If

    'handle first endpoint
    Dim xend, yend, xgap, xpxl1, ypxl1
    xend = _Round(x0)
    yend = y0 + gradient * (xend - x0)
    xgap = (1 - ((x0 + .5) - Int(x0 + .5)))
    xpxl1 = xend 'this will be used in the main loop
    ypxl1 = Int(yend)
    If steep Then
        plX = ypxl1
        plY = xpxl1
        plI = (1 - (yend - Int(yend))) * xgap
        GoSub plot

        plX = ypxl1 + 1
        plY = xpxl1
        plI = (yend - Int(yend)) * xgap
        GoSub plot
    Else
        plX = xpxl1
        plY = ypxl1
        plI = (1 - (yend - Int(yend))) * xgap
        GoSub plot

        plX = xpxl1
        plY = ypxl1 + 1
        plI = (yend - Int(yend)) * xgap
        GoSub plot
    End If

    Dim intery
    intery = yend + gradient 'first y-intersection for the main loop

    'handle second endpoint
    Dim xpxl2, ypxl2
    xend = _Round(x1)
    yend = y1 + gradient * (xend - x1)
    xgap = ((x1 + .5) - Int(x1 + .5))
    xpxl2 = xend 'this will be used in the main loop
    ypxl2 = Int(yend)
    If steep Then
        plX = ypxl2
        plY = xpxl2
        plI = (1 - (yend - Int(yend))) * xgap
        GoSub plot

        plX = ypxl2 + 1
        plY = xpxl2
        plI = (yend - Int(yend)) * xgap
        GoSub plot
    Else
        plX = xpxl2
        plY = ypxl2
        plI = (1 - (yend - Int(yend))) * xgap
        GoSub plot

        plX = xpxl2
        plY = ypxl2 + 1
        plI = (yend - Int(yend)) * xgap
        GoSub plot
    End If

    'main loop
    Dim x
    If steep Then
        For x = xpxl1 + 1 To xpxl2 - 1
            plX = Int(intery)
            plY = x
            plI = (1 - (intery - Int(intery)))
            GoSub plot

            plX = Int(intery) + 1
            plY = x
            plI = (intery - Int(intery))
            GoSub plot

            intery = intery + gradient
        Next
    Else
        For x = xpxl1 + 1 To xpxl2 - 1
            plX = x
            plY = Int(intery)
            plI = (1 - (intery - Int(intery)))
            GoSub plot

            plX = x
            plY = Int(intery) + 1
            plI = (intery - Int(intery))
            GoSub plot

            intery = intery + gradient
        Next
    End If

    Exit Sub

    plot:
    ' Change to regular PSET for standard coordinate orientation.
    Call CPset(plX, plY, _RGB32(_Red32(c), _Green32(c), _Blue32(c), plI * _Alpha32(c)))
    Return
End Sub

Function GatherMousePoints& (arr() As Vector, res As Double)
    Dim As Long i
    Dim As Double mx, my, xx, yy, delta, xold, yold
    xold = 0
    yold = 0
    i = 0
    Do
        Do While _MouseInput
            mx = _MouseX
            my = _MouseY
            If _MouseButton(1) Then
                xx = mx - (_Width / 2)
                yy = -my + (_Height / 2)
                delta = Sqr((xx - xold) ^ 2 + (yy - yold) ^ 2)
                If (delta > res) Then
                    Call CCircleF(xx, yy, 3, _RGB(0, 0, 0))
                    _Display
                    arr(i).x = xx
                    arr(i).y = yy
                    xold = xx
                    yold = yy
                    i = i + 1
                End If
            End If
        Loop
        If ((i > 2) And (Not _MouseButton(1))) Then Exit Do
        If (i > 999) Then Exit Do
    Loop
    GatherMousePoints& = i
End Function

Sub QuickSort (LowLimit As Long, HighLimit As Long, rad() As Double, phase() As Double, omega() As Long)
    Dim As Long piv
    If (LowLimit < HighLimit) Then
        piv = Partition(LowLimit, HighLimit, rad(), phase(), omega())
        Call QuickSort(LowLimit, piv - 1, rad(), phase(), omega())
        Call QuickSort(piv + 1, HighLimit, rad(), phase(), omega())
    End If
End Sub

Function Partition (LowLimit As Long, HighLimit As Long, rad() As Double, phase() As Double, omega() As Long)
    Dim As Long i, j
    Dim As Double pivot, tmp
    pivot = rad(HighLimit)
    i = LowLimit - 1
    For j = LowLimit To HighLimit - 1
        tmp = rad(j) - pivot
        If (tmp >= 0) Then
            i = i + 1
            Swap rad(i), rad(j)
            Swap phase(i), phase(j)
            Swap omega(i), omega(j)
        End If
    Next
    Swap rad(i + 1), rad(HighLimit)
    Swap phase(i + 1), phase(HighLimit)
    Swap omega(i + 1), omega(HighLimit)
    Partition = i + 1
End Function

[vince]

so impressive, Sheldon

[bplus]

Thumbs up to both bobalooie and triggered! Nice topic and code demos both.