%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%     ELABORATO 2       %%%
%%% Geometria delle masse %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear all
format short g

% matricola
MATRICOLA=[5 1 7 1 1 3 0]

SUM=sum(MATRICOLA)           % somma le cifre della matricola

% dati
muAl=1.00e-6                 %(Kg/mm^2) densità Alluminio (piattine)
muFe=3.28e-6                 %(Kg/mm^2) densità Ferro (struttura a Z)
s=mod(SUM,3)+1               % spessore piattine (mm)
b=2*(mod(SUM,7)+7)           % lunghezza base (mm)
d=(mod(SUM,4)+1)*s           % spessore centrale (mm)
r=2*s                        % raggio di curvatura dei raccordi (mm)
h=2*(b+d+s)                  % altezza profilato (mm)  

% calcolo delle masse (consultare figure elaborato)
m(1)=muFe/4*pi*r^2;
m(2)=muFe*(b-r+d)*r;
m(3)=muAl*s*(h-2*r);
m(4)=muFe*d*(h-2*r);
m(5)=m(3);
m(6)=m(2);
m(7)=m(1);
m(8)=muFe*r^2;
m(9)=-m(1);                % densità negativa
m(10)=m(8);
m(11)=-m(1);               % densità negativa

% massa intera struttura
M=0;                       % inizializzo la variabile M per il ciclo
for i=1:11
    M=M+m(i);
end
M                          % mostra il valore della massa totale

% coordinate baricentri
G_set=4*sqrt(2)*r/(pi*3)   % baricentro del settore circolare
                           % sull'asse di simmetria del medesimo
xG(1)=r-G_set/sqrt(2);
yG(1)=h-G_set/sqrt(2);

xG(2)=r+(b+d-r)/2;
yG(2)=h-r/2;

xG(3)=b-s/2;
yG(3)=(h-2*r)/2;

xG(4)=b+d/2;
yG(4)=r+(h-2*r)/2;

xG(5)=b+d+s/2;
yG(5)=(h-2*r)/2+2*r;

xG(6)=(d+b-r)/2+b;
yG(6)=r/2;

xG(7)=2*b+d-r+G_set/sqrt(2);
yG(7)=G_set/sqrt(2);

xG(8)=b-r/2;
yG(8)=h-3*r/2;

xG(9)=b-r+G_set/sqrt(2);
yG(9)=h-2*r+G_set/sqrt(2);

xG(10)=b+d+r/2;
yG(10)=3*r/2;

xG(11)=b+d+r-G_set/sqrt(2);
yG(11)=2*r-G_set/sqrt(2);

xG
yG

% momenti statici
Sx=0;
Sy=0;
for i=1:11
    Sx=Sx+m(i)*yG(i);
    Sy=Sy+m(i)*xG(i);
end
Sx
Sy

%baricentro struttura
XG=Sy/M;
YG=Sx/M;
G=[XG YG]

% momenti e prodotti d'inerzia delle singole figure
Iro_set=muFe*r^4/4*(pi/4-1/2)          
% ^ momento d'inerzia del settore circolare sull'asse di simmetria ^   

Is_set=muFe*r^4/4*(pi/4+1/2)                      
% ^ momento d'inerzia del settore circolare ^
% ^ sull'asse ortogonale a quello di simmetria ^

Iso_set=Is_set-m(1)*G_set^2                        
% ^ teorema di Huyghens (trasporto del momento sull'asse baricentrale)^

Ixo(1)=Iro_set*(cos(-pi/4))^2+Iso_set*(sin(-pi/4))^2;
% ^ dal teorema delle rotazioni ^

Iyo(1)=Ixo(1);
% Iyo(1)=Iso_set*(cos(-pi/4))^2+Iro_set*(sin(-pi/4))^2

Jxoyo(1)=(Iro_set-Iso_set)/2;
% ^ teorema delle rotazioni del prodotto d'inerzia ^

Ixo(2)=muFe*(b-r+d)*r^3/12;
Iyo(2)=muFe*(b-r+d)^3*r/12;
Jxoyo(2)=0;

Ixo(3)=muAl*s*(h-2*r)^3/12;
Iyo(3)=muAl*s^3*(h-2*r)/12;
Jxoyo(3)=0;

Ixo(4)=muFe*d*(h-2*r)^3/12;
Iyo(4)=muFe*d^3*(h-2*r)/12;
Jxoyo(4)=0;

Ixo(5)=Ixo(3);
Iyo(5)=Iyo(3);
Jxoyo(5)=0;

Ixo(6)=Ixo(2);
Iyo(6)=Iyo(2);
Jxoyo(6)=0;

Ixo(7)=Ixo(1);
Iyo(7)=Iyo(1);
Jxoyo(7)=Jxoyo(1);

% per queste ultime 4 figure si esegue 'quadrato (meno) settore circolare'
Ixo(8)=muFe*r^4/12;        
Iyo(8)=Ixo(8);
Jxoyo(8)=0;

Ixo(9)=-Ixo(1);
Iyo(9)=-Iyo(1);
Jxoyo(9)=-Jxoyo(1);

Ixo(10)=Ixo(8);
Iyo(10)=Iyo(8);
Jxoyo(10)=0;

Ixo(11)=-Ixo(1);
Iyo(11)=-Iyo(1);
Jxoyo(11)=-Jxoyo(1);

Ixo
Iyo
Jxoyo

% applico un isomorfismo e passo in coordinate baricentriche
xo=0;
yo=0;
for i=1:11
    xo(i)=xG(i)-XG;
    yo(i)=yG(i)-YG;
end
xo
yo

% momenti e prodotti d'inerzia della struttura rispetto agli
% assi baricentrali (uso ripetutamente il teorema di Huyghens)
for i=1:11
    IXo(i)=Ixo(i)+m(i)*yo(i)^2;
    IYo(i)=Iyo(i)+m(i)*xo(i)^2;
    JXoYo(i)=Jxoyo(i)+m(i)*xo(i)*yo(i);
end
IXo
IYo
JXoYo

IxG=sum(IXo)
IyG=sum(IYo)
JxGyG=sum(JXoYo)

% ricerca delle direzioni centrali d'inerzia
fistar=0.5*atan(2*JxGyG/(IyG-IxG))
fistar_deg=fistar*180/pi

% momenti rispetto alle direzioni centrali d'inerzia
Icsi=IxG*(cos(fistar))^2+IyG*(sin(fistar))^2-JxGyG*sin(2*fistar)
Ieta=IxG*(sin(fistar))^2+IyG*(cos(fistar))^2+JxGyG*sin(2*fistar)

% raggi d'inerzia
rhocsi=sqrt((Icsi)/M)
rhoeta=sqrt((Ieta)/M)

% calcolo del momento d'inerzia della struttura rispetto alla retta
% baricentrale r0 inclinata di 45° rispetto all'orizzontale 
Iro=IxG*(cos(pi/4))^2+IyG*(sin(pi/4))^2-JxGyG*sin(pi/2)
rho=sqrt(Iro/M)

% tensore d'inerzia
ve=[1/sqrt(2) 1/sqrt(2) 0];                     % versore della retta r0
e=ve'
yo=[IxG -JxGyG 0; -JxGyG IyG 0; 0 0 IxG+IyG]         % matrice d'inerzia
yoe=yo*e                         % applicazione tensore yo al versore 'e'
Irot=dot(yoe,e)

%% FINE %%