sigma_x = sum(x);
sigma_y = sum(y);
sigma_z = sum(z);

sigma_x2 = sum(x.^2);
sigma_y2 = sum(y.^2);
sigma_z2 = sum(z.^2);

sigma_x3 = sum(x.^3);
sigma_y3 = sum(y.^3);
sigma_z3 = sum(z.^3);

sigma_x4 = sum(x.^4);
sigma_y4 = sum(y.^4);
sigma_z4 = sum(z.^4);

sigma_x_y = dot(x,y);
sigma_x_z = dot(x,z);
sigma_y_z = dot(y,z);

sigma_x2_y2 = dot(x.^2,y.^2);
sigma_x2_z2 = dot(x.^2,z.^2);
sigma_y2_z2 = dot(y.^2,z.^2);

sigma_x2_y = dot(x.^2,y);
sigma_x2_z = dot(x.^2,z);

sigma_y2_z = dot(y.^2,z);
sigma_y2_x = dot(y.^2,x);

sigma_z2_x = dot(z.^2,x);
sigma_z2_y = dot(z.^2,y);

jacov = [sigma_x4 sigma_x2_y2 sigma_x2_z2 sigma_x3 sigma_x2_y sigma_x2_z ;
     sigma_x2_y2 sigma_y4 sigma_y2_z2 sigma_y2_x sigma_y3 sigma_y2_z ; 
     sigma_x2_z2 sigma_y2_z2 sigma_z4 sigma_z2_x sigma_z2_y sigma_z3 ;
     sigma_x3 sigma_y2_x sigma_z2_x sigma_x2 sigma_x_y sigma_x_z ;
     sigma_x2_y sigma_y3 sigma_z2_y sigma_x_y sigma_y2 sigma_y_z ;
     sigma_x2_z sigma_y2_z sigma_z3 sigma_x_z sigma_y_z sigma_z2];

b = [sigma_x2; sigma_y2; sigma_z2; sigma_x; sigma_y; sigma_z];

kai = inv(jacov) * b;

f = @(X,Y,Z) kai(1) * X.^2 + kai(2) * Y.^2 + kai(3) * Z.^2 + kai(4) * X + kai(5) * Y + kai(6) * Z - 1;
% fimplicit3(f);

cen_left = [2 * kai(1) 0 0; 0 2 * kai(2) 0; 0 0 2 * kai(3)];
cen_right = [-kai(4); -kai(5); -kai(6)];
a
xc = -kai(4) / (2 * kai(1));
yc = -kai(5) / (2 * kai(2));
zc = -kai(6) / (2 * kai(3));

Z = 1 + (kai(4)^2 / kai(1) + kai(5)^2 / kai(2) + kai(6)^2 / kai(3)) / 4;

xr = sqrt(Z / kai(1));
yr = sqrt(Z / kai(2));
zr = sqrt(Z / kai(3));

[x_ellip,y_ellip,z_ellip] = ellipsoid(xc,yc,zc,xr,yr,zr);
figure
surf(x_ellip,y_ellip,z_ellip);
hold on
scatter3(x,y,z)
scatter3(x - xc,(y - yc) * xr / yr, (z - zc) * xr / zr);

\begin{bmatrix}
\phi \\
\theta \\
\psi
\end{bmatrix}
=
\begin{bmatrix}
\tan^{-1}\frac{a_y}{a_z} \\
-\tan^{-1}\frac{a_x}{\sqrt{a_y^2+a_z^2}} \\
\tan^{-1}(- \frac{\cos\phi d_y - \sin\phi d_z}{\cos\theta d_x + \sin\theta \sin\phi d_y + \sin\theta \cos\phi d_z})
\end{bmatrix}