% calculo da expressão do binário

syms g t
syms l1 l2 l3
syms r1 r2 r3
syms m1 m2 m3
teta1 = sym('teta1(t)');
teta2 = sym('teta2(t)');
teta3 = sym('teta3(t)');


    % calculo dos centros de massa no espaço cartesiano
    
    x1 = r1*cos(teta1);
    y1 = r1*sin(teta1);
    
    x2 = l1*cos(teta1) + r2*cos(teta1+teta2);
    y2 = l1*sin(teta1) + r2*sin(teta1+teta2);

    x3 = l1*cos(teta1) + l2*cos(teta1+teta2) + r3*cos(teta1+teta2+teta3);
    y3 = l1*sin(teta1) + l2*sin(teta1+teta2) + r3*sin(teta1+teta2+teta3);
    
    
    
    % cálculo da velocidade linear dos centros de massa
    
    vx1 = diff(x1);
    vy1 = diff(y1);
    
    vx2 = diff(x2);
    vy2 = diff(y2);
    
    vx3 = diff(x3);
    vy3 = diff(y3);
    
    
    
    % cálculo das energias potencial e cinérica
    
    Ep1 = m1*g*y1;
    Ep2 = m2*g*y2;
    Ep3 = m3*g*y3;
    Ep = Ep1 + Ep2 + Ep3;
    
    Ec1 = 1/2*m1*(vx1^2+vy1^2);
    Ec2 = 1/2*m2*(vx2^2+vy2^2);
    Ec3 = 1/2*m3*(vx3^2+vy3^2);
    Ec = Ec1 + Ec2 + Ec3;
    
    
    
    % Função de Lagrange
    U = Ec-Ep;
    U_aux = subs(U, {diff('teta1(t)', t), diff('teta2(t)', t), diff('teta3(t)', t) },{'omega1', 'omega2', 'omega3'})
    U_aux = subs(U_aux, {'teta1(t)', 'teta2(t)', 'teta3(t)'}, {'teta1', 'teta2', 'teta3'});

    T1 = diff(U_aux, 'omega1'); subs(T1, 'omega1', 'teta1(t)\, 'omega1(t)');
    T1 = diff(T1)-diff(U,'teta1(t)');
    bin2 = diff(diff(U, 'omega2'))-diff(U,'teta2');
    subs(bin1, {diff(teta1), diff(teta1, 2)}, {'omega1', 'alfa1'})
 