Three swinging sticks kinetic energy

[neil]

I found this code online for a Freebasic swinging pendulum.
It uses the Runge-Kutta method. It's more complicated than mine.

Code: Select all

' FreeBasic swinging pendulum simulation using the Runge-Kutta method

screenres 800, 600, 32
dim shared as double g = 9.81 ' acceleration due to gravity
dim shared as double l = 200 ' length of the pendulum

function f(theta as double, omega as double) as double
    return -g/l * sin(theta) ' differential equation for pendulum motion
end function

sub RungeKutta(theta as double, omega as double, dt as double, ByRef theta_new as double, ByRef omega_new as double)
    dim as double k1_theta, k2_theta, k3_theta, k4_theta
    dim as double k1_omega, k2_omega, k3_omega, k4_omega

    k1_theta = omega
    k1_omega = f(theta, omega)

    k2_theta = omega + 0.5 * dt * k1_omega
    k2_omega = f(theta + 0.5 * dt * k1_theta, omega + 0.5 * dt * k1_omega)

    k3_theta = omega + 0.5 * dt * k2_omega
    k3_omega = f(theta + 0.5 * dt * k2_theta, omega + 0.5 * dt * k2_omega)

    k4_theta = omega + dt * k3_omega
    k4_omega = f(theta + dt * k3_theta, omega + dt * k3_omega)

    theta_new = theta + (dt/6) * (k1_theta + 2*k2_theta + 2*k3_theta + k4_theta)
    omega_new = omega + (dt/6) * (k1_omega + 2*k2_omega + 2*k3_omega + k4_omega)
end sub

dim as double theta = 1 ' initial angle in radians
dim as double omega = 0 ' initial angular velocity

dim as double dt = 0.08 ' time step
dim as double theta_new, omega_new

do
    screenlock
    cls

    circle(400 + l*sin(theta), 300 + l*cos(theta)), 20, rgb(255,255,0),,,,F ' draw pendulum bob
    line(400, 300) - (400 + l*sin(theta), 300 + l*cos(theta)),rgb(255,255,0) ' draw pendulum rod

    RungeKutta(theta, omega, dt, theta_new, omega_new)

    theta = theta_new
    omega = omega_new
   
    screenunlock
    sleep 10
loop until len(inkey)

[neil]

Double pendulum chaotic motion. This double pendulum was quite tricky to get working in Freebasic.
https://en.wikipedia.org/wiki/Double_pendulum

Code: Select all

' double pendulum chaotic motion by neil

const AS Integer W = 660
const AS Integer H = 660
const As Double pi = 3.14159265

screenres W, H, 32

dim as double l1 = 190 ' length of long pendulum
dim as double l2 = 110 ' length of short pendulum
dim as double m1 = 20 ' mass of long pendulum
dim as double m2 = 15 ' mass of short pendulum
dim as double g = 9.81 ' acceleration due to gravity
dim as double th1 = 3 * pi/4 ' initial position of long pendulum
dim as double th2 = pi ' initial position of short pendulum
dim as double th1_v = 0 ' initial angular velocity of long pendulum
dim as double th2_v = 0 ' initial angular velocity of short pendulum

dim as double th1_a = 0 ' acceleratiom of long pendulum
dim as double th2_a = 0 ' acceleratiom of short pendulum
dim as double x1, y1, x2, y2
dim as double dt = 0.20 ' time step
dim as long fps
dim as integer x0, y0

x0 = W / 2
y0 = H / 2

function Regulate(byval MyFps as long, byref fps as long) as long
  static as double timervalue, _lastsleeptime, t3, frames
  frames += 1
  if (timer - t3) >= 1 then t3 = timer : fps = frames : frames = 0
  var sleeptime = _lastsleeptime + ((1 / myfps) - timer + timervalue) * 1000
  if sleeptime < 1 then sleeptime = 1
  _lastsleeptime = sleeptime
  timervalue = timer
  return sleeptime
end function

do
    screenlock
    cls
    
     ' Calculate accelerations
    th1_a = (-g*(2*m1+m2)*sin(th1) - m2*g*sin(th1 - 2*th2) - 2*sin(th1-th2)*m2*(th2_v*th2_v*l2 + th1_v*th1_v*l1*cos(th1-th2))) / (l1*(2*m1+m2-m2*cos(2*th1-2*th2)))
    th2_a = (2*sin(th1-th2)*(th1_v*th1_v*l1*(m1+m2) + g*(m1+m2)*cos(th1) + th2_v*th2_v*l2*m2*cos(th1-th2))) / (l2*(2*m1+m2-m2*cos(2*th1-2*th2)))
    
    ' update angles
    th1_v += th1_a * dt
    th2_v += th2_a * dt
    
    ' Update velocities
    th1 += th1_v * dt
    th2 += th2_v * dt
    
    'update long pendulum position
    x1 = l1 * sin(th1)
    y1 = l1 * cos(th1)
  
   'update short pendulum position
    x2 = x1 + l2 * sin(th2)
    y2 = y1 + l2 * cos(th2)
    
    ' draw pendulums
    circle(X0+x1, Y0+y1), m1, rgb(0,255,0),,,,F
    circle(X0+x2, Y0+y2), m2, rgb(255,255,0),,,,F
    line(X0, Y0) - (X0+x1, Y0+y1), rgb(0,255,0)
    line(X0+x1, Y0+y1) - (X0+x2, Y0+y2), rgb(255,255,0)
    
    'circle for pivot
    circle(X0, Y0), 5, rgb(0,255,255),,,,F
  
    screenunlock
    sleep regulate(60, fps), 1
    
loop until len(inkey)