0%% DIMENSIONAMENTO STRUTTURALE
clear all
close all
clc

L=1.700; %m %lunghezza travi centrali (bracci)
EG=0.625; %m %distanza
EP_f=0.3; %m %distanza fissa imposta
EP_max=0.7; %m %distanza masssima
EP_min=0.1; %m %0.3 distanza minima
AAa=0.45; %m distanza tra la cerniera A della base e la cerniera A' (Aa) dell'attuatore %arbitraria
%AAa=0.3;
m=200; %kg %massa
g=10; %m/s^2 %accelerazione di gravità
P=m*g; %N %peso
H_f=1; %m %altezza fissa imposta
H_max=1.2; %m %altezza massima dal vincolo
H_min=0.1; %m altezza minima dal vincolo
v=100; %lunghezza vettori
p=6; %bar %pressione di alimentazione
p_a=1; %bar %pressione ambiente
k=0.1; %sovrastima della forza di spinta attuatore
%N.B.
%Studio1: PRIMO TENTATIVO (H=IMPOSTO; EP=IMPOSTO);
%Studio2: PARAMETRIZZAZIONE RISPETTO A H (H=VARIABILE; EP=IMPOSTO)->TROVO H_min EFFETTIVA
%Studio3: PARAMETRIZZAZIONE RISPETTO A EP (H=IMPOSTO; EP=VARIABILE)->TROVO EP EFFETTIVA

%% 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

%% PARAMETRIZZAZIONE RISPETTO A H
%Obiettivo: determinare l'altezza minima da terra al fine di non avere
%tensioni eccessivamente elevate sulle travi

%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=linspace(H_min,H_max,v); %m %altezza tra i vincoli
AD=sqrt(L^2-H.^2); %m %base
alpha=atan(H./AD); %radianti
%alpha_gradi=alpha*180/pi;
EP=EP*ones(1,length(H));

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 in funzione di H

figure
subplot(3,1,1)
plot(H, Nx1,'b',H,Nx3,'r') %Nx1=Nx2, Nx3=Nx4
title('Diagramma delle sollecitazioni Trave AB in funzione di H')
xlabel('Altezza H [m]')
ylabel('Sforzo normale N_x [N]')
legend('Campata 1','Campata 2')
grid on
subplot(3,1,2)
plot(H, Ty1,'b',H,Ty3,'r') %Ty1=Ty2, Ty3=Ty4
xlabel('Altezza H [m]')
ylabel('Sforzo tangenziale T_y [N]')
legend('Campata 1','Campata 2')
grid on
subplot(3,1,3)
plot(H, Mz3) %Mz1=Mz2=0, Mz3=Mz4
xlabel('Altezza H [m]')
ylabel('Momento flettente M_z [Nm]')
legend('Sezione C')
axis([H_min H_max Mz3(1)-100 Mz3(1)+100])
grid on

%TRAVE DE
figure
subplot(3,1,1)
plot(H, Nx5,'b',H,Nx7,'r',H,Nx9,'g') %Nx5=Nx6, Nx7=Nx8, Nx9=Nx10
title('Diagramma delle sollecitazioni Trave DE in funzione di H')
xlabel('Altezza H [m]')
ylabel('Sforzo normale N_x [N]')
legend('Campata 1','Campata 2','Campata 3')
grid on
subplot(3,1,2)
plot(H, Ty5,'b',H,Ty7,'r',H,Ty9,'g') %Ty5=Ty6, Ty7=Ty8, Ty9=Ty10
xlabel('Altezza H [m]')
ylabel('Sforzo tangenziale Ty [N]')
legend('Campata 1','Campata 2','Campata 3')
grid on
subplot(3,1,3)
plot(H, Mz7,'b',H,Mz9,'r') %Mz5=Mz6=0, Mz7=Mz8, Mz9=Mz10
xlabel('Altezza H [m]')
ylabel('Momento flettente [Nm]')
legend('Sezione P','Sezione C')
grid on

%Dai grafici ottenuti si osserva che al di sotto di un certo valore di H,
%gli sforzi normali tendono ad infinito. Pertanto, al fine di evitare
%condizioni critiche, è oppertuno mantenersi ad un'altezza minima di 0.3 m
%(300 mm)

H_min_eff=input('Valore di H_min effettivo, letto dai diagrammi precedenti, tale da non determinare sforzi eccessivi (in metri): ');
%0.3

%% PARAMETRIZZAZIONE RISPETTO A EP
%Obiettivo: determinare la posizione del punto di applicazione
%dell'attutaore pneumatico P

%Costanti costruttive (variabili manualmente)
EP=linspace(EP_min,EP_max,v); %m %distanza

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

%Variabili operative
H=H_f*ones(1,v); %m %alteazza %imponiamo un valore arbitrario di H in modo da
%vedere, per la generica altezza H, cosa accade in termini di EP
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 in funzione di EP

figure
subplot(3,1,1)
plot(EP, Nx1,'b',EP,Nx3,'r') %Nx1=Nx2, Nx3=Nx4
title('Diagramma delle sollecitazioni Trave AB in funzione di EP')
xlabel('Posizione punto di applicazione attuatore P [m]')
ylabel('Sforzo normale N_x [N]')
legend('Campata 1','Campata 2')
grid on
subplot(3,1,2)
plot(EP, Ty1,'b',EP,Ty3,'r') %Ty1=Ty2, Ty3=Ty4
xlabel('Posizione punto di applicazione attuatore P [m]')
ylabel('Sforzo tangenziale T_y [N]')
legend('Campata 1','Campata 2')
grid on
subplot(3,1,3)
plot(EP, Mz3) %Mz1=Mz2=0, Mz3=Mz4
xlabel('Posizione punto di applicazione attuatore P [m]')
ylabel('Momento flettente M_z [Nm]')
legend('Sezione C')
axis([EP_min EP_max Mz3(1)-100 Mz3(1)+100])
grid on

%TRAVE DE
figure
subplot(3,1,1)
plot(EP, Nx5,'b',EP,Nx7,'r',EP,Nx9,'g') %Nx5=Nx6, Nx7=Nx8, Nx9=Nx10
title('Diagramma delle sollecitazioni Trave DE in funzione di EP')
xlabel('Posizione punto di applicazione attuatore P [m]')
ylabel('Sforzo normale N_x [N]')
legend('Campata 1','Campata 2','Campata 3')
grid on
subplot(3,1,2)
plot(EP, Ty5,'b',EP,Ty7,'r',EP,Ty9,'g') %Ty5=Ty6, Ty7=Ty8, Ty9=Ty10
xlabel('Posizione punto di applicazione attuatore P [m]')
ylabel('Sforzo tangenziale Ty [N]')
legend('Campata 1','Campata 2','Campata 3')
grid on
subplot(3,1,3)
plot(EP, Mz7,'b',EP,Mz9,'r') %Mz5=Mz6=0, Mz7=Mz8, Mz9=Mz10
xlabel('Posizione punto di applicazione attuatore P [m]')
ylabel('Momento flettente [Nm]')
legend('Sezione P','Sezione C')
grid on

%Dai grafici ottenuti si osserva che al di sopra di un certo valore di EP,
%gli sforzi normali tendono ad infinito (in particolare per EP tendente a
%EC). Pertanto, al fine di evitare condizioni critiche, è oppertuno
%mantenersi a valori di EP bassi


%% Limite geometria attuatore pneumatico
%H_min_eff=input('Valore di H_min effettivo, letto dal diagramma precedente, tale da non determinare sforzi eccessivi')
H_min=H_min_eff; %non più 0.1 come in precedenza dato che abbiamo visto esserci
%problemi legati a sforzi normali eccessivamente elevati
AD_max=sqrt(L^2-H_max^2);
AD_min=sqrt(L^2-H_min^2);
alpha_max=atan(H_max/AD_max);
alpha_min=atan(H_min/AD_min);
alpha_a_max=atan(H_max/AAa);
alpha_a_min=atan(H_min/AAa);
gamma_max=pi-alpha_max-alpha_a_max;
gamma_min=pi-alpha_min-alpha_a_min;
AaE_max=sqrt(H_max^2+AAa^2);
AaE_min=sqrt(H_min^2+AAa^2);
%AP_max=sqrt(H_max^2+EP.^2-2*H_max*EP.*sin(alpha_max));
AaP_max=sqrt(EP.^2+AaE_max^2-2*AaE_max.*EP*cos(gamma_max));
%AP_min=sqrt(H_min^2+EP.^2-2*H_min*EP.*sin(alpha_min));
AaP_min=sqrt(EP.^2+AaE_min^2-2*AaE_min.*EP*cos(gamma_min));
AaP_max_min=AaP_max./AaP_min;

figure
plot(EP,AaP_max_min)
title('Andamento del rapporto AP_{max}/AP_{min} in funzione del punto di applicazione P')
xlabel('Posizione punto di applicazione attuatore P [m]')
ylabel('Rapporto AaP_{max}/AaP_{min}')
axis([EP_min EP_max AaP_max_min(v)-1 AaP_max_min(1)+1])
grid on

EP_eff=input('Valore di EP effettivo, letto dal diagramma precedente, tale che AP_max/AP_min>2 (in metri): ');
%0.43

%Anche se avevamo visto che, sarebbe conveniente avere un valore di EP
%più basso possibile così da avere gli sforzi più bassi, abbiamo un limite
%dettato dalla geometria dell'attuatore. Infatti, il rapporto AP_max/AP_min
%non può essere minore di 2 per ragioni costruttive. Quindi, nel nostro
%caso EP>0.40m


%% Forza, corsa e alesaggio attuatore
%In questa sezione H_min=H_min_eff dato che il dimensionamento deve tenere
%conto della condizione di H_min effettiva

%Costanti costruttive (variabili manualmente)
EP=EP_eff; %m %distanza

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

%Variabili operative
H=linspace(H_min,H_max,v); %m %altezza tra i vincoli
AD=sqrt(L^2-H.^2); %m %base
alpha=atan(H./AD); %radianti
%alpha_gradi=alpha*180/pi;
EP=EP*ones(1,length(H));

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;

%Forza massima e alesaggio minimo 
F=sqrt(Op.^2+Vp.^2);
F_max_t=max(F); %N %forza teorica di spinta
F_max_eff=F_max_t+F_max_t*k; %N %forza effettiva di spinta
d_min=sqrt((4*F_max_eff)/(pi*(p-p_a)*10^5)); %m %alesaggio minimo attuatore

figure
plot(H,F)
title('Andamento forza F di spinta attuatore rispetto ad altezza H')
xlabel('Altezza H [m]')
ylabel('Forza di spinta attuatore F [N]')
axis([H_min H_max F(end)-100 F(1)+100])
grid on

%Corsa dell'attuatore pneumatico
EP=EP_eff;
%AP_max=sqrt(H_max^2+EP^2-2*H_max*EP*sin(alpha_max));
AaP_max=sqrt(EP^2+AaE_max^2-2*EP*AaE_max*cos(gamma_max));
%AP_min=sqrt(H_min^2+EP^2-2*H_min*EP*sin(alpha_min));
AaP_min=sqrt(EP^2+AaE_min^2-2*EP*AaE_min*cos(gamma_min));
c=AaP_max-AaP_min; %m %corsa attuatore

%output
fprintf('Forza teorica spinta attuatore[N]= %.1f \n', F_max_t);
fprintf('Alesaggio minimo attuatore [mm]= %.1f \n', d_min*10^3);
fprintf('Corsa attuatore [mm]= %.1f \n', c*10^3);

%Forza nel punto di cerniera C
F_c=sqrt(Oc.^2+Vc.^2);
Fc_max_t=max(F_c);
fprintf('\nDati utili in "Dimensionamento_strutturale_2":\n');
fprintf('Forza teorica spinta attuatore F_max_t [N]= %.1f \n', F_max_t);
fprintf('Forza teorica cerniera C Fc_max_t [N]= %.1f \n', Fc_max_t);

figure
plot(H,F_c)
title('Andamento forza Fc_{tot} rispetto ad altezza H')
xlabel('Altezza H [m]')
ylabel('Forza sulla cerniera C      Fc_{tot} [N]')
axis([H_min H_max F(end)-100 F(1)+100])
grid on