update naming

1 files changed, 0 insertions, 767 deletions

diff --git a/t/plz.jl b/t/plz.jl
deleted file mode 100644
index 962398a..0000000
--- a/t/plz.jl
+++ /dev/null
@@ -1,767 +0,0 @@
-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)