%Eccentricity Calculator
function ecc=eccentricity(Yt)
l=0;
m=0;
n=0;

for i=1:8000
    %Io orbits
    if mod(i,1000)==0
        for a=1:1000
            dist_io(a)=norm(Yt(l+a+1,1:2));
            raio=max(dist_io);%norm(Yt(l+a+1,1:2)));
            rpio=min(dist_io);%norm(Yt(l+a+1,1:2)));
            io(i/1000)=(raio-rpio)/(raio+rpio);
        end
        l=l+1000;
        dist_io=0;
    end
    %Eu Orbits
    if mod(i,2000)==0
        for b=1:2000
            dist_eu(b)=norm(Yt(m+b+1,5:6));
            raeu=max(dist_eu);
            rpeu=min(dist_eu);
            eu(i/2000)=(raeu-rpeu)/(raeu+rpeu);
        end
        m=m+2000;
        dist_eu=0;
    end
    %Gy Orbits
    if mod(i,4000)==0
        for c=1:4000
            dist_gy(c)=norm(Yt(n+c+1,9:10));
            ragy=max(dist_gy);
            rpgy=min(dist_gy);
            gy(i/4000)=(ragy-rpgy)/(ragy+rpgy);
        end
        n=n+4000;
        dist_gy=0;
    end
end
%Building ecc
temp=[io,eu,gy];
ecc=zeros(3,8);
for i=1:14
    if i<=8
        ecc(1,i)=temp(i);
    elseif i>8 && i<=12
        x=2*(i-8);
        ecc(2,x-1)=temp(i);
    else
        z=4*(i-12);
        ecc(3,z-3)=temp(i);
    end
end