clear all;
close all;
clc;
clf;
%načtení s datového souboru
%fid=fopen('posuvy_5_2.dat');
    
%css1=fscanf(fid,'%g %g %g %g %g %g %g %g',[8 inf]);
css1=load('posuvy_5_2_20_10_k4_zbyt.dat');
%css1=css1';
%fclose(fid);

%css1(:,5:7)=uxyz(:,2:4);

poc_n=length(css1);      %počet uzlů 


%********************************************************
css=sortrows(css1,2);    %seřadí podle souřadnice x

css(:,1)=css(:,2);       %redukování o čísla uzlů
css(:,2)=css(:,3);
css(:,3)=css(:,4);
css(:,4)=css(:,5);
css(:,5)=css(:,6);
css(:,6)=css(:,7);
css(:,7)=0;

css(:,1:3)=css(:,1:3)+css(:,4:6);          %s posuvama

cs=css;
clear css;
n_bodu=length(cs);       %celkový počet uzlů   



%nalezení těžiště křivky
teziste=sum(cs,1)/n_bodu;
xT=teziste(1);
yT=teziste(2);
zT=teziste(3);

%vykreslení težiště křivky do grafu
plot3(cs(:,1),cs(:,2),cs(:,3));
hold on
plot3(xT,yT,zT,'*');
axis ([-3000 3000 -3000 3000 -3000 3000]);

%translace těžiště křivky do poč. souřad. soustavy
Tcs(:,1)=cs(:,1)-xT;
Tcs(:,2)=cs(:,2)-yT;
Tcs(:,3)=cs(:,3)-zT;

plot3(Tcs(:,1),Tcs(:,2),Tcs(:,3),'r');
plot3(0,0,0,'*r');
xlabel('x');
ylabel('y');
zlabel('z');


%sestavení tenzoru setrvačnosti
Jxx=0;
Jyy=0;
Jzz=0;
Dxy=0;
Dyz=0;
Dzx=0;
for i=1:n_bodu
Jxx=Jxx+Tcs(i,2)^2+Tcs(i,3)^2;
Jyy=Jyy+Tcs(i,1)^2+Tcs(i,3)^2;
Jzz=Jzz+Tcs(i,1)^2+Tcs(i,2)^2;
Dxy=Dxy+Tcs(i,1)*Tcs(i,2);
Dyz=Dyz+Tcs(i,2)*Tcs(i,3);
Dzx=Dzx+Tcs(i,3)*Tcs(i,1);
end

TS=[Jxx -Dxy -Dzx; -Dxy Jyy -Dyz; -Dzx -Dyz Jzz]; %tenzor setrvačnosti


[smer_kos,hl_m]=eig(TS);  %nalezne vlastní vektory a vlastní čísla (hl. směry hl. momenty)

for j=1:n_bodu
    HRcs(:,j)=smer_kos^-1*[Tcs(j,1);Tcs(j,2);Tcs(j,3)];
    Hcs=HRcs';
end;


figure
plot3(Hcs(:,1),Hcs(:,2),Hcs(:,3),'black');
xlabel('x');
ylabel('y');
zlabel('z');
axis([-500,500,-5,5,-5,5]);
hold on
%plot3(Tcs(:,1),Tcs(:,2),Tcs(:,3),'b');
obde(1:3)=min(Hcs);
obde(4:6)=max(Hcs);

p_hran=4; %počet hran kvádru

x=linspace(obde(1),obde(4),p_hran);
y=linspace(obde(2),obde(5),p_hran);
z=linspace(obde(3),obde(6),p_hran);
for i=1:p_hran
    xx(i,:)=x(:);
    yy(i,:)=y(:);
    zz(i,:)=z(:);
end
minx(1:p_hran)=obde(1);
plot3(minx,y,zz,'r');
maxx(1:p_hran)=obde(4);
plot3(maxx,y,zz,'r');
miny(1:p_hran)=obde(2);
plot3(x,miny,zz,'g');            
maxy(1:p_hran)=obde(5);
plot3(x,maxy,zz,'g');           
minz(1:p_hran)=obde(3);
plot3(xx,y,minz,'b');
maxz(1:p_hran)=obde(6);
plot3(xx,y,maxz,'b');

%souřadný systém
sx1=[0,200];
sy1=[0,0];
sz1=[0,0];
plot3(sx1,sy1,sz1,'g');
sx2=[0,0];
sy2=[0,-0.2];
sz2=[0,0];
plot3(sx2,sy2,sz2);
sx3=[0,0];
sy3=[0,0];
sz3=[0,0.2];
plot3(sx3,sy3,sz3,'r');

%krivost
krivost(1,1)=abs(maxy(1))+abs(miny(1));
krivost(1,2)=abs(maxz(1))+abs(minz(1));
krivost(1,:)

