A falling Slinky in vpython

falling-slinky-vpython.py

from vpython import *
import numpy as np

# Parameters
n = 200  # Increased number of segments for a longer slinky (smoother turns)
radius = 5.0  # Radius of the helical path (width of the slinky)

dt = 0.005  # Time step in seconds
endtime = 1  # Total simulation time in seconds
iters = round(endtime / dt) # Number of iterations to run the simulation, rounded to the nearest integer

k = 1 * n # Spring constant (stiffness)
m = 1 / n # Mass of each segment
g = 9.8 # acceleration due to gravity

# Construct the tridiagonal matrix A
A = np.diag(np.full(n, -1), -1) + np.diag(np.full(n + 1, 2), 0) + np.diag(np.full(n, -1), 1)
A[0, 0] = 1 # Fix the first element of the diagonal
A[-1, -1] = 1 # Fix the last element of the diagonal
e = np.ones(n + 1) # Vector of ones

# Initial conditions
j = np.arange(n + 1) # Array of integers from 0 to n
y = -(m * g / k) * (n * j - 0.5 * j * (j - 1)) # start from center of the canvas
# y = -(m * g / k) * (n * j - 0.5 * j * (j - 1)) + 10  # start from the top, add 10 to all y values

ydot = np.zeros(n + 1) # Initial velocity is zero
ydot2 = -(k / m) * A @ y - g * e # Second derivative of y with respect to time (acceleration)

# Initialize arrays to store the data for plotting
positions = np.zeros((n + 1, iters)) # Array to store the positions of each segment at each time step

# Simulation loop
for iter in range(iters):
    ydot += ydot2 * dt / 2 # Update the velocity
    y += ydot * dt / 2 # Update the position

    for i in range(n):
        for j in range(i + 1, n + 1):
            if y[j] < y[i]:
                break
        j -= 1
        ybar = np.mean(y[i:j + 1]) # Average height of the segments
        ydotbar = np.mean(ydot[i:j + 1]) # Average velocity of the segments
        y[i:j + 1] = ybar # Set the height of the segments to the average
        ydot[i:j + 1] = ydotbar # Set the velocity of the segments to the average

    y += ydot * dt / 2 # Update the position again
    ydot2 = -(k / m) * A @ y - g * e # Calculate the acceleration at the new positions
    ydot += ydot2 * dt / 2 # Update the velocity

    positions[:, iter] = y # Store the positions of the segments at this time step

# VPython setup
scene = canvas(title='Falling Slinky Animation', width=800, height=600, center=vector(0, 0, 0), background=color.white)
floor = box(pos=vector(0, -10, 0), size=vector(20, 0.1, 20), color=color.green)  # Create a floor
backwall = box(pos=vector(0, 0, -10), size=vector(20, 20, 0.1), color=color.orange)  # Create a back wall
# # mark the origin for easy reference
# curve(pos=[vector(-1, 0, -1), vector(1, 0, 1)], color=color.red)
# curve(pos=[vector(-1, 0, 1), vector(1, 0, -1)], color=color.red)
scene.camera.pos = vector(0, 0, 20) 
scene.camera.axis = vector(0, 0, -20)
scene.camera.fov = pi/3

# Create the helix segments
theta = np.linspace(0, 8 * np.pi, n + 1)  # Adjusted for more turns

# Animation loop
slinky = None
while True:
    for t in range(iters):
        rate(15)  # Control the speed of the animation
        if slinky:  # Make the previous slinky invisible
            slinky.visible = False
        curve_points = [vector(radius * np.cos(theta[i]), positions[i, t], radius * np.sin(theta[i])) for i in range(n + 1)] # Update the positions of the segments
        slinky = curve(pos=curve_points, color=color.blue)