%-----------------------------------------------------------
% arquivo: exemplo_9_1.m
% Cálculo das forças internas em um compressor de ar.
% Ver Doughty, Cap.9-Reações e Forças Internas
%-----------------------------------------------------------

clear all
close all
clc

%% GRÁFICOS ATIVOS
% 's' - ativo
% 'n' - inativo
grafico_pressao='s';
grafico_volume='n';
grafico_an_geral='n';
grafico_det_J='n';
grafico_x_y_PI='n';
grafico_veloc_acel_PI='n';
grafico_reacoes_externas='s';
grafico_F126='s';
%% DIMENSIONAMENTO DAS VARIÁVEIS
NP=361;
A=zeros(NP,1);
X=zeros(NP,1);
detJ=zeros(NP,1);
Ka=zeros(NP,1);
Kx=zeros(NP,1);
La=zeros(NP,1);
Lx=zeros(NP,1);
Xp1=zeros(NP,1);
Yp1=zeros(NP,1);
Kxp1=zeros(NP,1);
Kyp1=zeros(NP,1);
Lxp1=zeros(NP,1);
Lyp1=zeros(NP,1);
Xp2=zeros(NP,1);
Yp2=zeros(NP,1);
Kxp2=zeros(NP,1);
Kyp2=zeros(NP,1);
Lxp2=zeros(NP,1);
Lyp2=zeros(NP,1);
V=zeros(NP,1);
Mo=zeros(NP,1);
F1=zeros(NP,1);
F2=zeros(NP,1);
F3=zeros(NP,1);
F4=zeros(NP,1);
F5=zeros(NP,1);
F6=zeros(NP,1);
F7=zeros(NP,1);
Rx=zeros(NP,1);
Ry=zeros(NP,1);
Rz=zeros(NP,1);
pressao=zeros(NP,1);

% syms M1 M2 M3 P P1 Ac q qp qpp
% syms Kx1 Ky1 Kx2 Ky2 Kx3 Ky3 Lx1 Ly1 Lx2 Ly2 Lx3 Ly3
% syms U1 U2 L R X3
% syms A X
% syms Ka La Kx Lx
% syms I1 I2
% syms Co C1
% syms F1 F2 F3 F4 F5 F6 F7 Rx Ry Rz Mo

%% DADOS NUMÉRICOS (SI)
% D=0.1397;   %   m
% R=0.1016;   %   m
% L=0.5080;   %   m
% V1=0.0033;  %   m^3
% P1=101352;  %   Pa
% P2=517104;  %   Pa
% U1=-0.0024; %   m
% W1=1777;    %   N
% I1=9.3054;  %   kg.m^2
% W2=152;     %   N
% I2=0.3338;  %   kg.m^2
% U2=0.1326;  %   m
% W3=214;     %   N
% W4=58;      %   N
% W5=81;      %   N
% Co=16041.5; %   in-lb   coeficiente linear da curva de torque
% C1=-680.82; %   in-lb-s coeficiente angular da curva de torque
%=============================================
%% DADOS NUMÉRICOS (US)
D=5.5;      %   in 
R=4;        %   in
L=20;       %   in 
V1=202;     %   in^3
P1=14.7;    %   psia 
P4=P1;
P2=89.7;    %   psia 
P3=P2;
X1=16;
X2=22.166;
X3=24;
X4=22.674;
U1=-0.09373;%	in 
W1=399.6;   %   lbf
I1=82.36;   %   lbf-s2-in
W2=34.2;    %   lbf
I2=2.9543;  %   lbf-s2-in
U2=5.22;    %   in
W3=48.2;    %   lbf
W4=13;      %   lbf
W5=18.3;    %   lbf
Co=16041.5; %   in-lb
C1=-680.82; %   in-lb-s
%=============================================
gama=1.4;   %           expoente da politrópica
%=============================================
g=386.088;  %   386.088 in/s^2 
M1=W1/g;
M2=W2/g;
M3=W3/g;
%% NÚMERO DE PONTOS
qmin=0*pi/180;
qmax=360*pi/180;
q=linspace(qmin,qmax,NP);
%% EQUAÇÕES DE EQUILÍBRIO
% somente para apresentação
% f=[F1-F3-M1*(-qpp*U1*sin(q)-qp^2*U1*cos(q));
%    F2-F4-M1*(qpp*U1*cos(q)-qp^2*U1*sin(q));
%    Mo+F3*R*sin(q)-F4*R*cos(q)-I1*qpp;
%    F3-F5-M2*(qpp*Kx2+qp^2*Lx2);
%    F4-F6-M2*(qpp*Ky2+qp^2*Ly2);
%    F3*U2*sin(A)+F4*U2*cos(A)+F5*(L-U2)*sin(A)+F6*(L-U2)*cos(A)-I2*(qpp*Ka+qp^2*La);
%    F5-(P-P1)*Ac-M3*(qpp*Kx3+qp^2*Lx3);
%    F6-F7;
%    Rx-F1+(P-P1)*Ac;
%    Ry-F2+F6;
%    Rz+X3*F6;
%    Mo-(P-P1)*Ac*Kx3];
for i=1:NP
    %% CINEMÁTICA
    A(i)=asin((R*sin(q(i)))/L);
    X(i)=R*cos(q(i))+L*cos(A(i));
    J=[-L*sin(A(i))  -1;
       -L*cos(A(i))   0];
    detJ(i)=det(J);
    Q=[-R*sin(q(i));
        R*cos(q(i))];
    k=-inv(J)*Q;
    Ka(i)=k(1);
    Kx(i)=k(2);
    dJdq=[-Ka(i)*L*cos(A(i)) 0;
        Ka(i)*L*sin(A(i)) 0];
    dQdq=[-R*cos(q(i));
        -R*sin(q(i))];
    l=-inv(J)*(dJdq*k+dQdq);
    La(i)=l(1);
    Lx(i)=l(2);
    %% PONTO DE INTERESSE 1
    % posição
    Ori=[0;
        0];
    Rot=[cos(q(i)) -sin(q(i))
         sin(q(i))  cos(q(i))];
    Ploc=[U1;
           0];
    Pglob=Ori+Rot*Ploc;
    Xp1(i)=Pglob(1);
    Yp1(i)=Pglob(2);
    % velocidade
    dOdq=[0;
        0];
    dRotdq=[sin(q(i)) -cos(q(i));
            cos(q(i)) -sin(q(i))];
    K1 = dOdq+dRotdq*Ploc;
    Kxp1(i)=K1(1);
    Kyp1(i)=K1(2);
    %aceleração
    d2Odq2 = [0;
              0];
    d2Rotdq2 = [ cos(q(i))   sin(q(i));
                -sin(q(i))  -cos(q(i))];
    L1 = d2Odq2+d2Rotdq2*Ploc;
    Lxp1(i)=L1(1);
    Lyp1(i)=L1(2);
    %% PONTO DE INTERESSE 2
    % posição
    Ploc=[U2;
           0];
    Ori=[R*cos(q(i));
         R*sin(q(i))];
    Rot = [ cos(A(i)) sin(A(i));
           -sin(A(i)) cos(A(i))];
    Pglob=Ori+Rot*Ploc;
    Xp2(i)=Pglob(1);
    Yp2(i)=Pglob(2);
    % velocidade
    dOdq = [-R*sin(q(i));
             R*cos(q(i))];
    dRotdq = [-sin(A(i))  cos(A(i));
              -cos(A(i)) -sin(A(i))]*Ka(i);
    K2 = dOdq+dRotdq*Ploc;
    Kxp2(i)=K2(1);
    Kyp2(i)=K2(2);
    % aceleração
    d2Odq2 = [-R*cos(q(i));
              -R*sin(q(i))];
    d2Rotdq2 = [-cos(A(i))  -sin(A(i));
                 sin(A(i))  -cos(A(i))]*Ka(i)^2+...
        [-sin(A(i))   cos(A(i));
         -cos(A(i))  -sin(A(i))]*La(i);
    L1 = d2Odq2+d2Rotdq2*Ploc;
    Lxp2(i)=L1(1);
    Lyp2(i)=L1(2);
    %% PRESSÃO DO CILINDRO
    Ac=pi*D^2/4;
    Dc=V1/Ac-2*R;
    V2=Ac*(Dc+R+L-X2);
    V3=Ac*(Dc+R+L-X3);
    V4=Ac*(Dc+R+L-X4);
    V(i)=Ac*(Dc+R+L-X(i));
    if q(i)<pi
        if X(i)>X4
            pressao(i)=P3*(V3/V(i))^gama;
        else
            pressao(i)=P1;
        end
    else
        if X(i)<X2
            pressao(i)=P1*(V1/V(i))^gama;
        else
            pressao(i)=P2;
        end
    end
    %% SOLUÇÃO ESTÁTICA
%     qp=0;
%     qpp=0;
%     Mo(i)=(pressao(i)-P1)*Ac*Kx(i);
    %% SOLUÇÃO CINETOESTÁTICA
    qp=208*2*pi()/60;   % velocidade nominal
    qpp=0;
    Mo(i)=Co+C1*qp;
    
    F5(i)=M3*(qpp*Kx(i)+qp^2*Lx(i))+(pressao(i)-P1)*Ac;
    F3(i)=F5(i)+M2*(qpp*Kxp2(i)+qp^2*Lxp2(i));
    F1(i)=F3(i)-M1*(qpp*U1*sin(q(i))+qp^2*U1*cos(q(i)));
    F6(i)=(-F3(i)*U2*sin(A(i))-F5(i)*(L-U2)*sin(A(i))+I2*(qpp*Ka(i)+qp^2*La(i))-M2*(qpp*Kyp2(i)+qp^2*Lyp2(i)))/...
        (cos(A(i))*(L-U2));
    F4(i)=F6(i)+M2*(qpp*Kyp2(i)+qp^2*Lyp2(i));
    F2(i)=F4(i)+M1*(qpp*U1*cos(q(i))-qp^2*U1*sin(q(i)));
    F7(i)=F6(i);
    Rx(i)=F1(i)-(pressao(i)-P1)*Ac;
    Ry(i)=F2(i)-F6(i);
    Rz(i)=-X(i)*F6(i);
    

   
    
    
end
%=============================================================
if(grafico_pressao=='s')
%-------------------------------------------------------------
%% GRÁFICO DA PRESSÃO 
%-------------------------------------------------------------
figure('Name','Pressão x deslocamento')
plot(X,pressao,'LineWidth',2)
grid on
xlabel('X (m)')
ylabel('P (psia)')
title('Pressão')
end
if(grafico_volume=='s')
%-------------------------------------------------------------
%% GRÁFICO DO VOLUME
%-------------------------------------------------------------
figure('Name','Volume x deslocamento')
plot(X,V,'LineWidth',2)
grid on
xlabel('X (m)')
ylabel('V (in^3)')
title('Volume')
end
if(grafico_an_geral=='s')
%% GRÁFICO DA ANÁLISE GERAL
figure('Name','Análise Geral')  
% gráficos geminados
% primeiro gráfico
subplot(3,2,1) 
plot(q*180/pi,A*180/pi,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('A(q) (^o)')
title('Posição da biela')
xlim([qmin*180/pi qmax*180/pi])
% segundo gráfico
subplot(3,2,2)
plot(q*180/pi,X,'r','LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('X(q) (mm)')
title('Posição do pistão')
xlim([qmin*180/pi qmax*180/pi])
% terceiro gráfico
subplot(3,2,3)
plot(q*180/pi,Ka,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Ka(q) ')
title('Coeficiente de velocidade Ka')
xlim([qmin*180/pi qmax*180/pi])
% quarto gráfico
subplot(3,2,4)
plot(q*180/pi,Kx,'r','LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Kx(q) ')
title('Coeficiente de velocidade Kx')
xlim([qmin*180/pi qmax*180/pi])
% quinto gráfico
subplot(3,2,5)
plot(q*180/pi,La,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('La(q) ')
title('Derivada do coeficiente de velocidade La')
xlim([qmin*180/pi qmax*180/pi])
% sexto gráfico
subplot(3,2,6)
plot(q*180/pi,Lx,'r','LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Lx(q) ')
title('Derivada do coeficiente de velocidade Lx')
xlim([qmin*180/pi qmax*180/pi])
end
if(grafico_det_J=='s')
%-------------------------------------------------------------
%% GRÁFICO DO DETERMINANTE DA MATRIZ JACOBIANA
%-------------------------------------------------------------
figure('Name','Determinante da Matriz Jacobiana')% título na moldura 
plot(q*180/pi,detJ,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Det J(q)')
title('Determinante da Matriz Jacobiana')
xlim([qmin*180/pi qmax*180/pi])
ylim([-inf 0])
end
if(grafico_x_y_PI=='s')
%-------------------------------------------------------------
%% GRÁFICO X-Y DOS PONTOS DE INTERESSE
%-------------------------------------------------------------
figure('Name','Análise de Ponto de Interesse')
plot(Xp1,Yp1,Xp2,Yp2,'LineWidth',2)
grid on
xlabel('X (m)')
ylabel('Y (m)')
title('Posição')
axis equal
legend('Ponto 1','Ponto 2')
end
if(grafico_veloc_acel_PI=='s')
%-------------------------------------------------------------
%% GRÁFICO VELOCIDADE E ACELERAÇÃO DOS PONTOS DE INTERESSE
%-------------------------------------------------------------
figure('Name','Análise de Ponto de Interesse')
%------------------------------------
% gráfico 1
subplot(3,2,1)
plot(q*180/pi,Xp1,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Xp(q)')
title('Coordenada Xp(q)')
xlim([qmin*180/pi qmax*180/pi])
%------------------------------------
% gráfico 2
subplot(3,2,2)
plot(q*180/pi,Yp1,'r','LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Yp(q) ')
title('Coordenada Yp(q)')
xlim([qmin*180/pi qmax*180/pi])
%------------------------------------
% gráfico 3
subplot(3,2,3)
plot(q*180/pi,Kxp1,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Kpx(q)')
title('Coeficiente de velocidade Kpx')
xlim([qmin*180/pi qmax*180/pi])
%------------------------------------
% gráfico 4
subplot(3,2,4)
plot(q*180/pi,Kyp1,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Kpy(q) ')
title('Coeficiente de velocidade Kpy')
xlim([qmin*180/pi qmax*180/pi])
%------------------------------------
% gráfico 5
subplot(3,2,5)
plot(q*180/pi,Lxp1,'r','LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Lpx(q) ')
title('Derivada do coeficiente de velocidade Lpx')
xlim([qmin*180/pi qmax*180/pi])
%------------------------------------
% gráfico 6
subplot(3,2,6)
plot(q*180/pi,Lyp1,'r','LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('Lpy(q) ')
title('Derivada do coeficiente de velocidade Lpy')
xlim([qmin*180/pi qmax*180/pi])
end
if(grafico_reacoes_externas=='s')
%% GRÁFICO DAS REAÇÕES EXTERNAS
figure('Name','Reações Externas')
plot(q*180/pi,Rx,q*180/pi,Ry,q*180/pi,Rz/10,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('R (lbf)')
legend('R_x','R_y','R_z')
title('Reações Externas')
xlim([qmin*180/pi qmax*180/pi])
end
if(grafico_F126=='s')
%% GRÁFICO DAS REAÇÕES F1, F2 E F6
figure('Name','Reações F1, F2 e F6')
plot(q*180/pi,F1,q*180/pi,F2,q*180/pi,F6,'LineWidth',2)
grid on
xlabel('q (^o)')
ylabel('F (lbf)')
legend('F_1','F_2','F_6')
title('Reações F1, F2 e F6')
xlim([qmin*180/pi qmax*180/pi])
end