% Aplicaçao do formalismo de Lagrange para a deduçao das equaçoes
% diferenciais descritivas do sistema em causa.

% Inicializaçao do array dos resultados finais (binarios + factores
% aditivos usados para o calculo do binario)
% T=[T1   d/dDtheta1(EL)   d/dt(d/dDtheta1(EL))   d/dtheta1(EL)   Atrito1;
%    T2   d/dDtheta2(EL)   d/dt(d/dDtheta2(EL))   d/dtheta2(EL)   Atrito2;
%    ... ]
T=[];
% Binarios externos aplicados nas juntas -> T_ext=[Text1; Text2; ...]
T_ext=zeros(N,1);

% Deduçao de binario para cada junta i=1..N
for i=1:N

    % Calculo da derivada de EL em ordem a Dtheta_i -> d/dDthetai(EL))
    dEL_dDtheta=simple(simplify(diff(EL,theta(i,2))));
    
    % Atraves da expansao pela regra da cadeia calcular a derivada em ordem
    % ao tempo do resultado anterior -> d/dt (d/dDthetai (EL))
    dt_dEL_dDtheta=0;
    for j=1:N
        dt_dEL_dDtheta=dt_dEL_dDtheta + diff(dEL_dDtheta,theta(j,1))*theta(j,2)+diff(dEL_dDtheta,theta(j,2))*theta(j,3);
    end
    % Simplificaçao do resultado
    dt_dEL_dDtheta=simple(simplify(dt_dEL_dDtheta));
    
    % Calculo da derivada de EL em ordem a theta_i -> d/dthetai (EL)
    dEL_dtheta=simple(simplify(diff(EL,theta(i,1))));

    % Calculo dos atritos (viscoso+flexivel)
    atrito=simple(simplify(atr(i,1)*theta(i,1) + atr(i,2)*theta(i,2)));
    
    % Determinaçao da expressao final do binario para esta junta (i)
    bin=simple(simplify(dt_dEL_dDtheta-dEL_dtheta+atrito));

    % Actualizaçao do array T com os diversos resultados
    T=[T; bin dEL_dDtheta dt_dEL_dDtheta dEL_dtheta atrito];
    
    % Jacobiano transposto para a determinaçao do binario resultante da força
    % externa:
    J_trans=jacobian_trans(xyf(i,:),theta);
    % Soma do binario externo aplicado nas juntas a partir do Jacobiano
    T_ext=T_ext+J_trans*Fxy(i,:).';
end

% Simplificaçao das expressoes do binario externo
T_ext=simple(simplify(T_ext));

% Actualizaçao do binario total
T=[T T_ext];
T(:,1)=simple(simplify(T(:,1)-T_ext));