% Inicializaçao das variaveis simbolicas para o estudo do pendulo composto

syms g            % Variaveis simbolicas t (tempo) e g (aceleraçao da gravidade)

% === INICIALIZAÇAO DOS ARRAYS DE VARIAVEIS SIMBOLICAS =====================
theta=[];           % Variaveis simbolicas de theta/Dtheta/D2theta -> [theta1 Dtheta1 D2theta1; theta2 Dtheta2 D2theta2; ...]
m=[];               % Massas dos elos -> [m1; m2; ...]
r=[];               % Distancia da origem de cada elo ao seu centro de massa -> [r1; r2; ...]
rf=[];              % Distancia da origem de cada elo ao ponto de aplicaçao de uma força externa -> [rf1; rf2; ...]
L=[];               % Comprimento total de cada elo -> [L1; L2; ...]
atr=[];             % Coeficientes de atrito -> [k1 c1; k2 c2; ...]
xy=[];              % Cinematica directa dos centros de massa de cada elo -> [x1 y1; x2 y2; ...]
xyf=[];             % Cinematica directa dos pontos de aplicaçao de uma força externa em cada elo -> [xf1 yf1; xf2 yf2; ...]
Fxy=[];             % Força externa aplicada a cada elo -> [Fx1 Fy1; Fx2 Fy2; ...]
% =========================================================================

% Inicializaçao das variaveis simbolicas para os N graus de liberdade
for i=1:N
    % === Conversao do id do elo (numero) para string =====================
    % --- Esta string sera usada para a criaçao das variaveis simbolicas --
    str_numb=[];
    numb=i;
    while (numb~=0)
        str_numb=[char('0'+rem(numb,10)) str_numb];
        numb=floor(numb/10);
    end
    % =====================================================================
    
    % === Criaçao das variaveis simbolicas ================================
    m=[m; sym(['m' str_numb])];                             % Massa
    r=[r; sym(['r' str_numb])];                             % Distancia da origem ao CM
    rf=[rf; sym(['rf' str_numb])];                          % Distancia da origem ao ponto de aplicaçao de uma força externa
    L=[L; sym(['L' str_numb])];                             % Comprimento do elo
    atr=[atr; sym(['k' str_numb]) sym(['c' str_numb])];     % Coeficientes de atrito
    
    % Variaveis simbolicas theta, Dtheta e D2theta para cada elo
    theta=[theta; sym(['theta' str_numb]) sym(['Dtheta' str_numb]) sym(['D2theta' str_numb])];
    
    % Cinematica directa generalizada para N elos
    % Estas expressoes tem de ser deduzidas manualmente e inseridas aqui:
    
    % Usando angulos absolutos
%     xy=[xy; r(i)*sin(theta(i,1)) -r(i)*cos(theta(i,1))];
%     for k=1:i-1
%         xy(i,:)=xy(i,:)+[L(k)*sin(theta(k,1)) -L(k)*cos(theta(k,1))];
%     end

    % Usando angulos relativos ao da junta anterior
%     xy=[xy; r(i)*sin(sum(theta(1:i))) -r(i)*cos(sum(theta(1:i)))];
%     xyf=[xyf; rf(i)*sin(sum(theta(1:i))) -rf(i)*cos(sum(theta(1:i)))];
%     for k=1:i-1
%         xy(i,:)=xy(i,:)+[L(k)*sin(sum(theta(1:k))) -L(k)*cos(sum(theta(1:k)))];
%         xyf(i,:)=xyf(i,:)+[L(k)*sin(sum(theta(1:k))) -L(k)*cos(sum(theta(1:k)))];
%     end

    % Angulos relativos a junta anterior
    % Aqui considera-se theta=0 no ponto (x>0 e y=0)
    xy=[xy; r(i)*cos(sum(theta(1:i))) r(i)*sin(sum(theta(1:i)))];
    xyf=[xyf; rf(i)*cos(sum(theta(1:i))) rf(i)*sin(sum(theta(1:i)))];
    for k=1:i-1
        xy(i,:)=xy(i,:)+[L(k)*cos(sum(theta(1:k))) L(k)*sin(sum(theta(1:k)))];
        xyf(i,:)=xyf(i,:)+[L(k)*cos(sum(theta(1:k))) L(k)*sin(sum(theta(1:k)))];
    end

    % Força externa aplicada a cada um dos CM
    Fxy=[Fxy; sym(['Fx' str_numb]) sym(['Fy' str_numb])];
    % =====================================================================
end
