Epicycles - bobalooie - 10-29-2022
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
RE: Epicycles - triggered - 10-29-2022
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?
RE: Epicycles - bplus - 10-29-2022
This, epicycles, is same as this:
https://staging.qb64phoenix.com/showthread.php?tid=700&pid=7499#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!
RE: Epicycles - Kernelpanic - 10-29-2022
(10-29-2022, 04:17 PM)triggered Wrote: 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
RE: Epicycles - bobalooie - 10-29-2022
(10-29-2022, 04:17 PM)triggered Wrote: 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.
RE: Epicycles - bobalooie - 10-29-2022
(10-29-2022, 04:22 PM)bplus Wrote: This, epicycles, is same as this:
https://staging.qb64phoenix.com/showthread.php?tid=700&pid=7499#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.
RE: Epicycles - triggered - 10-29-2022
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
RE: Epicycles - vince - 10-29-2022
so impressive, Sheldon
RE: Epicycles - bplus - 10-29-2022
Thumbs up to both bobalooie and triggered! Nice topic and code demos both.
|