1 files changed, 134 insertions, 0 deletions
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)