//=============================================================
// Mecanismos
// Prof. Jorge Luiz Erthal
// Análise geral e de ponto de interesse do exemplo 2.1
// Solução analítica 
//=============================================================
clear
xdel(winsid())
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);
end

//-------------------------------------------------------------
//// gráficos da análise geral
//-------------------------------------------------------------
scf(1)
// gráfico 1
subplot(3,2,1)
plot(q*180/%pi,A*180/%pi,'LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)')
ylabel('A(q) (^o)')
title('Posição da Biela')
// gráfico 2
subplot(3,2,2)
plot(q*180/%pi,B,'r','LineWidth',2)
set(gca(),"grid",[1 1])
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
// gráfico 3
subplot(3,2,3)
plot(q*180/%pi,Ka,'LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)') 
ylabel('Ka(q)') 
title('Coeficiente de velocidade Ka') 
// gráfico 4
subplot(3,2,4)
plot(q*180/%pi,Kb,'r','LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)')
ylabel('Kb(q)') 
title('Coeficiente de velocidade Kb') 
// gráfico 5
subplot(3,2,5)
plot(q*180/%pi,La,'LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)')
ylabel('La(q)') 
title('Derivada do coeficiente de velocidade La')
// gráfico 6
subplot(3,2,6)
plot(q*180/%pi,Lb,'r','LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)') 
ylabel('Lb(q)') 
title('Derivada do coeficiente de velocidade Lb')
//-------------------------------------------------------------
//// gráfico do determinante da matriz jacobiana
//-------------------------------------------------------------
scf(2)
plot(q*180/%pi,detJ,'LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)') 
ylabel('detJ(q)')
title('Determinante da matriz Jacobiana') 
//-------------------------------------------------------------
//// gráfico da posição do ponto de interesse
//-------------------------------------------------------------
scf(3)
plot(Xp,Yp,'LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('X (mm)') 
ylabel('Y (mm)')
title('Posição do ponto de interesse') 
a=get("current_axes")//obtendo o manipulador dos novos eixos criados
a.axes_visible="on"; // tornando os eixos visíveis
a.font_size=3; //ajustando o tamanho da fonte dos rótulos dos tiques
a.x_location="middle"; //ajustando a posição do eixo x
a.data_bounds=[0,-1,0;C1+R,1,0]; //ajustando os valores limites para as coordenadas x,y e z
a.sub_tics=[5,0];
a.labels_font_color=5;
a.grid=[2,2];
a.box="off";





//-------------------------------------------------------------
//// gráficos da análise do ponto de interesse
//-------------------------------------------------------------
scf(4)
// gráfico 1
subplot(3,2,1)
plot(q*180/%pi,Xp,'LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)') 
ylabel('Xp(q) (mm)')
title('Coordenada X')
// gráfico 2
subplot(3,2,2)
plot(q*180/%pi,Yp,'r','LineWidth',2) 
set(gca(),"grid",[1 1])
xlabel('q (^o)')
ylabel('Yp(q) (mm)')
title('Coordenada Y') 
// gráfico 3
subplot(3,2,3)
plot(q*180/%pi,Kpx,'LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)')
ylabel('Kpx(q)')
title('Coeficiente de velocidade Kpx')
// gráfico 4
subplot(3,2,4)
plot(q*180/%pi,Kpy,'r','LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)')
ylabel('Kpy(q)')
title('Coeficiente de velocidade Kpy')
// gráfico 5
subplot(3,2,5)
plot(q*180/%pi,Lpx,'LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)')
ylabel('Lpx(q)')
title('Derivada do coeficiente de velocidade Lpx') 
// gráfico 6
subplot(3,2,6)
plot(q*180/%pi,Lpy,'r','LineWidth',2)
set(gca(),"grid",[1 1])
xlabel('q (^o)')
ylabel('Lpy(q)')
title('Derivada do coeficiente de velocidade Lpy')
//-------------------------------------------------------------
// FIM
//-------------------------------------------------------------
