// title: Damped Driven Double Pendulum Sonification and Visualization
// author: xavieredavenport
// description:
// 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.
// code:
(

///////////////////////
// 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;
};

)