%% Bola en barra "catapulta": ecuaciones del movimiento (2GL, 1GL forzado)
% © 2021, Antonio Sala Piqueras. Universitat Politècnica de València. Todos los derechos reservados.
% Presentaciones en vídeo:
%  http://personales.upv.es/asala/YT/V/ballbar2GL.html , http://personales.upv.es/asala/YT/V/ballbar1GL.html .
% Este código funcionó correctamente con Matlab R2021b
% Objetivo: obtener las ecuaciones de movimiento de un sistema "bola agujereada ensartada en una barra que gira", con Euler-Lagrange.
% Tabla de Contenidos
% Cinematica forzada
% Newton: balance de fuerzas
% Euler-Lagrange 1GL
% Modelado como sistema 2GL con par externo
%     Cinemática
%     Energías y ecs. movimiento Euler-Lagrange
%
% *Nota: no vamos a considerar energía potencial de la barra. Suponemos simétrica alrededor del eje de rotación, de modo que su CdG no se mueva. Si no, deberíamos incluir la masa de la barra y la posición de su CdG en la energía potencial en el modelo 2GL. Esto NO es relevante en el caso 1GL con la velocidad de la barra forzada, obviamente, dado que no importa el detalle de peso, inercia, etc. de la barra en el movimiento de la bola ensartada si "alguien" hace las fuerzas necesarias para que gire con la velocidad prefijada \omega que queramos.

syms rho t real
syms M g omega phi_0 real
%g=0 %plano horizontal; g=9.81, plano vertical

% Analicemos una bola ensartada en una barra que gira a velocidad constante (1GL, el otro GL sería "forzado").
%% Cinematica forzada

x=rho*cos(omega*t+phi_0); y=rho*sin(omega*t+phi_0); %bola en barra sin sujección radial... 
syms v_rho a_rho real %velocidades y aceleraciones 
r=[x;y] %vector de posición

% El único grado de libertad será q\equiv \rho.

q=rho; v_q=v_rho; a_q=a_rho;
radial_vector=diff(r,q)

% Implementemos cálculo de velocidad \frac{dr}{dt}=\frac{\partial r}{\partial q}\dot q+\frac{\partial r}{\partial t} = [\frac{\partial r}{\partial q}\ \frac{\partial r}{\partial t} ]\cdot \left(\begin{array}{c}\dot q \\ 1\end{array}\right)

veloc=jacobian(r,[q; t])*[v_q; 1] 
veloc_radial=simplify(radial_vector'*veloc)

% y, del mismo modo, la regla de la cadena dice que, con v(q(t),\dot q(t),t), su derivada temporal es: \frac{dv}{dt}=\frac{\partial v}{\partial q}\dot q+\frac{\partial v}{\partial \dot q}\ddot q+\frac{\partial v}{\partial t}

acel=simplify(jacobian(veloc,[q; v_q; t])*[v_q; a_q; 1])
acel_radial=simplify(radial_vector'*acel)

%% Newton: balance de fuerzas
% Balance tangencial:

Gravedad_radial=radial_vector'*[0;-M*g]
a_qNewton=simplify(solve(radial_vector'*acel==Gravedad_radial/M,a_q))

%% Euler-Lagrange 1GL

T=simplify(0.5*M*veloc'*veloc) %cinética, no hace falta cinética de la barra fozada
V=M*g*y %potencial
Lagrangian=T-V;
p=diff(Lagrangian,v_q)
dpdt=jacobian(p,[q v_q t])*[v_q; a_q; 1]
EcmovEL1GL= dpdt - diff(Lagrangian,q) ==0
a_q_EulerLagrange=simplify(solve(EcmovEL1GL  , a_q))

%%
%
%% Modelado como sistema 2GL con par externo
%% Cinemática
% Barra de inercia "J" alrededor de eje de rotación.
% Entonces  \rho es un grado de libertad, y su derivada es la v_\rho; \phi es el segundo grado de libertad, y su derivada es la \omega:

syms phi J real
syms  a_phi real
x=rho*cos(phi); y=rho*sin(phi); 
r=[x;y] %vector de posición
q=[rho;phi]; v_q=[v_rho;omega]; a_q=[a_rho;a_phi];

% Regla de la cadena; dado r(q(t)), su derivada: \frac{dr}{dt}=\frac{\partial r}{\partial q}\dot q

veloc=jacobian(r,q)*v_q

% Regla de la cadena; dado v(q(t),\dot q(t)), su derivada: \frac{dv}{dt}=\frac{\partial v}{\partial q}\dot q+\frac{\partial v}{\partial \dot q}\ddot q=[\frac{\partial v}{\partial q} \ \frac{\partial v}{\partial \dot q}]\cdot \left[\begin{array}{c}\dot q\\ \ddot q\end{array}\right]

acel=simplify(jacobian(veloc,[q; v_q])*[v_q; a_q]) %no hace falta para Euler-Lagrange

%% Energías y ecs. movimiento Euler-Lagrange
%

T=simplify(0.5*M*(veloc'*veloc)+0.5*J*omega^2) %cinética
V=M*g*y %potencial [¡¡¡CDG barra centro rotación!!!!]
Lagrangian=T-V 
p=jacobian(Lagrangian,v_q) %momento cinético generalizado

% Regla de la cadena, dado p(q(t),\dot q(t)), su derivada: \frac{dp}{dt}=\frac{\partial p}{\partial q}\dot q+\frac{\partial p}{\partial \dot q}\ddot q=[\frac{\partial p}{\partial q} \ \frac{\partial p}{\partial \dot q}]\cdot \left[\begin{array}{c}\dot q\\ \ddot q\end{array}\right]

dpdt=jacobian(p,[q; v_q ])*[v_q;  a_q]
syms tau real %par sobre barra (externo), trabajo virtual tau*ðphi

% Ecs. Euler Lagrange del movimiento \frac{dp}{dt}- \frac{\partial L}{\partial q}= Q

EcsMovimiento = dpdt - jacobian(Lagrangian,q)' == [0;tau]
a_q_EulerLagrange2GL=solve( EcsMovimiento  , a_q);
a_q_2gl=simplify(a_q_EulerLagrange2GL.a_rho) 
a_phi_2gl=simplify(a_q_EulerLagrange2GL.a_phi)

% Con lo cual, el par \tau necesario para mantener el giro a velocidad constante de la barra guía es:

tau_rFL=simplify(solve(a_phi_2gl==0,tau))

%