«Damped Driven Double Pendulum Sonification and Visualization» by xavieredavenport

on 29 Jun'21 23:25 in sonificationphysicspendulum

ctrl+enter to run ctrl+. to stop

Version 1. A little toy physics simulation of a double pendulum with damping and drive. The math is a little inconsistent in SuperCollider when dealing with numbers super small and super large, especially when trig functions are involved. I've tried to mitigate the errors by using sinPi instead of regular sin, but you still have to be careful when setting system variables. A strong damping is a must, and you have to be very careful when adding a drive to the system. Otherwise, we see the expected chaotic behavior for a range of easier-to-calculate parameters.

Pitch is determined by the distance of the second pendulum bob from the center of the screen, while amplitude is determined by the velocity of the bob.

log fork json 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
(

///////////////////////
// Variable declaration
///////////////////////

var centerx = 410, centery = 410, //Coordinates used to offset the pendulum fixed point
r1 = 200, //Length of pendulum rod 1
r2 = 200, //Length of pendulum rod 2
m1 = 20, //Mass of pendulum bob 1
m2 = 20, //Mass of pendulum bob 2
a1 = pi+(rrand(0.0001,0.0009)), //Initial angle of pendulum 1
a2 = pi, //Initial angle of pendulum 2
a1_v = 0, //Initial velocity of pendulum 1
a2_v = 0, //Initial velocity of pendulum 2
g = 1.0, //Gravitational constant

time_step = 1.0,
time = 0,
k1 = 10.0, //Damping coefficient of pendulum 1
k2 = 10.0, //Damping coefficient of pendulum 2, usually same as pendulum 1
force1 = 0.0, //Magnitude of oscillation strength
force2 = 0.0, //This is usually kept at 0
//freq1 = sqrt(g/r1)*sqrt(2+sqrt(2))/100,
freq1 = 0.0, //Drive frequency for pendulum 1
freq2 = 0.0, //Drive frequency for pendulum 2
iterations = 0, //Total frames drawn
alpha, beta, deltheta, num1, num2, num3, num4, gamma, den, //mathematical functions
a1_a, a2_a,  // accelerations
x1, y1, x2, y2, // positions

pitch = 440,

///////////////////////
// Window setup
///////////////////////


width = 820, height = 820,
w = Window("Double Pendulum", Rect(99, 99, width, height), false),
u = UserView(w, Rect(0, 0, width, height));

u.clearOnRefresh = true;
u.background = Color.white;
w.front;
u.frameRate = 40;
u.animate = true;
CmdPeriod.doOnce({if(w.isClosed.not, {w.close})});




///////////////////////
// Synthdef(s)
///////////////////////


SynthDef(\gliss, {
	arg freq=440, gate=1, amp=0.0, out=0;
	var sig, env;
	freq = Lag.kr(freq, 1);
	amp = Lag.kr(amp, 1);
	env = EnvGen.kr(Env.asr, gate, doneAction:2);
	sig = SinOsc.ar(freq)!2;
	sig = sig * amp;
	sig = sig * env;
	Out.ar(out, sig);
}).add;


x = Synth(\gliss);

u.drawFunc = {




	///////////////////////
	// Equations
	///////////////////////

	alpha = (k1 * a1_v) - (force1 * cosPi(freq1 * time));
	beta = (k2 * a2_v) - (force2 * cosPi(freq2 * time));
	deltheta = a1-a2;

	num1 = m2 * r1 * a1_v * a1_v * sinPi(2 * deltheta);
	num2 = 2 * m2 * r2 * a2_v * a2_v * sinPi(deltheta);
	num3 = 2 * g * m2 * cosPi(a2)*sinPi(deltheta);
	num4 = 2 * g * m1 * sinPi(a1);
	gamma = (2*alpha) - (2*beta*cosPi(deltheta));
	den = (-2 * r1*m1) + (r1*m2*sinPi(deltheta)*sinPi(deltheta));
	a1_a = (num1 + num2 + num3 + num4 + gamma)/den;

	num1 = m2 * r2 * a2_v * a2_v * sinPi(2 * deltheta);
	num2 = 2 * (m1+m2) * r1 * a1_v*a1_v * sinPi(deltheta);
	num3 = 2 * g * (m1+m2) * cosPi(a1)*sinPi(deltheta);
	gamma = (2*alpha*cosPi(deltheta)) - (2*(m1+m2)*beta/m2);
	den = (2 * r2*m1) + (r2*m2*sinPi(deltheta)*sinPi(deltheta));
	a2_a = (num1 + num2 + num3 + gamma)/den;

	x1 = r1 * sinPi(a1);
	y1 = r1 * cosPi(a1);

	x2 = x1 + (r2 * sinPi(a2));
	y2 = y1 + (r2 * cosPi(a2));



	///////////////////////
	// Advance position
	///////////////////////

	a1_v = a1_v + a1_a;
	a2_v = a2_v + a2_a;
	a1 = a1 + a1_v;
	a2 = a2 + a2_v;


	///////////////////////
	// Create Window and draw pendulum
	///////////////////////


	Pen.strokeColor = Color.blue;
	Pen.width = 2;
	Pen.moveTo(centerx@centery);
	Pen.lineTo((centerx+x1)@(centery+y1));
	Pen.lineTo((centerx+x2)@(centery+y2));
	Pen.stroke;
	Pen.strokeOval(Rect(centerx+x1-5, centery+y1-5, 10, 10));
	Pen.strokeOval(Rect(centerx+x2-5, centery+y2-5, 10, 10));

	pitch = sqrt((x2*x2)+(y2*y2));

	x.set(\freq, pitch);
	x.set(\amp, a2_v*5);

	time = time + time_step;

	iterations = iterations + 1;
};

)
raw 3599 chars (focus & ctrl+a+c to copy)
reception
comments