Animate a pendulum

From Rosetta Code

Jump to: navigation, search
Capture of the Oz version.
Animate a pendulum is a programming task. Visitors like you are encouraged to solve it according to the task description, using any language they may happen to know.
Add to BlogMarksAdd to del.icio.usAdd to diggAdd to NewsvineAdd to redditAdd to Slashdot

One good way of making an animation is by simulating a physical system and illustrating the variables in that system using a dynamically changing graphical display. The classic such physical system is a simple gravity pendulum.

For this task, create a simple physical model of a pendulum and animate it.


Contents

[edit] C

Library: SDL

Library: Cairo

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <cairo.h>
#include <SDL.h>
 
#define WIDTH 320.0
#define HEIGHT 240.0
 
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
 
const double L = HEIGHT*150.0/200.0; // length
double theta = M_PI/4.0;
double dtheta = 0.0;
int homeX = WIDTH / 2;
int homeY = HEIGHT*25.0/200.0;
 
static uint32_t update_anim(uint32_t iv, void *p)
{
double scaling = 30.0 / (L*L); // this seems better for this impl.
 
double firstDDtheta = -sin(theta) * scaling;
double midDtheta = dtheta + firstDDtheta;
double midtheta = theta + (dtheta + midDtheta)/2.0;
 
double midDDtheta = -sin(midtheta) * scaling;
midDtheta = dtheta + (firstDDtheta + midDDtheta)/2.0;
midtheta = theta + (dtheta + midDtheta)/2.0;
 
midDDtheta = -sin(midtheta) * scaling;
double lastDtheta = midDtheta + midDDtheta;
double lasttheta = midtheta + (midDtheta + lastDtheta)/2.0;
 
double lastDDtheta = -sin(lasttheta) * scaling;
lastDtheta = midDtheta + (midDDtheta + lastDDtheta)/2.0;
lasttheta = midtheta + (midDtheta + lastDtheta)/2.0;
 
dtheta = lastDtheta;
theta = lasttheta;
 
return iv;
}
 
void show_pendolum(SDL_Surface *s)
{
cairo_surface_t *surf =
cairo_image_surface_create_for_data( s->pixels,
CAIRO_FORMAT_RGB24,
s->w, s->h, s->pitch);
 
cairo_t *ct = cairo_create(surf);
 
double x = homeX + L*sin(theta);
double y = homeY + L*cos(theta);
 
cairo_set_source_rgba(ct, 0, 0, 0, 1);
cairo_new_path(ct);
cairo_rectangle(ct, 0, 0, WIDTH, HEIGHT);
cairo_fill(ct);
 
cairo_set_source_rgba(ct, 1, 1, 1, 1);
 
// rod
cairo_new_path(ct);
cairo_move_to(ct, homeX, homeY);
cairo_line_to(ct, x, y);
cairo_stroke(ct);
 
cairo_set_source_rgba(ct, 1, 1, 0, 1);
 
// bob
cairo_new_path(ct);
cairo_move_to(ct, x, y);
cairo_arc(ct, x, y, 20.0, 0, 2*M_PI);
cairo_fill(ct);
 
cairo_set_source_rgba(ct, 0.7, 0.7, 0.7, 1);
 
// plate
cairo_new_path(ct);
cairo_move_to(ct, 0, homeY);
cairo_line_to(ct, WIDTH, homeY);
cairo_stroke(ct);
 
// pivot
cairo_new_path(ct);
cairo_move_to(ct, homeX, homeY);
cairo_arc(ct, homeX, homeY, 5.0, 0, 2*M_PI);
cairo_fill(ct);
 
cairo_surface_destroy(surf);
cairo_destroy(ct);
}
 
int main()
{
SDL_Surface *scr, *t;
SDL_Event event[1];
bool quit = false;
 
if ( SDL_Init(SDL_INIT_VIDEO | SDL_INIT_TIMER) >= 0 ) {
atexit(SDL_Quit);
if ( SDL_SetVideoMode(WIDTH, HEIGHT, 32, SDL_DOUBLEBUF) != NULL ) {
 
scr = SDL_GetVideoSurface();
SDL_AddTimer(30, update_anim, NULL);
 
event->type = SDL_VIDEOEXPOSE;
SDL_PushEvent(event);
 
while(SDL_WaitEvent(event) && !quit) {
switch(event->type) {
case SDL_VIDEOEXPOSE:
while(SDL_LockSurface(scr) != 0) SDL_Delay(1);
 
show_pendolum(scr);
 
SDL_UnlockSurface(scr);
SDL_Flip(scr);
event->type = SDL_VIDEOEXPOSE;
SDL_PushEvent(event);
break;
case SDL_KEYDOWN:
if (event->key.keysym.sym == SDLK_q) quit = true;
break;
}
}
 
SDL_FreeSurface(scr);
}
}
 
return 0;
}

[edit] E

Works with: E-on-Java (Uses Java Swing for GUI. The animation logic is independent, however.)

The angle θ of a pendulum with length L and acceleration due to gravity g with all its mass at the end and no friction/air resistance has an acceleration at any given moment of

\frac{d^2θ}{dt^2}\theta = -\frac{g}{L} \sin \theta

This simulation uses this formula directly, updating the velocity from the acceleration and the position from the velocity; inaccuracy results from the finite timestep.

The event flow works like this: The clock object created by the simulation steps the simulation on the specified in the interval. The simulation writes its output to angle, which is a Lamport slot which can notify of updates. The whenever set up by makeDisplayComponent listens for updates and triggers redrawing as long as interest has been expressed, which is done whenever the component actually redraws, which happens only if the component's window is still on screen. When the window is closed, additionally, the simulation itself is stopped and the application allowed to exit. (This logic is more general than necessary; it is designed to be suitable for a larger application as well.)

#!/usr/bin/env rune
pragma.syntax("0.9")
 
def pi := (-1.0).acos()
def makeEPainter := <unsafe:com.zooko.tray.makeEPainter>
def makeLamportSlot := <import:org.erights.e.elib.slot.makeLamportSlot>
def whenever := <import:org.erights.e.elib.slot.whenever>
def colors := <import:java.awt.makeColor>
 
# --------------------------------------------------------------
# --- Definitions
 
def makePendulumSim(length_m :float64,
gravity_mps2 :float64,
initialAngle_rad :float64,
timestep_ms :int) {
var velocity := 0
def &angle := makeLamportSlot(initialAngle_rad)
def k := -gravity_mps2/length_m
def timestep_s := timestep_ms / 1000
def clock := timer.every(timestep_ms, fn _ {
def acceleration := k * angle.sin()
velocity += acceleration * timestep_s
angle += velocity * timestep_s
})
return [clock, &angle]
}
 
def makeDisplayComponent(&angle) {
def c
def updater := whenever([&angle], fn { c.repaint() })
 
bind c := makeEPainter(def paintCallback {
to paintComponent(g) {
try {
def originX := c.getWidth() // 2
def originY := c.getHeight() // 2
def pendRadius := (originX.min(originY) * 0.95).round()
def ballRadius := (originX.min(originY) * 0.04).round()
def ballX := (originX + angle.sin() * pendRadius).round()
def ballY := (originY + angle.cos() * pendRadius).round()
 
g.setColor(colors.getWhite())
g.fillRect(0, 0, c.getWidth(), c.getHeight())
g.setColor(colors.getBlack())
 
g.fillOval(originX - 2, originY - 2, 4, 4)
g.drawLine(originX, originY, ballX, ballY)
g.fillOval(ballX - ballRadius, ballY - ballRadius, ballRadius * 2, ballRadius * 2)
 
updater[] # provoke interest provided that we did get drawn (window not closed)
} catch p {
stderr.println(`In paint callback: $p${p.eStack()}`)
}
}
})
 
c.setPreferredSize(<awt:makeDimension>(300, 300))
return c
}
 
# --------------------------------------------------------------
# --- Application setup
 
def [clock, &angle] := makePendulumSim(1, 9.80665, pi*99/100, 10)
 
# Initialize AWT, move to AWT event thread
when (currentVat.morphInto("awt")) -> {
 
# Create the window
def frame := <unsafe:javax.swing.makeJFrame>("Pendulum")
frame.setContentPane(def display := makeDisplayComponent(&angle))
frame.addWindowListener(def mainWindowListener {
to windowClosing(_) {
clock.stop()
interp.continueAtTop()
}
match _ {}
})
frame.setLocation(50, 50)
frame.pack()
 
# Start and become visible
frame.show()
clock.start()
}
 
interp.blockAtTop()

[edit] Factor

Approximation of the pendulum for small swings : theta = theta0 * cos(omega0 * t)

USING: accessors alarms arrays calendar colors.constants kernel
locals math math.constants math.functions math.rectangles
math.vectors opengl sequences system ui ui.gadgets ui.render ;
IN: pendulum
 
CONSTANT: g 9.81
CONSTANT: l 20
CONSTANT: theta0 0.5
 
: current-time ( -- time ) nano-count -9 10^ * ;
 
: T0 ( -- T0 ) 2 pi l g / sqrt * * ;
: omega0 ( -- omega0 ) 2 pi * T0 / ;
: theta ( -- theta ) current-time omega0 * cos theta0 * ;
 
: relative-xy ( theta l -- xy )
[ [ sin ] [ cos ] bi ]
[ [ * ] curry ] bi* bi@ 2array ;
: theta-to-xy ( origin theta l -- xy ) relative-xy v+ ;
 
TUPLE: pendulum-gadget < gadget alarm ;
 
: O ( gadget -- origin ) rect-bounds [ drop ] [ first 2 / ] bi* 0 2array ;
: window-l ( gadget -- l ) rect-bounds [ drop ] [ second ] bi* ;
: gadget-xy ( gadget -- xy ) [ O ] [ drop theta ] [ window-l ] tri theta-to-xy ;
 
M: pendulum-gadget draw-gadget*
COLOR: black gl-color
[ O ] [ gadget-xy ] bi gl-line ;
 
M:: pendulum-gadget graft* ( gadget -- )
[ gadget relayout-1 ]
20 milliseconds every gadget (>>alarm) ;
M: pendulum-gadget ungraft* alarm>> cancel-alarm ;
 
: <pendulum-gadget> ( -- gadget )
pendulum-gadget new
{ 500 500 } >>pref-dim ;
: pendulum-main ( -- )
[ <pendulum-gadget> "pendulum" open-window ] with-ui ;
MAIN: pendulum-main
 

[edit] J

require 'gl2 trig'
coinsert 'jgl2'
 
DT =: %30 NB. seconds
ANGLE=: 0.25p1 NB. radians
L =: 1 NB. metres
G =: 9.80665 NB. ms_2
VEL =: 0 NB. ms_1
 
PEND=: noun define
pc pend;pn "Pendulum";
xywh 0 0 320 200;cc isi isigraph rightmove bottommove;
pas 0 0;pcenter;
rem form end;
)
 
pend_run =: verb def ' wd PEND,'';pshow;timer '',":DT * 1000 '
pend_close =: verb def ' wd ''timer 0; pclose'' '
pend_isi_paint=: verb def ' drawPendulum ANGLE '
 
sys_timer_z_=: verb define
recalcAngle ''
wd 'psel pend; setinvalid isi'
)
 
recalcAngle=: verb define
accel=. - (G % L) * sin ANGLE
VEL =: VEL + accel * DT
ANGLE=: ANGLE + VEL * DT
)
 
drawPendulum=: verb define
width=. {. glqwh''
ps=. (-: width) , 40
pe=. ps + 280 <.@* (cos , sin) 0.5p1 + y NB. adjust orientation
glrgb 91 91 91
glbrush''
gllines ps , pe
glellipse (,~ ps - -:) 40 15
glellipse (,~ pe - -:) 20 20
glrect 0 0 ,width, 40
)
 
pend_run'' NB. run animation

[edit] Logo

Works with: UCB Logo

make "angle 45
make "L 1
make "bob 10
 
to draw.pendulum
clearscreen
seth :angle+180 ; down on screen is 180
forward :L*100-:bob
penup
forward :bob
pendown
arc 360 :bob
end
 
make "G 9.80665
make "dt 1/30
make "acc 0
make "vel 0
 
to step.pendulum
make "acc -:G / :L * sin :angle
make "vel  :vel + :acc * :dt
make "angle :angle + :vel * :dt
wait :dt*60
draw.pendulum
end
 
hideturtle
until [key?] [step.pendulum]

[edit] Oz

Inspired by the E and Ruby versions.

declare
[QTk] = {Link ['x-oz://system/wp/QTk.ozf']}
 
Pi = 3.14159265
 
class PendulumModel
feat
K
attr
angle
velocity
 
meth init(length:L <= 1.0 %% meters
gravity:G <= 9.81 %% m/s²
initialAngle:A <= Pi/2.) %% radians
self.K = ~G / L
angle := A
velocity := 0.0
end
 
meth nextAngle(deltaT:DeltaTMS %% milliseconds
 ?Angle) %% radians
DeltaT = {Int.toFloat DeltaTMS} / 1000.0 %% seconds
Acceleration = self.K * {Sin @angle}
in
velocity := @velocity + Acceleration * DeltaT
angle := @angle + @velocity * DeltaT
Angle = @angle
end
end
 
%% Animates a pendulum on a given canvas.
class PendulumAnimation from Time.repeat
feat
Pend
Rod
Bob
home:pos(x:160 y:50)
length:140.0
delay
 
meth init(Pendulum Canvas delay:Delay <= 25) %% milliseconds
self.Pend = Pendulum
self.delay = Delay
%% plate and pivot
{Canvas create(line 0 self.home.y 320 self.home.y width:2 fill:grey50)}
{Canvas create(oval 155 self.home.y-5 165 self.home.y+5 fill:grey50 outline:black)}
%% the pendulum itself
self.Rod = {Canvas create(line 1 1 1 1 width:3 fill:black handle:$)}
self.Bob = {Canvas create(oval 1 1 2 2 fill:yellow outline:black handle:$)}
%%
{self setRepAll(action:Animate delay:Delay)}
end
 
meth Animate
Theta = {self.Pend nextAngle(deltaT:self.delay $)}
%% calculate x and y from angle
X = self.home.x + {Float.toInt self.length * {Sin Theta}}
Y = self.home.y + {Float.toInt self.length * {Cos Theta}}
in
%% update canvas
try
{self.Rod setCoords(self.home.x self.home.y X Y)}
{self.Bob setCoords(X-15 Y-15 X+15 Y+15)}
catch system(tk(alreadyClosed ...) ...) then skip end
end
end
 
Pendulum = {New PendulumModel init}
 
Canvas
GUI = td(title:"Pendulum"
canvas(width:320 height:210 handle:?Canvas)
action:proc {$} {Animation stop} {Window close} end
)
Window = {QTk.build GUI}
 
Animation = {New PendulumAnimation init(Pendulum Canvas)}
in
{Window show}
{Animation go}

[edit] Ruby

Library: Ruby/Tk

Translation of: Tcl

This does not have the window resizing handling that Tcl does -- I did not spend enough time in the docs to figure out how to get the new window size out of the configuration event. Of interest when running this pendulum side-by-side with the Tcl one: the Tcl pendulum swings noticibly faster.

require 'tk'
 
$root = TkRoot.new("title" => "Pendulum Animation")
$canvas = TkCanvas.new($root) do
width 320
height 200
create TkcLine, 0,25,320,25, 'tags' => 'plate', 'width' => 2, 'fill' => 'grey50'
create TkcOval, 155,20,165,30, 'tags' => 'pivot', 'outline' => "", 'fill' => 'grey50'
create TkcLine, 1,1,1,1, 'tags' => 'rod', 'width' => 3, 'fill' => 'black'
create TkcOval, 1,1,2,2, 'tags' => 'bob', 'outline' => 'black', 'fill' => 'yellow'
end
$canvas.raise('pivot')
$canvas.pack('fill' => 'both', 'expand' => true)
 
$Theta = 45.0
$dTheta = 0.0
$length = 150
$homeX = 160
$homeY = 25
 
def show_pendulum
angle = $Theta * Math::PI / 180
x = $homeX + $length * Math.sin(angle)
y = $homeY + $length * Math.cos(angle)
$canvas.coords('rod', $homeX, $homeY, x, y)
$canvas.coords('bob', x-15, y-15, x+15, y+15)
end
 
def recompute_angle
scaling = 3000.0 / ($length ** 2)
# first estimate
firstDDTheta = -Math.sin($Theta * Math::PI / 180) * scaling
midDTheta = $dTheta + firstDDTheta
midTheta = $Theta + ($dTheta + midDTheta)/2
# second estimate
midDDTheta = -Math.sin(midTheta * Math::PI / 180) * scaling
midDTheta = $dTheta + (firstDDTheta + midDDTheta)/2
midTheta = $Theta + ($dTheta + midDTheta)/2
# again, first
midDDTheta = -Math.sin(midTheta * Math::PI / 180) * scaling
lastDTheta = midDTheta + midDDTheta
lastTheta = midTheta + (midDTheta + lastDTheta)/2
# again, second
lastDDTheta = -Math.sin(lastTheta * Math::PI/180) * scaling
lastDTheta = midDTheta + (midDDTheta + lastDDTheta)/2
lastTheta = midTheta + (midDTheta + lastDTheta)/2
# Now put the values back in our globals
$dTheta = lastDTheta
$Theta = lastTheta
end
 
def animate
recompute_angle
show_pendulum
$after_id = $root.after(15) {animate}
end
 
show_pendulum
$after_id = $root.after(500) {animate}
 
$canvas.bind('<Destroy>') {$root.after_cancel($after_id)}
 
Tk.mainloop

[edit] Scheme

Library: Scheme/PsTk

Translation of: Ruby

This is a direct translation of the Ruby/Tk example into Scheme + PS/Tk.

#!r6rs
 
;;; R6RS implementation of Pendulum Animation
 
(import (rnrs)
(lib pstk main) ; change this for your pstk installation
)
 
(define PI 3.14159)
(define *conv-radians* (/ PI 180))
(define *theta* 45.0)
(define *d-theta* 0.0)
(define *length* 150)
(define *home-x* 160)
(define *home-y* 25)
 
;;; estimates new angle of pendulum
(define (recompute-angle)
(define (avg a b) (/ (+ a b) 2))
(let* ((scaling (/ 3000.0 (* *length* *length*)))
; first estimate
(first-dd-theta (- (* (sin (* *theta* *conv-radians*)) scaling)))
(mid-d-theta (+ *d-theta* first-dd-theta))
(mid-theta (+ *theta* (avg *d-theta* mid-d-theta)))
; second estimate
(mid-dd-theta (- (* (sin (* mid-theta *conv-radians*)) scaling)))
(mid-d-theta-2 (+ *d-theta* (avg first-dd-theta mid-dd-theta)))
(mid-theta-2 (+ *theta* (avg *d-theta* mid-d-theta-2)))
; again first
(mid-dd-theta-2 (- (* (sin (* mid-theta-2 *conv-radians*)) scaling)))
(last-d-theta (+ mid-d-theta-2 mid-dd-theta-2))
(last-theta (+ mid-theta-2 (avg mid-d-theta-2 last-d-theta)))
; again second
(last-dd-theta (- (* (sin (* last-theta *conv-radians*)) scaling)))
(last-d-theta-2 (+ mid-d-theta-2 (avg mid-dd-theta-2 last-dd-theta)))
(last-theta-2 (+ mid-theta-2 (avg mid-d-theta-2 last-d-theta-2))))
; put values back in globals
(set! *d-theta* last-d-theta-2)
(set! *theta* last-theta-2)))
 
;;; The main event loop and graphics context
(let ((tk (tk-start)))
(tk/wm 'title tk "Pendulum Animation")
(let ((canvas (tk 'create-widget 'canvas)))
 
;;; redraw the pendulum on canvas
;;; - uses angle and length to compute new (x,y) position of bob
(define (show-pendulum canvas)
(let* ((pendulum-angle (* *conv-radians* *theta*))
(x (+ *home-x* (* *length* (sin pendulum-angle))))
(y (+ *home-y* (* *length* (cos pendulum-angle)))))
(canvas 'coords 'rod *home-x* *home-y* x y)
(canvas 'coords 'bob (- x 15) (- y 15) (+ x 15) (+ y 15))))
 
;;; move the pendulum and repeat after 20ms
(define (animate)
(recompute-angle)
(show-pendulum canvas)
(tk/after 20 animate))
 
;; layout the canvas
(tk/grid canvas 'column: 0 'row: 0)
(canvas 'create 'line 0 25 320 25 'tags: 'plate 'width: 2 'fill: 'grey50)
(canvas 'create 'oval 155 20 165 30 'tags: 'pivot 'outline: "" 'fill: 'grey50)
(canvas 'create 'line 1 1 1 1 'tags: 'rod 'width: 3 'fill: 'black)
(canvas 'create 'oval 1 1 2 2 'tags: 'bob 'outline: 'black 'fill: 'yellow)
 
;; get everything started
(show-pendulum canvas)
(tk/after 500 animate)
(tk-event-loop tk)))
 

[edit] Tcl

Works with: Tcl version 8.5 and Library: Tk

package require Tcl 8.5
package require Tk
 
# Make the graphical entities
pack [canvas .c -width 320 -height 200] -fill both -expand 1
.c create line 0 25 320 25 -width 2 -fill grey50 -tags plate
.c create line 1 1 1 1 -tags rod -width 3 -fill black
.c create oval 1 1 2 2 -tags bob -fill yellow -outline black
.c create oval 155 20 165 30 -fill grey50 -outline {} -tags pivot
 
# Set some vars
set points {}
set Theta 45.0
set dTheta 0.0
set pi 3.1415926535897933
set length 150
set homeX 160
 
# How to respond to a changing in size of the window
proc resized {width} {
global homeX
.c coords plate 0 25 $width 25
set homeX [expr {$width / 2}]
.c coords pivot [expr {$homeX-5}] 20 [expr {$homeX+5}] 30
showPendulum
}
 
# How to actually arrange the pendulum, mapping the model to the display
proc showPendulum {} {
global Theta dTheta pi length homeX
set angle [expr {$Theta * $pi/180}]
set x [expr {$homeX + $length*sin($angle)}]
set y [expr {25 + $length*cos($angle)}]
.c coords rod $homeX 25 $x $y
.c coords bob [expr {$x-15}] [expr {$y-15}] [expr {$x+15}] [expr {$y+15}]
}
 
# The dynamic part of the display
proc recomputeAngle {} {
global Theta dTheta pi length
set scaling [expr {3000.0/$length**2}]
 
# first estimate
set firstDDTheta [expr {-sin($Theta * $pi/180)*$scaling}]
set midDTheta [expr {$dTheta + $firstDDTheta}]
set midTheta [expr {$Theta + ($dTheta + $midDTheta)/2}]
# second estimate
set midDDTheta [expr {-sin($midTheta * $pi/180)*$scaling}]
set midDTheta [expr {$dTheta + ($firstDDTheta + $midDDTheta)/2}]
set midTheta [expr {$Theta + ($dTheta + $midDTheta)/2}]
# Now we do a double-estimate approach for getting the final value
# first estimate
set midDDTheta [expr {-sin($midTheta * $pi/180)*$scaling}]
set lastDTheta [expr {$midDTheta + $midDDTheta}]
set lastTheta [expr {$midTheta + ($midDTheta + $lastDTheta)/2}]
# second estimate
set lastDDTheta [expr {-sin($lastTheta * $pi/180)*$scaling}]
set lastDTheta [expr {$midDTheta + ($midDDTheta + $lastDDTheta)/2}]
set lastTheta [expr {$midTheta + ($midDTheta + $lastDTheta)/2}]
# Now put the values back in our globals
set dTheta $lastDTheta
set Theta $lastTheta
}
 
# Run the animation by updating the physical model then the display
proc animate {} {
global animation
 
recomputeAngle
showPendulum
 
# Reschedule
set animation [after 15 animate]
}
set animation [after 500 animate]; # Extra initial delay is visually pleasing
 
# Callback to handle resizing of the canvas
bind .c <Configure> {resized %w}
# Callback to stop the animation cleanly when the GUI goes away
bind .c <Destroy> {after cancel $animation}
Personal tools
Google AdSense