update naming

25 files changed, 0 insertions, 1270 deletions

diff --git a/t/old/1d.jl b/t/old/1d.jl
deleted file mode 100644
index 203466c..0000000
--- a/t/old/1d.jl
+++ /dev/null
@@ -1,109 +0,0 @@
-using Plots
-
-N = 10  # number of masses in the 1D lattice
-K = 100   # elastic force constant
-m = 1   # mass of particles
-M = 1000 # big mass (one only)
-t_final = 2 # seconds
-dt = 0.001 # timestep
-
-initial_distance_between_masses = 10 # meters
-L = initial_distance_between_masses * N # length of the 1D lattice
-
-# 2D array of current positions and velocities (Hamiltonian state)
-q = zeros(N)
-p = zeros(N)
-
-# initialize the positions and velocities
-for i in 1:N
-	q[i] = initial_distance_between_masses * (i - 1)
-	p[i] = 0
-end
-
-
-# plot positions
-function plot_positions(q, t)
-	# only 1D, set y to 0
-	y = zeros(N)
-
-	# plot the masses
-	plot(
-		q, y,
-		seriestype = :scatter, label = "Masses @ t = $t",
-		# xlims = (-1, N + 1), 
-		ylims = (-1, 1),
-		markerstrokewidth = 0,
-	)
-
-	# plot the middle one as red
-	plot!(
-		[q[Int(N / 2)]], [0],
-		seriestype = :scatter, label = "Large mass @ t = $t",
-		color = :red,
-		markerstrokewidth = 0, markersize = 10,
-	)
-
-	return plot!()
-end
-
-# update the state
-function update_state!(q, p, dt)
-	new_q = copy(q)
-	new_p = copy(p)
-
-	# update the small masses state
-	for i in 2:N-1
-		dx_right = q[i+1] - q[i]
-		dx_left = q[i] - q[i-1]
-
-		new_q[i] += dt * p[i] / m
-		new_p[i] += dt * (K * dx_right - K * dx_left)
-	end
-
-	# handle the ends, since our 1D system is cyclic
-	# case where i = 1, first particle
-	dx_right = q[2] - q[1]
-	distance_from_L = L - q[N]
-	dy_left = q[1] - distance_from_L
-	new_q[1] += dt * p[1] / m
-	new_p[1] += dt * (K * dx_right - K * dy_left)
-
-	# case where i = N, last particle
-	distance_from_0 = q[1]
-	dx_right = L + q[1] - q[N]
-	dx_left = q[N] - q[N-1]
-	new_q[N] += dt * p[N] / m
-	new_p[N] += dt * (K * dx_right - K * dx_left)
-
-
-	# update the large mass in middle (different mass, difference case)
-	middle_index = Int(N / 2)
-	dx_right = q[middle_index+1] - q[middle_index]
-	dx_left = q[middle_index] - q[middle_index-1]
-	new_q[middle_index] += dt * p[middle_index] / M
-	new_p[Int(N / 2)] += dt * (K * dx_right - K * dx_left)
-
-	# update the state
-	for i in 1:N
-		q[i] = new_q[i]
-		p[i] = new_p[i]
-	end
-end
-
-display(plot_positions(q, 0))
-
-function progress_system(q, p, dt, t_final)
-	t = 0
-	while t < t_final
-		update_state!(q, p, dt)
-		t += dt
-	end
-end
-
-progress_system(q, p, dt, t_final)
-
-println("Final state:")
-println(q)
-println(p)
-
-display(plot_positions(q, t_final))
diff --git a/t/old/animate-positions-.1stringa.mp4 b/t/old/animate-positions-.1stringa.mp4
deleted file mode 100644
index 0964313..0000000
--- a/t/old/animate-positions-.1stringa.mp4
+++ /dev/null
diff --git a/t/old/animate-positions-.5vel.mp4 b/t/old/animate-positions-.5vel.mp4
deleted file mode 100644
index b56e372..0000000
--- a/t/old/animate-positions-.5vel.mp4
+++ /dev/null
diff --git a/t/old/animate-positions-10.mp4 b/t/old/animate-positions-10.mp4
deleted file mode 100644
index aa00547..0000000
--- a/t/old/animate-positions-10.mp4
+++ /dev/null
diff --git a/t/old/animate-positions-80k.mp4 b/t/old/animate-positions-80k.mp4
deleted file mode 100644
index 3551bf4..0000000
--- a/t/old/animate-positions-80k.mp4
+++ /dev/null
diff --git a/t/old/animate-positions-chaos.mp4 b/t/old/animate-positions-chaos.mp4
deleted file mode 100644
index 7edbc75..0000000
--- a/t/old/animate-positions-chaos.mp4
+++ /dev/null
diff --git a/t/old/animate-positions-stable.mp4 b/t/old/animate-positions-stable.mp4
deleted file mode 100644
index 0b858a2..0000000
--- a/t/old/animate-positions-stable.mp4
+++ /dev/null
diff --git a/t/old/animate-positions-t.mp4 b/t/old/animate-positions-t.mp4
deleted file mode 100644
index a4a34e0..0000000
--- a/t/old/animate-positions-t.mp4
+++ /dev/null
diff --git a/t/old/disp.jl b/t/old/disp.jl
deleted file mode 100644
index 18e0a1b..0000000
--- a/t/old/disp.jl
+++ /dev/null
@@ -1,272 +0,0 @@
-M = 10
-m = 1 # mass of particles
-
-function calculate_force(
-	left_pos,
-	middle_pos,
-	right_pos,
-	K,
-	alpha = 0.0,
-	beta = 1000.0,
-)
-	linear_force = K * (left_pos + right_pos - 2 * middle_pos)
-	quadratic_force = alpha * (left_pos - middle_pos)^2 + alpha * (right_pos - middle_pos)^2
-	cubic_force = beta * (left_pos - middle_pos)^3 + beta * (right_pos - middle_pos)^3
-
-	return linear_force + quadratic_force + cubic_force
-end
-
-function tendency!(du, u, p, t)
-	# unpack the params
-	N, K, m = p
-
-	# get the positions and momenta
-	qs = u[1:2:end]
-	ps = u[2:2:end]
-
-	# go over the points in the lattice and update the state
-	for i in 1:N-1
-		mass = m
-		if i == Int(N / 2)
-			mass = M
-		end
-
-		left_index = max(1, i - 1)
-		right_index = min(N, i + 1)
-
-		du[i*2-1] = ps[i] / mass
-		force = calculate_force(qs[left_index], qs[i], qs[right_index], K)
-		du[i*2] = force / mass
-	end
-
-	# make last point same as first
-	du[N*2-1] = du[1] = 0 # set to 0
-	du[N*2] = du[2] = 0
-
-
-	if t % 100000 == 0
-		println("TIME UPDATE: ", t)
-	end
-
-	# if ps[Int(N / 2)] / M >= 1
-	# 	println("(in sim!) Time: ", t, " Vel: ", ps[Int(N / 2)] / M)
-	# 	# println("Other Positions: ", qs)
-	# 	println("Other Velocities: ", ps, "\n")
-	# end
-end
-
-function get_initial_state(
-	N,
-	initial_displacement = 2,
-	initial_velocity = 0,
-)
-	state = zeros(2 * N)
-
-	middle_index = 2 * Int(N / 2) - 1 # middle mass
-	state[middle_index] = initial_displacement
-	state[middle_index+1] = initial_velocity * M
-	return state
-end
-
-using DifferentialEquations
-function run_simulation(
-	N,
-	K,
-	m,
-	final_time,
-	initial_displacement = 2,
-	initial_velocity = 0,
-)
-	println("Running simulation with N = $N, K = $K, m = $m, final_time = $final_time, initial_displacement = $initial_displacement, initial_velocity = $initial_velocity\n")
-	s_0 = get_initial_state(N, initial_displacement, initial_velocity)
-
-	calculate_energy(s_0)
-
-	# pack the params
-	p = N, K, m
-	t_span = (0.0, final_time)
-	prob = ODEProblem(tendency!, s_0, t_span, p)
-	sol = solve(prob, Tsit5(), reltol = 1e-5, abstol = 1e-5, maxiters = 1e10) # control simulation
-
-	calculate_energy(sol.u[end])
-
-	println("Done Running Sim!\n\n")
-	return sol
-end
-
-using Plots
-function animate_positions(
-	states,
-	time_steps,
-	time_min = 0,
-	time_max = 30,
-	red_threshold = 2,
-	shift = true,
-)
-	println("Animating positions")
-	anim = @animate for i in 1:length(time_steps)
-		t = time_steps[i]
-		if t < time_min
-			continue
-		end
-		if t > time_max
-			break
-		end
-		positions = states[i][1:2:end]
-		v_middle = states[i][Int(length(states[1]) / 2)] / M
-		p_middle = states[i][Int(length(states[1]) / 2)-1]
-		y_lims = shift ? (-3 + p_middle, 3 + p_middle) : (-3, 3)
-		# plot(positions, label = "t = $(round(t, digits = 3)), v_middle=$(round(v_middle, digits=3))", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", ylim = (-3, 3))
-		if v_middle >= red_threshold
-			plot(positions, label = "t = $(round(t, digits = 7)), v_middle=$(round(v_middle, digits=7))", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", ylim = y_lims,
-				color = :red, legend = :topright,
-			)
-		else
-			plot(positions, label = "t = $(round(t, digits = 7)), v_middle=$(round(v_middle, digits=7))", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", ylim = y_lims,
-				color = :blue, legend = :topright,
-			)
-		end
-	end
-	mp4(anim, "t/animate-positions.mp4", fps = 30)
-	println("Done animating positions")
-end
-
-function plot_starting_and_final_positions(
-	states,
-	time_steps,
-)
-	# plot the positions
-	middle_index = Int(length(states[1]) / 2) - 1
-	pos_init = [x - states[1][middle_index] for x in states[1][1:2:end]]
-	pos_final = [x - states[end][middle_index] for x in states[end][1:2:end]]
-	p1 = plot(pos_init, label = "Initial", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", title = "Positions Over Time")
-	plot!(p1, pos_final, label = "Final", marker = :circle)
-
-	# plot the vels
-	vels_init = [x / M for x in states[1][2:2:end]]
-	vels_final = [x / M for x in states[end][2:2:end]]
-	p2 = plot(states[1][2:2:end], label = "Initial", marker = :circle, xlabel = "Mass Number", ylabel = "Velocity", title = "Velocities Over Time")
-	plot!(p2, states[end][2:2:end], label = "Final t = $(time_steps[end])", marker = :circle)
-
-	# save the plots
-	savefig(p1, "t/initial-final-positions.png")
-	savefig(p2, "t/initial-final-velocities.png")
-end
-
-function analyize_vels(
-	states,
-	time_steps,
-	threshold = 1.975,
-)
-	println("Analyzing velocities:\n")
-	output = []
-	for i in 1:length(states)
-		v = states[i][Int(length(states[i]) / 2)] / M
-		if v >= threshold - 10e-6
-			push!(output, i)
-			println("Time: ", time_steps[i], " Vel: ", v)
-		end
-	end
-
-	data = []
-	for i in 1:length(states)
-		push!(data, states[i][Int(length(states[i]) / 2)])
-	end
-	p = plot(data, label = "Velocity Over Time", xlabel = "Time", ylabel = "Velocity")
-	savefig(p, "t/velocity-over-time.png")
-
-	println("\nDone!\n\n")
-	return output
-end
-
-function analyize_pos(
-	states,
-	time_steps,
-	threshold = 1.975,
-)
-	println("Analyzing positions:\n")
-	output = []
-	for i in 1:length(states)
-		if states[i][Int(length(states[i]) / 2)-1] >= threshold
-			push!(output, i)
-			println("Time: ", time_steps[i], " Position: ", states[i][Int(length(states[i]) / 2)])
-		end
-	end
-
-	# plot the first 10 seconds of Velocity
-	data = []
-	for i in 1:length(states)
-		if time_steps[i] > 10
-			break
-		end
-		push!(data, states[i][Int(length(states[i]) / 2)] - 1)
-	end
-	p = plot(data, label = "Position Over Time", xlabel = "Time", ylabel = "Position")
-	savefig(p, "t/pos-over-time.png")
-
-	println("\nDone!\n\n")
-	return output
-end
-
-function calculate_energy(state)
-	# calculate the kinetic energy
-	kinetic_energy = 0
-	vels = state[2:2:end]
-	for i in 1:N
-		mass = i == Int(N / 2) - 1 ? M : m
-		# calculate the kinetic energy
-		kinetic_energy += 0.5 * vels[i] * vels[i] / mass
-	end
-
-	# calcaute the potential energy
-	potential_energy = 0
-	pos = state[1:2:end]
-	for i in 1:N-1
-		left_index = max(1, i - 1)
-		right_index = min(N, i + 1)
-		potential_energy += 0.5 * K * (pos[left_index] - pos[i])^2
-		potential_energy += 0.5 * K * (pos[right_index] - pos[i])^2
-	end
-
-	# print the energy
-	println("Kinetic Energy: ", kinetic_energy)
-	println("Potential Energy: ", potential_energy)
-	println("Total Energy: ", kinetic_energy + potential_energy, "\n")
-end
-
-function plot_middle_mass_phase_space(states)
-	# get the index of the middle mass
-	middle_index = Int(length(states[1]) / 2) - 1
-	# build an array of the pos and vel over time
-	pos = []
-	vel = []
-	for i in 1:length(states)
-		push!(pos, states[i][middle_index])
-		push!(vel, states[i][middle_index+1])
-	end
-
-	# plot the phase space
-	p = plot(pos, vel, xlabel = "Position", ylabel = "Momentum", title = "Phase Space of Middle Mass in FPU", label = "Beta = 10, K = 1, N = 64, M = $M")
-
-	# save the plot
-	savefig(p, "t/phase-space.png")
-end
-
-
-# Run the simulation
-N = 64 # number of masses
-K = 1 # spring constant
-final_time = 1000 # seconds
-plot_data = []
-
-my_vel = 100
-
-sol = run_simulation(N, K, m, final_time, 0, my_vel)
-
-println("final time: ", sol.t[end])
-# s = sol.u[1:2:end]
-# analyize_vels(sol.u, sol.t, my_vel)
-# analyize_pos(sol.u, sol.t, 1.4)
-plot_starting_and_final_positions(sol.u, sol.t)
-# animate_positions(sol.u, sol.t, 0, 100, my_vel) # expect 80913.35854226245 for k=10?? rip
-plot_middle_mass_phase_space(sol.u)
diff --git a/t/old/initial-final-positions.png b/t/old/initial-final-positions.png
deleted file mode 100644
index e76793d..0000000
--- a/t/old/initial-final-positions.png
+++ /dev/null
diff --git a/t/old/initial-final-velocities.png b/t/old/initial-final-velocities.png
deleted file mode 100644
index b7cc859..0000000
--- a/t/old/initial-final-velocities.png
+++ /dev/null
diff --git a/t/old/modes-beta15.png b/t/old/modes-beta15.png
deleted file mode 100644
index dbb632d..0000000
--- a/t/old/modes-beta15.png
+++ /dev/null
diff --git a/t/old/phase-space-.5vel.png b/t/old/phase-space-.5vel.png
deleted file mode 100644
index 100cd2e..0000000
--- a/t/old/phase-space-.5vel.png
+++ /dev/null
diff --git a/t/old/phase-space-chos.png b/t/old/phase-space-chos.png
deleted file mode 100644
index cbbbf09..0000000
--- a/t/old/phase-space-chos.png
+++ /dev/null
diff --git a/t/old/phase-space-stable.png b/t/old/phase-space-stable.png
deleted file mode 100644
index 2ffa7c6..0000000
--- a/t/old/phase-space-stable.png
+++ /dev/null
diff --git a/t/old/phase-space.01strings.png b/t/old/phase-space.01strings.png
deleted file mode 100644
index 5e78901..0000000
--- a/t/old/phase-space.01strings.png
+++ /dev/null
diff --git a/t/old/phase-space.png b/t/old/phase-space.png
deleted file mode 100644
index 48e5368..0000000
--- a/t/old/phase-space.png
+++ /dev/null
diff --git a/t/old/plot_data.png b/t/old/plot_data.png
deleted file mode 100644
index d052983..0000000
--- a/t/old/plot_data.png
+++ /dev/null
diff --git a/t/old/plz.gif b/t/old/plz.gif
deleted file mode 100644
index febf678..0000000
--- a/t/old/plz.gif
+++ /dev/null
diff --git a/t/old/plz.mp4 b/t/old/plz.mp4
deleted file mode 100644
index 6a0d7b4..0000000
--- a/t/old/plz.mp4
+++ /dev/null
diff --git a/t/old/pos-over-time.png b/t/old/pos-over-time.png
deleted file mode 100644
index e3a732d..0000000
--- a/t/old/pos-over-time.png
+++ /dev/null
diff --git a/t/old/r.jl b/t/old/r.jl
deleted file mode 100644
index 516ab09..0000000
--- a/t/old/r.jl
+++ /dev/null
@@ -1,134 +0,0 @@
-function dynamics!(
-	state,
-	prev_state,
-	params,
-	states,
-	plot_data,
-)
-	# Unpack the parameters
-	N, K, beta, A, dt, num_steps = params
-
-	for j in 1:num_steps
-		next_state = zeros(N)
-
-		# update the left particle (set to rest)
-		state[1] = 0
-
-		# update the right particle (set to rest)
-		state[N] = 0
-
-		# update the middle particles
-		for i in 2:N-1
-			a = 2 * state[i] - prev_state[i]
-			b = dt * dt * K * (state[i+1] - 2 * state[i] + state[i-1])
-			c = dt * dt * beta * ((state[i+1] - state[i])^3 + (state[i-1] - state[i])^3)
-
-			this_mass = 1
-			if i == Int(N / 2)
-				# time = i * dt
-				# push!(plot_data, (state[i]))
-				this_mass = 1
-			end
-			next_state[i] = (a + b + c) / this_mass
-		end
-
-		# push!(states, next_state)
-
-		# update the previous state
-		for i in 1:N
-			prev_state[i] = state[i]
-		end
-		# update the current state
-		for i in 1:N
-			state[i] = next_state[i]
-		end
-
-		# if the middle mass is close to 2, print the time and position
-		if state[Int(N / 2)] > 1.995
-			println("Time: ", j * dt, " Positions: ", state[Int(N / 2)])
-			println("Other Positions: ", state, "\n")
-			push!(states, (j * dt, state))
-		end
-
-		if (j * dt % 100000) == 0
-			println("TIME UPDATE: ", j * dt)
-		end
-	end
-end
-
-function get_initial_state(
-	N,
-	modes,
-	beta,
-	A,
-)
-	state = zeros(N)
-	# amp = A * sqrt(2 / (N - 1))
-	# for i in 2:N-1
-	# 	state[i] = amp * sin((modes * π * (i - 1)) / (N - 1))
-	# end
-	# set the middle to be the largest
-	state[Int(N / 2)] = 2
-	return state
-end
-
-function run_simulation(
-	N,
-	modes,
-	beta,
-	A,
-	dt,
-	num_steps,
-	plot_data,
-)
-	states = []
-	state = get_initial_state(N, 1, beta, A)
-	prev_state = zeros(N)
-	for i in 1:N
-		prev_state[i] = state[i]
-	end
-	params = (N, modes, beta, A, dt, num_steps)
-	dynamics!(state, prev_state, params, states, plot_data)
-	return states
-end
-
-# Run the simulation
-N = 10 # number of masses
-beta = 0 # cubic string spring
-K = 100 # spring constant
-A = 10 # amplitude
-modes = 3 # number of modes to plot
-final_time = 10000000 # seconds
-dt = 0.001 # seconds
-num_steps = Int(final_time / dt)
-params = (N, K, beta, A, dt, num_steps)
-plot_data = []
-s = run_simulation(N, K, beta, A, dt, num_steps, plot_data)
-println("\nStates of interest:", s)
-
-# plot the inital positions and the final positions
-# using Plots
-# plot(get_initial_state(N, 1, beta, A), label = "Initial", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", title = "First Three Modes")
-# plot!(s[end], label = "Final", marker = :circle)
-
-# plot the middle of the lattice over time
-# time_steps = [i * dt for i in 1:num_steps]
-# plot(time_steps, plot_data, label = "Middle Mass Over Time", xlabel = "Time", ylabel = "Displacement")
-# savefig("t/plot_data.png")
-
-# print the points close to the origional displacement
-# output = []
-# for i in 1:length(plot_data)
-# 	if plot_data[i] > 1.975
-# 		push!(output, i)
-# 		println("Time: ", i * dt, " Position: ", plot_data[i])
-# 	end
-# end
-
-
-# animate the s array of positions
-# anim = @animate for i in 1:length(s)
-# 	plot(s[i], label = "t = $(round(i * dt, digits = 3))", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", ylim = (-3, 3))
-# end
-# mp4(anim, "t/plz.mp4", fps = 30)
-# println("Output: ", output)
diff --git a/t/old/semi-old b/t/old/semi-old
deleted file mode 100644
index e69de29..0000000
--- a/t/old/semi-old
+++ /dev/null
diff --git a/t/old/t.jl b/t/old/t.jl
deleted file mode 100644
index 6a06777..0000000
--- a/t/old/t.jl
+++ /dev/null
@@ -1,755 +0,0 @@
-# molecular dynamics 2d.
-
-# usage:
-# at the end of this script, under the header "DEMOS", 
-# you'll see some functions which implement demos from GN chapter 9.
-# simply load the script in your development environment 
-# (I strongly recommend not using jupiter)
-# and in the console/REPL run 
-#	demo_0()
-# etc.
-
-# demos 0,1,3 can optionally make an animated gif
-# if you call it with the optional argument demo_3(gif=1)
-
-# lmk if this script is giving you grief or if you find any bugs
-# kian@brown.edu
-
-using Statistics
-using StatsPlots
-
-mutable struct ParticleSystem
-	N::Int64# number of particles
-	L::Float64# square box side length
-	T₀::Float64# initial temperature
-	t::Float64# system time
-	dt::Float64# time step
-	state::Vector{Float64}# state space array
-	steps::Int64# number of steps
-	sampleInterval::Int64# interval for sampling data
-	timeData::Vector{Float64}# array of sampled time points
-	energyData::Vector{Float64}# array of sampled energy values
-	tempData::Vector{Float64}# array of sampled temperature values
-	tempAccumulator::Float64# temperature accumulator
-	squareTempAccumulator::Float64# T^2 accumulator
-	virialAccumulator::Float64# virial accumulator
-	xData::Vector{Vector{Float64}} # array of sampled position data
-	vData::Vector{Vector{Float64}} # array of sampled velocity data
-	forceType::String# force
-end
-
-function ParticleSystem(N::Int64 = 64, L::Float64 = 10.0, T₀::Float64 = 1.0)
-	t = 0.0
-	dt = 0.001
-	state = zeros(4N) # state space array, [x1,y1,x2,y2,...,vx1,vy1,...]
-	steps = 0
-	timeData = Float64[]
-	energyData = Float64[]
-	sampleInterval = 100
-	tempData = Float64[]
-	tempAccumulator = 0.0
-	squareTempAccumulator = 0.0
-	virialAccumulator = 0.0
-	xData = Vector{Float64}[]
-	vData = Vector{Float64}[]
-	forceType = "lennardJones"
-
-	return ParticleSystem(
-		N,
-		L,
-		T₀,
-		t,
-		dt,
-		state,
-		steps,
-		sampleInterval,
-		timeData,
-		energyData,
-		tempData,
-		tempAccumulator,
-		squareTempAccumulator,
-		virialAccumulator,
-		xData,
-		vData,
-		forceType,
-	)
-end
-
-# some useful "views" of the state array
-# (read the performance tips chapter of the julia manual)
-@views positions(state) = state[1:Int64(length(state) / 2)]
-@views velocities(state) = state[(Int64(length(state) / 2)+1):end]
-@views xcomponent(vector) = vector[1:2:end]
-@views ycomponent(vector) = vector[2:2:end]
-@views particle(n, vector) = [vector[2n-1], vector[2n]]
-
-# INITIALIZATION
-################################################################################
-
-function set_random_positions!(sys::ParticleSystem)
-	println("\tposition configuration: random")
-	positions(sys.state) .= rand(2 * sys.N) .* sys.L
-	cool!(sys)
-end
-
-function set_square_lattice_positions!(sys::ParticleSystem)
-	println("\tposition configuration: square lattice")
-
-	n = Int64(floor(sqrt(sys.N))) # num lattice points per axis
-	latticeSpacing = sys.L / n
-
-	if sys.N != n^2
-		println("\t\toops... your chosen N=$(sys.N) is not a square number")
-		println("\t\t-> resetting N to $(n^2).")
-		sys.N = n^2
-		sys.state = zeros(4 * sys.N)
-	end
-
-	for i in 0:(n-1)
-		for j in 0:(n-1)
-			sys.state[2*(i*n+j)+1] = (i + 0.5) * latticeSpacing
-			sys.state[2*(i*n+j)+2] = (j + 0.5) * latticeSpacing
-		end
-	end
-end
-
-function set_triangular_lattice_positions!(sys::ParticleSystem)
-end
-
-function add_position_jitter!(sys::ParticleSystem, jitter::Float64 = 0.5)
-	println("\tadding a wee bit of random jitter to particle positions...")
-
-	for i ∈ 1:length(positions(sys.state))
-		sys.state[i] += rand() - jitter
-	end
-end
-
-function set_random_velocities!(sys::ParticleSystem)
-	println("\tvelocity configuration: random")
-
-	velocities(sys.state) .= rand(2 * sys.N) .- 0.5
-	zero_total_momentum!(sys)
-	velocities(sys.state) .*= sqrt(sys.T₀ / temperature(sys))
-end
-
-function zero_total_momentum!(sys::ParticleSystem)
-	xcomponent(velocities(sys.state)) .-=
-		mean(xcomponent(velocities(sys.state)))
-	ycomponent(velocities(sys.state)) .-=
-		mean(ycomponent(velocities(sys.state)))
-end
-
-
-# FORCES / POTENTIALS
-################################################################################
-
-function force(sys::ParticleSystem)
-	if sys.forceType == "lennardJones"
-		force, virial = lennard_jones_force(sys)
-	elseif sys.forceType == "powerLaw"
-		force, virial = power_law_force(sys)
-	end
-
-	sys.virialAccumulator += virial
-
-	return force
-end
-
-# the minimum image approximation
-# (periodic boundary conditions)
-function minimum_image(xij::Float64, L::Float64)
-	if xij > (L / 2)
-		xij -= L
-	elseif xij < -(L / 2)
-		xij += L
-	end
-	return xij
-end
-
-function lennard_jones_force(sys::ParticleSystem)
-	x = xcomponent(positions(sys.state))
-	y = ycomponent(positions(sys.state))
-	virial = 0.0
-	force = zeros(2 * sys.N)
-
-	Threads.@threads for i ∈ 1:(sys.N-1)
-		for j ∈ (i+1):sys.N
-			dx = minimum_image(x[i] - x[j], sys.L)
-			dy = minimum_image(y[i] - y[j], sys.L)
-
-			r2inv = 1.0 / (dx^2 + dy^2)
-			f = 48.0 * r2inv^7 - 24.0 * r2inv^4
-			fx = dx * f
-			fy = dy * f
-
-			force[2*i-1] += fx
-			force[2*i] += fy
-			force[2*j-1] -= fx
-			force[2*j] -= fy
-
-			virial += fx * dx + fy * dy
-		end
-	end
-
-	return force, 0.5 * virial
-end
-
-function lennard_jones_potential(sys::ParticleSystem)
-	x = xcomponent(positions(sys.state))
-	y = ycomponent(positions(sys.state))
-	U = 0.0
-
-	Threads.@threads for i in 1:(sys.N-1)
-		for j in (i+1):sys.N
-			dx = minimum_image(x[i] - x[j], sys.L)
-			dy = minimum_image(y[i] - y[j], sys.L)
-
-			r2inv = 1.0 / (dx^2 + dy^2)
-			U += r2inv^6 - r2inv^3
-		end
-	end
-	return 4.0 * U
-end
-
-function power_law_force(sys::ParticleSystem)
-end
-
-function power_law_potential(sys::ParticleSystem)
-end
-
-# TIME EVOLUTION
-################################################################################
-
-function keep_particles_in_box!(sys::ParticleSystem)
-	for i in 1:(2*sys.N)
-		if positions(sys.state)[i] > sys.L
-			positions(sys.state)[i] -= sys.L
-		elseif positions(sys.state)[i] < 0.0
-			positions(sys.state)[i] += sys.L
-		end
-	end
-
-	#	# another way of doing this: use the ternary operator
-	#	for i in 1:(2 * sys.N)
-	#		positions(sys.state)[i] < 0.0 ?  
-	#		positions(sys.state)[i] % sys.L + sys.L : 
-	#		positions(sys.state)[i] % sys.L
-	#	end
-end
-
-function verlet_step!(sys::ParticleSystem)
-	# compute acceleration at current time
-	acceleration = force(sys)
-
-	# compute positions at t + dt
-	positions(sys.state) .+=
-		velocities(sys.state) .* sys.dt .+
-		0.5 .* acceleration .* (sys.dt)^2
-
-	# enforce boundary conditions
-	# (basically check if any particles left the box and put them back)
-	# see function implementation for deets
-	keep_particles_in_box!(sys)
-
-	# compute velocities at t + dt
-	velocities(sys.state) .+=
-		0.5 * sys.dt .* (acceleration + force(sys))
-end
-
-function evolve!(sys::ParticleSystem, runtime::Float64 = 10.0)
-	numsteps = Int64(abs(runtime / sys.dt) + 1)
-
-	print_evolution_message(runtime, numsteps)
-
-	@time for step in 1:numsteps
-		verlet_step!(sys)
-		zero_total_momentum!(sys)
-
-		if (step % sys.sampleInterval == 1)
-			push!(sys.timeData, sys.t)
-			push!(sys.energyData, energy(sys))
-			push!(sys.xData, positions(sys.state))
-			push!(sys.vData, velocities(sys.state))
-
-			T = temperature(sys)
-			push!(sys.tempData, T)
-			sys.tempAccumulator += T
-			sys.squareTempAccumulator += T^2
-		end
-
-		sys.t += sys.dt
-		sys.steps += 1
-	end
-	println("done.")
-end
-
-function reverse_time!(sys)
-	sys.dt *= -1
-	println("\ntime reversed! dt = $(sys.dt)")
-end
-
-function cool!(sys::ParticleSystem, cooltime::Float64 = 1.0)
-	numsteps = Int64(cooltime / sys.dt)
-	for step in 1:numsteps
-		verlet_step!(sys)
-		velocities(sys.state) .*= (1.0 - sys.dt)
-	end
-	reset_statistics!(sys)
-end
-
-# MEASUREMENTS
-################################################################################
-
-function kinetic_energy(sys::ParticleSystem)
-	return 0.5 * sum(velocities(sys.state) .* velocities(sys.state))
-end
-
-function potential_energy(sys::ParticleSystem)
-	return lennard_jones_potential(sys)
-end
-
-function temperature(sys::ParticleSystem)
-	return kinetic_energy(sys) / sys.N
-end
-
-function energy(sys::ParticleSystem)
-	return potential_energy(sys) + kinetic_energy(sys)
-end
-
-# STATISTICS
-################################################################################
-
-function reset_statistics!(sys::ParticleSystem)
-	sys.steps = 0
-	sys.tempAccumulator = 0.0
-	sys.squareTempAccumulator = 0.0
-	sys.virialAccumulator = 0.0
-	sys.xData = []
-	sys.vData = []
-end
-
-function mean_temperature(sys::ParticleSystem)
-	return sys.tempAccumulator / sys.steps
-end
-
-function mean_square_temperature(sys::ParticleSystem)
-	return sys.squareTempAccumulator / sys.steps
-end
-
-function mean_pressure(sys::ParticleSystem)
-	# factor of half because force is calculated twice each step
-	meanVirial = 0.5 * sys.virialAccumulator / sys.steps
-	return 1.0 + 0.5 * meanVirial / (sys.N * mean_temperature(sys))
-end
-
-function heat_capacity(sys::ParticleSystem)
-	meanTemperature = mean_temperature(sys)
-	meanSquareTemperature = mean_square_temperature(sys)
-	σ2 = meanSquareTemperature - meanTemperature^2
-	denom = 1.0 - σ2 * sys.N / meanTemperature^2
-	return sys.N / denom
-end
-
-function mean_energy(sys::ParticleSystem)
-	return mean(sys.energyData)
-end
-
-function std_energy(sys::ParticleSystem)
-	return std(sys.energyData)
-end
-
-# MATH / ADDITIONAL FUNCTIONS
-################################################################################
-
-function dot(v1::Vector{Float64}, v2::Vector{Float64})
-	return sum(v1 .* v2)
-end
-
-# GRAPHS
-################################################################################
-
-function initialize_plot()
-	plot(
-		size = (800, 800),
-		titlefontsize = 12,
-		guidefontsize = 12,
-	)
-end
-
-function plot_positions_t(sys::ParticleSystem, t::Int64)
-	initialize_plot()
-	for n ∈ 1:sys.N
-		scatter!(
-			[sys.xData[t][2n-1]],
-			[sys.xData[t][2n]],
-			markersize = 4.0,
-			markercolor = n,
-			markerstrokewidth = 0.4,
-			grid = true,
-			framestyle = :box,
-			legend = false,
-		)
-	end
-end
-
-function animate(sys::ParticleSystem, interval::Int64 = 1)
-	println("\ngenerating gif...")
-
-	scatter!()
-	animation = @animate for t in 1:length(sys.xData)
-		scatter()
-		for n ∈ 1:sys.N
-			scatter!(
-				[sys.xData[t][2n-1]],
-				[sys.xData[t][2n]],
-				#markersize = 4.0,
-				markercolor = n,
-				#markerstrokewidth = 0.4,
-				grid = true,
-				framestyle = :box,
-				legend = false,
-			)
-		end
-		xlims!(0, sys.L)
-		ylims!(0, sys.L)
-		xlabel!("x")
-		ylabel!("y")
-	end every interval
-
-	gif(animation, "./animation.gif")
-	println("done.")
-end
-
-function plot_positions(sys::ParticleSystem)
-	initialize_plot()
-	for n ∈ 1:sys.N
-		scatter!(
-			[xcomponent(positions(sys.state))[n]],
-			[ycomponent(positions(sys.state))[n]],
-			markersize = 4.0,
-			markercolor = n,
-			markerstrokewidth = 0.4,
-			grid = true,
-			framestyle = :box,
-			legend = false,
-		)
-	end
-	xlims!(0, sys.L)
-	ylims!(0, sys.L)
-	xlabel!("x")
-	ylabel!("y")
-	title!("positions at t=$(round(sys.t, digits=4))")
-end
-
-function plot_trajectories(sys::ParticleSystem, particles::Vector{Int64} = [1])
-	initialize_plot()
-	for n ∈ 1:sys.N
-		scatter!(
-			[xcomponent(positions(sys.state))[n]],
-			[ycomponent(positions(sys.state))[n]],
-			markersize = 4.0,
-			markercolor = n,
-			markerstrokewidth = 0.4,
-			grid = true,
-			framestyle = :box,
-			legend = false,
-		)
-	end
-
-	for n in collect(particles)
-		xdata = [sys.xData[i][2n-1] for i in 1:length(sys.xData)]
-		ydata = [sys.xData[i][2n] for i in 1:length(sys.xData)]
-
-		# plot trajectory line for nth particle
-		scatter!(
-			xdata,
-			ydata,
-			color = n,
-			#markerstrokewidth = 0,
-			markerstrokecolor = n,
-			markersize = 0.7,
-			markeralpha = 0.5,
-			label = false,
-			widen = false,
-		)
-
-		# plot initial position for nth particle
-		scatter!(
-			[sys.xData[1][2n-1]],
-			[sys.xData[1][2n]],
-			markersize = 4.0,
-			markercolor = n,
-			markerstrokewidth = 0.4,
-			markeralpha = 0.3,
-			#label = "pcl. $n @t=t₀",
-			widen = false,
-		)
-
-		# plot final position for nth particle
-		scatter!(
-			[sys.xData[end][2n-1]],
-			[sys.xData[end][2n]],
-			markersize = 4.0,
-			markercolor = n,
-			markerstrokewidth = 0.4,
-			markeralpha = 1.0,
-			#label = "pcl $n @t=t",
-			widen = false,
-		)
-	end
-	title!("positions & trajectories at time t=$(round(sys.t, digits=2))")
-	plot!()
-end
-
-function plot_temperature(sys::ParticleSystem)
-	initialize_plot()
-	plot!(
-		sys.timeData,
-		sys.tempData,
-		#widen = true,
-	)
-	ylims!(
-		mean(sys.tempData) - std(sys.tempData) * 20,
-		mean(sys.tempData) + std(sys.tempData) * 20,
-	)
-	xlabel!("t")
-	ylabel!("T(t)")
-	title!("temperature vs time")
-
-end
-
-function plot_energy(sys::ParticleSystem, ylimit::Float64 = 1.0)
-	initialize_plot()
-	plot!(
-		sys.timeData,
-		sys.energyData,
-		#widen = true,
-	)
-	ylims!(
-		#ylimit * (mean(sys.energyData) - 1), 
-		#ylimit * (mean(sys.energyData) + 1)
-		mean(sys.energyData) - std(sys.energyData) * 10,
-		mean(sys.energyData) + std(sys.energyData) * 10,
-	)
-	xlabel!("t")
-	ylabel!("E(t)")
-	title!("energy vs time")
-end
-
-function plot_speed_distribution(sys::ParticleSystem, numSamples::Int64 = 5)
-	initialize_plot()
-
-	numDataPoints = Int64(length(sys.vData))
-	interval = Int64(floor(numDataPoints / numSamples))
-
-	samples = collect(1:interval:numDataPoints)
-	for s in samples
-		speed = sqrt.(
-			xcomponent(sys.vData[s]) .^ 2 .*
-			ycomponent(sys.vData[s]) .^ 2
-		)
-		density!(
-			sys.vData[s],
-			normalize = :pdf,
-			label = "t = $(round(sys.timeData[s], digits=2))",
-		)
-	end
-	xlabel!("speed")
-	title!("speed distribution")
-end
-
-# CONSOLE PRINT DATA
-################################################################################
-
-function print_hello()
-	println("\nmolecular dynamics!")
-	println("number of threads: ", Threads.nthreads())
-end
-
-function print_bonjour()
-	println("\nbonjour")
-end
-
-function print_system_parameters(sys::ParticleSystem)
-	println("\nsystem parameters:")
-	println("\tN =  $(sys.N)   (number of particles)")
-	println("\tL =  $(sys.L)   (side length of square box)")
-	println("\tDT = $(sys.dt)  (time step)")
-end
-
-function print_system_data(sys::ParticleSystem)
-	println("\nsystem data at time t=$(round(sys.t, digits=4))")
-
-	if sys.steps == 0
-		println("\ttemperature:     $(temperature(sys))")
-		println("\tenergy:          $(energy(sys))")
-	else
-		println("\tsteps evolved:   $(sys.steps)")
-		println("\ttemperature:     $(temperature(sys))")
-		println("\tenergy:          $(energy(sys))")
-		println("\tmean energy:     $(mean_energy(sys))")
-		println("\tstd energy:      $(std_energy(sys))")
-		println("\theat capacity:   $(heat_capacity(sys))")
-		println("\tPV/NkT:          $(mean_pressure(sys))")
-	end
-end
-
-function print_evolution_message(runtime, numsteps)
-	println("\nevolving...")
-end
-
-# DEMOS
-################################################################################
-
-
-# DEMO 0: APPROACH TO EQUILIBRIUM
-function demo_0(; gif = 0)
-	println("\nDEMO 0: APPROACH TO EQUILIBRIUM")
-	println("----------------------------------------")
-
-	sys = ParticleSystem(64, 120.0, 1.0)
-	print_system_parameters(sys)
-
-	set_square_lattice_positions!(sys)
-	set_random_velocities!(sys)
-	print_system_data(sys)
-	p1 = plot_positions(sys)
-
-	evolve!(sys, 20.0)
-	print_system_data(sys)
-
-	p2 = plot_trajectories(sys, collect(1:64))
-	p3 = plot_energy(sys)
-	p4 = plot_temperature(sys)
-
-	# make gif
-	if gif == 1
-		animate(sys, 1)
-	end
-
-	plot(
-		p1, p2, p3, p4,
-		layout = grid(2, 2, heights = [0.7, 0.3]),
-		size = (1280, 720),
-	)
-end
-
-# DEMO 1: TIME REVERSAL TEST
-function demo_1(; gif = 0)
-	println("\nDEMO 1: TIME REVERSAL TEST")
-	println("----------------------------------------")
-
-	sys = ParticleSystem(64, 120.0, 1.0)
-	print_system_parameters(sys)
-
-	set_square_lattice_positions!(sys)
-	set_random_velocities!(sys)
-	print_system_data(sys)
-	p1 = plot_positions(sys)
-
-	evolve!(sys, 50.0)
-	#p2 = plot_trajectories(sys, collect(1:64))
-	p2 = plot_positions(sys)
-
-	reverse_time!(sys)
-	evolve!(sys, 50.0)
-	print_system_data(sys)
-	#p3 = plot_trajectories(sys, collect(1:64))
-	p3 = plot_positions(sys)
-
-	# make gif
-	if gif == 1
-		animate(sys, 4)
-	end
-
-	plot(
-		p1, p2, p3,
-		layout = (1, 3),
-		size = (1200, 400),
-	)
-end
-
-# DEMO 2: SPEED DISTRIBUTION
-function demo_2()
-	println("\nDEMO 2: SPEED DISTRIBUTION")
-	println("----------------------------------------")
-
-	sys = ParticleSystem[]
-
-	# array for speed distribution plots
-	ps = Plots.Plot{Plots.GRBackend}[]
-
-	# array for trajectory plots
-	pt = Plots.Plot{Plots.GRBackend}[]
-
-	# initialize three systems with different initial conditions 
-	# but same KE and PE, evolve, and save plots
-	for i ∈ 1:3
-		push!(sys, ParticleSystem(64, 120.0, 1.0))
-
-		println("\nSYSTEM $i")
-		print_system_parameters(sys[i])
-
-		set_square_lattice_positions!(sys[i])
-		add_position_jitter!(sys[i])
-		set_random_velocities!(sys[i])
-		print_system_data(sys[i])
-
-		evolve!(sys[i], 40.0)
-		print_system_data(sys[i])
-		push!(ps, plot_speed_distribution(sys[i], 5))
-		push!(pt, plot_trajectories(sys[i], collect(1:64)))
-	end
-
-
-	# plot speed distribution and trajectory plots
-	plot(
-		ps[1], ps[2], ps[3],
-		pt[1], pt[2], pt[3],
-		layout = (2, 3),
-		size = (1920, 1080),
-	)
-end
-
-# DEMO 3: MELTING TRANSITION
-function demo_3(; gif = 0)
-	println("\nDEMO 3: MELTING TRANSITION")
-	println("----------------------------------------")
-
-	# initialize system of particles on square lattice with zero velocity
-	sys = ParticleSystem(100, 10.0, 5.0)
-	set_square_lattice_positions!(sys)
-	print_system_data(sys)
-	p1 = plot_positions(sys)
-
-	# evolve the system and watch them "crystallize" 
-	# into a triangular lattice formation
-	evolve!(sys, 20.0)
-	print_system_data(sys)
-	p2 = plot_trajectories(sys, collect(1:100))
-
-	# now, increase the temperature of the system by giving the particles
-	# some velocity. evolve the system and plot the trajectories.
-	set_random_velocities!(sys)
-	evolve!(sys, 60.0)
-	print_system_data(sys)
-	p3 = plot_trajectories(sys, collect(1:100))
-
-	# some more plots
-	p4 = plot_energy(sys, 0.0)
-	p5 = plot_temperature(sys)
-	p6 = plot_speed_distribution(sys, 20)
-
-	# make gif
-	if gif == 1
-		animate(sys, 1)
-	end
-
-	plot(
-		p1, p2, p3, p4, p5, p6,
-		layout = (2, 3),
-		size = (1280, 720),
-	)
-end
-
-demo_0()
diff --git a/t/old/velocity-over-time.png b/t/old/velocity-over-time.png
deleted file mode 100644
index 3cab01f..0000000
--- a/t/old/velocity-over-time.png
+++ /dev/null