%% DIMENSIONAMENTO STRUTTURALE
clear all
close all
clc

L=1.700; %m %lunghezza travi centrali
EG=0.625; %m %distanza
EP_f=0.43; %m %distanza fissa imposta
AAa=0.45; %m distanza tra la cerniera A della base e la cerniera A' (Aa) dell'attuatore %arbitraria
m=200; %kg %massa
g=10; %m/s^2 %accelerazione di gravità
P=m*g; %N %peso
H_min_eff=input('Immettere valore di H_min effettivo, usato nel foglio "Dimensionamento_strutturale_1" (in metri): ');
%H_min_eff=0.3
EP_f=input('Immettere valore di EP effettivo, usato nel foglio "Dimensionamento_strutturale_1" (in metri): ');
%EP_f=0.43
H_f=H_min_eff; %m %altezza fissa imposta

d_min_L1=input('Immettere valore diametro esterno boccola su albero spinta pistone (in millimetri): ');
%d_min_L1=42 mm
d_min_L2=input('Immettere valore diametro estreno boccola su cerniera interna pantografo (in millimetri): ');
%d_min_L2=42 mm

%s_b=2.5; %mm %spessore boccola su braccio

Rp02=235; %MPa %carico di snervamento
CS=1.5; %coefficente di sicurezza
CS_m=2.5; %coefficente di sicurezza maggiorato dato che in un primo momento
%consideriamo solo il momento flettente e non lo sforzo di taglio
r=2; %rapporto geometrio h/b
s_t=3; %mm %spessore tubolare

%% PRIMO TENTATIVO (H=COST)
%Costanti costruttive (variabili manualmente)
EP=EP_f; %m %distanza

%Distanze geometriche derivanti
PC=L/2-EP; %m %distanza
PD=L-EP; %m %distanza

%Variabili operative
H=H_f; %m %altezza tra i vincoli
AD=sqrt(L^2-H^2); %m %base
alpha=atan(H/AD); %radianti
%alpha_gradi=alpha*180/pi;

alpha_a=atan(H/AAa);
alpha_aa=pi/2-alpha_a;
AaE=sqrt(H^2+AAa^2);
gamma=pi-alpha-alpha_a;
AaP=sqrt(EP^2+AaE^2-2*EP*AaE*cos(gamma));
%AP=sqrt(H^2+EP^2-2*H*EP*sin(alpha)); %m %lunghezza pistone
gamma_aa=asin((EP/AaP)*sin(gamma));
beta=pi/2-alpha_aa-gamma_aa;

%Reazioni vincolari
%Vd=P*EG/AD;
Vb=P*EG/AD;
Ve=P-Vb;
Oe=0;
%Oc=(Ve*EP-2*Vb*PC+Vd*PD)/(2*PC*tan(alpha));
Oc=(2*Vb*((PC/PD)-1)-Ve*(1+(EP/PD)))/(tan(alpha)*(1-2*(PC/PD))-tan(beta));
Vc=-2*Vb-Oc*tan(alpha);
Oa=-Oc;
Oaa=-Oa;
Op=Oc;
Vp=Op*tan(beta);
Vd=Ve-Vc-Vp;
Vaa=Vp;
%Va=Vb+Vc;
Va=P-Vd-Vaa;
%Vaa=P-Vd-Va;
%Vp=Vaa;

F=sqrt(Op^2+Vp^2);


%Studio delle campate
%TRAVE AB
%Campata1 (0<x<L/2)
%Inizio
x1=0;
Nx1=-Oa*cos(alpha)-Va*sin(alpha);
Ty1=-Va*cos(alpha)+Oa*sin(alpha);
Mz1=(Va*cos(alpha)-Oa*sin(alpha))*x1;
%Fine
x2=L/2;
Nx2=Nx1;
Ty2=Ty1;
Mz2=(Va*cos(alpha)-Oa*sin(alpha))*x2;
%Campata2 (L/2<x<L)
%Inizio
x3=L/2;
Nx3=-(Oa+Oc)*cos(alpha)+(Vc-Va)*sin(alpha);
Ty3=(Vc-Va)*cos(alpha)+(Oa+Oc)*sin(alpha);
Mz3=(Va*cos(alpha)-Oa*sin(alpha))*x3-(Vc*cos(alpha)+Oc*sin(alpha))*(x3-L/2);
%Fine
x4=L;
Nx4=Nx3;
Ty4=Ty3;
Mz4=(Va*cos(alpha)-Oa*sin(alpha))*x4-(Vc*cos(alpha)+Oc*sin(alpha))*(x4-L/2);


%TRAVE DE
%Campata1 (0<x<EP)
%Inizio
x5=0;
Nx5=-Ve*sin(alpha);
Ty5=Ve*cos(alpha);
Mz5=-Ve*cos(alpha)*x5;
%Fine
x6=EP;
Nx6=Nx5;
Ty6=Ty5;
Mz6=-Ve*cos(alpha)*x6;
%Campata2 (EP<x<L/2)
%Inizio
x7=EP;
Nx7=(Vp-Ve)*sin(alpha)-Op*cos(alpha);
Ty7=(Ve-Vp)*cos(alpha)-Op*sin(alpha);
Mz7=-Ve*cos(alpha)*x7+(Vp*cos(alpha)+Op*sin(alpha))*(x7-EP);
%Fine
x8=L/2;
Nx8=Nx7;
Ty8=Ty7;
Mz8=-Ve*cos(alpha)*x8+(Vp*cos(alpha)+Op*sin(alpha))*(x8-EP);
%Campata3 (L/2<x<L)
%Inizio
x9=L/2;
Nx9=(Vp+Vc-Ve)*sin(alpha)+(Oc-Op)*cos(alpha);
Ty9=(Ve-Vp-Vc)*cos(alpha)+(Oc-Op)*sin(alpha);
Mz9=-Ve*cos(alpha)*x9+(Vp*cos(alpha)+Op*sin(alpha))*(x9-EP)+(Vc*cos(alpha)-Oc*sin(alpha))*(x9-L/2);
%Fine
x10=L;
Nx10=Nx9;
Ty10=Ty9;
Mz10=-Ve*cos(alpha)*x10+(Vp*cos(alpha)+Op*sin(alpha))*(x10-EP)+(Vc*cos(alpha)-Oc*sin(alpha))*(x10-L/2);

%Diagrammi delle sollecitazioni
%TRAVE AB
x_ab=[0,0,L/2,L/2,L,L];
Nx_ab=[0,Nx1,Nx2,Nx3,Nx4,0];
Ty_ab=[0,Ty1,Ty2,Ty3,Ty4,0];
Mz_ab=[0,Mz1,Mz2,Mz3,Mz4,0];

figure

subplot(3,1,1)
plot(x_ab, Nx_ab)
title('Diagramma delle sollecitazioni Trave AB')
xlabel('Coordinata longitudinale x [m]')
ylabel('Sforzo normale N_x [N]')
grid on
subplot(3,1,2)
plot(x_ab, Ty_ab,'r')
xlabel('Coordinata longitudinale x [m]')
ylabel('Sforzo tangenziale T_y [N]')
grid on
subplot(3,1,3)
plot(x_ab, Mz_ab,'g')
xlabel('Coordinata longitudinale x [m]')
ylabel('Momento flettente M_z [Nm]')
grid on

%TRAVE DE
x_de=[0,0,EP,EP,L/2,L/2,L,L];
Nx_de=[0,Nx5,Nx6,Nx7,Nx8,Nx9,Nx10,0];
Ty_de=[0,Ty5,Ty6,Ty7,Ty8,Ty9,Ty10,0];
Mz_de=[0,Mz5,Mz6,Mz7,Mz8,Mz9,Mz10,0];

figure
subplot(3,1,1)
plot(x_de, Nx_de)
title('Diagramma delle sollecitazioni Trave DE')
xlabel('Coordinata longitudinale x [m]')
ylabel('Sforzo normale N_x [N]')
grid on
subplot(3,1,2)
plot(x_de, Ty_de,'r')
xlabel('Coordinata longitudinale x [m]')
ylabel('Sforzo tangenziale T_y [N]')
grid on
subplot(3,1,3)
plot(x_de, Mz_de,'g')
xlabel('Coordinata longitudinale x [m]')
ylabel('Momento flettente M_z [Nm]')
grid on

%A seguire i valori in modulo
%TRAVE AB
%Sezione C
N_ab_c=abs(Nx2);
T_ab_c=abs(Ty2);
Mf_ab_c=abs(Mz2);
%d_ab_c=d_min_L2+2*s_b;
d_ab_c=d_min_L2;

%TRAVE ED
%Sezione C
N_ed_c=abs(Nx8);
T_ed_c=abs(Ty8);
Mf_ed_c=abs(Mz8);
%d_ed_c=d_min_L2+2*s_b;
d_ed_c=d_min_L2;

%Sezione P
N_ed_p=abs(Nx7);
T_ed_p=abs(Ty7);
Mf_ed_p=abs(Mz7);
%d_ed_p=d_min_L1+2*s_b;
d_ed_p=d_min_L1;

%fprintf
fprintf('\nTRAVE AB:\nSezione C: N=%.1f N T=%.1f N M_f=%.1f Nm d=%.1f mm', N_ab_c,T_ab_c,Mf_ab_c,d_ab_c);
fprintf('\nTRAVE EP:\nSezione C: N=%.1f N T=%.1f N M_f=%.1f Nm d=%.1f mm', N_ed_c,T_ed_c,Mf_ed_c,d_ed_c);
fprintf('\nSezione P: N=%.1f N T=%.1f N M_f=%.1f Nm d=%.1f mm\n', N_ed_p,T_ed_p,Mf_ed_p,d_ed_p);

%% Dimensionamento e verifica del punto C dei 2 bracci

Sig_id_max_V_c=Rp02/CS_m;
Wz_ab_c=(Mf_ab_c*10^3)/Sig_id_max_V_c; %mm^3
h_c=(6*r*Wz_ab_c+d_ab_c^3)^(1/3); %usiamo la formula della sezione forata piena
%per trovare una h_c e una b_c i primo tentativo
h_c=round(h_c); %mm
b_c=h_c/r;
b_c=round(b_c); %mm
A_c=b_c*(h_c-d_ab_c); %mm^2
fprintf('\nTRAVE AB-EC:\nSezione C:\nValori di primo tentativo (considerando una sezione forata piena): h=%.1f mm b=%.1f mm\n', h_c,b_c);

%Stato di tensione
%Bordo esterno superiore e inferiore (y=+/-h/2)
%Trave AB
A_ab_c=b_c*h_c-(b_c-2*s_t)*(h_c-2*s_t)-2*d_ab_c*s_t; %mm^2
Wz_ab_c=(1/(6*h_c))*(b_c*(h_c^3-d_ab_c^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ab_c^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_N_ab_c=N_ab_c/A_ab_c; %MPa
Sig_Mf_ab_c=(Mf_ab_c*10^3)/Wz_ab_c; %MPa
Sig_ab_c=Sig_N_ab_c+Sig_Mf_ab_c; %Mpa
Tau_ab_c=0;
Sig_id_max_V_ab_c=Sig_ab_c; %Mpa
CS_ab_c=Rp02/Sig_id_max_V_ab_c;
%Trave ED
A_ed_c=b_c*h_c-(b_c-2*s_t)*(h_c-2*s_t)-2*d_ed_c*s_t; %mm^2
Wz_ed_c=(1/(6*h_c))*(b_c*(h_c^3-d_ed_c^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ed_c^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_N_ed_c=N_ed_c/A_ed_c;
Sig_Mf_ed_c=(Mf_ed_c*10^3)/Wz_ed_c;
Sig_ed_c=Sig_N_ed_c+Sig_Mf_ed_c;
Tau_ed_c=0;
Sig_id_max_V_ed_c=Sig_ed_c;
CS_ed_c=Rp02/Sig_id_max_V_ed_c;


%Bordo foro ( (y=+/-d/2))
S_c_f=b_c/2*(h_c^2/4-d_ab_c^2/4); %(in y=d/2)
%Trave AB
Jz_ab_c=(1/12)*(b_c*(h_c^3-d_ab_c^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ab_c^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_Mf_ab_c_f=(Mf_ab_c*10^3*d_ab_c)/(Jz_ab_c*2); %MPa %f=foro
Sig_ab_c_f=Sig_N_ab_c+Sig_Mf_ab_c_f; %Mpa
Tau_ab_c_f=(T_ab_c*S_c_f)/(Jz_ab_c*b_c); %MPa
Sig_id_max_V_ab_c_f=sqrt(Sig_ab_c_f^2+3*Tau_ab_c_f^2); %Mpa
CS_ab_c_f=Rp02/Sig_id_max_V_ab_c_f;
%Trave ED
Jz_ed_c=(1/12)*(b_c*(h_c^3-d_ed_c^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ed_c^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_Mf_ed_c_f=(Mf_ed_c*10^3*d_ed_c)/(Jz_ed_c*2); %MPa %f=foro
Sig_ed_c_f=Sig_N_ed_c+Sig_Mf_ed_c_f; %Mpa
Tau_ed_c_f=(T_ed_c*S_c_f)/(Jz_ed_c*b_c); %MPa
Sig_id_max_V_ed_c_f=sqrt(Sig_ed_c_f^2+3*Tau_ed_c_f^2); %Mpa
CS_ed_c_f=Rp02/Sig_id_max_V_ed_c_f;


fprintf('\nSEZIONE C:\nTrave AB:\nBordo Esterno: Sigma=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ab_c,Tau_ab_c,Sig_id_max_V_ab_c,CS_ab_c);
fprintf('Bordo Interno: Sigma=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ab_c_f,Tau_ab_c_f,Sig_id_max_V_ab_c_f,CS_ab_c_f);
fprintf('Trave ED:\nBordo Esterno: Sigma_N=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ed_c,Tau_ed_c,Sig_id_max_V_ed_c,CS_ed_c);
fprintf('Bordo Interno: Sigma=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ed_c_f,Tau_ed_c_f,Sig_id_max_V_ed_c_f,CS_ed_c_f);

%% Verifica del punto P del braccio ED
%Trave ED
%Consideriamo gli stessi valori di a e b
%Stato di tensione
A_ed_p=b_c*h_c-(b_c-2*s_t)*(h_c-2*s_t)-2*d_ed_p*s_t; %mm^2
Wz_ed_p=(1/(6*h_c))*(b_c*(h_c^3-d_ed_p^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ed_p^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
%Bordo esterno superiore e inferiore (y=+/-h/2)
Sig_N_ed_p=N_ed_p/A_ed_p; %MPa
Sig_Mf_ed_p=(Mf_ed_p*10^3)/Wz_ed_p; %MPa
Sig_ed_p=Sig_N_ed_p+Sig_Mf_ed_p; %Mpa
Tau_ed_p=0;
Sig_id_max_V_ed_p=Sig_ed_p; %Mpa
CS_ed_p=Rp02/Sig_id_max_V_ed_p;

%Bordo foro ( (y=+/-d/2))
S_p_f=b_c/2*(h_c^2/4-d_ed_p^2/4); %(in y=d/2)
Jz_ed_p=(1/12)*(b_c*(h_c^3-d_ed_p^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ed_p^3)); %mm^3 %valore ottenuto considerando la sezione forata cava

Sig_Mf_ed_p_f=(Mf_ed_p*10^3*d_ed_p)/(Jz_ed_p*2); %MPa %f=foro
Sig_ed_p_f=Sig_N_ed_p+Sig_Mf_ed_p_f; %Mpa
Tau_ed_p_f=(T_ed_p*S_p_f)/(Jz_ed_p*b_c); %MPa
Sig_id_max_V_ed_p_f=sqrt(Sig_ed_p_f^2+3*Tau_ed_p_f^2); %Mpa
CS_ed_p_f=Rp02/Sig_id_max_V_ed_p_f;
fprintf('\nSEZIONE P:\nTrave ED:\nBordo Esterno Sigma=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ed_p,Tau_ed_p,Sig_id_max_V_ed_p,CS_ed_p);
fprintf('Bordo Interno Sigma_N=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ed_p_f,Tau_ed_p_f,Sig_id_max_V_ed_p_f,CS_ed_p_f);

%Se CS>1.5 vuol dire che i valori di b e h individuati veificano anche il
%punto P

%% Valori di h e b che soddisfano la condizione CS>1.5
b_c=input('\nImmettere valore di b (larghezza tubolare) desiderata (in millimetri): ');
%b_c=40 
h_c=input('Immettere valore di h (altezza tubolare) desiderata (in millimetri): ');
%h_c=80 
s_t=input('Immettere valore di s(spessore tubolare) desiderato (in millimetri): ');
%s_t=3 

%Sezione C
%Bordo esterno superiore e inferiore (y=+/-h/2)
%Trave AB
A_ab_c=b_c*h_c-(b_c-2*s_t)*(h_c-2*s_t)-2*d_ab_c*s_t; %mm^2
Wz_ab_c=(1/(6*h_c))*(b_c*(h_c^3-d_ab_c^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ab_c^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_N_ab_c=N_ab_c/A_ab_c; %MPa
Sig_Mf_ab_c=(Mf_ab_c*10^3)/Wz_ab_c; %MPa
Sig_ab_c=Sig_N_ab_c+Sig_Mf_ab_c; %Mpa
Tau_ab_c=0;
Sig_id_max_V_ab_c=Sig_ab_c; %Mpa
CS_ab_c=Rp02/Sig_id_max_V_ab_c;
%Trave ED
A_ed_c=b_c*h_c-(b_c-2*s_t)*(h_c-2*s_t)-2*d_ed_c*s_t; %mm^2
Wz_ed_c=(1/(6*h_c))*(b_c*(h_c^3-d_ed_c^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ed_c^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_N_ed_c=N_ed_c/A_ed_c;
Sig_Mf_ed_c=(Mf_ed_c*10^3)/Wz_ed_c;
Sig_ed_c=Sig_N_ed_c+Sig_Mf_ed_c;
Tau_ed_c=0;
Sig_id_max_V_ed_c=Sig_ed_c;
CS_ed_c=Rp02/Sig_id_max_V_ed_c;

%Bordo foro ( (y=+/-d/2))
S_c_f=b_c/2*(h_c^2/4-d_ab_c^2/4); %(in y=d/2)
%Trave AB
Jz_ab_c=(1/12)*(b_c*(h_c^3-d_ab_c^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ab_c^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_Mf_ab_c_f=(Mf_ab_c*10^3*d_ab_c)/(Jz_ab_c*2); %MPa %f=foro
Sig_ab_c_f=Sig_N_ab_c+Sig_Mf_ab_c_f; %Mpa
Tau_ab_c_f=(T_ab_c*S_c_f)/(Jz_ab_c*b_c); %MPa
Sig_id_max_V_ab_c_f=sqrt(Sig_ab_c_f^2+3*Tau_ab_c_f^2); %Mpa
CS_ab_c_f=Rp02/Sig_id_max_V_ab_c_f;
%Trave ED
Jz_ed_c=(1/12)*(b_c*(h_c^3-d_ed_c^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ed_c^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_Mf_ed_c_f=(Mf_ed_c*10^3*d_ed_c)/(Jz_ed_c*2); %MPa %f=foro
Sig_ed_c_f=Sig_N_ed_c+Sig_Mf_ed_c_f; %Mpa
Tau_ed_c_f=(T_ed_c*S_c_f)/(Jz_ed_c*b_c); %MPa
Sig_id_max_V_ed_c_f=sqrt(Sig_ed_c_f^2+3*Tau_ed_c_f^2); %Mpa
CS_ed_c_f=Rp02/Sig_id_max_V_ed_c_f;

fprintf('\nSEZIONE C:\nTrave AB:\nBordo Esterno: Sigma=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ab_c,Tau_ab_c,Sig_id_max_V_ab_c,CS_ab_c);
fprintf('Bordo Interno: Sigma=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ab_c_f,Tau_ab_c_f,Sig_id_max_V_ab_c_f,CS_ab_c_f);
fprintf('Trave ED:\nBordo Esterno: Sigma_N=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ed_c,Tau_ed_c,Sig_id_max_V_ed_c,CS_ed_c);
fprintf('Bordo Interno: Sigma=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ed_c_f,Tau_ed_c_f,Sig_id_max_V_ed_c_f,CS_ed_c_f);

%Sezione P
%Trave ED
%Consideriamo gli stessi valori di a e b
%Stato di tensione
A_ed_p=b_c*h_c-(b_c-2*s_t)*(h_c-2*s_t)-2*d_ed_p*s_t; %mm^2
Wz_ed_p=(1/(6*h_c))*(b_c*(h_c^3-d_ed_p^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ed_p^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
%Bordo esterno superiore e inferiore (y=+/-h/2)
Sig_N_ed_p=N_ed_p/A_ed_p; %MPa
Sig_Mf_ed_p=(Mf_ed_p*10^3)/Wz_ed_p; %MPa
Sig_ed_p=Sig_N_ed_p+Sig_Mf_ed_p; %Mpa
Tau_ed_p=0;
Sig_id_max_V_ed_p=Sig_ed_p; %Mpa
CS_ed_p=Rp02/Sig_id_max_V_ed_p;

%Bordo foro ( (y=+/-d/2))
S_p_f=b_c/2*(h_c^2/4-d_ed_p^2/4); %(in y=d/2)
Jz_ed_p=(1/12)*(b_c*(h_c^3-d_ed_p^3)-(b_c-2*s_t)*((h_c-2*s_t)^3-d_ed_p^3)); %mm^3 %valore ottenuto considerando la sezione forata cava
Sig_Mf_ed_p_f=(Mf_ed_p*10^3*d_ed_p)/(Jz_ed_p*2); %MPa %f=foro
Sig_ed_p_f=Sig_N_ed_p+Sig_Mf_ed_p_f; %Mpa
Tau_ed_p_f=(T_ed_p*S_p_f)/(Jz_ed_p*b_c); %MPa
Sig_id_max_V_ed_p_f=sqrt(Sig_ed_p_f^2+3*Tau_ed_p_f^2); %Mpa
CS_ed_p_f=Rp02/Sig_id_max_V_ed_p_f;
fprintf('\nSEZIONE P:\nTrave ED:\nBordo Esterno Sigma=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ed_p,Tau_ed_p,Sig_id_max_V_ed_p,CS_ed_p);
fprintf('Bordo Interno Sigma_N=%.1f MPa Tau=%.1f MPa Sigma_id_max_VonMises=%.1f MPa CS=%.1f\n', Sig_ed_p_f,Tau_ed_p_f,Sig_id_max_V_ed_p_f,CS_ed_p_f);

%Se CS>1.5 vuol dire che i valori di b e h individuati veificano anche il
%punto P
