function dydt=func_bielamanivela(t,y)
%-----------------------------------------------------------------------
% mecanismo baseado na Figura 7.2
%-------------------------------------------------------------
%% dados cinemática
%-------------------------------------------------------------
R  = 0.25;      % Comprimento da manivela (m)                       <**>
C  = 0.00;      % Distancia do pistao ao centro da manivela (m)     <**>
L  = 0.45;      % comprimento da biela (m)                          <**>
C1 = 1.00;      % Extremidade fixa da mola (m)                      <**>
%-------------------------------------------------------------
% pontos de interesse
up2=R/2;             % centro de massa da manivela                   <**>
vp2=0;               
up3=L/2;             % centro de massa da biela                      <**>
vp3=0;               
%-------------------------------------------------------------
% número de pontos e distribuição no intervalo
NP = 361;           % numero de pontos calculados
qmin = 0*pi/180;    % valor de q no inicio do intervalo             <**>
qmax = 360*pi/180;  % valor de q no fim do intervalo                <**>
%-----------------------------------------------------------------------
%% dados estática
%-------------------------------------------------------------
m2=0.01;    % massa da manivela (kg)
m3=0.01;    % massa da biela (kg)
m4=0.4;     % massa do pistao (kg)
%-------------------------------------------------------------
km = 2;     % constante da mola (N/m)
S0 = 0.2;   % comprimento livre da mola (m)
%-------------------------------------------------------------
F = 1;      % Força externa aplicada (N)
%-----------------------------------------------------------------------
grav = 9.81;        % aceleração da gravidade (m/s^2)
%-----------------------------------------------------------------------
%% Dados Dinâmica
%-----------------------------------------------------------------------
T=0.303926;      % torque aplicado para equilíbrio em q = 35 graus (N.m)
%-----------------------------------------------------------------------
B=0.2;        % coeficiente de amortecimento do amortecedor   (N.s/m)
%-----------------------------------------------------------------------
% inércias
I2=m2*R^2/12;   % momento de inercia da manivela                (kg.m^2)
I3=m3*L^2/12;   % momento de inercia da biela                   (kg.m^2)
%-----------------------------------------------------------------------
%% posição
q=y(1);             % posição q
ef=10^(-6);         % tolerancia para o metodo de Newton Raphson
a0 = 0*pi/180;      % valores iniciais das variaveis secundarias
x0 = (R+L);
s=[a0; x0];
f=[1;1];            %funçao f cuja norma e maior que ef

while norm(f) > ef
    a=s(1);
    x=s(2);
    f = [R*cos(q)+L*cos(a)-x;
        R*sin(q)-L*sin(a)-C];
    J = [-L*sin(a) -1;
        -L*cos(a)  0];
    s=s-inv(J)*f;
end
% distribuição do resultado nas variáveis secundárias
a=s(1);
x=s(2);
%-----------------------------------------------------------------------
%% coeficientes de velocidade
% matriz jacobiana
J=[-L*sin(a)   -1;
    -L*cos(a)    0];
% derivada das equacoes cinematicas em relacao a q
dfdq=[-R*sin(q);
    R*cos(q)];
% coeficientes de velocidade
k=-inv(J)*dfdq;
ka=k(1);
kx=k(2);
%-----------------------------------------------------------------------
%% coeficientes de aceleração
% derivada da matriz
dJdq=[-ka*L*cos(a) 0;
       ka*L*sin(a) 0];
d2fdq2=[-R*cos(q);
        -R*sin(q)];
l=-inv(J)*(dJdq*k+d2fdq2);
la=l(1);
lx=l(2);
%-----------------------------------------------------------------------
%% análise de posição do CM da manivela                                <**>
Loc2=[up2;
      vp2];
Ori2=[0;
      0];
Rot2=[cos(q) -sin(q)
      sin(q)  cos(q)];
Glob2=Ori2+Rot2*Loc2;
xp2=Glob2(1);
yp2=Glob2(2);
%-------------------------------------------------------------------
%% coeficientes de velocidade do CM da manivela         <**>
dOri2dq = [0;
           0];
dRot2dq = [-sin(q)  -cos(q);
            cos(q)  -sin(q)];
Kp2=dOri2dq + dRot2dq*Loc2;
kpx2=Kp2(1);
kpy2=Kp2(2);
%-------------------------------------------------------------------
%% derivadas dos coeficientes de velocidade do CM da manivela         <**>
ddOri2ddq = [0;
             0];
ddRot2ddq = [-cos(q)  sin(q);
             -sin(q) -cos(q)];
Lp2=ddOri2ddq + ddRot2ddq*Loc2;
lpx2=Lp2(1);
lpy2=Lp2(2);
%-------------------------------------------------------------------
%% análise de posição do CM da biela                                <**>
Loc3=[up3;
      vp3];
Ori3=[R*cos(q);
      R*sin(q)];
Rot3=[ cos(a) sin(a)
      -sin(a) cos(a)];
Glob3=Ori3+Rot3*Loc3;
xp3=Glob3(1);
yp3=Glob3(2);
%-------------------------------------------------------------------
%% coeficientes de velocidade do CM da biela         <**>
dOri3dq = [-R*sin(q);
            R*cos(q)];
dRot3dq = [-sin(a)*ka  cos(a)*ka;
           -cos(a)*ka -sin(a)*ka];
Kp3=dOri3dq + dRot3dq*Loc3;
kpx3=Kp3(1);
kpy3=Kp3(2);
%-------------------------------------------------------------------
%% derivadas dos coeficientes de velocidade do CM da biela         <**>
ddOri3ddq = [-R*cos(q);
             -R*sin(q)];
ddRot3ddq = [-cos(a)*ka^2-sin(a)*la -sin(a)*ka^2+cos(a)*la;
              sin(a)*ka^2-cos(a)*la -cos(a)*ka^2-sin(a)*la];
Lp3=ddOri3ddq + ddRot3ddq*Loc3;
lpx3=Lp3(1);
lpy3=Lp3(2);
%-------------------------------------------------------------------
%% inércia generalizada
Iq=m2*(kpx2^2+kpy2^2) + I2 + m3*(kpx3^2+kpy3^2) + I3*ka^2 + m4*kx^2;
%-----------------------------------------------------------------------
%% coeficiente centrípeto
Cq=m2*(kpx2*lpx2+kpy2*lpy2)+m3*(kpx3*lpx3+kpy3*lpy3)+I3*ka*la+m4*kx*lx;
%-----------------------------------------------------------------------
%% forças conservativas (mola e gravidade)
dVdq = -km*kx*(C1-x-S0)+m2*grav*kpy2+m3*grav*kpy3;
%-----------------------------------------------------------------------
%% forças não-conservativas (amortecedor e forças externas)
Qnc = F*kx + T - B*kx^2*y(2);
%-----------------------------------------------------------------------
%% montagem da equação do movimento
dydt=[y(2); (Qnc - Cq*y(2)^2 - dVdq)/Iq];
%-----------------------------------------------------------------------

