///// double pendulum model (second order ODE with Runge Kutta) ////////// ///// Emily Graber ///// 220a final project ///////////////////////////equations of motion for pendulum//////////////////////////////// 19.81 => float g; // gravity 1.0 => float L; // pendulum arm in meters fun float[] f( float X[], float t ){ X[0] => float theta1; X[1] => float omega1 => float f_theta1; X[2] => float theta2; X[3] => float omega2 => float f_theta2; -(Math.pow(omega1, 2)*Math.sin(2*theta1-2*theta2)+2*Math.pow(omega2, 2)*Math.sin(theta1-theta2)+(g/L)*Math.sin(theta1-2*theta2)+(g/L)*3*Math.sin(theta1))/(3-Math.cos(2*theta1-2*theta2)) => float f_omega1; (4*Math.pow(omega1, 2)*Math.sin(theta1-theta2)+Math.pow(omega2, 2)*Math.sin(2*theta1-2*theta2)+2*(g/L)*Math.sin(2*theta1-theta2)-2*(g/L)*Math.sin(theta2))/(3-Math.cos(2*theta1-2*theta2)) => float f_omega2; [f_theta1, f_omega1, f_theta2, f_omega2] @=> float X_temp[]; return X_temp; } ///////////// set initial conditions of pendulum /////////////////////////////// 0.0 => float omega1; // initial speed of pendulum 0.0 => float omega2; 5*pi/6 => float theta1; // initial position of pendulum pi/2 => float theta2; 1.0 => float m; // bob mass (kg) [ theta1, omega1, theta2, omega2 ] @=> float X[]; //array of angular position and speed //////////////////////////////////////////////////////////////////////////////////////// float X_temp[4]; float bob1[2]; float bob2[2]; SndBuf s1; SndBuf s2; "/user/e/emgraber/Documents/220a_final/perlman_partita.wav" => s1.read; "/user/e/emgraber/Documents/220a_final/perlman_partita.wav" => s2.read; DBAP4e b1; DBAP4e b2; //0.2 => s1.gain; //3.0 => s2.gain; float duration; 0.0 => float t; // starting time 100.0 => float b; // ending time 10000 => float n; // number of steps (b-t)/n => float h; // step size while( true ){ h +=> t; Math.fabs( X[3] ) => duration; duration::day => now ; // this is not working?????????????????????????????????? s1 => b1; s2 => b2; /////////// set position of masses ////////////////////////////////////////////// L*Math.cos(X[0] - pi/2) => bob1[0]; // x position of bob 1 L*Math.sin(X[0] - pi/2) => bob1[1]; // y position of bob 1 L*Math.cos(X[0] - pi/2)+L*Math.cos(X[2] - pi/2) => bob2[0]; // x position of bob 2 L*Math.sin(X[0] - pi/2)+L*Math.sin(X[2] - pi/2) => bob2[1]; // y position of bob 2 b1.setPosition( bob1[0], bob1[1] ); b2.setPosition( bob2[0], bob2[1] ); b1.setReverb(0); b2.setReverb(0); //////////// OSCsend //////////////////////////////////////////////////////////////// //////////// update anglular position and speed array, X /////////////////////////// f(X,t) @=> X_temp; float k1[4]; float k1over2[4]; float Xplusk1over2[4]; for ( 0 => int i; i <= 3; i++ ){ h*X_temp[i] => k1[i]; //h*f(X[],t) => foat k1; k1[i]/2 => k1over2[i]; //k1/2 X[i] + k1over2[i] => Xplusk1over2[i]; //X + k1/2 } f(Xplusk1over2, t+h/2) @=> X_temp; float k2[4]; float k2over2[4]; float Xplusk2over2[4]; for ( 0 => int i; i <= 3; i++ ){ h*X_temp[i] => k2[i]; //h*f(X+k1/2, t+h/2) => k2; k2[i]/2 => k2over2[i]; //k2/2 X[i] + k2over2[i] => Xplusk2over2[i]; } f(Xplusk2over2, t+h/2) @=> X_temp; float k3[4]; float Xplusk3[4]; for ( 0 => int i; i <= 3; i++ ){ h*X_temp[i] => k3[i]; //h*f(X+k2/2, t+h/2) => k3; X[i] + k3[i] => Xplusk3[i]; } f(Xplusk3, t+h) @=> X_temp; float k4[4]; for ( 0 => int i; i <= 3; i++ ){ h*X_temp[i] => k4[i]; //h*f(X+k3, t+h) => k4 } for ( 0 => int i; i <= 3; i++ ){ X[i] + k1[i]/6.0 + 2*k2[i]/6.0 + 2*k3[i]/6.0 + k4[i]/6.0 => X[i]; //updates X } }