function dydt=func_ex7_01_DINAMICA(t,y)
%-------------------------------------------------------------
%% dados numéricos
%-------------------------------------------------------------
% dimensões dos elos
C1=0.3;
C2=0.2;
%-------------------------------------------------------------
% pontos de interesse
% ponto de aplicação das forças
UP1=C1;
VP1=C2;
% centro de massa do elo triangular
UP2=C1*2/3;
VP2=C2/3;
%-------------------------------------------------------------
% ESTÁTICA
grav = 9.81;    % m/s^2
m1=2;       % kg
m2=2;       % kg
m3=5;       % kg
%-------------------------------------------------------------
kmola = 50; % N/m
S0 = C1/2;  % m
%-------------------------------------------------------------
Fx = 60;    % N
Fy = 40;    % N
%-------------------------------------------------------------
% DINÂMICA
Torque=-5.773;  % N.m, na posição q = 0.1 m
B=10;           % N.s/m
I3=0.05;        % kg.m^2
%-------------------------------------------------------------
%% Análise
% posição
q=y(1);
%-------------------------------------------------------------
ZERO=10^(-6);     % tolerancia para o metodo de Newton Raphson
a0 = 0*pi/180;   % valores iniciais das variaveis secundarias
x0 = 0;
s=[a0; x0];
f=[1;1];          %funçao f cuja norma e maior que ZERO
while norm(f) > ZERO
    A = s(1);
    X = s(2);
    f = [ C1*cos(A)-X;
        C1*sin(A)-q];
    jaco=[ -C1*sin(A), -1;
        C1*cos(A),  0];
    s=s-jaco\f;
end
% distribuição do resultado nas variáveis secundárias
A=s(1);
X=s(2);
%-------------------------------------------------------------
%% velocidade
%-------------------------------------------------------------
Q = [ 0;
    -1];
K = -jaco\Q;
Ka = K(1);
Kx = K(2);
%-------------------------------------------------------------
%% aceleração
%-------------------------------------------------------------
dJdq = [-C1*cos(A) 0;
    -C1*sin(A) 0];
dQdq = [0;
    0];
L = -jaco\(dJdq*K+dQdq);
La = L(1);
Lx = L(2);
%-------------------------------------------------------------
%% Ponto de interesse 1 - ponto de aplicação das forças
%-------------------------------------------------------------
UP=UP1;
VP=VP1;
%-------------------------------------------------------------
% origem do sistema local
O = [0;
    q];
% matriz de rotação
rot = [ cos(A)  sin(A)
    -sin(A)  cos(A)];
% coordenadas locais do ponto de interesse
Ploc = [UP; VP];
% coordenadas globais do ponto de interesse
Pglob = O+rot*Ploc;
Xp1 = Pglob(1);
Yp1 = Pglob(2);
%----------------------------------------------------
% Velocidade
% derivada da posição da origem do sistema local
dOdq = [0;
    1];
% derivada da matriz de rotação
dRdq = [-sin(A)*Ka  cos(A)*Ka;
    -cos(A)*Ka -sin(A)*Ka];
% coeficientes de velocidade
Kp = dOdq + dRdq*Ploc;
Kpx1 = Kp(1);
Kpy1 = Kp(2);
%-------------------------------------------------------------
% Aceleração
%-------------------------------------------------------------
% derivada segunda da posição da origem do sistema local
ddOddq = [0;
    0];
% derivada segunda da matriz de rotação
ddRddq = [-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];
% derivada dos coeficientes de velocidade
Lp = ddOddq + ddRddq*Ploc;
Lpx1 = Lp(1);
Lpy1 = Lp(2);
%-------------------------------------------------------------
%% Ponto de interesse 2
%-------------------------------------------------------------
UP=UP2;
VP=VP2;
%-------------------------------------------------------------
% origem do sistema local
O = [0;
    q];
% matriz de rotação
rot = [ cos(A)  sin(A)
    -sin(A)  cos(A)];
% coordenadas locais do ponto de interesse
Ploc = [UP; VP];
% coordenadas globais do ponto de interesse
Pglob = O+rot*Ploc;
Xp2 = Pglob(1);
Yp2 = Pglob(2);
%----------------------------------------------------
% Velocidade
% derivada da posição da origem do sistema local
dOdq = [0;
    1];
% derivada da matriz de rotação
dRdq = [-sin(A)*Ka  cos(A)*Ka;
    -cos(A)*Ka -sin(A)*Ka];
% coeficientes de velocidade
Kp = dOdq + dRdq*Ploc;
Kpx2 = Kp(1);
Kpy2 = Kp(2);
%-------------------------------------------------------------
% Aceleração
%-------------------------------------------------------------
% derivada segunda da posição da origem do sistema local
ddOddq = [0;
    0];
% derivada segunda da matriz de rotação
ddRddq = [-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];
% derivada dos coeficientes de velocidade
Lp = ddOddq + ddRddq*Ploc;
Lpx2 = Lp(1);
Lpy2 = Lp(2);
%-------------------------------------------------------------------
%% inércia generalizada
Iq=m1+m2*Kx^2+m3*(Kpx2^2+Kpy2^2)+I3*Ka^2;
%-----------------------------------------------------------------------
%% coeficiente centrípeto
Cq=m2*Kx*Lx+m3*(Kpx2*Lpx2+Kpy2*Lpy2)+I3*Ka*La;
%-----------------------------------------------------------------------
%% forças conservativas (mola e gravidade)
dVdq = kmola*(X-S0)*Kx+m1*grav+m3*grav*Kpy2;
%-----------------------------------------------------------------------
%% forças não-conservativas (amortecedor e forças externas)
Qnc = Fx*Kpx1 + Fy*Kpy1 - B*y(2) - Torque*Ka;
%-----------------------------------------------------------------------
%% montagem da equação do movimento
dydt=[y(2); (Qnc - Cq*y(2)^2 - dVdq)/Iq];
%-----------------------------------------------------------------------

