/*Jimmy Corno *September 17, 2003 *Stat Mech 2 * *Runge-Kutta 5th order integrator for the magnetic spinner toy *("magnetron"). The motion is all in terms of the angles of the first magnet *on each rotor (theta1 and theta2). */ class Magnetron3 { static double theta1, theta2, omega1, omega2; //Angles and angular //velocities static double theta1new, theta2new; //Updated variables for use static double omega1new, omega2new; //by Poincare function. static double x11, x12, x13, x21, x22, x23; //Cartesian positions of static double y11, y12, y13, y21, y22, y23; //magnets static double theta1tmp, theta2tmp, omega1tmp, //Temporary variables for omega2tmp; //integrator use static double[] omegaPrime = new double[2]; static double R, D, m, I; //Constants //R=rotor arm length //D=rotor spacing as a //fraction of R //m=dipole moment //I=rotor moment of inertia static double C = 1.0e-7; //mu0/4Pi static double d11to21, d11to22, d11to23; //magnet distances static double d12to21, d12to22, d12to23; static double d13to21, d13to22, d13to23; static double[] F_11on21 = new double[2]; //Interaction forces static double[] F_11on22 = new double[2]; static double[] F_11on23 = new double[2]; static double[] F_12on21 = new double[2]; static double[] F_12on22 = new double[2]; static double[] F_12on23 = new double[2]; static double[] F_13on21 = new double[2]; static double[] F_13on22 = new double[2]; static double[] F_13on23 = new double[2]; static double[] B_11at21 = new double[2]; //Interaction forces static double[] B_11at22 = new double[2]; static double[] B_11at23 = new double[2]; static double[] B_12at21 = new double[2]; static double[] B_12at22 = new double[2]; static double[] B_12at23 = new double[2]; static double[] B_13at21 = new double[2]; static double[] B_13at22 = new double[2]; static double[] B_13at23 = new double[2]; static double h, hnew, hmax, hmin; //time step and limits static double delta, err; //error controls static double errTheta1, errTheta2; static double errOmega1, errOmega2; static double f0, f1, f2, f3, f4, f5; //derivatives static double g0, g1, g2, g3, g4, g5; static double j0, j1, j2, j3, j4, j5; static double k0, k1, k2, k3, k4, k5; /*These doubles will be used to avoid the use of all but 4 sines and *cosines with every step, since the magnets have a fixed angular *separation.*/ static double sin120 = Math.sin(Math.PI*2.0/3.0); static double cos120 = Math.cos(Math.PI*2.0/3.0); static double sin240 = -sin120; static double cos240 = cos120; static boolean retry = false; static boolean stepProblem = false; static double time, stopTime; static int n = 0; //step count static int stopN; static int i = 0; //step count excluding retries public static void main( String[] Args ) { delta = .0001; //Maximum allowable error for a step h = .001; //Initial step size hmax = .05; hmin = .000001; hnew = h; theta1 = 0.0; theta2 = Math.PI + 0.1; omega1 = Math.PI*2.0; omega2 = 0; R = 1.0; D=2.1; m=100.0; I=0.002; time = 0.0; stopTime = 20.0; stopN = 60000; System.out.println( "Ei = " + energy(theta1, theta2, omega1, omega2) ); //Integrates until escape or a problem. do { /* if ( !retry && (n%1000 == 0) ) { System.out.print( time + "," + theta1 % Math.PI ); System.out.print( "," + theta2 % Math.PI + "," ); System.out.println( omega1 + "," + omega2 ); } */ advance(); /* if( (omega2 < 0) && (omega2new > 0) ) { Poincare( theta1new, omega1new, omega2new, theta1, omega1, omega2 ); } **/ theta1 = theta1new; theta2 = theta2new; omega1 = omega1new; omega2 = omega2new; } while ( (time < stopTime) && !(stepProblem) ); if( stepProblem ) { System.out.println( "Integration error: step size below minimum" ); } System.out.println( "n = " + n ); System.out.println( "t = " + time ); System.out.println( "Ef = " + energy(theta1, theta2, omega1, omega2) ); } //Advance is the integrator method, uses RK4-5 adaptive step method public static void advance() { theta1tmp = theta1; theta2tmp = theta2; omega1tmp = omega1; omega2tmp = omega2; f0 = omega1tmp; g0 = omega2tmp; omegaPrime = alpha( theta1tmp, theta2tmp ); j0 = omegaPrime[0]; k0 = omegaPrime[1]; theta1tmp = theta1 + b10*h*f0; theta2tmp = theta2 + b10*h*g0; omega1tmp = omega1 + b10*h*j0; omega2tmp = omega2 + b10*h*k0; f1 = omega1tmp; g1 = omega2tmp; omegaPrime = alpha( theta1tmp, theta2tmp ); j1 = omegaPrime[0]; k1 = omegaPrime[1]; theta1tmp = theta1 + b20*h*f0 + b21*h*f1; theta2tmp = theta2 + b20*h*g0 + b21*h*g1; omega1tmp = omega1 + b20*h*j0 + b21*h*j1; omega2tmp = omega2 + b20*h*k0 + b21*h*k1; f2 = omega1tmp; g2 = omega2tmp; omegaPrime = alpha( theta1tmp, theta2tmp ); j2 = omegaPrime[0]; k2 = omegaPrime[1]; theta1tmp = theta1 + b30*h*f0 + b31*h*f1 + b32*h*f2; theta2tmp = theta2 + b30*h*g0 + b31*h*g1 + b32*h*g2; omega1tmp = omega1 + b30*h*j0 + b31*h*j1 + b32*h*j2; omega2tmp = omega2 + b30*h*k0 + b31*h*k1 + b32*h*k2; f3 = omega1tmp; g3 = omega2tmp; omegaPrime = alpha( theta1tmp, theta2tmp ); j3 = omegaPrime[0]; k3 = omegaPrime[1]; theta1tmp = theta1 + b40*h*f0 + b41*h*f1 + b42*h*f2 + b43*h*f3; theta2tmp = theta2 + b40*h*g0 + b41*h*g1 + b42*h*g2 + b43*h*g3; omega1tmp = omega1 + b40*h*j0 + b41*h*j1 + b42*h*j2 + b43*h*j3; omega2tmp = omega2 + b40*h*k0 + b41*h*k1 + b42*h*k2 + b43*h*k3; f4 = omega1tmp; g4 = omega2tmp; omegaPrime = alpha( theta1tmp, theta2tmp ); j4 = omegaPrime[0]; k4 = omegaPrime[1]; theta1tmp = theta1 + b50*h*f0 + b51*h*f1 + b52*h*f2 + b53*h*f3 + b54*h*f4; theta2tmp = theta2 + b50*h*g0 + b51*h*g1 + b52*h*g2 + b53*h*g3 + b54*h*g4; omega1tmp = omega1 + b50*h*j0 + b51*h*j1 + b52*h*j2 + b53*h*j3 + b54*h*j4; omega1tmp = omega1 + b50*h*k0 + b51*h*k1 + b52*h*k2 + b53*h*k3 + b54*h*k4; f5 = omega1tmp; g5 = omega2tmp; omegaPrime = alpha( theta1tmp, theta2tmp ); j5 = omegaPrime[0]; k5 = omegaPrime[1]; theta1new = theta1 + h * ( c0*f0 + c2*f2 + c3*f3 + c4*f4 + c5*f5 ); theta2new = theta2 + h * ( c0*g0 + c2*g2 + c3*g3 + c4*g4 + c5*g5 ); omega1new = omega1 + h * ( c0*j0 + c2*j2 + c3*j3 + c4*j4 + c5*j5 ); omega2new = omega2 + h * ( c0*k0 + c2*k2 + c3*k3 + c4*k4 + c5*k5 ); time += h; i++; n++; } //This doubles calculates interaction forces. public static double[] F_AonB (double xA, double yA, double xB, double yB) { double[] F = new double[2]; double dx = D + xB - xA; double dy = yB - yA; double dist = Math.sqrt(dx*dx + dy*dy); F[0] = 3.0*C*m*m*( (2.0*xA*xA*dx + xA*(dx*xA + dy*yA) + yA*yA*dx + xA*yA*dy)/Math.pow(dist, 5) - 5*(dx*xA + dy*yA) *(Math.pow(dx, 3)+dx*dy*dy)/Math.pow(dist, 7)); F[1] = 3.0*C*m*m*( (2.0*yA*yA*dy + yA*(dx*xA + dy*yA) + xA*xA*dy + xA*yA*dx)/Math.pow(dist, 5) - 5*(dx*xA + dy*yA) *(Math.pow(dy, 3)+dx*dx*dy)/Math.pow(dist, 7)); return F; } //This double calculates angular acceleration. public static double[] alpha ( double Theta1, double Theta2 ) { double[] Alpha = new double [2]; //Angular acceleration //returned to the advance() //function. x11 = Math.cos( Theta1 ); y11 = Math.sin( Theta1 ); x21 = Math.cos( Theta2 ); y21 = Math.sin( Theta2 ); x12 = x11*cos120 - y11*sin120; x13 = x11*cos240 - y11*sin240; y12 = y11*cos120 + x11*sin120; y13 = y11*cos240 + x11*sin240; x22 = x21*cos120 - y21*sin120; x23 = x21*cos240 - y21*sin240; y22 = y21*cos120 + x21*sin120; y23 = y21*cos240 + x21*sin240; F_11on21 = F_AonB( x11, y11, x21, y21 ); F_11on22 = F_AonB( x11, y11, x22, y22 ); F_11on23 = F_AonB( x11, y11, x23, y23 ); F_12on21 = F_AonB( x12, y12, x21, y21 ); F_12on22 = F_AonB( x12, y12, x22, y22 ); F_12on23 = F_AonB( x12, y12, x23, y23 ); F_13on21 = F_AonB( x13, y13, x21, y21 ); F_13on22 = F_AonB( x13, y13, x22, y22 ); F_13on23 = F_AonB( x13, y13, x23, y23 ); Alpha[0] = -x11 * ( F_11on21[1] + F_11on22[1] + F_11on23[1] ) + y11 * ( F_11on21[0] + F_11on22[0] + F_11on23[0] ) - x12 * ( F_12on21[1] + F_12on22[1] + F_12on23[1] ) + y12 * ( F_12on21[0] + F_12on22[0] + F_12on23[0] ) - x13 * ( F_13on21[1] + F_13on22[1] + F_13on23[1] ) + y13 * ( F_13on21[0] + F_13on22[0] + F_13on23[0] ); Alpha[1] = x21 * ( F_11on21[1] + F_12on21[1] + F_13on21[1] ) - y21 * ( F_11on21[0] + F_12on21[0] + F_13on21[0] ) + x22 * ( F_11on22[1] + F_12on22[1] + F_13on22[1] ) - y22 * ( F_11on22[0] + F_12on22[0] + F_13on22[0] ) + x23 * ( F_11on23[1] + F_12on23[1] + F_13on23[1] ) - y23 * ( F_11on23[0] + F_12on23[0] + F_13on23[0] ); Alpha[0] = Alpha[0]/I; //Alpha calculations were actually torques. Alpha[1] = Alpha[1]/I; return Alpha; } //B_1at2 calculates the magnetic field due to one of the magnets on the //first rotor at the position of one of the magnets on the second rotor. //Its only purpose is to be called by the energy function. public static double[] B_1at2 ( double xA, double yA, double xB, double yB ) { double[] B = new double[2]; double dx = D + xB - xA; double dy = yB - yA; double dist = Math.sqrt(dx*dx + dy*dy); B[0] = C*m*( 3 * dx * (dx*xB + dy*yB)/(dist*dist*dist) - xB/dist ); B[1] = C*m*( 3 * dy * (dx*xB + dy*yB)/(dist*dist*dist) - yB/dist ); return B; } public static double energy ( double t1, double t2, double o1, double o2 ) { double U21, U22, U23, KE; x11 = Math.cos( t1 ); y11 = Math.sin( t1 ); x21 = Math.cos( t2 ); y21 = Math.sin( t2 ); x12 = x11*cos120 - y11*sin120; x13 = x11*cos240 - y11*sin240; y12 = y11*cos120 + x11*sin120; y13 = y11*cos240 + x11*sin240; x22 = x21*cos120 - y21*sin120; x23 = x21*cos240 - y21*sin240; y22 = y21*cos120 + x21*sin120; y23 = y21*cos240 + x21*sin240; B_11at21 = B_1at2( x11, y11, x21, y21 ); B_11at22 = B_1at2( x11, y11, x22, y22 ); B_11at23 = B_1at2( x11, y11, x23, y23 ); B_12at21 = B_1at2( x12, y12, x21, y21 ); B_12at22 = B_1at2( x12, y12, x22, y22 ); B_12at23 = B_1at2( x12, y12, x23, y23 ); B_13at21 = B_1at2( x13, y13, x21, y21 ); B_13at22 = B_1at2( x13, y13, x22, y22 ); B_13at23 = B_1at2( x13, y13, x23, y23 ); U21 = -m * x21 * ( B_11at21[0] + B_12at21[0] + B_13at21[0] ) - m * y21 * ( B_11at21[1] + B_12at21[1] + B_13at21[1] ); U22 = -m * x22 * ( B_11at22[0] + B_12at22[0] + B_13at22[0] ) - m * y22 * ( B_11at22[1] + B_12at22[1] + B_13at22[1] ); U23 = -m * x23 * ( B_11at23[0] + B_12at23[0] + B_13at23[0] ) - m * y23 * ( B_11at23[1] + B_12at23[1] + B_13at23[1] ); KE = I * ( o1*o1 + o2*o2 ); return U21 + U22 + U23 + KE; } //Poincare does a linear interpolation of the position and angular velocity //of the first rotor as the second rotor changes from negative to positive //angular velocity. public static void Poincare ( double t1, double o1, double o2, double t1old, double o1old, double o2old ) { double change = o2old / ( o2old - o2 ); double T1 = change * ( t1 - t1old ); double O1 = change * ( o1 - o1old ); System.out.println( T1 + "," + O1 ); } //These are the coefficients for the integrator. static final double a1 = 0.25; static final double a2 = 0.375; static final double a3 = 12.0 / 13.0; static final double a4 = 1.0; static final double a5 = 0.5; static final double b10 = 0.25; static final double b20 = 3 / 32.0; static final double b21 = 9 / 32.0; static final double b30 = 1932.0 / 2197.0; static final double b31 = -7200.0 / 2197.0; static final double b32 = 7296.0 / 2197.0; static final double b40 = 439.0 / 216.0; static final double b41 = -8; static final double b42 = 3680.0 / 513.0; static final double b43 = -845.0 / 4104.0; static final double b50 = -8.0 / 27.0; static final double b51 = 2; static final double b52 = -3544.0 / 2565.0; static final double b53 = 1859.0 / 4104.0; static final double b54 = -11.0 / 40.0; static final double c0 = 16.0 / 135.0; static final double c2 = 6656.0 / 12825.0; static final double c3 = 28561.0 / 56430.0; static final double c4 = -0.18; static final double c5 = 2.0 / 55.0; static final double d0 = 1.0 / 360.0; static final double d2 = -128.0 / 4275.0; static final double d3 = -2197.0 / 75240.0; static final double d4 = 0.02; static final double d5 = 2.0 / 55.0; }