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)