# This is my first attempt at a Billiards Simulator.  Later versions will be more user-friendly.
# Inputting the number of sides " N ", initial point " x0 " and angle " theta ", 
# the program generates mathematica commands that draw the picture.
import math, cmath
	
def lineA(x,y):
	return [y[1]-x[1],x[0]-y[0],y[1]*x[0]-x[1]*y[0]]

def lineB(x, theta):
	return [-1*math.sin(theta),math.cos(theta),-1*x[0]*math.sin(theta) + x[1]*math.cos(theta)]

def intersect(L1,L2):
	#print "L1:", L1
	#print "L2:", L2
	det = L1[0]*L2[1] - L1[1]*L2[0]
	if(det):
		return [(L1[2]*L2[1] - L1[1]*L2[2])/det, (L1[0]*L2[2] - L1[2]*L2[0])/det]
	else:
		return None

def step(x0, theta, N):
	for x in range(0,N):
		pt = intersect( lineA(vertices[x],vertices[x+1] ),lineB(x0, theta) )
		print x, "check"
		# how to pick out the correct side?
		frac = [(pt[0] - vertices[x+1][0])/(vertices[x][0] - vertices[x+1][0]), (pt[1] - vertices[x+1][1])/(vertices[x][1] - vertices[x+1][1])]
		t = [(pt[0] - x0[0])/math.cos(theta), (pt[1] - x0[1])/math.sin(theta)]
		if( frac[0] > 0.0001 and frac[0] < 0.9999 and t[0] > 0.0001):
			out = "{{ {0},{1} }}, \n".format(pt[0], pt[1])
			f.write(out)
			return [pt, 2*(2*math.pi*(x+0.5)/N) + math.pi - theta]
		if( frac[1] > 0.0001 and frac[1] < 0.9999 and t[1] > 0.0001):
			out = "{{ {0},{1}  }}, \n".format(pt[0], pt[1])
			f.write(out)
			return [pt,  2*(2*math.pi*(x+0.5)/N) + math.pi - theta]

f = open("billiardData.txt", "w")

# number of sides
N = 11
vertices = [[cmath.exp( 2*cmath.pi*1j*k/N).real, cmath.exp( 2*cmath.pi*1j*k/N).imag] for k in range(0,N+1) ]

out = "AA = ListPlot[{"
for x in vertices[:-1]:
	if abs(x[0]) < 0.0001: x[0] = 0
	if abs(x[1]) < 0.0001: x[1] = 0
	out += "{{{0}, {1}}},".format(x[0],x[1])
if abs(vertices[-1][0]) < 0.0001: vertices[-1][0] = 0
if abs(vertices[-1][1]) < 0.0001: vertices[-1][1] = 0
out += "{{{0}, {1}}}}}, Joined -> True, AspectRatio -> 1, Axes -> False];\n".format(vertices[-1][0],vertices[-1][1])

f.write(out)

f.write("BB = ListPlot[{")
x0 = [0, 0.25]
theta = math.pi/4
z = step(x0, theta, N)
try:
	for x in range(0,10):
		z = step(z[0] ,z[1], N )
except:
	pass
f.write("}, PlotStyle -> PointSize[Medium], Joined -> True, Axes -> False];\n")
f.write("Show[AA,BB,AspectRatio-> 1]")