/* sim1.java contains several CSimulation subclasses */ /* Part of the www.MyPhysicsLab.com physics simulation applet. Copyright (c) 2001 Erik Neumann This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Contact Erik Neumann at erikn@MyPhysicsLab.com or 610 N. 65th St. Seattle WA 98103 */ import java.applet.*; ///////////////////////////////////////////////////////////////////////////// // CSpringMass1 class // // One spring connected to one mass // // /* origin = connection of spring to wall vars[0] = position (x) with origin as above vars[1] = velocity (v=x') R = rest length len = current length of spring = x - origin.x L = how much spring is stretched from rest length L = len - R = x - origin.x - R k = spring constant b = damping constant F = -kL -bv = -k(x - origin.x - R) -bv = m v' so diffeq's are: x' = v v' = -(k/m)(x - origin.x - R) -(b/m)v */ class CSpringMass1 extends CSimulation { private CSpring m_Spring; private CMass m_Mass; private CWall m_Wall; //private CText tx; public double diffeq0(double t, double[] x) // time & array of variables { return x[1]; // x' = v } public double diffeq1(double t, double[] x) // time & array of variables { // v' = -(k/m)(x - R) - (b/m) v double r = -m_Spring.m_SpringConst*(x[0] - m_Spring.m_X1 - m_Spring.m_RestLength)- m_Mass.m_Damping*x[1]; return r/m_Mass.m_Mass; } public CSpringMass1(Applet app) { super(app); var_names = new String[] { "position", "velocity", "acceleration", // calculated in collisionTime "kinetic energy", "spring energy", "total energy" }; //add(tx = new CText(3, 1, "energy ")); numVars = 2; // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y map = new CoordMap(CoordMap.INCREASE_DOWN, -0.5, 6.0, -1.5, 1.5, CoordMap.ALIGN_LEFT, CoordMap.ALIGN_MIDDLE); m_Wall = new CWall (-0.4, -1.5, 0, 1.5, 0); //x1, y1, x2, y2, angle add(m_Wall); m_Spring = new CSpring (0, 0, 2.5, 0.5); // x1, y1, restLen, thickness m_Spring.m_SpringConst=3; add(m_Spring); double w = 0.3; /* x1, y1, width, height, drawmode */ m_Mass = new CMass(0.5, -w, 2*w, 2*w, CElement.MODE_RECT); m_Spring.setX2(m_Mass.m_X1); m_Mass.m_Mass = 0.5; m_Mass.m_Damping = 0; add(m_Mass); resetVariables(null); params = new String[] {"mass of block","damping","spring stiffness","spring rest length"}; } public double get(int param) { switch(param) { case 0: return m_Mass.m_Mass; case 1: return m_Mass.m_Damping; case 2: return m_Spring.m_SpringConst; case 3: return m_Spring.m_RestLength; default: return 0; } } public void set(int param, double value) { switch(param) { case 0: m_Mass.m_Mass = value; break; case 1: m_Mass.m_Damping = value; break; case 2: m_Spring.m_SpringConst = value; break; case 3: m_Spring.m_RestLength = value; break; } } public double[] getRange(int var) { double[] r = new double[2]; switch(var) { case 0: r[0]=-1; r[1]=5; break; // position case 1: r[0]=-8; r[1]=8; break; // velocity case 2: r[0]=-12; r[1]=12; break; // acceleration case 3: r[0]=-1; r[1]=10; break; // energy case 4: r[0]=-1; r[1]=10; break; // energy case 5: r[0]=-1; r[1]=10; break; // energy } return r; } public void resetVariables(CElement p) { super.resetVariables(p); //starting displacement vars[0] = m_Mass.m_X1; vars[1] = 0; // velocity is zero at start } public void doModifyObjects() { m_Mass.setX1(vars[0]); m_Spring.setX2(m_Mass.m_X1); } public void startDrag(CElement e) { if (e==m_Mass) { calc[0] = false; calc[1] = false; } } public void constrainedSet(CElement e, double x, double y) { if (e==m_Mass) { // don't allow vertical dragging, so only set horizontal component //m_Mass.setX1(x); //m_Spring.setX2(x); // force spring to follow along vars[0] = x; vars[1] = 0; } // objects other than mass are not allowed to be dragged } public double collisionTime(double[] ov, double h) { // use this to set derived variables /* here is a way to get acceleration numerically from the last two values. Problem is its accel for the time in between the two values, so to plot against it you also have to figure out the average of other values. // acceleration: vars[2] = (vars[1] - ov[1])/h; // position at time t-1/2 (average of current and old position) vars[3] = (vars[0] + ov[0])/2; */ // instead, we use the diffeq to get the value of acceleration // at the current time... works great! vars[2] = diffeq1(0, vars); // acceleration // kinetic energy = 0.5*m*v^2 vars[3] = 0.5*m_Mass.m_Mass*vars[1]*vars[1]; // potential = 0.5*k*x^2 where k = spring const, x = spring stretch double x = vars[0] - m_Spring.m_X1 - m_Spring.m_RestLength; vars[4] = 0.5*m_Spring.m_SpringConst*x*x; // total energy vars[5] = vars[3] + vars[4]; //tx.setNumber(vars[5]); return 99999999; } } // end class CSpringMass1 ///////////////////////////////////////////////////////////////////////////// // CSpringMass2 class // Wall, 2 Springs & 2 Masses connected as W-S1-M1-S2-M2 /* origin = connection of spring1 to wall vars[0] = u1 = x1 position of mass1 with origin as above vars[1] = u2 = x1 position of mass2 with origin as above vars[2] = v1 = velocity of mass1 vars[3] = v2 = velocity of mass2 Abbreviations for diff eqs R = rest length of spring w = width of mass k = spring constant m = mass u = mass.x1 position L = how much spring is stretched F = force len = current length of spring L1 = u1 - R1 len2 = u2 - (u1 + w1) L2 = len2 - R2 = u2 - (u1 + w1) - R2 F1 = -k1 L1 + k2 L2 = m1 v1' F2 = -k2 L2 = m2 v2' so the diff eq's are u1' = v1 u2' = v2 v1' = -(k1/m1) (u1 - R1) + (k2/m1) (u2 - (u1 + w1) - R2) v2' = -(k2/m2) (u2 - (u1 + w1) - R2) */ class CSpringMass2 extends CSimulation { private CSpring m_Spring1, m_Spring2; private CMass m_Mass1, m_Mass2; private CWall m_Wall; // dx1 public double diffeq0(double t, double[] x) { return x[2]; // x1' = v1 } // dx2 public double diffeq1(double t, double[] x) { return x[3]; // x2' = v2 } // dv1 public double diffeq2(double t, double[] x) { // v1' = -(k1/m1) (u1 - R1) + (k2/m1) (u2 - (u1 + w1) - R2) double r1 = -(m_Spring1.m_SpringConst / m_Mass1.m_Mass)* (x[0] - m_Spring1.m_RestLength); double r2 = (m_Spring2.m_SpringConst / m_Mass1.m_Mass)* (x[1] - x[0] - m_Mass1.m_Width - m_Spring2.m_RestLength); return r1+r2; } // dv2 public double diffeq3(double t, double[] x) { // v2' = -(k2/m2) (u2 - (u1 + w1) - R2) return -(m_Spring2.m_SpringConst / m_Mass2.m_Mass)* (x[1] - x[0] - m_Mass1.m_Width - m_Spring2.m_RestLength); } public CSpringMass2(Applet app) { super(app); numVars = 4; var_names = new String[] { "position1", "position2", "velocity1", "velocity2", "acceleration1", "acceleration2", "total energy" }; // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y map = new CoordMap(CoordMap.INCREASE_DOWN, -0.5, 9.5, -1.5, 1.5, CoordMap.ALIGN_LEFT, CoordMap.ALIGN_MIDDLE); m_Wall = new CWall (-0.3, -1.5, 0, 1.5, 0); //x1, y1, x2, y2, angle add(m_Wall); m_Spring1 = new CSpring (0, 0, 2, 0.5); // x1, y1, restLen, thickness m_Spring1.m_SpringConst=6; add(m_Spring1); double w = 0.75; /* x1, y1, width, height, drawmode */ m_Mass1 = new CMass(m_Spring1.m_X2+0.5, -w/2, w, w, CElement.MODE_RECT); m_Mass1.m_Mass = 1; add(m_Mass1); m_Spring2 = new CSpring (m_Mass1.m_X2, 0, 2.0, 0.5); // x1, y1, restLen, height m_Spring2.m_SpringConst=6; add(m_Spring2); m_Mass2 = new CMass(m_Spring2.m_X2+1.5, -w/2, w, w, CElement.MODE_RECT); // x1, y1, height, width m_Mass2.m_Mass = 1; add(m_Mass2); resetVariables(null); params = new String[] {"block 1 mass", "spring 1 stiffness", "spring 1 rest length","block 2 mass", "spring 2 stiffness", "spring 2 rest length" }; } public double get(int param) { switch(param) { case 0: return m_Mass1.m_Mass; case 1: return m_Spring1.m_SpringConst; case 2: return m_Spring1.m_RestLength; case 3: return m_Mass2.m_Mass; case 4: return m_Spring2.m_SpringConst; case 5: return m_Spring2.m_RestLength; default: return 0; } } public void set(int param, double value) { switch(param) { case 0: m_Mass1.m_Mass = value; break; case 1: m_Spring1.m_SpringConst = value; break; case 2: m_Spring1.m_RestLength = value; break; case 3: m_Mass2.m_Mass = value; break; case 4: m_Spring2.m_SpringConst = value; break; case 5: m_Spring2.m_RestLength = value; break; } } public double[] getRange(int var) { double[] r = new double[2]; switch(var) { case 0: r[0]=0; r[1]=5; break; // position of mass 1 case 1: r[0]=2; r[1]=7; break; // position of mass 2 case 2: r[0]=-8; r[1]=8; break; // velocity of mass 1 case 3: r[0]=-8; r[1]=8; break; // velocity of mass 2 case 4: r[0]=-8; r[1]=8; break; // acceleration of mass 1 case 5: r[0]=-8; r[1]=8; break; // acceleration of mass 2 case 6: r[0]=-1; r[1]=10; break; // total energy } return r; } public void stop() { vars[0] = m_Spring1.m_RestLength; // position 1 vars[1] = vars[0] + m_Mass1.m_Width + m_Spring2.m_RestLength; // position 2 vars[2] = 0; // velocity 1 vars[3] = 0; // velocity 2 } public void resetVariables(CElement p) { super.resetVariables(p); //starting displacement vars[0] = m_Mass1.m_X1; vars[1] = m_Mass2.m_X1; vars[2] = 0; // v1 = 0 at start vars[3] = 0; // v2 = 0 } public void doModifyObjects() { m_Mass1.setX1(vars[0]); m_Spring1.setX2(m_Mass1.m_X1); m_Mass2.setX1(vars[1]); m_Spring2.setX1(m_Mass1.m_X2); m_Spring2.setX2(m_Mass2.m_X1); } public void startDrag(CElement e) { if (e==m_Mass1) { calc[0] = false; calc[2] = false; } else if (e==m_Mass2) { calc[1] = false; calc[3] = false; } } public void constrainedSet(CElement e, double x, double y) { if (e==m_Mass1) { // don't allow vertical dragging, so only set horizontal component //m_Mass1.setX1(x); //m_Spring1.setX2(x); // force spring to follow along //m_Spring2.setX1(x + m_Mass1.m_Width); vars[0] = x; vars[2] = 0; // velocity } else if (e==m_Mass2) { //m_Mass2.setX1(x); // only horizontal dragging allowed //m_Spring2.setX2(x); // force spring to follow vars[1] = x; vars[3] = 0; // velocity } // objects other than mass are not allowed to be dragged } public double collisionTime(double[] ov, double h) { // use this to set derived variables vars[4] = diffeq2(0, vars); // acceleration1 vars[5] = diffeq3(0, vars); // acceleration2 double k = 0.5*(m_Mass1.m_Mass*vars[2]*vars[2] + m_Mass2.m_Mass*vars[3]*vars[3]); double x1 = vars[0] - m_Spring1.m_RestLength; double x2 = vars[1] - vars[0] - m_Mass1.m_Width - m_Spring2.m_RestLength; double p = 0.5*(m_Spring1.m_SpringConst*x1*x1 + m_Spring2.m_SpringConst*x2*x2); vars[6] = k + p; // total energy return 9999999; } } // end class CSpringMass2 ///////////////////////////////////////////////////////////////////////////// // CPendulum1 class // // 2D pendulum // /* mass is suspended from ceiling on a stick origin = connection point of stick to ceiling, with y increasing downwards th = angle formed with vertical, positive is counter clockwise U = position of CENTER of mass v = velocity of angle = d(th)/dt m = mass of mass g = gravity constant L = length of rope b = friction constant A = amplitude of driving force k = related to frequency of driving force Note: we use 2*radius of arc as driving force Regard th as the only degree of freedom of the system (one-dimensional). there is no acceleration in the direction of the rope the only acceleration is perpendicular to the rope, Use the rotational analog of Newton's second law: Sum(torques) = I a where I = rotational inertia, and a = angular acceleration. Rotational inertia = mL^2 Torque due to gravity is -Lmg sin(th) Torque due to friction is -b v Torque due to driving force is A cos(w) where A is constant amplitude and w = k t is a linear function of time. Then we have -Lmg sin(th) -b v +A cos(w) = mL^2 a If our variables are as above, we get the equations th' = v v' = -(g/L) sin(th) -(b/mL^2) v + (A/mL^2) cos(k t) Compare to equation 3.1 in Chaotic Dynamics by Baker/Gollub. (I've translated that equation to equivalent variables here...) v' = - sin(th) - v/q + A cos(k t) The range of chaos is: q=2, 0.5 Math.PI) vars[0] = vars[0] - 2*Math.PI*Math.floor(vars[0]/Math.PI); else if (vars[0] < -Math.PI) vars[0] = vars[0] - 2*Math.PI*Math.ceil(vars[0]/Math.PI); // set the position of the pendulum according to the angle double len = m_Spring.m_RestLength; double w = m_Mass.m_Width/2; m_Mass.setX1(len*Math.sin(vars[0]) - w); m_Mass.setY1(len*Math.cos(vars[0]) - w); m_Spring.setX2(m_Mass.m_X1 + w); m_Spring.setY2(m_Mass.m_Y1 + w); // show the driving torque as a line circling about origin double t = m_DriveFrequency*sim_time; // note: sim_time is actually the old time... slightly inaccurate // angle is the angle from the startAngle, so from -90 to 90 degrees t = 180*t/Math.PI; // convert to degrees, starting at 0 t = t - 360 *Math.floor(t/360); // mod 360, range is 0 to 360 // here we generate a ramp that works as follows: // we want to represent cos(k t) // 0 90 180 270 360 // 90 0 -90 0 90 if ((t>0) && (t<=180)) // 0 to 180 is reversed and offset t = 90 - t; else t = t - 270; //if ((t>90) && (t<270)) // t = 180 - t; // 90 to 270 is mapped to 90 to -90 //else if (t>270) // t = t - 360; // 270 to 360 mapped to -90 to 0 m_Drive.m_Angle = t; } public void startDrag(CElement e) { if (e==m_Mass) { calc[0] = false; calc[1] = false; } } public void constrainedSet(CElement e, double x, double y) { if (e==m_Mass) { // only allow movement along circular arc // use line from origin to mouse position double len = m_Spring.m_RestLength; double w = m_Mass.m_Width/2; // calculate angle theta given current mass position & width & origin setting, etc. double xx = x + w; double yy = y + w; double th = Math.atan2(xx, yy); vars[0] = th; vars[1] = 0; } // objects other than mass are not allowed to be dragged } public double collisionTime(double[] ov, double h) { // use this to set derived variables vars[2] = diffeq1(0, vars); // acceleration return 99999999; } } // end class CPendulum1 ///////////////////////////////////////////////////////////////////////////// // CSpring2D class // // An immoveable but draggable mass with a spring and mass hanging below // and swinging in 2D. // /* 2-D spring simulation with gravity spring is suspended from top mass origin = topleft corner th = angle formed with vertical, positive is counter clockwise L = displacement of spring from rest length R = rest length U = position of CENTER of mass S = position of spring's X1,Y1 point V = velocity of mass k = spring constant b = damping constant m = mass of mass w = width (radius) of mass Fx = -kL sin(th) -bVx = m Vx' Fy = mg -kL cos(th) -bVy = m Vy' xx = Ux - Sx yy = Uy - Sy len = Sqrt(xx^2+yy^2) L = len - R th = atan(xx/yy) cos(th) = yy / len sin(th) = xx / len so here are the four variables and their diffeq's vars[0] = Ux vars[1] = Uy vars[2] = Vx vars[3] = Vy Ux' = Vx Uy' = Vy Vx' = -(k/m)L sin(th) -(b/m)Vx Vy' = g - (k/m)L cos(th) -(b/m)Vy */ class CSpring2D extends CSimulation { private CSpring m_Spring; private CMass m_Mass, m_TopMass; private double m_Gravity = 9.8; // dUx public double diffeq0(double t, double[] x) { return x[2]; //Ux' = Vx } //dUy public double diffeq1(double t, double[] x) { return x[3]; //Uy' = Vy } // dVx public double diffeq2(double t, double[] x) { //xx = Ux - Sx //yy = Uy - Sy //len = Sqrt(xx^2+yy^2) double xx = x[0] - m_Spring.m_X1; double yy = x[1] - m_Spring.m_Y1; double len = Math.sqrt(xx*xx + yy*yy); double m = m_Mass.m_Mass; //L = len - R //sin(th) = xx / len //Vx' = -(k/m)L sin(th) double r = -(m_Spring.m_SpringConst/m)* (len - m_Spring.m_RestLength) * xx / len; // damping: - (b/m) Vx double b = m_Mass.m_Damping; if (b!=0) r -= (b/m)*x[2]; return r; } // dVy public double diffeq3(double t, double[] x) { //xx = Ux - Sx //yy = Uy - Sy //len = Sqrt(xx^2+yy^2) double xx = x[0] - m_Spring.m_X1; double yy = x[1] - m_Spring.m_Y1; double len = Math.sqrt(xx*xx + yy*yy); double m = m_Mass.m_Mass; //L = len - R //cos(th) = yy / len //Vy' = g - (k/m)L cos(th) double r = m_Gravity - (m_Spring.m_SpringConst/m)* (len - m_Spring.m_RestLength) * yy / len; // damping: - (b/m) Vy double b = m_Mass.m_Damping; if (b!=0) r -= (b/m)*x[3]; return r; } public CSpring2D(Applet app) { super(app); numVars = 4; var_names = new String[] { "x position", "y position", "x velocity", "y velocity" }; // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y map = new CoordMap(CoordMap.INCREASE_DOWN, -6, 6, -6, 6, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE); double xx = 0; double yy = -2; m_Spring = new CSpring (xx, yy, 2.5, 0.6); // x1, y1, restLen, height m_Spring.setX2(xx); double s = yy + m_Spring.m_RestLength; m_Spring.setY2(s); m_Spring.m_SpringConst=6; add(m_Spring); double w = 1; m_Mass = new CMass(xx-w/2, s-w/2, w, w, CElement.MODE_CIRCLE); m_Mass.m_Mass = .5; // set the position of the mass to be where spring stretch balances weight s = s + m_Mass.m_Mass*m_Gravity/m_Spring.m_SpringConst; m_Mass.setY1(s-w/2); m_Mass.setX1(-2); // so it starts springing m_Mass.setY1(s+1); m_Mass.m_Damping = 0; add(m_Mass); m_TopMass = new CMass(xx-w/2, yy-w, w, w, CElement.MODE_RECT); add(m_TopMass); resetVariables(null); params = new String[] {"mass", "damping", "spring stiffness", "spring rest length"}; } public double get(int param) { switch(param) { case 0: return m_Mass.m_Mass; case 1: return m_Mass.m_Damping; case 2: return m_Spring.m_SpringConst; case 3: return m_Spring.m_RestLength; default: return 0; } } public void set(int param, double value) { switch(param) { case 0: m_Mass.m_Mass = value; break; case 1: m_Mass.m_Damping = value; break; case 2: m_Spring.m_SpringConst = value; break; case 3: m_Spring.m_RestLength = value; break; } } public double[] getRange(int var) { double[] r = new double[2]; switch(var) { case 0: r[0]=-5; r[1]=5; break; // x position case 1: r[0]=-5; r[1]=7; break; // y position case 2: r[0]=-8; r[1]=8; break; // x velocity case 3: r[0]=-8; r[1]=8; break; // y velocity } return r; } public void stop() { vars[0] = m_TopMass.m_X1 + m_TopMass.m_Width/2; vars[1] = m_TopMass.m_Y2 + m_Spring.m_RestLength + m_Mass.m_Mass*m_Gravity/m_Spring.m_SpringConst; vars[2] = 0; vars[3] = 0; } public void resetVariables(CElement p) { super.resetVariables(p); vars[0] = m_Mass.m_X1 + m_Mass.m_Width/2; vars[1] = m_Mass.m_Y1 + m_Mass.m_Height/2; } public void doModifyObjects() { m_Mass.setX1(vars[0] - m_Mass.m_Width/2); m_Mass.setY1(vars[1] - m_Mass.m_Height/2); m_Spring.setX2(m_Mass.m_X1 + m_Mass.m_Width/2); m_Spring.setY2(m_Mass.m_Y1 + m_Mass.m_Height/2); } public void startDrag(CElement e) { if (e==m_Mass) { calc[0] = false; calc[1] = false; calc[2] = false; calc[3] = false; } } public void constrainedSet(CElement e, double x, double y) { if (e==m_TopMass) { e.setX1(x); e.setY1(y); m_Spring.setX1(x + m_TopMass.m_Width/2); // force spring to follow along m_Spring.setY1(y + m_TopMass.m_Height); } else if (e==m_Mass) { vars[0] = x + m_Mass.m_Width/2; vars[1] = y + m_Mass.m_Height/2; vars[2] = 0; vars[3] = 0; } } } // end class CSpring2D ///////////////////////////////////////////////////////////////////////////// // CDoubleSpring2D class // // An immoveable but draggable mass with a 2 springs and 2 masses hanging below // and swinging in 2D. // /* 2-D spring simulation with gravity spring is suspended from top mass origin = topleft corner th = angle formed with vertical, positive is counter clockwise L = displacement of spring from rest length R = rest length U = position of CENTER of mass S1 = position of spring1 X1,Y1 point V = velocity of mass k = spring constant b = damping constant m = mass of mass w = width (radius) of mass F1x = -k1 L1 sin(th1) +k2 L2 sin(th2) -b1 V1x = m1 V1x' F1y = m1 g -k1 L1 cos(th1) +k2 L2 cos(th2) -b1 V1y = m1 V1y' F2x = -k2 L2 sin(th2) -b2 V2x = m2 V2x' F2y = m2 g -k2 L2 cos(th2) -b2 V2y = m2 V2y' xx = U1x - S1x yy = U1y - S1y len1 = Sqrt(xx^2+yy^2) L1 = len1 - R1 th1 = atan(xx/yy) cos(th1) = yy / len1 sin(th1) = xx / len1 xx2 = U2x - U1x yy2 = U2y - U1y len2 = sqrt(xx2^2+yy2^2) L2 = len2 - R2 cos(th2) = yy2 / len2 sin(th2) = xx2 / len2 here are the variables vars[0] = U1x vars[1] = U1y vars[2] = U2x vars[3] = U2y vars[4] = V1x vars[5] = V1y vars[6] = V2x vars[7] = V2y */ class CDoubleSpring2D extends CSimulation { private CSpring m_Spring1, m_Spring2; private CMass m_Mass1, m_Mass2, m_TopMass; private double m_Gravity = 9.8; // dU1x public double diffeq0(double t, double[] x) { return x[4]; //U1x' = V1x } // dU1y public double diffeq1(double t, double[] x) { return x[5]; //U1y' = V1y } // dU2x public double diffeq2(double t, double[] x) { return x[6]; //U2x' = V2x } // dU2y public double diffeq3(double t, double[] x) { return x[7]; //U2y' = V2y } // dV1x public double diffeq4(double t, double[] x) { // F1x = -k1 L1 sin(th1) +k2 L2 sin(th2) -b1 V1x = m1 V1x' //xx = U1x - S1x //yy = U1y - S1y //len1 = Sqrt(xx^2+yy^2) double xx = x[0] - m_Spring1.m_X1; double yy = x[1] - m_Spring1.m_Y1; double len1 = Math.sqrt(xx*xx + yy*yy); double m1 = m_Mass1.m_Mass; //L1 = len1 - R1 //sin(th1) = xx / len1 // -(k1/m1)L1 sin(th1) double r = -(m_Spring1.m_SpringConst/m1)* (len1 - m_Spring1.m_RestLength) * xx / len1; //xx2 = U2x - U1x //yy2 = U2y - U1y //len2 = sqrt(xx2^2+yy2^2) double xx2 = x[2] - x[0]; double yy2 = x[3] - x[1]; double len2 = Math.sqrt(xx2*xx2 + yy2*yy2); //L2 = len2 - R2 //sin(th2) = xx2 / len2 //+(k2/m1) L2 sin(th2) r += (m_Spring2.m_SpringConst/m1)* (len2 - m_Spring2.m_RestLength) * xx2 / len2; // damping: - (b1/m1) V1x double b1 = m_Mass1.m_Damping; if (b1!=0) r -= (b1/m1)*x[4]; return r; } // dV1y public double diffeq5(double t, double[] x) { //F1y = m1 g -k1 L1 cos(th1) +k2 L2 cos(th2) -b1 V1y = m1 V1y' //xx = U1x - S1x //yy = U1y - S1y //len1 = Sqrt(xx^2+yy^2) double xx = x[0] - m_Spring1.m_X1; double yy = x[1] - m_Spring1.m_Y1; double len1 = Math.sqrt(xx*xx + yy*yy); double m1 = m_Mass1.m_Mass; //L1 = len1 - R1 //cos(th1) = yy / len1 // g -(k1/m1)L1 cos(th1) double r = m_Gravity -(m_Spring1.m_SpringConst/m1)* (len1 - m_Spring1.m_RestLength) * yy / len1; //xx2 = U2x - U1x //yy2 = U2y - U1y //len2 = sqrt(xx2^2+yy2^2) double xx2 = x[2] - x[0]; double yy2 = x[3] - x[1]; double len2 = Math.sqrt(xx2*xx2 + yy2*yy2); //L2 = len2 - R2 //cos(th2) = yy2 / len2 //+(k2/m1) L2 cos(th2) r += (m_Spring2.m_SpringConst/m1)* (len2 - m_Spring2.m_RestLength) * yy2 / len2; // damping: - (b1/m1) V1y double b1 = m_Mass1.m_Damping; if (b1!=0) r -= (b1/m1)*x[5]; return r; } // dV2x public double diffeq6(double t, double[] x) { // F2x = -k2 L2 sin(th2) -b2 V2x = m2 V2x' //xx2 = U2x - U1x //yy2 = U2y - U1y //len2 = sqrt(xx2^2+yy2^2) double xx2 = x[2] - x[0]; double yy2 = x[3] - x[1]; double len2 = Math.sqrt(xx2*xx2 + yy2*yy2); double m2 = m_Mass2.m_Mass; //L2 = len2 - R2 //sin(th2) = xx2 / len2 //-(k2/m2) L2 sin(th2) double r = -(m_Spring2.m_SpringConst/m2)* (len2 - m_Spring2.m_RestLength) * xx2 / len2; // damping: - (b2/m2) V2x double b2 = m_Mass2.m_Damping; if (b2!=0) r -= (b2/m2)*x[6]; return r; } // dV2y public double diffeq7(double t, double[] x) { // F2y = m2 g -k2 L2 cos(th2) -b2 V2y = m2 V2y' //xx2 = U2x - U1x //yy2 = U2y - U1y //len2 = sqrt(xx2^2+yy2^2) double xx2 = x[2] - x[0]; double yy2 = x[3] - x[1]; double len2 = Math.sqrt(xx2*xx2 + yy2*yy2); double m2 = m_Mass2.m_Mass; //L2 = len2 - R2 //cos(th2) = yy2 / len2 //g -(k2/m2) L2 cos(th2) double r = m_Gravity -(m_Spring2.m_SpringConst/m2)* (len2 - m_Spring2.m_RestLength) * yy2 / len2; // damping: - (b2/m2) V2y double b2 = m_Mass2.m_Damping; if (b2!=0) r -= (b2/m2)*x[7]; return r; } public CDoubleSpring2D(Applet app) { super(app); numVars = 8; var_names = new String[] { "x1 position", "y1 position", "x2 position", "y2 position", "x1 velocity", "y1 velocity", "x2 velocity", "y2 velocity" }; // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y map = new CoordMap(CoordMap.INCREASE_DOWN, -6, 6, -6, 6, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE); double xx = 0; double yy = -2; double w = 0.5; m_TopMass = new CMass(xx-w/2, yy-w, w, w, CElement.MODE_RECT); add(m_TopMass); m_Spring1 = new CSpring (xx, yy, 1.0, 0.3); // x1, y1, restLen, thick m_Spring1.setX2(xx); m_Spring1.m_SpringConst=6; add(m_Spring1); m_Mass1 = new CMass(xx-w/2, 0, w, w, CElement.MODE_CIRCLE); m_Mass1.m_Mass = .5; m_Mass2 = new CMass(xx-w/2, 0, w, w, CElement.MODE_CIRCLE); m_Mass2.m_Mass = .5; // set the position of the mass to be where spring stretch balances weight double yy2 = yy + m_Spring1.m_RestLength + (m_Mass1.m_Mass+m_Mass2.m_Mass)*m_Gravity/m_Spring1.m_SpringConst; m_Mass1.setY1(yy2-w/2); m_Spring1.setY2(yy2); m_Mass1.m_Damping = 0; add(m_Mass1); m_Spring2 = new CSpring (xx, yy2, 1.0, 0.3); // x1, y1, restLen, thick m_Spring2.setX2(xx); m_Spring2.m_SpringConst=6; add(m_Spring2); // set the position of the mass to be where spring stretch balances weight double yy3 = yy2 + m_Spring2.m_RestLength + m_Mass2.m_Mass*m_Gravity/m_Spring2.m_SpringConst; m_Mass2.setY1(yy3-w/2); m_Spring2.setY2(yy3); m_Mass2.m_Damping = 0; add(m_Mass2); resetVariables(null); params = new String[] {"mass 1", "damping 1", "spring 1 stiffness", "spring 1 rest length", "mass 2", "damping 2", "spring 2 stiffness", "spring 2 rest length","gravity"}; } public double get(int param) { switch(param) { case 0: return m_Mass1.m_Mass; case 1: return m_Mass1.m_Damping; case 2: return m_Spring1.m_SpringConst; case 3: return m_Spring1.m_RestLength; case 4: return m_Mass2.m_Mass; case 5: return m_Mass2.m_Damping; case 6: return m_Spring2.m_SpringConst; case 7: return m_Spring2.m_RestLength; case 8: return m_Gravity; default: return 0; } } public void set(int param, double value) { switch(param) { case 0: m_Mass1.m_Mass = value; break; case 1: m_Mass1.m_Damping = value; break; case 2: m_Spring1.m_SpringConst = value; break; case 3: m_Spring1.m_RestLength = value; break; case 4: m_Mass2.m_Mass = value; break; case 5: m_Mass2.m_Damping = value; break; case 6: m_Spring2.m_SpringConst = value; break; case 7: m_Spring2.m_RestLength = value; break; case 8: m_Gravity = value; break; } } public double[] getRange(int var) { double[] r = new double[2]; switch(var) { case 0: r[0]=-5; r[1]=5; break; // x1 position case 1: r[0]=-5; r[1]=7; break; // y1 position case 2: r[0]=-5; r[1]=5; break; // x2 position case 3: r[0]=-5; r[1]=7; break; // y2 position case 4: r[0]=-8; r[1]=8; break; // x1 velocity case 5: r[0]=-8; r[1]=8; break; // y1 velocity case 6: r[0]=-8; r[1]=8; break; // x2 velocity case 7: r[0]=-8; r[1]=8; break; // y2 velocity } return r; } public void stop() { double m1 = m_Mass1.m_Mass; double m2 = m_Mass2.m_Mass; double g = m_Gravity; double k1 = m_Spring1.m_SpringConst; double k2 = m_Spring2.m_SpringConst; double r1 = m_Spring1.m_RestLength; double r2 = m_Spring2.m_RestLength; double T = m_TopMass.m_Y2; //x1 & x2 position vars[0] = vars[2] = m_TopMass.m_X1 + m_TopMass.m_Width/2; // derive these by writing the force equations to yield zero accel // when everything is lined up vertically. //y1 position vars[1] = g*(m1+m2)/k1 + r1 + T; //y2 position vars[3] = g*(m2/k2 + (m1+m2)/k1) + r1 + r2 + T; //velocities are all zero vars[4] = 0; vars[5] = 0; vars[6] = 0; vars[7] = 0; } public void resetVariables(CElement p) { super.resetVariables(p); // assume all masses are same width & height double w = m_Mass1.m_Width/2; // calculate the position of the center of the mass vars[0] = m_Mass1.m_X1 + w; vars[1] = m_Mass1.m_Y1 + w; vars[2] = m_Mass2.m_X1 + w; vars[3] = m_Mass2.m_Y1 + w; //vars[4] = 0; // velocity is zero at start //vars[5] = 0; //vars[6] = 0; //vars[7] = 0; } public void doModifyObjects() { // assume all masses are same width & height double w = m_Mass1.m_Width/2; m_Mass1.setX1(vars[0] - w); m_Mass1.setY1(vars[1] - w); m_Mass2.setX1(vars[2] - w); m_Mass2.setY1(vars[3] - w); m_Spring1.setX2(m_Mass1.m_X1 + w); m_Spring1.setY2(m_Mass1.m_Y1 + w); m_Spring2.setX1(m_Mass1.m_X1 + w); m_Spring2.setY1(m_Mass1.m_Y1 + w); m_Spring2.setX2(m_Mass2.m_X1 + w); m_Spring2.setY2(m_Mass2.m_Y1 + w); } public void startDrag(CElement e) { if (e==m_Mass1) { calc[0] = calc[1] = calc[4] = calc[5] = false; } else if (e==m_Mass2) { calc[2] = calc[3] = calc[6] = calc[7] = false; } } public void constrainedSet(CElement e, double x, double y) { // assume all masses are same width & height double w = m_Mass1.m_Width/2; if (e==m_TopMass) { e.setX1(x); e.setY1(y); m_Spring1.setX1(x + m_TopMass.m_Width/2); // force spring to follow along m_Spring1.setY1(y + m_TopMass.m_Height); } else if (e==m_Mass1) { vars[0] = x + w; vars[1] = y + w; vars[4] = vars[5] = 0; /* e.setX1(x); e.setY1(y); m_Spring1.setX2(x + w); // force spring1 to follow along m_Spring1.setY2(y + w); m_Spring2.setX1(x + w); // force spring2 to follow along m_Spring2.setY1(y + w); */ } else if (e==m_Mass2) { vars[2] = x + w; vars[3] = y + w; vars[6] = vars[7] = 0; /* e.setX1(x); e.setY1(y); m_Spring2.setX2(x + w); // force spring to follow along m_Spring2.setY2(y + w); */ } // objects other than mass are not allowed to be dragged } } // end class CDoubleSpring2D ///////////////////////////////////////////////////////////////////////////// // CCollideSpring1 class // // One spring connected to one mass, with another free moving mass, in 1D // /* xxx origin = connection of spring to wall origin = topleft vars[0] = position (x) with origin as above vars[1] = velocity (v=x') vars[2] = x1 position of mass 2 vars[3] = velocity of mass 2 R = rest length S1 = left end of spring S2 = right end of spring (same as x?) len = current length of spring = x - S1.x L = how much spring is stretched from rest length L = len - R = x - S1.x - R k = spring constant b = damping constant F = -kL -bv = -k(x - S1.x - R) -bv = m v' so diffeq's are: x' = v v' = -(k/m)(x - S1.x - R) -(b/m)v */ class CCollideSpring1 extends CSimulation { private CSpring m_Spring; private CMass m_Mass, m_Mass2; private CWall m_Wall, m_Wall2; // dx public double diffeq0(double t, double[] x) // time & array of variables { return x[1]; // x' = v } // dv public double diffeq1(double t, double[] x) // time & array of variables { // v' = -(k/m)(x - S1.x - R) - (b/m) v double r = -m_Spring.m_SpringConst*(x[0] - m_Spring.m_X1 - m_Spring.m_RestLength) - m_Mass.m_Damping*x[1]; return r/m_Mass.m_Mass; } // dx2 public double diffeq2(double t, double[] x) // time & array of variables { return x[3]; // x' = v } // dv2 public double diffeq3(double t, double[] x) // time & array of variables { return 0; // v' = 0 because constant velocity } public CCollideSpring1(Applet app) { super(app); numVars = 4; var_names = new String[] { "position 1", "velocity 1", "position 2", "velocity 2" }; // CoordMap inputs are direction, x1, x2, y1, y2, align_x, align_y map = new CoordMap(CoordMap.INCREASE_DOWN, -0.5, 7.5, -2, 2, CoordMap.ALIGN_LEFT, CoordMap.ALIGN_MIDDLE); double xx = 0.0; double yy = 0.0; m_Wall = new CWall (xx-0.3, yy-2, xx, yy+2, 0); //x1, y1, x2, y2 add(m_Wall); m_Spring = new CSpring (xx, yy, 2.5, 0.4); // x1, y1, restLen, thickness m_Spring.m_SpringConst=6; add(m_Spring); double w = 0.3; m_Mass = new CMass(xx+m_Spring.m_RestLength-2.0, yy-w, 2*w, 2*w, CElement.MODE_RECT); m_Spring.setX2(m_Mass.m_X1); m_Mass.m_Mass = 0.5; m_Mass.m_Damping = 0; m_Mass.ix = 0; // index of the x position variable in vars[] m_Mass.ivx = 1; add(m_Mass); m_Mass2 = new CMass(m_Mass.m_X2+1.0, yy-w, 2*w, 2*w, CElement.MODE_RECT); m_Mass2.m_Mass = 1.5; m_Mass2.ix = 2; m_Mass2.ivx = 3; add(m_Mass2); m_Wall2 = new CWall(7, yy-2, 7.4, yy+2, -90); add(m_Wall2); resetVariables(null); params = new String[] {"mass 1","damping 1","spring stiffness","spring rest length", "mass 2"}; } public double get(int param) { switch(param) { case 0: return m_Mass.m_Mass; case 1: return m_Mass.m_Damping; case 2: return m_Spring.m_SpringConst; case 3: return m_Spring.m_RestLength; case 4: return m_Mass2.m_Mass; default: return 0; } } public void set(int param, double value) { switch(param) { case 0: m_Mass.m_Mass = value; break; case 1: m_Mass.m_Damping = value; break; case 2: m_Spring.m_SpringConst = value; break; case 3: m_Spring.m_RestLength = value; break; case 4: m_Mass2.m_Mass = value; break; } } public double[] getRange(int var) { double[] r = new double[2]; switch(var) { case 0: r[0]=-1; r[1]=5; break; // position 1 case 1: r[0]=-8; r[1]=8; break; // velocity 1 case 2: r[0]=-1; r[1]=7; break; // position 2 case 3: r[0]=-8; r[1]=8; break; // velocity 2 } return r; } public void resetVariables(CElement p) { super.resetVariables(p); //starting displacement vars[0] = m_Mass.m_X1; vars[1] = 0; // velocity is zero at start vars[2] = m_Mass2.m_X1; vars[3] = 0; } public void doModifyObjects() { m_Mass.setX1(vars[0]); m_Spring.setX2(m_Mass.m_X1); m_Mass2.setX1(vars[2]); } public void startDrag(CElement e) { m_Animating = false; } public void constrainedSet(CElement e, double x, double y) { if (e==m_Mass) { if (x < m_Wall.m_X2) // don't allow drag past wall x = m_Wall.m_X2; // don't drag past wall2 if (x + m_Mass.m_Width + m_Mass2.m_Width > m_Wall2.m_X1) x = m_Wall2.m_X1 - m_Mass2.m_Width - m_Mass.m_Width; if (x+m_Mass.m_Width > m_Mass2.m_X1) // move other block m_Mass2.setX1(x+m_Mass.m_Width); m_Mass.setX1(x); // only horizontal dragging allowed m_Spring.setX2(x); // force spring to follow along vars[0] = m_Mass.m_X1; vars[1] = 0; // velocity is zero at start vars[2] = m_Mass2.m_X1; vars[3] = 0; } else if (e==m_Mass2) { if (x+m_Mass2.m_Width > m_Wall2.m_X1) // don't allow drag past wall2 x = m_Wall2.m_X1 - m_Mass2.m_Width; if (x - m_Mass.m_Width < m_Wall.m_X2) // don't allow drag past wall1 x = m_Wall.m_X2 + m_Mass.m_Width; if (x < m_Mass.m_X2) // move other block { m_Mass.setX1(x - m_Mass.m_Width); m_Spring.setX2(x - m_Mass.m_Width); // force spring to follow along } m_Mass2.setX1(x); // only horizontal dragging allowed vars[0] = m_Mass.m_X1; vars[1] = 0; // velocity is zero at start vars[2] = m_Mass2.m_X1; vars[3] = 0; } // objects other than mass 1 and mass 2 are not allowed to be dragged } public double collisionTime(double[] ov, double h) { /* (assume only 2 things can collide at same time) get access to the pre and post RK vectors. assume no intersections in the pre vector check whether there are any intersections in the post vectors (this is n-squared!) for each intersection: Find average of pre & post velocity (or just use pre?) Find velocity of the center of mass. Find velocities of each particle in center of mass frame Assuming constant velocity, find time at which collision would occur Find collision point Reflect velocities across line joining centers of disks (for 2D) (for 2D rects, need torque considerations here). Advance from collision point to current time, using reflected velocity Set dependent variables like springs. */ // vars contains the new values as computed by RK routine and diffeqs // for sim_time + h // ov contains the old values at sim_time // assume no collisions in old values // if we detect a collision in new values, we abandon diffeqs and // assume constant velocity for this time period (is this a problem?) // so start by looking for intersections in new values // for now, be specific for this simulation, generalize later? // note that bounds are not good enough to detect collision, // for example, 2 discs you need the diameter and centers. // (bounds may be helpful as a check beforehand) // want a comparison pattern like this: //5-4, 5-3, 5-2, 5-1, 5-0 //4-3, 4-2, 4-1, 4-0 //3-2, 3-1, 3-0 //2-1, 2-0 //1-0 doModifyObjects(); // ensures that the objects are up-to-date boolean collisionFound = false; int ind = m_oba.size(); // index of last element (origin zero) CElement p1, p2; while (--ind >= 1) // ind is now origin zero, and don't do for the last element { p1 = (CElement)m_oba.elementAt(ind); if (p1.m_Mass > 0) // only for objects with mass { int ind2 = ind; // start with previous object in list while (--ind2 >= 0) { p2 = (CElement)m_oba.elementAt(ind2); if (p2.m_Mass > 0) // only for objects with mass { // first do a bounds intersection test if (p2.m_X2 <= p1.m_X1) continue; if (p2.m_X1 >= p1.m_X2) continue; if (p2.m_Y2 <= p1.m_Y1) continue; if (p2.m_Y1 >= p1.m_Y2) continue; /* A potential collision has been found! Pass the positions & velocities of the two particles to a function which modifies the positions & velocities of the corresponding vars[] vector entries. */ if (DoCollide(ov, h, p1, p2)) collisionFound = true; /* then we can update the vector with new velocities and call doModifyObjects to modify them accordingly. */ } } } } return 9999999; } public boolean DoCollide(double[] ov, double h, CElement p1, CElement p2) { /* (Eventually need to know the types of the objects ...typeid?) (so that we know how to detect a collision between them.) (Could maybe use an overloaded C++ method, with specific types.) */ //Find average of pre & post velocity (or just use pre?) double vx1 = 0; double vx2 = 0; if (p1.ivx != -1) vx1 = (ov[p1.ivx] + vars[p1.ivx])/2; if (p2.ivx != -1) vx2 = (ov[p2.ivx] + vars[p2.ivx])/2; /* Find velocity of the center of mass. For collsion with things that don't move (walls) assume they are so massive that the center of mass isn't moving. */ double vxcm = 0; if ((p1.ivx != -1) && (p2.ivx != -1)) vxcm = ((p1.m_Mass*vx1)+(p2.m_Mass*vx2))/(p1.m_Mass+p2.m_Mass); // Find amount of overlap between the particles // NOTE: this is dependent on types of particles!!! double xgap = 0; if (p1.m_X1 < p2.m_X2) xgap = p2.m_X2 - p1.m_X1; else if (p2.m_X1 < p1.m_X2) xgap = p1.m_X2 - p2.m_X1; else System.out.println( "bad collision data: objects should be overlapping but are not"); // The particles approach each other with velocity = abs(v1-v2) // So the time it took after collision occurred was xgap/abs(v1-v2) double dt = xgap/Math.abs(vx1-vx2); if (dt > h) System.out.println( "bad collision data: calculated time of collsion is out of range"); // adjust the velocities of particles // To find new velocity, find the velocity in the center of mass frame // and reflect it. This works out to -v + 2*vcm. // Here's the derivation: // Velocity in cm frame is v-vcm. // In cm frame, total momentum is zero, after collision momentum is // preserved, so we just reverse signs of each velocity in cm frame. // Reflection of velocity is -(v-vcm) = vcm-v // Add back vcm to get to laboratory frame: = vcm + (vcm-v) = 2*vcm - v if (p1.ivx != -1) vars[p1.ivx] = -vx1 + 2*vxcm; if (p2.ivx != -1) vars[p2.ivx] = -vx2 + 2*vxcm; // adjust the positions of particles // Start with old position, // add distance travelled to collision point = (h-dt) * vx1 // add distance away from collision = dt(-v + 2*vcm) if (p1.ix != -1) vars[p1.ix] = ov[p1.ix] + (h-dt)*vx1 + dt*(-vx1 + 2*vxcm); if (p2.ix != -1) vars[p2.ix] = ov[p2.ix] + (h-dt)*vx2 + dt*(-vx2 + 2*vxcm); return true; // for now assume a collision always happens } } // end class CollideSpring1