Text Only | Text with Attachments

QB64.org Forum

Active Forums => Programs => Topic started by: FellippeHeitor on September 27, 2020, 04:04:40 pm

Title: Pendula
Post by: FellippeHeitor on September 27, 2020, 04:04:40 pm
Pendula of varying mass, hanging by strings of varying lengths.

Inspired by watching this post on reddit (https://www.reddit.com/r/BetterEveryLoop/comments/j0sdgz/slightly_different_string_lengths_cause_patterns/?utm_source=share&utm_medium=ios_app&utm_name=iossmf).


Code: QB64: [Select]
  2.  
  3. SCREEN _NEWIMAGE(1000, 600, 32)
  4.  
  5. CONST maxPendulum = 15
  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)
  18.  
  19. px = _WIDTH / 2
  20. py = 0
  21.  
  23. _TITLE "Pendula"
  24.  
  26.     CLS
  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
  45.  
  46. FUNCTION map! (value!, minRange!, maxRange!, newMinRange!, newMaxRange!)
  47.     map! = ((value! - minRange!) / (maxRange! - minRange!)) * (newMaxRange! - newMinRange!) + newMinRange!
  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
  78.  

PS: You can use your mouse.
Title: Re: Pendula
Post by: bplus on September 27, 2020, 04:07:27 pm
Nice effect with many pendulum.
Title: Re: Pendula
Post by: SpriggsySpriggs on September 27, 2020, 04:07:38 pm
@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!
Title: Re: Pendula
Post by: FellippeHeitor on September 27, 2020, 04:11:10 pm
@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!
Title: Re: Pendula
Post by: SpriggsySpriggs on September 27, 2020, 04:17:47 pm
@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!
Title: Re: Pendula
Post by: FellippeHeitor on September 27, 2020, 04:19:48 pm
Map() is really useful. Highly recommended video on the matter:

Title: Re: Pendula
Post by: OldMoses on September 27, 2020, 04:37:07 pm
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.
Title: Re: Pendula
Post by: bplus on September 27, 2020, 04:51:57 pm
Yeah that map function is probably particularly handy when one or the other ranges don't start at 0.
Title: Re: Pendula
Post by: STxAxTIC on September 27, 2020, 05:20:45 pm
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.
Title: Re: Pendula
Post by: FellippeHeitor on September 27, 2020, 05:23:04 pm
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.
Title: Re: Pendula
Post by: STxAxTIC on September 27, 2020, 05:27:32 pm
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.
Title: Re: Pendula
Post by: OldMoses on September 27, 2020, 05:29:52 pm
Quote from: bplus on September 27, 2020, 04:51:57 pm
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.
Title: Re: Pendula
Post by: STxAxTIC on September 27, 2020, 05:53:29 pm
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:
  14. COLOR White
  15. PRINT " Type a problem number and press ENTER."
  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"; ""
  43. INPUT " Enter a problem (ex. 1a): ", problem$
  44. problem$ = LTRIM$(RTRIM$(LCASE$(problem$)))
  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
  161.  
  162. GOSUB drawaxes
  163.  
  164. ' Main loop.
  165. iterations = 0
  166. q1 = q10
  167. p1 = p10
  168. q2 = q20
  169. p2 = p20
  170.  
  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.  
  288. GOTO showmenu
  289.  
  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)
  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
  329.  
  330. SUB ccircle (x1, y1, rad, col AS _UNSIGNED LONG)
  331.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), rad, col
  333.  
  334. SUB cpset (x1, y1, col AS _UNSIGNED LONG)
  335.     PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
  337.  

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

http://wfbarnes.net/homedata/Mechalculus.pdf (http://wfbarnes.net/homedata/Mechalculus.pdf)
Title: Re: Pendula
Post by: _vince on September 28, 2020, 04:51:23 pm
I'm now extremely curious to see a nth order generalization to the chaotic double pendulum -- too much work for me
Title: Re: Pendula
Post by: FellippeHeitor on September 28, 2020, 05:26:08 pm
Quote from: STxAxTIC on September 27, 2020, 05:53:29 pm
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 (http://wfbarnes.net/homedata/Mechalculus.pdf)

That is always a fun render to watch, STx.
Title: Re: Pendula
Post by: DANILIN on September 29, 2020, 04:07:06 pm
more  of pendulums: long read text

http://rosettacode.org/wiki/Animate_a_pendulum
rosettacode.org/wiki/Animate_a_pendulum

and me learn another language graphic by this page
Title: Re: Pendula
Post by: Ashish on September 30, 2020, 09:06:40 am
My attempt to recreate this in 3D -

Code: QB64: [Select]
  1. 'Pendula in 3D : By Ashish in QB64 using OpenGL
  2. 'idea from here - https://www.qb64.org/forum/index.php?topic=3056.msg123288#msg123288
  3. _TITLE "Pendula in 3D"
  4. SCREEN _NEWIMAGE(800, 600, 32)
  5.  
  6. TYPE pendulum_type
  7.     ang AS SINGLE 'angle
  8.     ang_factor AS SINGLE
  9.     min_ang AS SINGLE
  10.     max_ang AS SINGLE
  11.     mass AS SINGLE
  12.     clr AS _UNSIGNED LONG 'color of the sphere
  13.     rad AS SINGLE 'radius of the sphere
  14.     rope_length AS SINGLE 'length of the rope on which the sphere is hanging
  16.  
  18.     SUB glutSolidSphere (BYVAL radius AS DOUBLE, BYVAL slices AS LONG, BYVAL stack AS LONG) 'use to draw a sphere
  19.     SUB gluLookAt (BYVAL eyeX#, BYVAL eyeY#, BYVAL eyeZ#, BYVAL centerX#, BYVAL centerY#, BYVAL centerZ#, BYVAL upX#, BYVAL upY#, BYVAL upZ#) 'for camera purpose
  21.  
  22. CONST MAX_PENDULUM = 15
  23.  
  24. DIM SHARED pendulum(1 TO MAX_PENDULUM) AS pendulum_type, init AS _BYTE
  25. DIM SHARED ix, iy, iz
  26.  
  27. ix = 0: iy = 1.5: iz = 2
  28.  
  29. FOR i = 1 TO MAX_PENDULUM
  30.     pendulum(i).rad = 0.25
  31.     pendulum(i).min_ang = _PI(1.25)
  32.     pendulum(i).max_ang = _PI(1.75)
  33.     pendulum(i).clr = _RGB(RND * 255, RND * 255, RND * 255)
  34.     pendulum(i).rope_length = 2
  35.     pendulum(i).mass = 3 + i * 0.5
  37.  
  38. init = 1
  40.     FOR i = 1 TO MAX_PENDULUM
  41.         pendulum(i).ang_factor = pendulum(i).ang_factor + 0.01 * pendulum(i).mass
  42.         pendulum(i).ang = map(SIN(pendulum(i).ang_factor), -1, 1, pendulum(i).min_ang, pendulum(i).max_ang)
  43.     NEXT
  44.     _LIMIT 60
  46.  
  47. SUB _GL ()
  48.     STATIC t
  49.     t = t + 0.01
  50.     IF init = 0 THEN EXIT SUB
  52.    
  53.     ' _glClearColor 0.3,0.3,0.3,1
  54.     _glClear _GL_COLOR_BUFFER_BIT OR _GL_DEPTH_BUFFER_BIT
  55.    
  56.     _glEnable _GL_LIGHTING
  57.     _glEnable _GL_LIGHT0
  58.     _glEnable _GL_DEPTH_TEST
  59.     _glEnable _GL_COLOR_MATERIAL
  60.    
  61.     DIM v(3) AS SINGLE
  62.     v(0) = 0: v(1) = 0: v(2) = 1: v(3) = 0
  63.    
  64.     _glLightfv _GL_LIGHT0, _GL_POSITION, _OFFSET(v())
  65.    
  66.     _glMatrixMode _GL_PROJECTION
  68.     _gluPerspective 50, _WIDTH / _HEIGHT, 0.1, 25
  69.    
  70.     _glMatrixMode _GL_MODELVIEW
  72.    
  73.     gluLookAt 0, 2, 5, 0, 0, 0, 0, 1, 0
  74.    
  75.     drawXZPlane 0, -2, 0, 8, 16, -1
  76.    
  77.     DIM z
  78.     z = iz
  79.     _glLineWidth 5.0
  80.     FOR i = 1 TO MAX_PENDULUM
  81.         _glColor3ub 255, 255, 0
  82.         _glBegin _GL_LINES
  83.         _glVertex3f ix, iy, z
  84.         _glVertex3f ix + pendulum(i).rope_length * COS(pendulum(i).ang), iy + pendulum(i).rope_length * SIN(pendulum(i).ang), z
  85.         _glEnd
  86.        
  87.         _glPushMatrix
  88.         _glTranslatef ix + (pendulum(i).rope_length + pendulum(i).rad) * COS(pendulum(i).ang), iy + (pendulum(i).rope_length + pendulum(i).rad) * SIN(pendulum(i).ang), z
  89.         _glColor3ub _RED(pendulum(i).clr), _GREEN(pendulum(i).clr), _BLUE(pendulum(i).clr)
  90.         glutSolidSphere pendulum(i).rad, 50, 50
  91.         _glPopMatrix
  92.        
  93.         z = z - 2.3 * pendulum(i).rad
  94.     NEXT
  95.     _glColor3ub 255, 255, 0
  96.     _glBegin _GL_LINES
  97.     _glVertex3f ix, iy, iz
  98.     _glVertex3f ix, iy, z
  99.     _glEnd
  100.    
  101.     _glFlush
  103.  
  104. SUB drawXZPlane (x, y, z, w, d, n) '(x,z)->location to draw, w->width, d->depth of the plane, n->normal direction in Y direction
  105.     _glBegin _GL_QUADS
  106.     _glTexCoord2f 1, 0
  107.     _glNormal3f 0, n, 0
  108.     _glVertex3f x - (w / 2), y, z - d / 2
  109.     _glTexCoord2f 1, 1
  110.     _glNormal3f 0, n, 0
  111.     _glVertex3f x + (w / 2), y, z - d / 2
  112.     _glTexCoord2f 0, 1
  113.     _glNormal3f 0, n, 0
  114.     _glVertex3f x + (w / 2), y, z + d / 2
  115.     _glTexCoord2f 0, 0
  116.     _glNormal3f 0, n, 0
  117.     _glVertex3f x - (w / 2), y, z + d / 2
  118.     _glEnd
  120.  
  121. 'from p5js.bas
  122. FUNCTION map! (value!, minRange!, maxRange!, newMinRange!, newMaxRange!)
  123.     map! = ((value! - minRange!) / (maxRange! - minRange!)) * (newMaxRange! - newMinRange!) + newMinRange!
  125.  
  126.  
Title: Re: Pendula
Post by: FellippeHeitor on September 30, 2020, 09:20:35 am
Quote from: Ashish on September 30, 2020, 09:06:40 am
My attempt to recreate this in 3D -

Whoooooooooa!
🤩
Title: Re: Pendula
Post by: DANILIN on January 06, 2022, 12:38:38 am
improvement:

for colorize:
Randomize Timer

thik lines:
_glLineWidth 0.0

or symbol ':
'_glBegin _GL_LINES 
Title: Re: Pendula
Post by: _vince on January 07, 2022, 05:13:58 am
We've all seen the chaotic double pendulum.  I'll leave my simple Runge-Kutta one done with canned formulas from wikipedia at the end.

Challenge:  Generalize it to n-pendulum!  Triple, quadruple, whatever the user wants.  The next ball is to be attached to the previous.  It's going to be a fairly tedious system of differential equations based on the lagrangian of the system (https://en.wikipedia.org/wiki/Double_pendulum#Analysis_and_interpretation (https://en.wikipedia.org/wiki/Double_pendulum#Analysis_and_interpretation)).  Probably so tedious you're going to need a large matrix and some linear algebra.  But goddamn I'll be so impressed if you do it, you'll come out of it totally competent to program everything DEs.

Here's the trivial case:
Code: QB64: [Select]
  1. defsng a-z
  2. pi = 3.141593
  3. r = 100
  4. a1 = -pi/2
  5. a2 = -pi/2
  6. h=0.001
  7. m=1000
  8. g=1000
  9. sw = 640
  10. sh = 480
  11.  
  13. p1 = _newimage(sw,sh,12)
  14. screen _newimage(sw,sh,12)
  15.  
  17.         a1p = (6/(m*r^2))*(2*m*v1-3*cos(a1-a2)*m*v2)/(16-9*(cos(a1-a2))^2)
  18.         a1k1 = a1p
  19.         a1k2 = a1p + h*a1k1/2
  20.         a1k3 = a1p + h*a1k2/2
  21.         a1k4 = a1p + h*a1k3
  22.         a2p = (6/(m*r^2))*(8*m*v2-3*cos(a1-a2)*m*v1)/(16-9*(cos(a1-a2))^2)
  23.         a2k1 = a2p
  24.         a2k2 = a2p + h*a2k1/2
  25.         a2k3 = a2p + h*a2k2/2
  26.         a2k4 = a2p + h*a2k3
  27.         na1=a1+h*(a1k1+2*a1k2+2*a1k3+a1k4)/6
  28.         na2=a2+h*(a2k1+2*a2k2+2*a2k3+a2k4)/6
  29.         v1p = -0.5*r^2*(a1p*a2p*sin(a1-a2)+3*g*sin(a1)/r)
  30.         v1k1 = v1p
  31.         v1k2 = v1p + h*v1k1/2
  32.         v1k3 = v1p + h*v1k2/2
  33.         v1k4 = v1p + h*v1k3
  34.         v2p = -0.5*r^2*(-a1p*a2p*sin(a1-a2)+g*sin(a2)/r)
  35.         v2k1 = v2p
  36.         v2k2 = v2p + h*v2k1/2
  37.         v2k3 = v2p + h*v2k2/2
  38.         v2k4 = v2p + h*v2k3
  39.         nv1=v1+h*(v1k1+2*v1k2+2*v1k3+v1k4)/6
  40.         nv2=v2+h*(v2k1+2*v2k2+2*v2k3+v2k4)/6
  41.         a1=na1
  42.         a2=na2
  43.         v1=nv1
  44.         v2=nv2
  45.  
  46.         _dest p1
  47.         pset (sw/2+r*cos(a1-pi/2)+r*cos(a2-pi/2), sh/2-r*sin(a1-pi/2)-r*sin(a2-pi/2)),12
  48.         _dest 0
  49.         _putimage, p1
  50.  
  51.         line (sw/2, sh/2)-step(r*cos(a1-pi/2), -r*sin(a1-pi/2))
  52.         circle step(0,0),5
  53.         line -step(r*cos(a2-pi/2), -r*sin(a2-pi/2))
  54.         circle step(0,0),5
  55.  
  56.         _limit 200
  57.         _display
  60.  
Title: Re: Pendula
Post by: STxAxTIC on January 07, 2022, 10:14:12 am
That's an interesting challenge. Bonus points if all pendula are constrained to never self-intersect, so that the resulting object should be like a hanging chain and exhibit the same physics. I wonder if this is even *that* terrible to do, becau - ya know what? I'll shut up.... Nice challenge!
Text Only | Text with Attachments