%% GAUSSIANA 3D
close all;
clear all;
clc;


u = dlmread('samples_50.txt');
u = u(:,2:end);

lambda_1 = 0.3e-03;
lambda_2 = 0.3e-03;
lambda_3 = 1.7e-03;

angle = [0,0];

D = create_diffusion_tensor([lambda_1,lambda_2,lambda_3],angle);

l = length(u);

E = zeros(l,1);
b = 1000;
for i=1:l
    E(i) = exp(-b* u(i,:)*D*u(i,:)'); 
end

draw_elipsoid(D)


A = zeros(length(u), 6);
for i=1:length(u)
    A(i,:) = [u(i,1)^2, 2*u(i,1)*u(i,2), 2*u(i,1)*u(i,3), 2*u(i,2)*u(i,3), u(i,2)^2, u(i,3)^2]; 
end

c = pinv(A'*A)*A'*(-log(E)/b);

D_res= [[c(1),c(2), c(3)];[c(2),c(5),c(4)];[c(3),c(4),c(6)]];


draw_elipsoid(D_res)

[U,L]=eig(D_res);
[m,i] = max(diag(L));
vec = U(:,i);

angle = deg2rad(angle);
[v(1),v(2),v(3)] = sph2cart(angle(1),angle(2),1);

v*vec

[th,ph,r] = cart2sph(vec(1),vec(2),vec(3));

angle2 = [rad2deg(th), rad2deg(ph)]