% Coment out the comet not to compute
load('./C-2017-K2.mat') % First comet
% load('./P-2011-NO1.mat') % Second comet

x = xy(:,1);
y = xy(:,2);

X = [x.^2, x.*y, y.^2, x, y];
Y = -ones(size(X, 1), 1);

theta = X.' * X \ X.' * Y

eccentricity(theta(1), theta(2), theta(3), theta(4), theta(5), 1)

syms xp yp
fimplicit(theta(1)*xp^2 + theta(2)*xp*yp + theta(3)*yp^2 + theta(4)*xp + theta(5)*yp == 1)
hold on

scatter(x, y, 'r');