A 2D simulation of perfectly elastic collisions of a set of circles with two fixed line segments, with constant vertical acceleration. The circles don't mutually interact or have rotational inertia.
A full cycle of the periodic trajectories of the circles on the left is shown. The other circles bounce out of the trough after several cycles.
Same physics and geometry, except with an initial vertical arrangement of circles. A full cycle of the circle with the longest period (white, top) is shown.
For this second initial configuration, here's a representation of the full range of positions the circles will occupy:
Processing code:
//
// Circle_vs_line_v15
//
// Physics simulation, composed of 2-D elastic collisions of a
// mobile circle with fixed lines.
//
// Mark Dow, April 13, 2016 (v12, multiple non-interacting circles, equal timesteps)
// Mark Dow, April 13, 2016 (v13, tic-toc-tuc geometry)
// Mark Dow, April 14, 2016 (v15, multiple particle logic corrections)
import processing.sound.*;
float TAU = 2*PI;
Line[] Lines = new Line[2];
Circle[] Circs = new Circle[25];
// Circle[] Circs = new Circle[22];
boolean bPaused = false;
boolean bTrails = false;
boolean bStep = false;
boolean bShowOrnaments = false;
float tSim = 0.0; // Simulation time.
// Drawing of all elements occurs at integral boundaries of the time,
// but updates occur if collisions (strong interactions) occur within
// one of these timesteps.
TriOsc triOsc;
Env env;
float attackTime = 0.001;
float sustainTime = 0.0;
float sustainLevel = 0.5;
float releaseTime = 0.0;
void setup() {
// size( 750, 750 );
size( 540, 540 );
//Circs[0] = new Circle( new PVector( 51, -280 ), new PVector( 0, 0 ), color(255,0,0) );
//Circs[1] = new Circle( new PVector( 77, -280 ), new PVector( 0, 0 ), color(255,255,0) );
//Circs[2] = new Circle( new PVector( 103, -280 ), new PVector( 0, 0 ), color(0,255,0) );
//Circs[3] = new Circle( new PVector( 129, -280 ), new PVector( 0, 0 ), color(0,128,255) );
//Circs[4] = new Circle( new PVector( 155, -280 ), new PVector( 0, 0 ), color(128,0,255) );
for ( int iC = 0; iC < Circs.length; iC++ ){
// Circs[iC] = new Circle( new PVector( 51 + 16*iC, -280 ), new PVector( 0, 0 ),
// Circs[iC] = new Circle( new PVector( 50 + 18*iC, -80 ), new PVector( 0, 0 ),
Circs[iC] = new Circle( new PVector( 160, -20 - 12*iC ), new PVector( 0, 0 ),
color( max(0,255-iC*5), max(0,255-iC*10), max(0,255-iC*15) ) );
}
// Corner-well lines.
int xVertex = 320;
int yVertex = -440;
float lengthCWLine0 = 410;
float lengthCWLine1 = 340;
float alpha = 3*TAU/8;
// float beta = TAU/8;
float beta = TAU/8 +.2318240;
float dx1 = lengthCWLine0 * cos( alpha );
float dy1 = lengthCWLine0 * sin( alpha );
float dx2 = lengthCWLine1 * cos( beta );
float dy2 = lengthCWLine1 * sin( beta );
Lines[0] = new Line( new PVector( xVertex, yVertex ), new PVector( xVertex+dx1, yVertex+dy1 ) );
Lines[1] = new Line( new PVector( xVertex, yVertex ), new PVector( xVertex+dx2, yVertex+dy2 ) );
background( 140, 200, 255 );
// Create triangle wave and envelope.
triOsc = new TriOsc(this);
env = new Env(this);
}
void draw() {
if ( !bPaused || bStep ) {
// Find if a collision will occur within the next time step, and if so the
// minimum time to collision.
float tcLineMin = -1;
float tcLineNext;
int iLMin = -1;
int iCMin = -1;
// For each circle:
for ( int iC = 0; iC < Circs.length; iC++ ){
// For each line:
for ( int iL = 0; iL < Lines.length; iL++ ){
// Time to collision with line.
tcLineNext = Circs[iC].TimeToCollisionWithLine( Lines[iL], 1.0, bShowOrnaments );
// If collision will occur within the next time step:
if ( tcLineNext != -1.0 && tcLineNext < 1.0 ) {
// println( " " );
// println( tcLineNext );
// Time to collision with line with better time to collision estimate.
tcLineNext = Circs[iC].TimeToCollisionWithLine( Lines[iL], tcLineNext, bShowOrnaments );
// println( tcLineNext );
// If collision will occur within the next time step and is earliest so far:
if ( ( tcLineNext != -1.0 && tcLineNext < 1.0 )
&& ( tcLineNext < tcLineMin && tcLineMin > 0 ) || tcLineMin < 0 ) {
// New earliest time to collision.
tcLineMin = tcLineNext;
iLMin = iL;
iCMin = iC;
// To Do: Handle simultaneous collisions.
}
}
}
}
float tNext = (float) Math.ceil( tSim ) - tSim;
if ( tNext == 0 ) { tNext = 1.0; }
// If collision will occur within this time step:
if ( iLMin != -1 && tcLineMin <= tNext ) {
//println( tSim - (float) Math.floor( tSim ) );
// For each circle:
for ( int iC = 0; iC < Circs.length; iC++ ){
// Update all cicles to collision timepoint, pre-collision.
Circs[iC].update( tcLineMin );
}
// Update collision, post-collision.
Circs[iCMin].UpdateCollisionWithLine( Lines[iLMin], bShowOrnaments );
tSim += tcLineMin;
// Draw at collision timepoint, not at timestep boundaries.
Circs[iCMin].display( false );
// println( tcLineMin );
// println( tSim );
// Collision sound effect.
triOsc.play();
env.play(triOsc, attackTime, sustainTime, sustainLevel, releaseTime);
}
// Else no collision before timestep boundary:
else {
// For each circle:
for ( int iC = 0; iC < Circs.length; iC++ ){
// Update to timestep boundary.
Circs[iC].update( tNext );
}
tSim += tNext;
if ( bTrails ) {
// Draw fading background.
fill( 140, 200, 255, 5 ); noStroke();
rect( 0, 0, width, height );
}
else {
background( 140, 200, 255 );
}
for( Circle c : Circs ) { c.display( true ); }
// To Do: Do static lines need to be displayed every frame?
for( Line l : Lines ) { l.display(); }
//// Saves each frame as line-000001.png, line-000002.png, etc.
//saveFrame("test-###.png");
}
}
bStep = false;
showKeyCommands();
// println( tSim );
}
void keyPressed() {
if ( key == 't' ) { bTrails = !bTrails; }
else if ( key == 's' ) { bStep = true; }
else if ( key == 'p' ) { bPaused = !bPaused; }
else if ( key == 'o' ) { bShowOrnaments = !bShowOrnaments; }
else {
background( 140, 200, 255 );
setup();
bPaused = true;
// To Do: Update to next timestep boundary before display.
for( Circle c : Circs ) { c.display( true ); }
for( Line l : Lines ) { l.display(); }
}
}
//void mousePressed() {
//
// bPaused = !bPaused;
//}
void showKeyCommands() {
color clrON = color( 255, 150, 150 );
color clrOFF = color( 150, 150, 255 );
color clrGRAY = color( 190, 190, 190 );
fill( clrOFF );
text( "n - new ball", 5, height - 20 );
if ( bPaused ) { fill( clrON ); } else { fill( clrOFF ); }
text( "p - pause", 75, height - 5 );
fill( clrGRAY );
if ( bPaused ) { fill( clrOFF ); }
text( "s - step", 145, height - 5 );
// Show state of trails, by color and/or button rectangle.
if ( bTrails ) { fill( clrON ); } else { fill( clrOFF ); }
text( "t - trails", 5, height - 5 );
// if ( bShowOrnaments ) { fill( clrON ); } else { fill( clrOFF ); }
// text( "o - ornaments", 200, height - 5 );
}
class Circle {
float r;
PVector g;
PVector position;
PVector velocity;
color clrFill;
float phi; // velocity direction angle, where CCW is positive relative to the x-axis
Circle( PVector pos_, PVector vel_, color clr_ ) {
r = 12;
g = new PVector( 0, -.4 );
position = new PVector( pos_.x, pos_.y );
velocity = new PVector( vel_.x, vel_.y );
clrFill = clr_;
}
void update( float tFraction ) {
// Increment velocity and position, using a discrete linear approximation.
// p( t + Dt ) = p ( t ) + Dt * Dp/Dt
// v( t + Dt ) = v ( t ) + Dt * Dv/Dt
// Insure temporal symmetry by treating the first and second half of the
// temporal increment equivalently.
// v( t + Dt/2 ) = v ( t ) + Dt/2 * Dv/Dt
// p( t + Dt ) = p ( t ) + Dt * v( t + Dt/2 )
// v( t + Dt ) = v ( t ) + Dt * Dv/Dt = v( t + Dt/2 ) + Dt/2 * Dv/Dt
// update velocity based on acceleration, first half increment
PVector dg = PVector.mult( g, .5*tFraction );
velocity.add( dg );
// update position based on velocity
PVector dv = PVector.mult( velocity, tFraction );
position.add( dv );
// update velocity based on acceleration, second half increment
dg = PVector.mult( g, .5*tFraction );
velocity.add( dg );
}
void display( boolean bTBoundary ) {
stroke(0); strokeWeight(1);
if ( bTBoundary ) {
fill( clrFill );
}
else {
// noFill();
fill( clrFill, 92 );
}
ellipse( position.x, -position.y, 2*r-3, 2*r-3 );
}
void UpdateCollisionWithLine( Line line, boolean bShowOrnaments ) {
phi = atan2( velocity.y, velocity.x );
// println( " " );
// println( "line angle: " + degrees(line.theta) );
// println( "ball angle in: " + degrees(b.phi) );
// Collision and rebound of ball with line.
// Update velocity change due to collision (but don't update
// position or velocity change due to gravity).
float vel_in_mag = velocity.mag();
velocity.x = vel_in_mag * cos( 2*line.theta - phi );
velocity.y = vel_in_mag * sin( 2*line.theta - phi );
if ( bShowOrnaments ) {
// Draw ball's Vout vector.
stroke( 255, 0, 0 );
line( position.x, -position.y, position.x + 10*velocity.x, -position.y - 10*velocity.y );
}
phi = atan2( velocity.y, velocity.x );
// println( "ball angle out: " + degrees(phi) );
// println( "ball speed in: " + vel_in_mag );
// println( "ball speed out: " + velocity.mag() );
}
float TimeToCollisionWithLine( Line line, float tFraction_est, boolean bShowOrnaments ) {
// Velocity at half estimated time increment.
PVector dg = PVector.mult( g, 0.5*tFraction_est ); //<>// //<>//
PVector velocity_next = new PVector( velocity.x, velocity.y );
velocity_next.add( dg );
float phi_in = atan2( velocity_next.y, velocity_next.x );
// // Side of line, using y-value:
// // a*x + b*y + c = 0;
// // y = ( -c - ax )/b
// // Note: Fails with b == 0, a vertical line.
// float signAbove = -1;
// float yl_bx = ( -line.C - line.A*b.position.x ) / line.B;
// if ( b.position.y > yl_bx ) {
// signAbove = 1;
// }
// Perpendicular distance change numerator:
float pv = line.A*position.x + line.B*position.y + line.C;
float p_dir = pv/abs(pv); // above or below line?
// float dpv = line.A*velocity_next.x + line.B*velocity_next.y + line.C;
// float dp_dir = dpv/abs(dpv);
// Find contact point on circle. This is the point where the tangent matches the
// line orientation, on the side of the circle nearest to the line.
float xposC_b = position.x + r*cos( line.theta + p_dir*PI/2 );
float yposC_b = position.y + r*sin( line.theta + p_dir*PI/2 );
// Perpendicular distance between contact point on circle and line,
// using line standard coefficients:
// distance( ax + by + c = 0, (x0,y0) ) = | ax0 + by0 + c | / sqrt( a^2 + b^2 )
float perp_dist = abs( line.A*xposC_b + line.B*yposC_b + line.C )
/ sqrt( line.A*line.A + line.B*line.B );
float a_incidence = PI/2 + phi_in - line.theta;
if ( p_dir > 0 ) a_incidence += PI;
if ( a_incidence > PI ) a_incidence -= 2*PI;
if ( a_incidence < -PI ) a_incidence += 2*PI;
// Point on line closest to contact point on circle (and position, center, of circle):
// x = ( b*( b*x0 - a*y0 ) - a*c )/ ( a^2 + b^2 )
// y = ( a*( -b*x0 + a*y0 ) + b*c )/ ( a^2 + b^2 )
float xcl = ( line.B*( line.B*position.x - line.A*position.y ) - line.A*line.C )
/ ( line.A*line.A + line.B*line.B );
float ycl = ( line.A*( -line.B*position.x + line.A*position.y ) - line.B*line.C )
/ ( line.A*line.A + line.B*line.B );
// If showing ornaments and near to line:
if ( bShowOrnaments && perp_dist <= 3*r ) {
line( xposC_b, -yposC_b, xcl, -ycl );
println( " " );
println( "perpendicular distance: " + perp_dist );
println( "angle of incidence: " + degrees( a_incidence ) );
}
// println( "perpendicular distance: " + perp_dist );
// If circle is moving toward the line:
if ( a_incidence > -PI/2
&& a_incidence < PI/2 ) {
float tIntersect = perp_dist / ( velocity_next.mag()*cos(a_incidence) );
//// Position at intersection with projection of line.
//float xposIntersect = position.x + velocity_next.x*tIntersect;
//float yposIntersect = position.y + velocity_next.y*tIntersect;
// If closest (contact) point is between the line endpoints:
if ( xcl >= line.end1.x
&& xcl <= line.end2.x
&& ycl >= min( line.end2.y, line.end1.y )
&& ycl <= max( line.end2.y, line.end1.y ) ) {
// To Do: Order these line endpoints on creation?
// time_to_collision = distance_to_collision / speed_of_approach
// println( " " );
// println( "time to collision with a line: " + tIntersect );
return tIntersect;
}
else {
return -1.0;
}
}
return -1.0; //<>//
}
}
class Line {
PVector end1;
PVector end2;
float A, B, C;
float theta;
Line( PVector end1_, PVector end2_ ) {
if ( end1_.x < end2_.x ) {
end1 = end1_;
end2 = end2_;
}
else {
end1 = end2_;
end2 = end1_;
}
// Standard form of line coefficients:
// Ax + By + C = 0
if ( end2.x - end1.x != 0 ) {
A = ( end2.y - end1.y ) / ( end2.x - end1.x ); // slope
B = -1;
C = end1.y - (end1.x * A); // intercept
}
else {
A = 1;
B = 0;
C = -end1.x;
}
theta = atan2( end2.y-end1.y, end2.x-end1.x );
}
void display() {
strokeWeight(2);
stroke( 0, 0, 0 );
line( end1.x, -end1.y, end2.x, -end2.y );
}
}
No comments:
Post a Comment