Pendula

[FellippeHeitor]

Pendula of varying mass, hanging by strings of varying lengths.

Inspired by watching this post on reddit.

Code: QB64: [Select]

1. OPTION _EXPLICIT
2.  
3. SCREEN _NEWIMAGE(1000, 600, 32)
4.  
5. CONST maxPendulum = 15
6. DIM g, i AS INTEGER
7. DIM theta(maxPendulum), accel(maxPendulum), speed(maxPendulum)
8. DIM length(maxPendulum), mass(maxPendulum), c(maxPendulum) AS _UNSIGNED LONG
9. DIM px, py, bx, by
10. g = 9.81
11. FOR i = 0 TO maxPendulum
12.     theta(i) = _PI / 1.5
13.     speed(i) = 0
14.     length(i) = 100 + i * 30
15.     mass(i) = 1 - i * .01
16.     c(i) = _RGB32(255 * RND, 255 * RND, 255 * RND)
17. NEXT
18.  
19. px = _WIDTH / 2
20. py = 0
21.  
22. DO UNTIL _SCREENEXISTS: _LIMIT 1: LOOP
23. _TITLE "Pendula"
24.  
25. DO
26.     CLS
27.     WHILE _MOUSEINPUT: WEND
28.     FOR i = 0 TO maxPendulum
29.         IF _MOUSEBUTTON(1) THEN
30.             theta(i) = map(_MOUSEX, 0, _WIDTH - 1, 4.7, 1.57)
31.             speed(i) = 0
32.         ELSE
33.             accel(i) = mass(i) * g * SIN(theta(i)) / 100
34.             speed(i) = speed(i) + accel(i) / 100
35.             theta(i) = theta(i) + speed(i)
36.         END IF
37.         bx = px + length(i) * SIN(theta(i))
38.         by = py - length(i) * COS(theta(i))
39.         LINE (px, py)-(bx, by), c(i)
40.         CircleFill bx, by, map(i, 0, maxPendulum, 20, 40), c(i)
41.     NEXT
42.     _DISPLAY
43.     _LIMIT 60
44. LOOP
45.  
46. FUNCTION map! (value!, minRange!, maxRange!, newMinRange!, newMaxRange!)
47.     map! = ((value! - minRange!) / (maxRange! - minRange!)) * (newMaxRange! - newMinRange!) + newMinRange!
48. END FUNCTION
49.  
50. SUB CircleFill (CX AS INTEGER, CY AS INTEGER, R AS INTEGER, C AS _UNSIGNED LONG)
51.     ' CX = center x coordinate
52.     ' CY = center y coordinate
53.     '  R = radius
54.     '  C = fill color
55.     DIM Radius AS INTEGER, RadiusError AS INTEGER
56.     DIM X AS INTEGER, Y AS INTEGER
57.     Radius = ABS(R)
58.     RadiusError = -Radius
59.     X = Radius
60.     Y = 0
61.     IF Radius = 0 THEN PSET (CX, CY), C: EXIT SUB
62.     LINE (CX - X, CY)-(CX + X, CY), C, BF
63.     WHILE X > Y
64.         RadiusError = RadiusError + Y * 2 + 1
65.         IF RadiusError >= 0 THEN
66.             IF X <> Y + 1 THEN
67.                 LINE (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
68.                 LINE (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
69.             END IF
70.             X = X - 1
71.             RadiusError = RadiusError - X * 2
72.         END IF
73.         Y = Y + 1
74.         LINE (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
75.         LINE (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
76.     WEND
77. END SUB
78.  

PS: You can use your mouse.

[bplus]

Nice effect with many pendulum.

[SpriggsySpriggs]

@FellippeHeitor This blows my mind away!!!!

EDIT: Gasp! You have the map function! That is the function I needed to make a clock in p5js!

[FellippeHeitor]

@bplus thank you!

@SpriggsySpriggs here's mine and Ashish's translation of p5.js into p5js.bas: https://github.com/FellippeHeitor/p5js.bas <-- grab what you need!

[SpriggsySpriggs]

@FellippeHeitor I do have p5js.bas but I didn't know what to do with the map function so it was good to see this example!

[FellippeHeitor]

Map() is really useful. Highly recommended video on the matter:

[OldMoses]

Great effect. I like that map function too. I seem to be doing things like it fairly often and have always been reinventing the wheel when I do.

[bplus]

Yeah that map function is probably particularly handy when one or the other ranges don't start at 0.

[STxAxTIC]

Are those pendulums connected in this simulation or do they hang separate, but are just plotted closely?

If not connected, meaning each swings free, then the varying mass should not matter.

Angular frequency equals square root of gravity constant divided by pendulum length, mass cancels out. Try changing your masses to random numbers, the output should be identical.

If they are connected then nevermind, different equation.

[FellippeHeitor]

They're... just... there... this is as close as I could get from the original post.

I tried using the pendulum code you'd provided for my pendulum game, but then there's the extra tweaks that got added for the game itself.

I copied the core calculations in this one from the rosetta stone page on pendulum animations. And tweaked it from there.

[STxAxTIC]

Say, remember my post on integrators and  that pdf on how to "compute anything"? Just before covid... I will dig that up if nobody links it sooner, but it fleshes this out pretty well.

[OldMoses]

Yeah that map function is probably particularly handy when one or the other ranges don't start at 0.

Indeed it is, I define a WINDOW that is 620 x 620 actual, but -1000 to 1000 left to right X and 1000 to -1000 top to bottom Y and starts at upper left corner 19,560 on the main screen. It was a PITA calculation to convert _MOUSEX and _MOUSEY, but map walked in and took over like a boss.

[STxAxTIC]

Stopped home in the middle of work just to dig out that code. Here is every right way and wrong way to do a pendulum:

Code: QB64: [Select]

1. CONST Black = _RGB32(0, 0, 0)
2. CONST Blue = _RGB32(0, 0, 255)
3. CONST Gray = _RGB32(128, 128, 128)
4. CONST Green = _RGB32(0, 128, 0)
5. CONST Red = _RGB32(255, 0, 0)
6. CONST White = _RGB32(255, 255, 255)
7. CONST Yellow = _RGB32(255, 255, 0)
8.  
9. SCREEN _NEWIMAGE(1280, 480, 32)
10.  
11. showmenu:
12. CLS
13. _KEYCLEAR
14. COLOR White
15. PRINT " Type a problem number and press ENTER."
16. PRINT
17. PRINT " PROBLEM"
18. PRINT " 1) Particle in r^-1 central potential"
19. PRINT "    a) Perturbed motion ......................... Stable, Symplectic"
20. PRINT "    b) Elliptical motion ........................ Stable, Symplectic"
21. PRINT "    c) Eccentric elliptical motion .............. Stable, Symplectic"
22. PRINT " 2) Particle in r^-2 central potential"
23. PRINT "    a) Attempted circular motion ................ Unstable, Symplectic"
24. PRINT " 3) Mass on a spring"
25. PRINT "    a) Initial Momentum ..........*Incorrect*.... Unstable, Euler"
26. PRINT "    b) Initial Displacement ......*Incorrect*.... Unstable, Euler"
27. PRINT " 4) Two equal masses connected by three springs"
28. PRINT "    a) Symmetric mode ........................... Stable, Symplectic"
29. PRINT "    b) Antisymmetric mode ....................... Stable, Symplectic"
30. PRINT "    c) Perturbed motion ......................... Stable, Symplectic"
31. PRINT " 5) Plane pendulum (Part I)"
32. PRINT "    a) Initial Momentum ..........*Incorrect*.... Unstable, Euler"
33. PRINT " 6) Plane pendulum (Part III)"
34. PRINT "    a) Sub-critical case ........................ Stable, RK4"
35. PRINT "    b) Over-critical case ....................... Stable, RK4"
36. PRINT " 7) Plane pendulum (Part III)"
37. PRINT "    a) Initial Displacement ..................... Stable, Symplectic"
38. PRINT "    b) Sub-critical case ........................ Stable, Symplectic"
39. PRINT "    c) Over-critical case ....................... Stable, Symplectic"
40. PRINT " 8) Plane pendulum (Part IV)"
41. PRINT "    a) Initial Displacement ..................... Stable, Mixed Euler"; ""
42. PRINT
43. INPUT " Enter a problem (ex. 1a): ", problem$
44. problem$ = LTRIM$(RTRIM$(LCASE$(problem$)))
45. CLS
46.  
47. ' Initial conditions
48. q10 = 0
49. p10 = 0
50. q20 = 0
51. p20 = 0
52. SELECT CASE VAL(LEFT$(problem$, 1))
53.     CASE 1
54.         dt = .003
55.         cyclic = 1
56.         delayfactor = 1000
57.         scalebig = 20
58.         scalesmall = dt * 5
59.         m1 = 1
60.         m2 = 1
61.         g = 1
62.         SELECT CASE (RIGHT$(problem$, 1))
63.             CASE "a"
64.                 q10 = 1: p10 = -.2: q20 = 0: p20 = 1
65.             CASE "b"
66.                 q10 = .5: p10 = 1: q20 = -.5: p20 = 1
67.             CASE "c"
68.                 q10 = 1: p10 = .5: q20 = 0: p20 = 1.15
69.         END SELECT
70.     CASE 2
71.         dt = .001
72.         cyclic = 0
73.         delayfactor = 1000
74.         scalebig = 20
75.         scalesmall = dt * 10
76.         m1 = 1
77.         m2 = 1
78.         g = 1
79.         SELECT CASE (RIGHT$(problem$, 1))
80.             CASE "a"
81.                 q10 = 1: p10 = 0: q20 = 0: p20 = 1.22437
82.         END SELECT
83.     CASE 3
84.         dt = .03
85.         cyclic = 0
86.         delayfactor = 1000
87.         scalebig = 10
88.         scalesmall = dt * 7.5
89.         m = 1
90.         k = 1
91.         SELECT CASE (RIGHT$(problem$, 1))
92.             CASE "a"
93.                 q10 = 0: p10 = 2
94.             CASE "b"
95.                 q10 = 2: p10 = 0
96.         END SELECT
97.     CASE 4
98.         dt = .003
99.         cyclic = 1
100.         delayfactor = 50000
101.         scalebig = 12
102.         scalesmall = dt * 7.5
103.         m = 1
104.         k = 1
105.         SELECT CASE (RIGHT$(problem$, 1))
106.             CASE "a"
107.                 q10 = 0: p10 = 2: q20 = 0: p20 = 2
108.             CASE "b"
109.                 q10 = 2: p10 = 0: q20 = -2: p20 = 0
110.             CASE "c"
111.                 q10 = 2: p10 = 0: q20 = 0: p20 = 0
112.         END SELECT
113.     CASE 5
114.         dt = .03
115.         cyclic = 1
116.         delayfactor = 500000
117.         scalebig = 10
118.         scalesmall = dt * 5
119.         m = 1
120.         g = 1
121.         l = 1
122.         SELECT CASE (RIGHT$(problem$, 1))
123.             CASE "a"
124.                 q10 = 0: p10 = 1.5: scalebig = 10
125.         END SELECT
126.     CASE 6
127.         dt = .03
128.         cyclic = 1
129.         delayfactor = 500000
130.         scalebig = 10
131.         scalesmall = dt * 5
132.         m = 1
133.         g = 1
134.         l = 1
135.         SELECT CASE (RIGHT$(problem$, 1))
136.             CASE "a"
137.                 q10 = 0: p10 = 2.19
138.             CASE "b"
139.                 q10 = 0: p10 = 2.25
140.         END SELECT
141.     CASE 7, 8
142.         dt = .03
143.         cyclic = 1
144.         delayfactor = 500000
145.         scalebig = 10
146.         scalesmall = dt * 5
147.         m = 1
148.         g = 1
149.         l = 1
150.         SELECT CASE (RIGHT$(problem$, 1))
151.             CASE "a"
152.                 q10 = 0: p10 = -1.5
153.             CASE "b"
154.                 q10 = 0: p10 = 1.999
155.             CASE "c"
156.                 q10 = 0: p10 = 2.001
157.         END SELECT
158.     CASE ELSE
159.         GOTO showmenu
160. END SELECT
161.  
162. GOSUB drawaxes
163.  
164. ' Main loop.
165. iterations = 0
166. q1 = q10
167. p1 = p10
168. q2 = q20
169. p2 = p20
170.  
171. DO
172.  
173.     FOR thedelay = 0 TO delayfactor: NEXT
174.  
175.     q1temp = q1
176.     p1temp = p1
177.     q2temp = q2
178.     p2temp = p2
179.  
180.     SELECT CASE VAL(LEFT$(problem$, 1))
181.         CASE 1
182.             ' Particle in r^-1 central potential - Symplectic integrator
183.             q1 = q1temp + dt * (p1temp / m2)
184.             q2 = q2temp + dt * (p2temp / m2)
185.             p1 = p1temp - dt * g * m1 * m2 * (q1 / ((q1 ^ 2 + q2 ^ 2) ^ (3 / 2)))
186.             p2 = p2temp - dt * g * m1 * m2 * (q2 / ((q1 ^ 2 + q2 ^ 2) ^ (3 / 2)))
187.         CASE 2
188.             ' Particle in r^-2 central potential - Symplectic integrator
189.             q1 = q1temp + dt * (p1temp / m2)
190.             q2 = q2temp + dt * (p2temp / m2)
191.             p1 = p1temp - dt * g * m1 * m2 * ((3 / 2) * q1 / ((q1 ^ 2 + q2 ^ 2) ^ (5 / 2)))
192.             p2 = p2temp - dt * g * m1 * m2 * ((3 / 2) * q2 / ((q1 ^ 2 + q2 ^ 2) ^ (5 / 2)))
193.         CASE 3
194.             ' Mass on a spring - Forward Euler method
195.             q1 = q1temp + (p1temp / m) * dt
196.             p1 = p1temp - (q1temp * k) * dt
197.             ' Mass on a spring - Backward Euler method
198.             q2 = (q2temp + (p2temp / m) * dt) / (1 + dt ^ 2)
199.             p2 = (p2temp - (q2temp * k) * dt) / (1 + dt ^ 2)
200.         CASE 4
201.             ' Two equal masses connected by three springs - Symplectic integrator
202.             q1 = q1temp + m * (p1temp) * dt
203.             p1 = p1temp - dt * k * (2 * (q1temp + m * (p1temp) * dt) - (q2temp + m * (p2temp) * dt))
204.             q2 = q2temp + m * (p2temp) * dt
205.             p2 = p2temp - dt * k * (2 * (q2temp + m * (p2temp) * dt) - (q1temp + m * (p1temp) * dt))
206.         CASE 5
207.             ' Plane pendulum - Forward Euler method
208.             q1 = q1temp + (p1temp / m) * dt
209.             p1 = p1temp - (g / l) * SIN(q1temp) * dt
210.         CASE 6
211.             ' Plane pendulum - Runge-Kutta 4th
212.             k1w = -(g / l) * SIN(q1temp)
213.             k1t = p1temp
214.             w2 = p1temp + k1w * dt / 2
215.             t2 = q1temp + k1t * dt / 2
216.             k2w = -(g / l) * SIN(t2)
217.             k2t = w2
218.             w3 = p1temp + k2w * dt / 2
219.             t3 = q1temp + k2t * dt / 2
220.             k3w = -(g / l) * SIN(t3)
221.             k3t = w3
222.             w4 = p1temp + k3w * dt
223.             t4 = q1temp + k3t * dt
224.             k4w = -(g / l) * SIN(t4)
225.             dwdt = (k1w + 2 * k2w + 2 * k3w + k4w) / 6
226.             dtdt = (k1t + 2 * k2t + 2 * k3t + k4t) / 6
227.             p1 = p1temp + dwdt * dt
228.             q1 = q1temp + dtdt * dt
229.         CASE 7
230.             ' Plane pendulum - Symplectic integrator
231.             q1 = q1temp + dt * (p1temp / m)
232.             p1 = p1temp - dt * (g / l) * (SIN(q1))
233.         CASE 8
234.             ' Plane pendulum - Modified Euler method
235.             q1 = q1temp + (p1temp / m) * dt
236.             p1 = p1temp - (g / l) * SIN(q1) * dt
237.     END SELECT
238.  
239.     x = (iterations * scalesmall - 180)
240.     IF (x <= 80) THEN
241.         x = x + (320 - _WIDTH / 2)
242.         CALL cpset(x, (q1 * scalebig + 160), Yellow) ' q1 plot
243.         CALL cpset(x, (p1 * scalebig + 60), Yellow) ' p1 plot
244.         CALL cpset(x, (q2 * scalebig - 60), Red) ' q2 plot
245.         CALL cpset(x, (p2 * scalebig - 160), Red) ' p2 plot
246.     ELSE
247.         IF (cyclic = 1) THEN
248.             iterations = 0
249.             PAINT (1, 1), _RGBA(0, 0, 0, 100)
250.             GOSUB drawaxes
251.         ELSE
252.             EXIT DO
253.         END IF
254.     END IF
255.  
256.     ' Phase portrait
257.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1temp * (2 * scalebig) + 100), Yellow)
258.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p1 * (2 * scalebig) + 100), Gray)
259.     CALL cpset((q2temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2temp * (2 * scalebig) + 100), Red)
260.     CALL cpset((q2 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (p2 * (2 * scalebig) + 100), White)
261.  
262.     ' Position portrait
263.     CALL cpset((q1temp * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2temp * (2 * scalebig) - 100), Blue)
264.     CALL cpset((q1 * (2 * scalebig) + (190 + (320 - _WIDTH / 2))), (q2 * (2 * scalebig) - 100), White)
265.  
266.     ' System portrait - Requires 2x wide window.
267.     SELECT CASE VAL(LEFT$(problem$, 1))
268.         CASE 4
269.             CALL cline(160, 0, 220 + 20 * q1, 0, Green)
270.             CALL cline(220 + 20 * q1, 0, 380 + 20 * q2, 0, Yellow)
271.             CALL cline(380 + 20 * q2, 0, 440, 0, Green)
272.             CALL ccircle(220 + 20 * q1temp, 0, 15, Black)
273.             CALL ccircle(220 + 20 * q1, 0, 15, Blue)
274.             CALL ccircle(380 + 20 * q2temp, 0, 15, Black)
275.             CALL ccircle(380 + 20 * q2, 0, 15, Red)
276.         CASE 5, 6, 7, 8
277.             CALL cline(200, 0, 400, 0, White)
278.             CALL cline(300, 0, 300 + 100 * SIN(q1temp), -100 * COS(q1temp), Black)
279.             CALL cline(300, 0, 300 + 100 * SIN(q1), -100 * COS(q1), Blue)
280.     END SELECT
281.  
282.     iterations = iterations + 1
283.  
284.     '_DISPLAY
285.  
286. LOOP UNTIL INKEY$ <> ""
287. DO: LOOP UNTIL INKEY$ <> ""
288. GOTO showmenu
289.  
290. END
291.  
292. drawaxes:
293. ' Axis for q1 plot
294. COLOR White
295. LOCATE 2, 3: PRINT "Generalized"
296. LOCATE 3, 3: PRINT "Coordinates"
297. CALL cline(-_WIDTH / 2 + 140, 160, -_WIDTH / 2 + 400, 160, Gray)
298. COLOR White: LOCATE 5, 3: PRINT "Position q1(t)"
299. COLOR Gray: LOCATE 6, 3: PRINT "q1(0) ="; q10
300. ' Axis for p1 plot
301. CALL cline(-_WIDTH / 2 + 140, 60, -_WIDTH / 2 + 400, 60, Gray)
302. COLOR White: LOCATE 11, 3: PRINT "Momentum p1(t)"
303. COLOR Gray: LOCATE 12, 3: PRINT "p1(0) ="; p10
304. ' Axis for q2 plot
305. CALL cline(-_WIDTH / 2 + 140, -60, -_WIDTH / 2 + 400, -60, Gray)
306. COLOR White: LOCATE 18, 3: PRINT "Position q2(t)"
307. COLOR Gray: LOCATE 19, 3: PRINT "q2(0) ="; q20
308. ' Axis for p2 plot
309. CALL cline(-_WIDTH / 2 + 140, -160, -_WIDTH / 2 + 400, -160, Gray)
310. COLOR White: LOCATE 25, 3: PRINT "Momentum p2(t)"
311. COLOR Gray: LOCATE 26, 3: PRINT "p2(0) ="; p20
312. ' Axes for q & p plots
313. COLOR White
314. LOCATE 2, 60: PRINT "Phase and"
315. LOCATE 3, 58: PRINT "Position Plots"
316. LOCATE 4, 66: PRINT "p"
317. LOCATE 10, 76: PRINT "q"
318. CALL cline((190 + (320 - _WIDTH / 2)), 80 + 100, (190 + (320 - _WIDTH / 2)), -80 + 100, Gray)
319. CALL cline((190 + (320 - _WIDTH / 2)) - 100, 100, (190 + (320 - _WIDTH / 2)) + 100, 100, Gray)
320. LOCATE 17, 66: PRINT "q2"
321. LOCATE 23, 75: PRINT "q1"
322. CALL cline((190 + (320 - _WIDTH / 2)), -80 - 100, (190 + (320 - _WIDTH / 2)), 80 - 100, Gray)
323. CALL cline((190 + (320 - _WIDTH / 2)) - 100, -100, (190 + (320 - _WIDTH / 2)) + 100, -100, Gray)
324. RETURN
325.  
326. SUB cline (x1, y1, x2, y2, col AS _UNSIGNED LONG)
327.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
328. END SUB
329.  
330. SUB ccircle (x1, y1, rad, col AS _UNSIGNED LONG)
331.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), rad, col
332. END SUB
333.  
334. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
335.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
336. END SUB
337.  

... And the accompanying document is a tad calculus heavy but justifies everything in the above:

http://wfbarnes.net/homedata/Mechalculus.pdf

[_vince]

I'm now extremely curious to see a nth order generalization to the chaotic double pendulum -- too much work for me

[FellippeHeitor]

Stopped home in the middle of work just to dig out that code. Here is every right way and wrong way to do a pendulum:

<code>

... And the accompanying document is a tad calculus heavy but justifies everything in the above:

http://wfbarnes.net/homedata/Mechalculus.pdf

That is always a fun render to watch, STx.