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))