1 files changed, 109 insertions, 0 deletions
diff --git a/final-project/old/1d.jl b/final-project/old/1d.jl
new file mode 100644
index 0000000..203466c
--- /dev/null
+++ b/final-project/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))