#!/usr/bin/python




import csv
import sys
import numpy as np
from numpy.linalg import eig, inv

def fitEllipse(x,y):
	x = x[:,np.newaxis]
	y = y[:,np.newaxis]
	D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
	S = np.dot(D.T,D)
	C = np.zeros([6,6])
	C[0,2] = C[2,0] = 2; C[1,1] = -1
	E, V =  eig(np.dot(inv(S), C))
	n = np.argmax(np.abs(E))
	a = V[:,n]
	return a


def ellipse_center(a):
	b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
	num = b*b-a*c
	x0=(c*d-b*f)/num
	y0=(a*f-b*d)/num
	return np.array([x0,y0])

def ellipse_angle_of_rotation( a ):
	b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
	return 0.5*np.arctan(2*b/(a-c))
    
    
    
def ellipse_angle_of_rotation2( a ):
	b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
	return (0.5*np.arctan(2*b/(a-c))) / np.pi * 180 % 180
    
def ellipse_angle_degrees(phi):
	yang = np.cos(phi)
	xang = np.sin(phi)
	
	return np.arctan2(yang ,xang ) * 180 / np.pi
    
def ellipse_axis_length( a ):
	b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
	up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
	down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
	down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
	res1=np.sqrt(up/down1)
	res2=np.sqrt(up/down2)
	return np.array([res1, res2])


arc = 2
R = np.arange(0,arc*np.pi, 0.01)
#x = 1.5*np.cos(R) + 2 + 0.1*np.random.rand(len(R))
#y = np.sin(R) + 1. + 0.1*np.random.rand(len(R))


import numpy as np
d = np.loadtxt(sys.argv[1], delimiter="\t")
x =  d[:,sys.argv[2]]
y =  d[:,sys.argv[3]]  


a = fitEllipse(x,y)
center = ellipse_center(a)
phi = ellipse_angle_of_rotation(a)
angle = ellipse_angle_of_rotation2(a) 
axes = ellipse_axis_length(a)


if axes[1] > axes[0]:
	angle = angle - 90
	ellipticity = 1 - (axes[0]/axes[1])

else:
	ellipticity = 1 - (axes[1]/axes[0])

print ellipticity, angle % 180


print "center = ", center
print "angle of rotation = ", angle % 180
print "axes = ", axes

a, b = axes
xx = center[0] + a*np.cos(R)*np.cos(phi) - b*np.sin(R)*np.sin(phi)
yy = center[1] + a*np.cos(R)*np.sin(phi) + b*np.sin(R)*np.cos(phi)

from pylab import *
import matplotlib.pyplot as plt
plt.plot(x,y, 'o')
plot(xx,yy, color = 'red')
axis('equal')
show()
