25 files changed, 1270 insertions, 0 deletions

diff --git a/t/old/1d.jl b/t/old/1d.jl
new file mode 100644
index 0000000..203466c
--- /dev/null
+++ b/t/old/1d.jl
@@ -0,0 +1,109 @@
+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
new file mode 100644
index 0000000..0964313
--- /dev/null
+++ b/t/old/animate-positions-.1stringa.mp4
diff --git a/t/old/animate-positions-.5vel.mp4 b/t/old/animate-positions-.5vel.mp4
new file mode 100644
index 0000000..b56e372
--- /dev/null
+++ b/t/old/animate-positions-.5vel.mp4
diff --git a/t/old/animate-positions-10.mp4 b/t/old/animate-positions-10.mp4
new file mode 100644
index 0000000..aa00547
--- /dev/null
+++ b/t/old/animate-positions-10.mp4
diff --git a/t/old/animate-positions-80k.mp4 b/t/old/animate-positions-80k.mp4
new file mode 100644
index 0000000..3551bf4
--- /dev/null
+++ b/t/old/animate-positions-80k.mp4
diff --git a/t/old/animate-positions-chaos.mp4 b/t/old/animate-positions-chaos.mp4
new file mode 100644
index 0000000..7edbc75
--- /dev/null
+++ b/t/old/animate-positions-chaos.mp4
diff --git a/t/old/animate-positions-stable.mp4 b/t/old/animate-positions-stable.mp4
new file mode 100644
index 0000000..0b858a2
--- /dev/null
+++ b/t/old/animate-positions-stable.mp4
diff --git a/t/old/animate-positions-t.mp4 b/t/old/animate-positions-t.mp4
new file mode 100644
index 0000000..a4a34e0
--- /dev/null
+++ b/t/old/animate-positions-t.mp4
diff --git a/t/old/disp.jl b/t/old/disp.jl
new file mode 100644
index 0000000..18e0a1b
--- /dev/null
+++ b/t/old/disp.jl
@@ -0,0 +1,272 @@
+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
new file mode 100644
index 0000000..e76793d
--- /dev/null
+++ b/t/old/initial-final-positions.png
diff --git a/t/old/initial-final-velocities.png b/t/old/initial-final-velocities.png
new file mode 100644
index 0000000..b7cc859
--- /dev/null
+++ b/t/old/initial-final-velocities.png
diff --git a/t/old/modes-beta15.png b/t/old/modes-beta15.png
new file mode 100644
index 0000000..dbb632d
--- /dev/null
+++ b/t/old/modes-beta15.png
diff --git a/t/old/phase-space-.5vel.png b/t/old/phase-space-.5vel.png
new file mode 100644
index 0000000..100cd2e
--- /dev/null
+++ b/t/old/phase-space-.5vel.png
diff --git a/t/old/phase-space-chos.png b/t/old/phase-space-chos.png
new file mode 100644
index 0000000..cbbbf09
--- /dev/null
+++ b/t/old/phase-space-chos.png
diff --git a/t/old/phase-space-stable.png b/t/old/phase-space-stable.png
new file mode 100644
index 0000000..2ffa7c6
--- /dev/null
+++ b/t/old/phase-space-stable.png
diff --git a/t/old/phase-space.01strings.png b/t/old/phase-space.01strings.png
new file mode 100644
index 0000000..5e78901
--- /dev/null
+++ b/t/old/phase-space.01strings.png
diff --git a/t/old/phase-space.png b/t/old/phase-space.png
new file mode 100644
index 0000000..48e5368
--- /dev/null
+++ b/t/old/phase-space.png
diff --git a/t/old/plot_data.png b/t/old/plot_data.png
new file mode 100644
index 0000000..d052983
--- /dev/null
+++ b/t/old/plot_data.png
diff --git a/t/old/plz.gif b/t/old/plz.gif
new file mode 100644
index 0000000..febf678
--- /dev/null
+++ b/t/old/plz.gif
diff --git a/t/old/plz.mp4 b/t/old/plz.mp4
new file mode 100644
index 0000000..6a0d7b4
--- /dev/null
+++ b/t/old/plz.mp4
diff --git a/t/old/pos-over-time.png b/t/old/pos-over-time.png
new file mode 100644
index 0000000..e3a732d
--- /dev/null
+++ b/t/old/pos-over-time.png
diff --git a/t/old/r.jl b/t/old/r.jl
new file mode 100644
index 0000000..516ab09
--- /dev/null
+++ b/t/old/r.jl
@@ -0,0 +1,134 @@
+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
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/t/old/semi-old
diff --git a/t/old/t.jl b/t/old/t.jl
new file mode 100644
index 0000000..6a06777
--- /dev/null
+++ b/t/old/t.jl
@@ -0,0 +1,755 @@
+# 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
new file mode 100644
index 0000000..3cab01f
--- /dev/null
+++ b/t/old/velocity-over-time.png