%=============================================================
% Mecanismos
% Prof. Jorge Luiz Erthal
% Análise geral e de ponto de interesse do exemplo 2.1
% Solução analítica 
%=============================================================
clear all
close all
clc
%-------------------------------------------------------------
%% dados numéricos
%-------------------------------------------------------------
R=0.5;         % Comprimento da manivela (m)0,5
C=2.5;         % Distancia do pistao ao centro da manivela (m)2,5
C1=3.5;
%-------------------------------------------------------------
% ponto de interesse
Up=C1;
Vp=0;
%-------------------------------------------------------------
%% número de pontos
%-------------------------------------------------------------
NP = 361;
qmin = 0*pi/180;
qmax = 360*pi/180;
for i = 1:NP
    q(i)=qmin+(i-1)*(qmax-qmin)/(NP-1);
    %-------------------------------------------------------------
    %% análise geral de posição
    %-------------------------------------------------------------
    A(i)=atan(R*sin(q(i))/(C+R*cos(q(i))));
    if q(i)==0
        B(i)=(C+R*cos(q(i))/cos(A(i)));
    else
        B(i)=R*sin(q(i))/sin(A(i));
    end
    %-------------------------------------------------------------
    %% determinante da matriz Jacobiana
    %-------------------------------------------------------------
    detJ(i) = -B(i);
    %-------------------------------------------------------------
    %% análise geral de velocidade
    %-------------------------------------------------------------
    Ka(i) = R/B(i)*cos(A(i)-q(i));
    Kb(i) = R*sin(A(i)-q(i));
    %-------------------------------------------------------------
    %% análise geral de aceleração
    %-------------------------------------------------------------
    La(i) = -2*Ka(i)*Kb(i)/B(i)+R/B(i)*sin(A(i)-q(i));
    Lb(i) = Ka(i)^2*B(i)-R*cos(A(i)-q(i));
    %-------------------------------------------------------------
    %% análise de pontos de interesse - posição
    %-------------------------------------------------------------
    Xp(i) = Up*cos(A(i)) - Vp*sin(A(i));
    Yp(i) = Up*sin(A(i)) + Vp*cos(A(i));
    %-------------------------------------------------------------
    %% análise de pontos de interesse - velocidade
    %-------------------------------------------------------------
    Kpx(i) = -Up*sin(A(i))*Ka(i) - Vp*cos(A(i))*Ka(i);
    Kpy(i) =  Up*cos(A(i))*Ka(i) - Vp*sin(A(i))*Ka(i);
    %-------------------------------------------------------------
    %% análise de pontos de interesse - aceleração
    %-------------------------------------------------------------
    Lpx(i) = -Up*cos(A(i))*Ka(i)^2-Up*sin(A(i))*La(i) +...
                       Vp*sin(A(i))*Ka(i)^2-Vp*cos(A(i))*La(i);
    Lpy(i) = -Up*sin(A(i))*Ka(i)^2+Up*cos(A(i))*La(i) -...
                       Vp*cos(A(i))*Ka(i)^2-Vp*sin(A(i))*La(i);
    %% velocidade 
end
%-------------------------------------------------------------
%% gráficos da análise geral
%-------------------------------------------------------------
figure('Name','ANÁLISE GERAL DO MECANISMO DO EXEMPLO 2.1')    % título na moldura superior
% gráfico 1
subplot(3,2,1)
plot(q*180/pi,A*180/pi,...                       
    'LineWidth',2)                % espessura da linha
grid on                           % linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('A(q) (^o)')               % rótulo do eixo y
title('Posição da Biela')         % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x

% gráfico 2
subplot(3,2,2)
plot(q*180/pi,B,...
    'r',...                       % cor da curva
    'LineWidth',2)                % espessura da linha
grid on                           % insere linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('B(q) (mm)')               % rótulo do eixo y
title('Posição do Pistão')        % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
ylim([0 inf])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 3
subplot(3,2,3)
plot(q*180/pi,Ka,...                       
    'LineWidth',2)                % espessura da linha
grid on                           % linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Ka(q)')               % rótulo do eixo y
title('Coeficiente de velocidade Ka')         % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 4
subplot(3,2,4)
plot(q*180/pi,Kb,...
    'r',...                       % cor da curva
    'LineWidth',2)                % espessura da linha
grid on                           % insere linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Kb(q)')               % rótulo do eixo y
title('Coeficiente de velocidade Kb')        % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 5
subplot(3,2,5)
plot(q*180/pi,La,...                       
    'LineWidth',2)                % espessura da linha
grid on                           % linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('La(q)')               % rótulo do eixo y
title('Derivada do coeficiente de velocidade La')         % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 6
subplot(3,2,6)
plot(q*180/pi,Lb,...
    'r',...                       % cor da curva
    'LineWidth',2)                % espessura da linha
grid on                           % insere linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Lb(q)')               % rótulo do eixo y
title('Derivada do coeficiente de velocidade Lb')        % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x

%-------------------------------------------------------------
%% gráfico do determinante da matriz jacobiana
%-------------------------------------------------------------
figure('Name','DETERMINANTE DA MATRIZ JACOBIANA')    % título na moldura superior
plot(q*180/pi,detJ,...                       
    'LineWidth',2)                % espessura da linha
grid on                           % linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('detJ(q)')               % rótulo do eixo y
title('Determinante da matriz Jacobiana')         % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
ylim([-inf 0])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x

%-------------------------------------------------------------
%% gráfico da posição do ponto de interesse
%-------------------------------------------------------------
figure('Name','POSIÇÃO DO PONTO DE INTERESSE')    % título na moldura superior
plot(Xp,Yp,...                       
    'LineWidth',2)                % espessura da linha
grid on                           % linhas de grade
xlabel('X (mm)')                  % rótulo do eixo x
ylabel('Y (mm)')               % rótulo do eixo y
title('Posição do ponto de interesse')         % título do gráfico
axis equal
xlim([0 C1*1.1])
%ylim([-3 3])
%-------------------------------------------------------------
%% gráficos da análise do ponto de interesse
%-------------------------------------------------------------
figure('Name','ANÁLISE DE PONTO DE INTERESSE DO MECANISMO DO EXEMPLO 2.1')  
% gráfico 1
subplot(3,2,1)
plot(q*180/pi,Xp,...                       
    'LineWidth',2)                % espessura da linha
grid on                           % linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Xp(q) (mm)')               % rótulo do eixo y
title('Coordenada X')         % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 2
subplot(3,2,2)
plot(q*180/pi,Yp,...
    'r',...                       % cor da curva
    'LineWidth',2)                % espessura da linha
grid on                           % insere linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Yp(q) (mm)')               % rótulo do eixo y
title('Coordenada Y')        % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 3
subplot(3,2,3)
plot(q*180/pi,Kpx,...                       
    'LineWidth',2)                % espessura da linha
grid on                           % linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Kpx(q)')               % rótulo do eixo y
title('Coeficiente de velocidade Kpx')         % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 4
subplot(3,2,4)
plot(q*180/pi,Kpy,...
    'r',...                       % cor da curva
    'LineWidth',2)                % espessura da linha
grid on                           % insere linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Kpy(q)')               % rótulo do eixo y
title('Coeficiente de velocidade Kpy')        % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 5
subplot(3,2,5)
plot(q*180/pi,Lpx,...                       
    'LineWidth',2)                % espessura da linha
grid on                           % linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Lpx(q)')               % rótulo do eixo y
title('Derivada do coeficiente de velocidade Lpx')         % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x
% gráfico 6
subplot(3,2,6)
plot(q*180/pi,Lpy,...
    'r',...                       % cor da curva
    'LineWidth',2)                % espessura da linha
grid on                           % insere linhas de grade
xlabel('q (^o)')                  % rótulo do eixo x
ylabel('Lpy(q)')               % rótulo do eixo y
title('Derivada do coeficiente de velocidade Lpy')        % título do gráfico
xlim([qmin*180/pi qmax*180/pi])
set(gca,'XTick',qmin*180/pi:45:qmax*180/pi) % marcas do eixo x



%-------------------------------------------------------------
% FIM
%-------------------------------------------------------------
