update naming

1 files changed, 767 insertions, 0 deletions

diff --git a/final-project/fpu.jl b/final-project/fpu.jl
new file mode 100644
index 0000000..962398a
--- /dev/null
+++ b/final-project/fpu.jl
@@ -0,0 +1,767 @@
+function calculate_force(
+	left_pos,
+	middle_pos,
+	right_pos,
+	K,
+	alpha = 0.0,
+	beta = 0.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, params, t)
+	# unpack the params
+	N, modes, beta, A, dt, num_steps = params
+
+	# get the positions and momenta
+	qs = u[1:2:end]
+	ps = u[2:2:end]
+
+	num_masses = length(qs)
+
+	# go over the points in the lattice and update the state
+	for i in 2:num_masses-1
+		# left_index = max(1, i - 1)
+		# right_index = min(N, i + 1)
+
+		du[i*2-1] = ps[i]
+		force = calculate_force(qs[i-1], qs[i], qs[i+1], 1, 0, beta)
+		du[i*2] = force
+	end
+
+	# make last point same as first at rest
+	du[num_masses*2-1] = du[1] = 0 # set to 0
+	du[num_masses*2] = du[2] = 0
+
+	# 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,
+	modes,
+	beta,
+	A,
+)
+	state = zeros(N + 2)
+	amp = A * sqrt(2 / (N + 1))
+	for i in 2:N+1
+		state[i] = amp * sin((modes * π * (i - 1)) / (N + 1))
+	end
+	return state
+end
+
+using DifferentialEquations
+function run_simulation(
+	N,
+	modes,
+	beta,
+	A,
+	dt,
+	num_steps,
+)
+	println("Running simulation with params: N=$N, modes=$modes, beta=$beta, A=$A, dt=$dt, num_steps=$num_steps")
+	# 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)
+
+	# state above is position, need to add in momentum
+	s_0 = zeros(2 * (N + 2))
+	for i in 1:length(state)
+		s_0[i*2-1] = state[i]
+	end
+	final_time = num_steps * dt
+	t_span = (0.0, final_time)
+	prob = ODEProblem(tendency!, s_0, t_span, params)
+	sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8, maxiters = 1e10) # control simulation
+
+	return sol
+end
+
+function get_pos_from_state(state)
+	return state[1:2:end]
+end
+
+# Run the simulation
+# N = 32 # number of masses
+# beta = 2.3 # cubic string spring
+# A = 10 # amplitude
+# modes = 3 # number of modes to plot
+# final_time = 100000 # seconds
+# dt = 0.05 # seconds
+# num_steps = Int(final_time / dt)
+# params = (N, modes, beta, A, dt, num_steps)
+# println("Running simulation with params: ", params)
+# sol = run_simulation(N, modes, beta, A, dt, num_steps)
+# states = sol.u
+# timesteps = sol.t
+# println("Final state: ", states[end])
+
+# plot the inital positions and the final positions
+using Plots
+# init_pos = get_initial_state(N, 1, beta, A)
+# final_pos = get_pos_from_state(states[end])
+# println("Initial: ", init_pos)
+# println("Final: ", final_pos)
+# plot(init_pos, label = "Initial", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", title = "First Three Modes")
+# plot!(final_pos, label = "Final", marker = :circle)
+# println("Saving initial-final-beta15.png")
+# savefig("t/initial-final-beta-$beta.png")
+
+function animate_states(states, timesteps, from = 1, to = 1000)
+	println("\nAnimating states")
+	# animate the s array of positions
+	# find the timestep that maps to the state index
+	i_start = 1
+	i_end = length(timesteps)
+	for i in 1:length(timesteps)
+		if timesteps[i] <= from
+			i_start = i
+			continue
+		end
+		if timesteps[i] >= to
+			i_end = i
+			break
+		end
+	end
+
+	anim = @animate for i in i_start:i_end
+		s = get_pos_from_state(states[i])
+		t = timesteps[i]
+		plot(s, label = "t = $t", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", ylim = (-3, 3))
+	end
+
+	mp4(anim, "t/animated-fpu.mp4", fps = 30)
+	println("Done animating, wrote to animated-fpu.mp4\n")
+end
+# animate_states(states, timesteps, 149500, 150000)
+
+function caluclate_energies_for_mode(states, mode)
+	N = length(states[1]) / 2 - 2
+
+	total = []
+	kinetic = []
+	potential = []
+
+	for i in 1:length(states)
+		total_energy = 0
+		kinetic_energy = 0
+		potential_energy = 0
+
+		# find a_k for this mode
+		poses = get_pos_from_state(states[i])
+		a_k = 0
+		for j in 1:length(poses)
+			a_k += poses[j] * sin((mode * (j - 1) * π) / (N + 1))
+		end
+		amp = sqrt(2 / (N + 1))
+		a_k *= amp
+
+		# get the potenital energy
+		omega_k = 2 * sin((mode * π) / (2 * (N + 1)))
+		p = 0.5 * omega_k^2 * a_k^2
+
+		# find the kinetic energy
+		v_k = 0
+		vels = states[i][2:2:end]
+		for j in 1:length(vels)
+			v_k += vels[j] * sin((mode * (j - 1) * π) / (N + 1))
+		end
+		v_k *= amp
+		k = 0.5 * v_k^2
+
+
+		total_energy += k + p
+		kinetic_energy += k
+		potential_energy += p
+		push!(total, total_energy)
+		push!(kinetic, kinetic_energy)
+		push!(potential, potential_energy)
+	end
+	return total, kinetic, potential
+end
+
+function calculate_total_energy_by_state(states, timesteps)
+	total = []
+	kinetic = []
+	potential = []
+	for i in 1:length(states)
+		vels = states[i][2:2:end]
+		k = 0.5 * sum(vels .^ 2)
+
+		poses = get_pos_from_state(states[i])
+		p = 0
+		for j in 2:length(poses)
+			p += beta * (poses[j] - poses[j-1])^2
+		end
+
+		push!(total, k + p)
+		push!(kinetic, k)
+		push!(potential, p)
+	end
+
+	# plot all three over time on the same plot
+	plot(timesteps, total, label = "Total Energy", xlabel = "Time", ylabel = "Energy", title = "Total Energy Over Time")
+	plot!(timesteps, kinetic, label = "Kinetic Energy")
+	plot!(timesteps, potential, label = "Potential Energy")
+
+	return total, kinetic, potential
+end
+
+function print_energies_for_modes(states, max_mode = 3, index = 1)
+	# get the sum of total energy over all mode
+	total = 0
+	kinetic = 0
+	potential = 0
+	for i in 1:max_mode
+		t, k, p = caluclate_energies_for_mode(states, i)
+		total += t[index]
+		kinetic += k[index]
+		potential += p[index]
+	end
+	println("Total Energy: ", total)
+	println("Kinetic Energy: ", kinetic)
+	println("Potential Energy: ", potential)
+end
+
+function plot_phase_space(states, mass_number, intersecting_mass = 4)
+	# get the positions and momenta
+	qs = []
+	ps = []
+
+	for i in 2:length(states)
+		# only store if the other mass is crossing 0
+		if (states[i-1][intersecting_mass*2-1] < 0 && states[i][intersecting_mass*2-1] > 0) || (states[i-1][intersecting_mass*2-1] > 0 && states[i][intersecting_mass*2-1] < 0)
+			s = states[i]
+			q = s[mass_number*2-1]
+			p = s[mass_number*2]
+			push!(qs, q)
+			push!(ps, p)
+		end
+	end
+
+	s = scatter(
+		qs,
+		ps,
+		label = "Mass $mass_number",
+		xlabel = "Displacement",
+		ylabel = "Momentum",
+		title = "Phase Space for Masses when Mass $intersecting_mass crosses 0",
+		legend = :topleft,
+		marker = :circle,
+		markersize = 2,
+		xlim = (-1.5, 1.5),
+		ylim = (-1.5, 1.5),
+	)
+
+	# savefig(s, "t/phase-space-mass$mass_number-beta15.png")
+
+	return s
+end
+
+function analyze_energies_of_n_modes(modes::Array{Int}, states, timesteps, beta)
+	println("Analyzing energies of modes: ", modes)
+
+	total_by_mode = Dict()
+	kinetic_by_mode = Dict()
+	potential_by_mode = Dict()
+	for j in modes
+		t, k, p = caluclate_energies_for_mode(states, j)
+
+		total_by_mode[j] = t
+		kinetic_by_mode[j] = k
+		potential_by_mode[j] = p
+	end
+
+	total_energy = []
+	kinetic_energy = []
+	potential_energy = []
+	for i in 1:length(states)
+		total = 0
+		kinetic = 0
+		potential = 0
+
+		for j in modes
+			total += total_by_mode[j][i]
+			kinetic += kinetic_by_mode[j][i]
+			potential += potential_by_mode[j][i]
+		end
+
+		push!(total_energy, total)
+		push!(kinetic_energy, kinetic)
+		push!(potential_energy, potential)
+	end
+
+	# plot all three over time on the same plot
+	plot(timesteps, total_energy, label = "Total Energy", xlabel = "Time", ylabel = "Energy", title = "Total Energy Over Time")
+	plot!(timesteps, kinetic_energy, label = "Kinetic Energy")
+	plot!(timesteps, potential_energy, label = "Potential Energy")
+	savefig("t/procs/energies-modes-sum-beta-$beta.png")
+	println("Saved energies-modes-sum-beta-$beta.png. Summary Below:\n")
+
+	# plot the energies for each mode over time
+	plt = plot(title = "Energy Over Time for Modes $(join(modes, ", "))", xlabel = "Time", ylabel = "Energy")
+	for mode in modes
+		plot!(plt, timesteps, total_by_mode[mode], label = "Mode $mode")
+		# plot!(plt, timesteps, kinetic_by_mode[mode], label = "Mode $mode Kinetic")
+		# plot!(plt, timesteps, potential_by_mode[mode], label = "Mode $mode Potential")
+	end
+	savefig(plt, "t/procs/energies-modes-separate-beta-$beta.png")
+
+	# print the energies
+	println("Total Energy Initial: ", total_energy[1])
+	println("Total Energy Final: ", total_energy[end])
+	println("Kinetic Energy Initial: ", kinetic_energy[1])
+	println("Kinetic Energy Final: ", kinetic_energy[end])
+	println("Potential Energy Initial: ", potential_energy[1])
+	println("Potential Energy Final: ", potential_energy[end], "\n")
+
+	println("Done analyzing energies of modes: ", modes, "\n")
+
+	return total_energy, kinetic_energy, potential_energy, total_by_mode, kinetic_by_mode, potential_by_mode
+end
+
+function animate_masses_phase_space(states, mass_numbers, intersecting_mass = 4)
+	println("\nAnimating phase space for mass $mass_numbers when mass $intersecting_mass crosses 0")
+	# animate the s array of positions
+	# find the timestep that maps to the state index
+	anim = @animate for i in mass_numbers
+		s = plot_phase_space(states, i, intersecting_mass)
+		s
+	end
+
+	mp4(anim, "t/animated-phase-space-masses-beta-$beta.mp4", fps = 24)
+end
+
+function plot_mode_phase_space(states, mode_num, intersecting_mode, beta)
+	println("\nPlotting phase space for mode $mode_num")
+
+	# convert the state to a mode
+	a_modes = []
+	vel_modes = []
+	a_inter_modes = []
+	vel_inter_modes = []
+
+	N = length(states[1]) / 2 - 2
+
+	for i in 1:length(states)
+		a_mode = 0
+		vel_mode = 0
+		a_inter_mode = 0
+		vel_inter_mode = 0
+
+		positions = states[i][1:2:end]
+		velocities = states[i][2:2:end]
+		for j in 2:length(positions)-1
+			a_mode += positions[j] * sin((mode_num * (j - 1) * π) / (N + 1))
+			vel_mode += velocities[j] * sin((mode_num * (j - 1) * π) / (N + 1))
+			a_inter_mode += positions[j] * sin((intersecting_mode * (j - 1) * π) / (N + 1))
+			vel_inter_mode += velocities[j] * sin((intersecting_mode * (j - 1) * π) / (N + 1))
+		end
+
+		amp = sqrt(2 / (N + 1))
+		a_mode *= amp
+		vel_mode *= amp
+		a_inter_mode *= amp
+		vel_inter_mode *= amp
+
+		if i == 1
+			push!(a_inter_modes, a_inter_mode)
+			push!(vel_inter_modes, vel_inter_mode)
+			continue
+		end
+
+		if ((a_inter_modes[end] < 0 && a_inter_mode > 0) || (a_inter_modes[end] > 0 && a_inter_mode < 0))
+			push!(a_modes, a_mode)
+			push!(vel_modes, vel_mode)
+		end
+		push!(a_inter_modes, a_inter_mode)
+		push!(vel_inter_modes, vel_inter_mode)
+	end
+
+
+	# plot the phase space
+	scatter(a_modes, vel_modes, label = "Beta = $beta", xlabel = "Displacement", ylabel = "Momentum",
+		title = "Phase Space for Mode $mode_num when a_$intersecting_mode crosses 0",
+		xlim = (-10, 10),
+		ylim = (-1, 1),
+		color = :red,
+		legend = :topright,
+	)
+	# savefig("t/frames/phase-space-mode-$mode_num-beta-$beta.png")
+	println("Saved phase-space-mode-$mode_num-beta-$beta.png")
+end
+
+function animate_phase_space_plots_over_beta(betas, t_max, desired_mode, intersecting_mode)
+	N = 32 # number of masses
+	A = 10 # amplitude
+	modes = 3 # number of modes to plot
+	final_time = t_max # seconds
+	dt = 0.05 # seconds
+	num_steps = Int(final_time / dt)
+
+	results = []
+	for b in betas
+		s = run_simulation(N, modes, b, A, dt, num_steps)
+		states = s.u
+		timesteps = s.t
+
+		push!(results, states)
+
+		println("Simualted beta: ", b, " to t=", timesteps[end])
+	end
+
+	anim = @animate for i in 1:length(betas)
+		plot_mode_phase_space(results[i], desired_mode, intersecting_mode, betas[i])
+	end
+
+	mp4(anim, "t/animated-phase-space-modes-over-beta.mp4", fps = 15)
+end
+
+# animate_phase_space_plots_over_beta(collect(0:0.025:4), 25000, 1, 3)
+
+# analyze_energies_of_n_modes([1, 3, 5, 7, 9, 11], timesteps)
+
+# animate_masses_phase_space(states, collect(2:N+1), 17)
+
+# plot_mode_phase_space(states, 1, 3)
+
+function estimate_lynapov_exponent(t_max, beta, jitter = 0.001)
+	println("Estimating Lynapov Exponent for beta: ", beta)
+
+	N = 32 # number of masses
+	A = 10 # amplitude
+	modes = 3 # number of modes to plot
+	println("t_msx: ", t_max)
+	final_time = t_max # seconds
+	dt = 0.05 # seconds
+	num_steps = Int(final_time / dt)
+
+	t_span = (0.0, final_time)
+	params = (N, modes, beta, A, dt, num_steps)
+
+	# run a simulation with normal parameters
+	is = get_initial_state(N, 1, beta, A)
+	# plot iniital state
+	plot(is, label = "Initial State", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", title = "Initial Positions of Masses", markersize = 4)
+	is_0 = zeros(2 * (N + 2))
+	for i in 1:length(is)
+		is_0[i*2-1] = is[i]
+	end
+	prob_s = ODEProblem(tendency!, is_0, t_span, params)
+	# solve so that the timestep is exactly set to dt
+	s = solve(prob_s, SSPRK33(), dt = dt)
+	states_control = s.u
+	timesteps = s.t
+
+
+	# run a simulation with a small perturbation 
+	# dynamics!(state, prev_state, params, states)
+	# state above is position, need to add in momentum
+	is2 = get_initial_state(N, 1, beta, A) .* (1 - jitter) # perturb the initial state
+	# is2 = get_initial_state(N, 1, beta, A) .- jitter # perturb the initial state
+	plot!(is2, label = "Perturbed Initial State", marker = :cross, xlabel = "Mass Number", ylabel = "Displacement", title = "Perturbed Initial Positions of Masses", msw = 2)
+	savefig("t/new/inits-jitter.png")
+	s_0 = zeros(2 * (N + 2))
+	for i in 1:length(is2)
+		s_0[i*2-1] = is2[i]
+	end
+
+	# add random jitter to the initial state
+	# middle_index = Int(N / 2) * 2
+	# s_0[middle_index] -= jitter # subtract a little from momenetum of middle particle
+
+	final_time = num_steps * dt
+	t_span = (0.0, final_time)
+	prob = ODEProblem(tendency!, s_0, t_span, params)
+	sol = solve(prob, SSPRK33(), dt = dt) # perturbed simulation
+	state_perturbed = sol.u
+	time_perturbed = sol.t
+
+	# println(time_perturbed .- timesteps, "\n")
+
+	println("Finished simulations")
+	println("Control len: ", length(states_control), " Perturbed len: ", length(state_perturbed), " Time len: ", length(time_perturbed))
+
+
+	# # print the difference in the two initial states
+	# println("Initial State Difference (pertrubed - control): ", state_perturbed[1] .- states_control[1])
+	# println("Final State Difference (pertrubed - control): ", state_perturbed[end] .- states_control[end])
+
+	# plot the difference in the two states over time
+	differences = []
+	for i in 1:length(state_perturbed)
+		# calculate mode 1 difference
+		a_k = 0
+		amp = sqrt(2 / (N + 1))
+		for j in 1:N
+			diff = abs(state_perturbed[i][j*2-1] - states_control[i][j*2-1])
+			a_k += diff * sin((1 * (j - 1) * π) / (N + 1))
+		end
+		a_k *= amp
+		push!(differences, a_k)
+	end
+	plot(time_perturbed, differences, label = "Difference in States", xlabel = "Time", ylabel = "Difference", title = "Absolute Difference in States Over Time", lw = 1.5, color = :blue)
+	savefig("t/frames/difference-in-states-beta-$beta.png")
+	# perform a linear regression on the log of the differences
+	# cut the data set into everythign after t=500
+	# find the index where t = 500
+	# index = 1
+	# for i in 1:length(time_perturbed)
+	# 	if time_perturbed[i] >= 500
+	# 		index = i
+	# 		break
+	# 	end
+	# end
+	# differences = differences[index:end]
+	# time_perturbed = time_perturbed[index:end]
+	ln_differeces = log.(differences)
+	(a, b) = linear_regression(time_perturbed, ln_differeces)
+	a = round(a, digits = 4)
+	b = round(b, digits = 4)
+	plot(time_perturbed, differences, label = "ln(difference in mode1)", xlabel = "Time", ylabel = "Difference", title = "Absolute Difference in States Over Time", msw = 0.0, yscale = :ln, color = :blue, lw = 1.5)
+	plot!(time_perturbed, exp.(a * time_perturbed .+ b), label = "exp($a*t + $b)", lw = 2, color = :red, linestyle = :dash)
+	savefig("t/frames/difference-in-states-beta-$beta-log.png")
+
+	println("Saved difference-in-states-beta-$beta.png and difference-in-states-beta-$beta-log.png\n")
+	println("Lynapov Exponent: ", a)
+
+	# check that solve doesn't explode
+	# animate_states(state_perturbed, time_perturbed, 0, 1)
+	return a
+end
+
+# linear regression of the log of the difference
+function linear_regression(x, y)
+	n = length(x)
+	x̄ = sum(x) / n
+	ȳ = sum(y) / n
+	a = sum((x .- x̄) .* (y .- ȳ)) / sum((x .- x̄) .^ 2)
+	b = ȳ - a * x̄
+	return (a, b)
+end
+
+function analyze_equiparition(t_max, beta)
+	# run a simulation with normal parameters
+	N = 32 # number of masses
+	A = 10 # amplitude
+	modes = 3 # number of modes to plot
+	final_time = t_max # seconds
+	dt = 0.05 # seconds
+	num_steps = Int(final_time / dt)
+
+	t_span = (0.0, final_time)
+	params = (N, modes, beta, A, dt, num_steps)
+	s = run_simulation(N, modes, beta, A, dt, num_steps)
+
+	states = s.u
+	timesteps = s.t
+
+	analyze_energies_of_n_modes([1, 3, 5, 7], states, timesteps, beta)
+end
+
+# my_beta = 0.3
+# estimate_lynapov_exponent(10000, my_beta, 0.001)
+# analyze_equiparition(10000, my_beta)
+
+# function plot_betas_versus_lynapov_exponent(betas, t_max, jitter = 0.001)
+# 	lynapovs = []
+# 	for b in betas
+# 		lynapov = estimate_lynapov_exponent(t_max, b, jitter)
+# 		push!(lynapovs, lynapov)
+# 	end
+
+# 	plot(betas, lynapovs, label = "Lynapov Exponent", xlabel = "Beta", ylabel = "Lynapov Exponent", title = "Lynapov Exponent Over Beta")
+# 	savefig("t/lynapov-over-beta.png")
+# end
+
+function plot_probability_distrubtion_of_mode(energy_of_mode, total_energy, timesteps, beta)
+	# get the probability distribution of the mode
+	prob = []
+	for i in 1:length(timesteps)
+		push!(prob, energy_of_mode[i] / total_energy[1]) # energy is conserved
+	end
+
+	plot(timesteps, prob, label = "Probability of Mode 1", xlabel = "Time", ylabel = "Probability", title = "Probability of Mode 1 Over Time with beta=$beta", ylim = (0, 1))
+	savefig("t/probability-mode1-beta-$beta.png")
+end
+
+function analyze_probability_distribution_of_mode(t_max, beta)
+	# run a simulation with normal parameters
+	N = 32 # number of masses
+	A = 10 # amplitude
+	modes = 3 # number of modes to plot
+	final_time = t_max # seconds
+	dt = 0.05 # seconds
+	num_steps = Int(final_time / dt)
+
+	t_span = (0.0, final_time)
+	params = (N, modes, beta, A, dt, num_steps)
+	s = run_simulation(N, modes, beta, A, dt, num_steps)
+
+	states = s.u
+	timesteps = s.t
+
+	# get the energies of the mode
+	total_energy, kinetic_energy, potential_energy, total_by_mode, kinetic_by_mode, potential_by_mode = analyze_energies_of_n_modes([1, 3, 5, 7, 9, 11, 13], states, timesteps, beta)
+
+	# get the probability distribution of the mode
+	# plot_probability_distrubtion_of_mode(total_by_mode[1], total_energy, timesteps, beta)
+
+	# print the time when the probability gets lower than .5
+	for i in 1:length(timesteps)
+		if (total_by_mode[1][i] / total_energy[1]) < 0.6
+			println("TIME when probability of mode 1 is less than 0.6 for beta=$beta: ", timesteps[i])
+			return timesteps[i]
+		end
+	end
+
+	return -1
+end
+
+# estimate_lynapov_exponent(1000, 0.3, 0.001)
+
+function plot_pos_at_time_t(t, beta)
+	# run a simulation with normal parameters
+	N = 32 # number of masses
+	A = 10 # amplitude
+	modes = 3 # number of modes to plot
+	final_time = t # seconds
+	dt = 0.05 # seconds
+	num_steps = Int(final_time / dt)
+
+	t_span = (0.0, final_time)
+	params = (N, modes, beta, A, dt, num_steps)
+	s = run_simulation(N, modes, beta, A, dt, num_steps)
+
+	states = s.u
+	timesteps = s.t
+	positions = s.u[end][1:2:end]
+	inital_positions = s.u[1][1:2:end]
+
+	plot(positions, label = "Positions at t=$t", xlabel = "Mass Number", ylabel = "Displacement", title = "Positions of Masses at t=$t, beta=$beta", lw = 2, marker = :circle)
+	plot!(inital_positions, label = "Initial Positions", lw = 2, marker = :circle)
+	savefig("t/new/positions-at-t-$t-beta-$beta.png")
+end
+
+function get_breakdown_time(beta, t_find = 3600)
+	t = analyze_probability_distribution_of_mode(t_find, beta)
+	if t == -1
+		println("Could not find time for beta: ", beta)
+		return -1
+	end
+	return t
+end
+
+
+# plot_pos_at_time_t(1000, 10)
+# plot_pos_at_time_t(1000, 0)
+
+
+
+# plot_betas_versus_lynapov_exponent(collect(0:1:50), 1000, 0.001)
+
+# analyze_probability_distribution_of_mode(100000, 2.0)
+
+# break_down_times = []
+# bs = collect(1:0.5:50)
+# for b_s in bs
+# 	t_find = 3600
+# 	t = analyze_probability_distribution_of_mode(t_find, b_s)
+# 	if t == -1
+# 		println("Could not find time for beta: ", b_s)
+# 		continue
+# 	end
+# 	push!(break_down_times, t)
+# end
+# # plot beta vs breakdown time
+# plot(bs, break_down_times, label = "Breakdown Time", xlabel = "Beta", ylabel = "Time", title = "Breakdown Time Over Beta")
+# savefig("t/breakdown-time-over-beta.png")
+
+# # plot the log of this
+# one_over_breakdown_time = break_down_times .^ (-1)
+# (a, b) = linear_regression(bs, one_over_breakdown_time)
+# a = round(a, digits = 4)
+# b = round(b, digits = 4)
+# plot(bs, one_over_breakdown_time, label = "1 / (breakdown time)", xlabel = "Beta", ylabel = "1 / (breakdown time)", title = "Breakdown Time ^ -1 Over Beta", msw = 0.0, color = :blue, lw = 1.5)
+# plot!(bs, a * bs .+ b, label = "1 / t = $a * beta + $b", lw = 2, color = :red, linestyle = :dash)
+# savefig("t/ln-breakdown-time-over-beta.png")
+
+# println("Done!")
+# # for a range of betas, find the breakdown time and the lynapov exponent
+function tmp(betas)
+	break_down_times = []
+	lynapovs = []
+	t_find = 350000
+	for b in betas
+		# if b < 0.5
+		# 	t_find = 10000
+		# end
+		t = get_breakdown_time(b, t_find)
+		if t == -1
+			continue
+		end
+		t_find = round(Int, t) + 100
+
+		l = estimate_lynapov_exponent(t_find, b, 0.001)
+		if l == 0
+			println("Lynapov for beta was zero: ", b)
+			continue
+		end
+		push!(break_down_times, t)
+		push!(lynapovs, l)
+	end
+
+	# plot lynapov vs breakdown time
+	scatter(break_down_times, lynapovs, label = "Lynapov Exponent", xlabel = "Breakdown Time", ylabel = "Lynapov Exponent", title = "Lynapov Exponent Over Breakdown Time")
+	savefig("t/new/new-lynapov-over-breakdown-time.png")
+
+	# plot this over log of lynapovs
+	ln_lynapovs = log.(lynapovs)
+	ln_break_down_times = log.(break_down_times)
+	(a, b) = linear_regression(ln_break_down_times, ln_lynapovs)
+	a = round(a, digits = 4)
+	b = round(b, digits = 4)
+	scatter(ln_break_down_times, ln_lynapovs, label = "ln(Lynapov Exponent)", xlabel = "ln(Breakdown Time)", ylabel = "ln(Lynapov Exponent)", title = "ln(Lynapov Exponent) Over ln(Breakdown Time)", msw = 0.0, color = :blue, lw = 1.5)
+	plot!(ln_break_down_times, a * ln_break_down_times .+ b, label = "ln(beta) = $a * ln(t) + $b", lw = 2, color = :red, linestyle = :dash)
+	savefig("t/new/new-ln-lynapov-over-breakdown-time.png")
+
+	# only one log
+	(a, b) = linear_regression(lynapovs, ln_break_down_times)
+	a = round(a, digits = 4)
+	b = round(b, digits = 4)
+	scatter(lynapovs, ln_break_down_times, label = "ln(Breakdown Time)", xlabel = "Lynapov Exponent", ylabel = "ln(Breakdown Time)", title = "ln(Breakdown Time) Over Lynapov Exponent", msw = 0.0, color = :blue, lw = 1.5)
+	plot!(lynapovs, a * lynapovs .+ b, label = "ln(t) = $a * beta + $b", lw = 2, color = :red, linestyle = :dash)
+	savefig("t/new/new-new-ln-breakdown-time-over-lynapov.png")
+
+
+
+	# plot the lynapoovs to the ^-1
+	one_over_breakdown_time = break_down_times .^ (-1)
+	(a, b) = linear_regression(lynapovs, one_over_breakdown_time)
+	a = round(a, digits = 4)
+	b = round(b, digits = 4)
+	scatter(lynapovs, one_over_breakdown_time, label = "1 / (breakdown time)", xlabel = "Lynapov Exponent", ylabel = "1 / (breakdown time)", title = "1 / (breakdown time) Over Lynapov Exponent", msw = 0.0, color = :blue, lw = 1.5)
+	plot!(lynapovs, a * lynapovs .+ b, label = "1 / t = $a * beta + $b", lw = 2, color = :red, linestyle = :dash)
+	savefig("t/new/new-one-over-lynapov-over-breakdown-time.png")
+
+	# plot the lynapov over betas
+	scatter(betas, lynapovs, label = "Lynapov Exponent", xlabel = "Beta", ylabel = "Lynapov Exponent", title = "Lynapov Exponent Over Beta")
+	savefig("t/new/new-lynapov-over-beta.png")
+end
+
+tmp(collect(0.3:0.025:10))
+
+# estimate_lynapov_exponent(500000, 0.3, 0.001)