#!/usr/bin/env pythonw
"""
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 <dbailey@physics.utoronto.ca>
                            (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.
    """

from __future__ import print_function   # Use Python 3 consistent print

# 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          = 50       # 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 = "horizontal line"
chain_shape = "hanging trefoil"
chain_shape = "horizontal trefoil"
#   Current choices are: "hanging trefoil", "horizontal trefoil",
#                       "horizontal line"

# Vibrating circular plate parameters
plate_diameter  = 0.30
plate_thickness = 0.002
amplitude       = 0.006      # Amplitude of vertical oscillations
frequency       = 15.      # Frequency of vertical oscillations in Hz

# Inelasticities (fractional speed change per bounce) of beads
plate_dv = 0.30 # Bead collisions with plate  
bead_dv  = 0.80 # 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("brg")
cmap=cm.get_cmap("gist_rainbow")

# Import vPython modules
from visual     import color, cylinder, distant_light, norm,\
                       rate, 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, orientation):
    """
    Define positions as either hanging vertical or horizontal laid out
        trefoil knot starting at pos
       If vertical, currently the knot hangs 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]
    if orientation == "horizontal" :
        spaces = vector(-spacing,0,0)
    else :
        # Default is vertical hanging knot
        spaces = vector(0,0,-spacing)
    for i in range(1,n_bead):
        positions[i] = positions[i-1] + spaces
    # 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 + 0.001

    start_bead = straight
    current_bead = start_bead + 1
    print("Current bead = ", current_bead)
    for i in range(0,n):
        t = i/float(n)
        if orientation == "horizontal" :
            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+r*cos(nu_1*2*pi*t))*cos(nu_2*2*pi*t)),
                positions[start_bead][2] + scale*   r*sin(nu_1*2*pi*t)
                )
        else :
            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] + spaces
       
    # Debugging section, to make sure bead spacings are okay
    print ("Bead diameter   = ", bead_diameter)
    print ("Rod length      = ", rod_length)
    print ("Spacing         = ", spacing)
    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
bead_spacing = bead_diameter+rod_length/2.0
if chain_shape == "horizontal line" :
    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,
            "vertical")
elif chain_shape == "horizontal trefoil":
    bead_pos = trefoil(
            vector((n_bead-min_loop_size)*bead_spacing/2, 0.0,
                    max_separation*((n_bead-2*min_loop_size)/2.+4)), 
            bead_diameter+rod_length/2.0,
            "horizontal")
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)<B(1)<A(2)<B(2).
   
   """

graph1 = gdisplay(xtitle='time (seconds)',ytitle='Location of crossings', ymax=n_bead,
                  x=0, y=scene.height)
crossings_plot = []
crossings_plot.append( gdots(gdisplay=graph1, color=color.red,     size=4) )
crossings_plot.append( gdots(gdisplay=graph1, color=color.green,   size=4) )
crossings_plot.append( gdots(gdisplay=graph1, color=color.blue,    size=4) )
crossings_plot.append( gdots(gdisplay=graph1, color=color.yellow,  size=4) )
crossings_plot.append( gdots(gdisplay=graph1, color=color.cyan,    size=4) )
crossings_plot.append( gdots(gdisplay=graph1, color=color.magenta, size=4) )
crossings_plot.append( gdots(gdisplay=graph1, color=color.white,   size=4) )
max_crossings = len(crossings_plot)-1

# Zero time that will track whether the chain is really unkotted.
unknotted_time = 0.0

# Calculated some useful constant values outside loop
plated = bead_radius + plate_thickness
edged  = plate.radius - bead_radius
edged2 = edged**2

# Create list all possible unique bead pairs
k_ij = []
k_plus1 = []
for i in range(n_bead-1) :
    k_plus1.append(int(i*(n_bead-0.5*(1+i))))
    for j in range(i+1,n_bead) :
        k_ij.append([i,j])

n = -1          # Initialize iteration step counter
t = 0.          # Initialize time
timer("Before Main Loop")

# Main simulation loop
while(True):
    timer("Start of Main Loop")

    # Check if time to exit loop
    if maximum_number_of_steps >= 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 :
                    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 j<n_bead-1 :
                        if intersect(x[i],x[i+1],x[j],x[j+1]) :
                            crossings.append([i+0.5,j+0.5])
                    if i<n_bead-1:
                        if intersect(x[i],x[i+1],x[j-1],x[j]) :
                            crossings.append([i+0.5,j-0.5])
                    if i>0:
                        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<n_bead-1 :
                        if intersect(x[i-1],x[i],x[j],x[j+1]) :
                            crossings.append([i-0.5,j+0.5])
            
        # Remove duplicates
        unique_crossings = []
        for crossing in crossings:
            if crossing not in unique_crossings:
                unique_crossings.append(crossing)
        return unique_crossings
    
    #   Only update scene every n_update_scene time steps
    if not n%n_update_scene :
        rate(30)
        chain_pos = vector(0,0,0)
        for i in range(n_bead) :
            bead[i].pos = bead_pos[i]
            chain_pos   = chain_pos + bead_pos[i]
        # Move plate
        plate.z = amplitude*sin(omega*t) - plate_thickness
        # View follows chain
        scene.center = chain_pos/n_bead
        # Track knot crossing points
        crossings = check_for_crossings(0)
        n_crossings = len(crossings)
        # Sort crossings, with entangled crossings first
        #   and other crossings last.
        entangled_crossings = []
        other_crossings = []
        i=0
        while i < n_crossings :
            try :
                # A set of 3 crossings is labelled entangled if
                #   they are consistent with being a trefoil knot lying
                #   flat without any extra bits of chain crossing the knot.
                if (    crossings[i  ][1] >= 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 chain_shape != "horizontal line" and 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 beaded_chain.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.
"""