Source code in the nim programming language

# Pendulum simulation.

import math
import times

import opengl
import opengl/glut

var
  # Simulation variables.
  lg: float         # Pendulum length.
  g: float          # Gravity (should be positive).
  currTime: Time    # Current time.
  theta0: float     # Initial angle.
  theta: float      # Current angle.
  omega: float      # Angular velocity = derivative of theta.
  accel: float      # Angular acceleration = derivative of omega.
  e: float          # Total energy.

#---------------------------------------------------------------------------------------------------

proc initSimulation(length, gravitation, start: float) =
  ## Initialize the simulation.

  lg = length
  g = gravitation
  currTime = getTime()
  theta0 = start                    # Initial angle for which omega = 0.
  theta = start
  omega = 0
  accel = -g / lg * sin(theta0)
  e = g * lg * (1 - cos(theta0))    # Total energy = potential energy when starting.

#---------------------------------------------------------------------------------------------------

proc elapsed(): float =
  ## Return the elapsed time since previous call, expressed in seconds.

  let nextTime = getTime()
  result = (nextTime - currTime).inMicroseconds.float / 1e6
  currTime = nextTime

#---------------------------------------------------------------------------------------------------

proc resize(w, h: GLsizei) =
  ## Resize the window.

  glViewport(0, 0, w, h)
  glMatrixMode(GL_PROJECTION)
  glLoadIdentity()

  glMatrixMode(GL_MODELVIEW)
  glLoadIdentity()
  glOrtho(0, GLdouble(w), GLdouble(h), 0, -1, 1)

#---------------------------------------------------------------------------------------------------

proc render() {.cdecl.} =
  ## Render the window.

  # Compute the position of the mass.
  var x = 320 + 300 * sin(theta)
  var y = 300 * cos(theta)

  resize(640, 320)
  glClear(GL_COLOR_BUFFER_BIT)

  # Draw the line from pivot to mass.
  glBegin(GL_LINES)
  glVertex2d(320, 0)
  glVertex2d(x, y)
  glEnd()
  glFlush()

  # Update theta and omega.
  let dt = elapsed()
  theta += (omega + dt * accel / 2) * dt
  omega += accel * dt

  # If, due to computation errors, potential energy is greater than total energy,
  # reset theta to ±theta0 and omega to 0.
  if lg * g * (1 - cos(theta)) >= e:
    theta = sgn(theta).toFloat * theta0
    omega = 0

  accel = -g / lg * sin(theta)

#---------------------------------------------------------------------------------------------------

proc initGfx(argc: ptr cint; argv: pointer) =
  ## Initialize OpenGL rendering.

  glutInit(argc, argv)
  glutInitDisplayMode(GLUT_RGB)
  glutInitWindowSize(640, 320)
  glutIdleFunc(render)
  discard glutCreateWindow("Pendulum")
  glutDisplayFunc(render)
  loadExtensions()

#———————————————————————————————————————————————————————————————————————————————————————————————————

initSimulation(length = 5, gravitation = 9.81, start = PI / 3)

var argc: cint = 0
initGfx(addr(argc), nil)
glutMainLoop()


# Pendulum simulation.

import math
import times

import gintro/[gobject, gdk, gtk, gio, cairo]
import gintro/glib except Pi

type

  # Description of the simulation.
  Simulation = ref object
    area: DrawingArea       # Drawing area.
    length: float           # Pendulum length.
    g: float                # Gravity (should be positive).
    time: Time              # Current time.
    theta0: float           # initial angle.
    theta: float            # Current angle.
    omega: float            # Angular velocity = derivative of theta.
    accel: float            # Angular acceleration = derivative of omega.
    e: float                # Total energy.

#---------------------------------------------------------------------------------------------------

proc newSimulation(area: DrawingArea; length, g, theta0: float): Simulation {.noInit.} =
  ## Allocate and initialize the simulation object.

  new(result)
  result.area = area
  result.length = length
  result.g = g
  result.time = getTime()
  result.theta0 = theta0
  result.theta = theta0
  result.omega = 0
  result.accel = -g / length * sin(theta0)
  result.e = g * length * (1 - cos(theta0))    # Total energy = potential energy when starting.

#---------------------------------------------------------------------------------------------------

template toFloat(dt: Duration): float = dt.inNanoseconds.float / 1e9

#---------------------------------------------------------------------------------------------------

const Origin = (x: 320.0, y: 100.0)   # Pivot coordinates.
const Scale = 300                     # Coordinates scaling constant.

proc draw(sim: Simulation; context: cairo.Context) =
  ## Draw the pendulum.

  # Compute coordinates in drawing area.
  let x = Origin.x + sin(sim.theta) * Scale
  let y = Origin.y + cos(sim.theta) * Scale

  # Clear the region.
  context.moveTo(0, 0)
  context.setSource(0.0, 0.0, 0.0)
  context.paint()

  # Draw pendulum.
  context.moveTo(Origin.x, Origin.y)
  context.setSource(0.3, 1.0, 0.3)
  context.lineTo(x, y)
  context.stroke()

  # Draw pivot.
  context.setSource(0.3, 0.3, 1.0)
  context.arc(Origin.x, Origin.y, 8, 0, 2 * Pi)
  context.fill()

  # Draw mass.
  context.setSource(1.0, 0.3, 0.3)
  context.arc(x, y, 8, 0, 2 * Pi)
  context.fill()

#---------------------------------------------------------------------------------------------------

proc update(sim: Simulation): gboolean =
  ## Update the simulation state.

  # compute time interval.
  let nextTime = getTime()
  let dt = (nextTime - sim.time).toFloat
  sim.time = nextTime

  # Update theta and omega.
  sim.theta += (sim.omega + dt * sim.accel / 2) * dt
  sim.omega += sim.accel * dt

  # If, due to computation errors, potential energy is greater than total energy,
  # reset theta to ±theta0 and omega to 0.
  if sim.length * sim.g * (1 - cos(sim.theta)) >= sim.e:
    sim.theta = sgn(sim.theta).toFloat * sim.theta0
    sim.omega = 0

  # Compute acceleration.
  sim.accel = -sim.g / sim.length * sin(sim.theta)

  result = gboolean(1)

  sim.draw(sim.area.window.cairoCreate())

#---------------------------------------------------------------------------------------------------

proc activate(app: Application) =
  ## Activate the application.

  let window = app.newApplicationWindow()
  window.setSizeRequest(640, 480)
  window.setTitle("Pendulum simulation")

  let area = newDrawingArea()
  window.add(area)

  let sim = newSimulation(area, length = 5, g = 9.81, theta0 = PI / 3)

  timeoutAdd(10, update, sim)

  window.showAll()

#———————————————————————————————————————————————————————————————————————————————————————————————————

let app = newApplication(Application, "Rosetta.pendulum")
discard app.connect("activate", activate)
discard app.run()