""" beaded_chain.py: A vPython simulation of a beaded chain on a vibrating plate. This version (27 March 2012) vibrates a chain with a trefoil knot, plots the position of any points where the chain crosses itself, and stops when the chain unknots. Copyright (c) 2012 University of Toronto Last Modification: 27 March 2012 by David Bailey Original Version: 23 March 2012 by David Bailey Contact: David Bailey (http://www.physics.utoronto.ca/~dbailey) License: Released under the MIT License; the full terms are this license are appended to the end of this file, and are also available at http://www.opensource.org/licenses/mit-license.php. Physics notes: All units in SI (i.e. metres, seconds, Newtons, ...) Software notes: This software is primarily pedagogical, so readability is more important than maximum speed. The simulation runs about 20 times slower than real time for a 60 bead chain running on a late 2010 MacBook Pro. Suggestions on improving the speed are welcome. Numpy is used in some places, but the complex loops are not easy to convert. Trying to switch back and forth between python lists and numpy arrays is very slow, as are the sometimes necessary operations on individual numpy array elements. It may be possible to speed it up using C (e.g. Cython) or Fortran (e.g. f2py) modules, but that would make the code less portable and less readable for students who only know python. A minor warning: Using "+=" and "-=" does not always work properly with visual vectors, so avoided here. """ # Imports follow these parameter definitions # Time step (in simulated seconds) between main loop iterations. dt = 0.0001 # If this is too large, the chain can pass through itself; # if too small, the simulation will be very slow. # Maximum number of steps in main loop, use -1 for infinite loop. maximum_number_of_steps = -1 # How many steps between vPython scene updates n_update_scene = 20 # Bead parameters n_bead = 60 # Number of beads in chain bead_diameter = 0.00388 # Diameter of each bead in metres rod_length = 0.0009 # Effective length of connecting rod # Using a realistic rod length (e.g. 0.0014 for bead diameter=0.00388) # makes it likely that the chain will pass through itself, since this # simulation does not yet check for collisions between connecting rods. # Minimum number of beads that can form a simple closed loop. min_loop_size = 9 # Initial chain shape chain_shape = "hanging trefoil" # Current choices are: "hanging trefoil" or "horizontal line" # Vibrating circular plate parameters plate_diameter = 0.30 plate_thickness = 0.002 amplitude = 0.006 # Amplitude of vertical oscillations frequency = 13.37 # Frequency of vertical oscillations # Inelasticities (fractional speed change per bounce) of beads plate_dv = 0.50 # Bead collisions with plate bead_dv = 0.94 # Bead collisions with each other # Maximum size of small random directional perturbations on bounces from plate plate_dd = 0.2 # Include some matplotlib color modules # and select colormap for coloring beads while debugging from matplotlib import cm, colors cmap=cm.get_cmap("gist_rainbow") # Import vPython modules from visual import color, cylinder, distant_light, norm,\ scene, sphere, vector # import graphing features from visual.graph import gdisplay, gdots # The following numpy modules are included in visual from visual import dot, nonzero, pi, random # Special distance finding routine from scipy, for speed from scipy.spatial.distance import pdist # Use math library functions, since corresponding numpy functions # are currently slower since operations are currently mostly made on # elements, not whole arrays. from math import cos, sin, sqrt import os class TimeProfiler : # Simple profiling tool that measures time # Typical use: # timer = TimeProfiler() # ... some code ... # timer(1) # ... code to be timed ... # timer(2) # ... some more code ... # timer("moose") # ... some more code to be timed ... # timer("elk") # ... some more code ... # timer.printer() def __init__(self): import time self.last_time = time.time() self.times = {} self.keylist = [] def __call__(self,key) : # The key can be either an integer or a string, i.e. # anything that works for a python dictionary. import time t = time.time() try : self.times[key] += t-self.last_time except : # If key does not already exist, create it self.times[key] = t-self.last_time self.keylist.append(key) self.last_time=t def printer(self) : # Output time profiling results for key in self.keylist : print( "{0:20s}: {1:8.3f} seconds, {2:6.1f}%". format( str(key), self.times[key], 100*self.times[key]/sum(self.times.values()) ) ) print( "{0:20s}: {1:8.3f} seconds, {2:6.1f}%". format( "Total Time", sum(self.times.values()),100 ) ) timer = TimeProfiler() timer("Start of Program") # Create scene scene.title = 'Bouncing Beaded Chain' scene.width = 600 scene.height = 600 scene.autoscale = False # No automatic scaling in scene # Set camera direction: (0,0,-1) is looking straight down, # (0,1,-0.03) looks sideways at chain scene.forward = (0,0,-1) scene.background = (0.95,0.95,0.95) # background color scene.lights = [ distant_light(direction=(0, 0, 1), # top lighting color=color.gray(1)), distant_light(direction=(0, -1, 0), # side lighting color=color.gray(0.55))] # Create visual (vPython) plate object omega = 2*pi*frequency # Angular Frequency of plate plate_axis = (0, 0, plate_thickness) plate = cylinder( pos = (0,0,-plate_thickness), axis=(0,0,plate_thickness), radius=0.5*plate_diameter, color=(0,0,0.1) ) # Set area viewed scene.range = 1.0*vector([plate.radius,plate.radius,plate.radius]) # Create visual (vPython) chain bead_radius = bead_diameter/2.0 # Maximum distance between centres of connected beads max_separation = bead_diameter+rod_length # Minimum angle formed by 3 consecutive beads in a chain, # determined from the minimum loop size. min_cos = cos(pi*(1.0-2.0/min_loop_size)) def line(pos, spacing, direction): # Define positions equally spaced along a line starting at pos # along direction vector positions = n_bead*[pos] delta_position = spacing*norm(direction) for i in range(1,n_bead): positions[i] = positions[i-1] + delta_position return positions def trefoil(pos, spacing): # Define positions as hanging trefoil knot starting at pos # Currently the knot hangs vertically, aligned along the y axis; # the size of the trefoil varies by 1 bead, depending # on whether the chain has an even or odd number of beads. # The sign of "spacing" defines whether the ends or the knot are at bottom # The parametric equation is from http://arxiv.org/pdf/physics/0103039.pdf # "The ideal trefoil knot", P. Pieranski and S. Przybyl. knot_size = 2*min_loop_size # Length of straight hanging sections of chain straight = (n_bead-knot_size)/2+1 # Initialize bead positions positions = n_bead*[pos] for i in range(1,n_bead): positions[i] = positions[i-1] + vector(0,0,-spacing) # Catch impossible knots if n_bead < knot_size : print "\n**** Chain is too short to tie trefoil knot! ****\n" return positions # Parameters for trefoil knot parametric equation nu_1 = 3 nu_2 = 2 r = 1.0 R = 2.0 n = 1000 # Measure the length of trefoil knot described by parameterization length = 0.0 previous_position = vector(0,0,0) for i in range(n) : t = i /float(n) current_position = vector( (R+r*cos(nu_1*2*pi*t))*sin(nu_2*2*pi*t), r*sin(nu_1*2*pi*t), (R+r*cos(nu_1*2*pi*t))*cos(nu_2*2*pi*t) ) length += (current_position-previous_position).mag previous_position = current_position # Calculate scale factor to map parameterization onto real chain length scale = knot_size*bead_diameter/length start_bead = straight current_bead = start_bead + 1 for i in range(1,n): t = i/float(n) positions[current_bead] = vector( positions[start_bead][0] + scale*(R+r*cos(nu_1*2*pi*t))*sin(nu_2*2*pi*t), positions[start_bead][1] + scale* r*sin(nu_1*2*pi*t), positions[start_bead][2] + scale*((R+r*cos(nu_1*2*pi*t))*cos(nu_2*2*pi*t)-R-r) ) # Check that the current bead does not overlap with any previous beads done = True for j in range(current_bead) : if (positions[current_bead]-positions[j]).mag < bead_diameter : done = False break if done : current_bead += 1 if current_bead >= n_bead : break for i in range(current_bead,n_bead): positions[i] = positions[i-1] - vector(0,0,-spacing) # Debugging section, to make sure bead spacings are okay for i in range(0,n_bead-1): if (positions[i+1]-positions[i]).mag > bead_diameter+rod_length : print i, "Too Far" print( i, (positions[i+1]-positions[i]).mag/spacing, (positions[i+1]-positions[i]).mag/bead_diameter, (positions[i+1]-positions[i]).mag/(bead_diameter+rod_length) ) print positions[i+1]-positions[i] print positions[i+1],positions[i] for j in range(i+1,n_bead): if (positions[j]-positions[i]).mag < bead_diameter : print i,j, "Too Close", n_bead, straight print((positions[j]-positions[i]).mag/spacing, (positions[j]-positions[i]).mag/bead_diameter, (positions[j]-positions[i]).mag/(bead_diameter+rod_length) ) print positions[j]-positions[i] return positions # Create and initialize chain if chain_shape == "horizontal line" : bead_spacing = bead_diameter+rod_length/2.0 bead_0_position = vector(-n_bead*bead_spacing/2,0,0.005) bead_pos = line( bead_0_position, bead_spacing, vector(1,0,0)) elif chain_shape == "hanging trefoil": bead_pos = trefoil( vector(-0.03,0.0,max_separation*((n_bead-2*min_loop_size)/2.+4)), bead_diameter+rod_length/2.0 ) else : print "Incorrect configuration specified." print print "Please check your code and try again" # Exit gracefully os._exit(1) # Create visual (vPython) bead objects bead = n_bead*[None] bead_vel = n_bead*[vector(0,0,0)] for i in range(n_bead) : bead[i] = sphere(pos = bead_pos[i], radius = bead_radius, color = colors.colorConverter.to_rgb(cmap(i/float(n_bead)))) """ Define graph for crossing positions A crossing is defined as any point where the chain crosses itself. A crossing, C, position is defined by the two beads, C(1) & C(2) where the crossing occurs, where beads are identified by their index counted from one end of the chain. The first 3 crossings curves (red, green, blue) are for entangled crossings. Two crossings, A & B, are described as entangled if A(1)= 0 and n > maximum_number_of_steps: break n += 1 #track iteration steps t += dt #track time # Initialize bead acceleration to normal gravity bead_acc = n_bead*[vector(0,0,-9.81)] def bounce_off_plate() : # Useful constant platedz = plated + amplitude*sin(omega*t) - plate_thickness plate_velz = amplitude*omega*cos(omega*t) # Bounce bead off plate if it is in contact with it. for i in range(n_bead) : # Bounce bead if it is on plate, otherwise let it fall if bead_pos[i].z < platedz : # Move vertical position to above plate bead_pos[i].z = 2*platedz - bead_pos[i].z # Calculate bead velocity in plate reference frame, with energy loss bead_vel[i].z = plate_dv*(bead_vel[i].z - plate_velz) # If bead falling, bounce up (with a bit of randomness) if bead_vel[i].z <= 0 : vz1 = bead_vel[i].z bead_vel[i] = bead_vel[i] + (bead_vel[i].z * 2 * norm( vector( plate_dd*(random.random()-0.5), plate_dd*(random.random()-0.5) , -1) ) ) # Revert to velocity in lab reference frame bead_vel[i].z = bead_vel[i].z + plate_velz # Bounce off walls at edge of plate if (bead_pos[i].x**2 + bead_pos[i].y**2) > edged2 : bead_pos[i] = 2*( norm((bead_pos[i].x,bead_pos[i].y,0)) * edged) +\ (-bead_pos[i].x, -bead_pos[i].y, bead_pos[i].z) if dot(norm((bead_pos[i].x,bead_pos[i].y,0)),bead_vel[i]) > 0 : bead_vel[i] = bead_vel[i] - ( (1+plate_dv)* dot(norm((bead_pos[i].x,bead_pos[i].y,0)), bead_vel[i]) * norm((bead_pos[i].x,bead_pos[i].y,0)) ) bounce_off_plate() timer("Bounced Off Plate") # Calculate distances between all pairs of beads # This would be about a factor of 10 faster if bead_pos was # already a numpy array, but other complex loops would be much # slower, more than offsetting the time saving here. bead_separation = pdist(bead_pos) timer("After pdist") k_bead=10 def bounce_beads(i, j, direction) : # Bounces two beads off each other; direction of # bounce can be either away (direction = +1) or # towards (-1), since bounce can occur either because # beads have touched or because they have reached # the end of their connecting chain. # Since all beads have same mass, just use velocity # conservation instead of momentum conservation mean_velocity = 0.5*(bead_vel[i] + bead_vel[j]) # transform velocities to cm frame of two beads, vcm_i = bead_vel[i] - mean_velocity vcm_j = bead_vel[j] - mean_velocity # Direction vector from centre of bead i to bead j n_ij = norm(bead_pos[j] - bead_pos[i]) # Change velocities if beads moving in wrong direction if dot(n_ij,vcm_i)*direction > 0 : # Bounce in cm frame, with inelasticity vcm_i = bead_dv*(vcm_i - 2*dot(vcm_i,n_ij)*n_ij) vcm_j = bead_dv*(vcm_j - 2*dot(vcm_j,n_ij)*n_ij) # transform velocities back to lab frame bead_vel[i] = vcm_i + mean_velocity bead_vel[j] = vcm_j + mean_velocity # Add a weak acceleration to the bounce # This is for cases where two beads overlap (or are too far) # apart, but which have zero relative velocity, in which # case reversing the velocity does nothing. bead_acc[i] = bead_acc[i]-0.5*direction*k_bead*n_ij bead_acc[j] = bead_acc[j]+0.5*direction*k_bead*n_ij # Bounce if beads too close # This takes advantage of numpy to do the if statement all at once # to save huge amount of time. for k in nonzero(bead_separation < bead_diameter)[0] : bounce_beads(k_ij[k][0],k_ij[k][1],1) # Bounce next bead if too far for k in k_plus1 : if bead_separation[k] > max_separation : bounce_beads(k_ij[k][0],k_ij[k][1],-1) timer("After Bounce Beads") def check_angle_between_beads() : k_bend = 10 b2 = (bead_pos[1]-bead_pos[0]).mag2 for i in range(n_bead-2) : a2 = b2 b2 = (bead_pos[i+2]-bead_pos[i+1]).mag2 c2 = (bead_pos[i+2]-bead_pos[i ]).mag2 cos_angle = (a2+b2-c2)/(2*sqrt(a2*b2)) if cos_angle > min_cos: # Get mean position and velocity of outer beads outer_pos = 0.5*(bead_pos[i] + bead_pos[i+2]) outer_vel = 0.5*(bead_vel[i] + bead_vel[i+2]) # Follow method of bounce, treating the outer beads as a double # mass bead # Mean velocity of all 3 beads mean_velocity = (bead_vel[i]+bead_vel[i+1]+bead_vel[i+2])/3.0 vcm_i = bead_vel[i+1] - mean_velocity vcm_o = outer_vel - mean_velocity n_io = norm(outer_pos - bead_pos[i+1]) # Change velocities if beads moving in wrong direction if dot(n_io,vcm_i) < 0 : # Bounce in cm frame, with inelasticity vcm_i = bead_dv*(vcm_i - 2*dot(vcm_i,n_io)*n_io) vcm_o = bead_dv*(vcm_o - 2*dot(vcm_o,n_io)*n_io) # transform velocities back to lab frame bead_vel[i+1] = vcm_i + mean_velocity vcm_io = bead_vel[i ] - outer_vel vcm_jo = bead_vel[i+2] - outer_vel outer_vel = vcm_o + mean_velocity bead_vel[i ] = vcm_io + outer_vel bead_vel[i+2] = vcm_jo + outer_vel # Add a weak acceleration to the bounce # Acceleration vector from bead i+1 to midpoint between # beads i and i+2 acceleration = k_bend*(min_cos-cos_angle)*\ norm(0.5*(bead_pos[i+2] + bead_pos[i]) - bead_pos[i+1]) bead_acc[i ] = bead_acc[i ] - 0.5*acceleration bead_acc[i+1] = bead_acc[i+1] + acceleration bead_acc[i+2] = bead_acc[i+2] - 0.5*acceleration check_angle_between_beads() timer("Angle test") # Increment velocity, position for i in range(n_bead) : bead_vel[i] = bead_vel[i] + bead_acc[i]*dt # Move bead to new position based on current velocity and time step bead_pos[i] = bead_pos[i] + bead_vel[i]*dt timer("Increment bead x, v") def check_for_crossings(axis) : # A chain crossing point is defined as # One bead above another bead, as long as # - there is not another crossing point within the # minimum loop size along the chain def ccw(a,b,c) : # Determine if three 2-Dimensional points are in counterclockwise order # # Based on method found in Jeff Erickson's notes at # http://compgeom.cs.uiuc.edu/~jeffe/teaching/compgeom/notes/01-convexhull.pdf return (c[1]-a[1])*(b[0]-a[0]) > (b[1]-a[1])*(c[0]-a[0]) def intersect(a,b,c,d) : # Determine if two 2-Dimensional line segments, (a,b) & (c,d), intersect # Based on method found in Jeff Erickson's notes at # http://compgeom.cs.uiuc.edu/~jeffe/teaching/373/notes/x06-sweepline.pdf return ccw(a,c,d) != ccw(b,c,d) and ccw(a,b,c) != ccw(a,b,d) crossings = [] # Calculate distance between beads in chosen plane. # The default plane is xy, i.e. looking from top. # Planes 1 and 2 are rotated a bit around x and y axes, # so that obscured beads can be seen x = [] if axis == 0 : for i in range(n_bead) : x.append([bead_pos[i][0], bead_pos[i][1]]) elif axis == 1 : for i in range(n_bead) : x.append( [ bead_pos[i][0] - 0.17*bead_pos[i][2], bead_pos[i][1] ] ) elif axis == 2 : for i in range(n_bead) : x.append( [ bead_pos[i][0], bead_pos[i][1] - 0.17*bead_pos[i][2] ] ) elif axis == 3 : for i in range(n_bead) : x.append( [ 0.707106781*(bead_pos[i][0]+bead_pos[i][1]), 0.707106781*(bead_pos[i][0]-bead_pos[i][1]) - 0.17*bead_pos[i][2] ] ) elif axis == 4 : for i in range(n_bead) : x.append( [ 0.707106781*(bead_pos[i][0]+bead_pos[i][1]) - 0.17*bead_pos[i][2], 0.707106781*(bead_pos[i][0]-bead_pos[i][1]) ]) else : print "Invalid check_for_crossings call" os._exit(1) # Exit gracefully. bead_xy_separation = pdist(x) # Look for beads that are close together # Take advantage of numpy to do the if statement. too_close = bead_xy_separation < bead_diameter min_sep = min_loop_size-1 for i in range(n_bead-min_sep) : for j in range(i+min_sep,n_bead) : k = (i*(2*n_bead-i-1))/2 + j-i -1 # Identify non-neighbour bead pairs that are close together. if too_close[k] : # Any crossing could be either side of the selected beads, # so four line segments must be checked. if j0: if intersect(x[i-1],x[i],x[j-1],x[j]) : crossings.append([i-0.5,j-0.5]) if i>0 and j= crossings[i+1][0] and crossings[i+1][1] >= crossings[i+2][0] and crossings[i ][1] <= crossings[i+1][1] and crossings[i+1][1] <= crossings[i+2][1] ) : entangled_crossings.append(crossings[i]) entangled_crossings.append(crossings[i+1]) entangled_crossings.append(crossings[i+2]) i=i+3 else : other_crossings.append(crossings[i]) i=i+1 except : other_crossings.append(crossings[i]) i=i+1 # Put crossings back together crossings = entangled_crossings + other_crossings # Plot crossing positions for i in range(n_crossings) : if i < max_crossings : crossings_plot[i].plot(pos=(t,crossings[i][0])) crossings_plot[i].plot(pos=(t,crossings[i][1])) else : crossings_plot[max_crossings].plot( pos=(t,crossings[i][0])) crossings_plot[max_crossings].plot( pos=(t,crossings[i][1])) if n_crossings < 3: unknotted_time = unknotted_time + n_update_scene*dt # Confirm that this is not a temporary fluctuation if unknotted_time > 0.05 : print "Unknotted at ", t, " seconds, after ",n," iterations" print "Crossings in xy plane: ", crossings # Check_for_crossings(True,0) print "Number of observed crossings: ", \ len(check_for_crossings(0)),\ len(check_for_crossings(1)),\ len(check_for_crossings(2)) break else: unknotted_time = 0.0 timer("End of Loop") print "\nTermination after ", n, " loops; simulated time: ",t," seconds" timer("End of Program") timer.printer() ##### End extended_curve_fit_to_data.py """ Full text of MIT License: Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """