#!/usr/bin/env python

# NOTE: This script uses 'matplotlib' and 'scipy', which can be installed via:
#   python3 -m pip install --user matplotlib scipy
# However, the script will attempt to install them automatically if they are not present.

from math import hypot, sqrt, pi

def derivatives_1_and_2(f, coords, h, i, j):
	coord = coords[i]
	x = coord[j]
	coord[j] = x + h; fp = f(coords)
	coord[j] = x - h; fn = f(coords)
	coord[j] = x    ; f0 = f(coords)
	deriv1 = (fp - fn) / (2 * h)
	deriv2 = ((fp - f0) - (f0 - fn)) / (h * h)
	return (deriv1, deriv2)

def create_newton_minimizer(f, coords, h):  # 'h' is the step size for differentiation: (f(x + h) - f(x - h)) / (2 * h)
	fixed_boundary = 1
	n = len(coords) - 2 * fixed_boundary
	m = len(coords[0])
	i_range = range(fixed_boundary, n + fixed_boundary)
	j_range = range(m)
	def iterate():
		delta = 0.0
		for i in i_range:
			coord = coords[i]
			for j in j_range:
				(deriv1, deriv2) = derivatives_1_and_2(f, coords, h, i, j)
				if deriv2 > h:
					step = -deriv1 / deriv2
				else:
					step = -deriv1 / n
				delta += step if step >= 0 else -step
				coord[j] += step
		return delta
	return iterate

def quadratic_roots(two_a, b, c):
	if not two_a:
		roots = [-c / b]
	else:
		D = b * b - 2 * two_a * c
		if D >= 0:
			left = -b
			if D > 0:
				right = sqrt(D)
				roots = [(left + right) / two_a, (left - right) / two_a]
			else:
				roots = [left / two_a]
		else:
			roots = []
	return roots

def falling_time(coords):
	gravity = -9.80665
	speed = 0.0
	t = 0
	prev = None
	for coord in coords:
		if prev is not None:
			# Assume we are always on the ground (NOT a safe assumption!)
			dx = coord[0] - prev[0]
			dy = coord[1] - prev[1]
			d = hypot(dx, dy)
			acceleration = gravity * dy / d
			for dt in quadratic_roots(acceleration, speed, -d):
				if dt > 0:
					speed += acceleration * dt
					t += dt
		prev = coord
	return t

def split_every(items, n):  # splits the given list every n items
	result = []
	current_list = None
	for item in items:
		if current_list is None:
			current_list = []
		current_list.append(item)
		if len(current_list) >= n:
			result.append(current_list)
			current_list = None
	if current_list is not None:
		result.append(current_list)
	return result

def generate_via_callback(f, capacity=1):
	from Queue import Queue
	from threading import Thread
	def target(f, finished, queue):
		try: return f(queue)
		finally: queue.put(finished); queue.join()
	queue = Queue(capacity)
	finished = object()
	thd = Thread(target=target, args=(f, finished, queue))
	thd.setDaemon(True)
	thd.start()
	def generator(queue, finished):
		while 1:
			item = queue.get()
			try:
				if item is finished:
					break
				yield item
			finally:
				queue.task_done()
	return generator(queue, finished)

def import_or_install_module(name, *args):  # imports a module, installing it if it doesn't exist
	for k in range(3):
		try:
			if k >= 2: import ensurepip; ensurepip._main()
			if k >= 1:
				import pip
				pip_main = getattr(pip, 'main', None)
				if pip_main is None: from pip.__main__ import _main as pip_main
				pip_main(['install', '--user'] + list(args) + [name])
			if k >= 0: return __import__(name)
			break
		except ImportError: pass

def create_scipy_minimizer(f, coords, h):
	fixed_boundary = 1
	from itertools import chain
	import_or_install_module('scipy')
	import scipy.optimize
	generator = generate_via_callback(lambda queue: scipy.optimize.minimize(
		lambda guess: f(coords[:+fixed_boundary] + split_every(guess, len(coords[0])) + coords[-fixed_boundary:]),
		# Join all x & y coordinates into one long vector (x1, y1, x2, y2, ...)
		list(chain.from_iterable(coords[+fixed_boundary:-fixed_boundary])),
		tol=4 * h ** 0.5,
		callback=queue.put))
	# Make sure imports happen before solver is returned
	it = iter(generator)
	def iterate():
		delta = 0.0
		try:
			guess = next(it)
			old_coords = coords[:]
			coords[+fixed_boundary:-fixed_boundary] = split_every(guess, len(coords[0]))
			delta += sum(map(lambda a, b: abs(a - b), chain.from_iterable(coords), chain.from_iterable(old_coords)))
		except StopIteration: pass
		return delta
	return iterate

def create_step_generator(stepper):
	while 1:
		delta = stepper()
		if not delta: break
		yield delta

def main(program, n=30, altitude=1.0, h=0.000001):
	n = int(n)
	h = float(h)
	altitude = float(altitude)
	visualize = True
	use_scipy_minimizer = False
	radius = 1.0
	slope = -altitude / radius
	(y1, y2) = (altitude, 0.0)
	(x1, x2) = (0.0, radius)
	coords = list(map(lambda i: [
		x1 + (x2 - x1) * i / n,
		y1 + (y2 - y1) * i / n
	], range(n + 1)))
	if use_scipy_minimizer:
		create_minimizer = create_scipy_minimizer
	else:
		create_minimizer = create_newton_minimizer
	stepper = create_minimizer(falling_time, coords, h)
	step_generator = create_step_generator(stepper)
	if visualize:
		import_or_install_module('matplotlib')
		import matplotlib.pyplot, matplotlib.animation
		# matplotlib.pyplot.rc('font', family='sans-serif', **{'sans-serif': 'Latin Modern Sans'})
		fig = matplotlib.pyplot.figure(frameon=False)
		ax = fig.add_subplot(1, 1, 1)
		ax.grid(True, which='both')
		ax.margins(x=0, y=0)
		artists = ax.plot(*zip(*coords), color='blue')
		fig.tight_layout()
		blit = True  # faster drawing, but doesn't move around plot or redraw axes labels
		def update_plot(delta):
			if not blit:
				ax.relim()
				ax.autoscale_view(True, True, True)
			artists[0].set_data(zip(*coords))
			if not blit:
				fig.canvas.draw()
			return artists
		matplotlib.pyplot.show(matplotlib.animation.FuncAnimation(fig, update_plot, step_generator, repeat=False, blit=blit, interval=0))
	else:
		from timeit import default_timer as timer
		tstart = timer()
		iterations = len(list(step_generator))  # iterates through steps too
		tend = timer()
		print("%.0f ms = %d iterations * %.2f ms/iteration" % (1000 * (tend - tstart), iterations, 1000 * (tend - tstart) / iterations))
		print(coords)

if __name__ == '__main__':
	import sys
	raise SystemExit(main(*sys.argv))