#!/usr/bin/env python
# coding: utf-8
# # Polar and Cylindrical Frame of Reference
#
# > Renato Naville Watanabe, Marcos Duarte
# > [Laboratory of Biomechanics and Motor Control](http://pesquisa.ufabc.edu.br/bmclab)
# > Federal University of ABC, Brazil
#
# ## Python setup
# In[39]:
import numpy as np
import sympy as sym
from sympy.plotting import plot_parametric,plot3d_parametric_line
from sympy.vector import CoordSys3D
import matplotlib.pyplot as plt
# from matplotlib import rc
# rc('text', usetex=True)
sym.init_printing()
# Consider that we have the position vector $\bf\vec{r}$ of a particle, moving in a circular path indicated in the figure below by a dashed line. This vector $\bf\vec{r(t)}$ is described in a fixed reference frame as:
#
#
#
# \begin{equation}
# {\bf\vec{r}}(t) = x(t){\bf\hat{i}} + y(t){\bf\hat{j}} + z(t){\bf\hat{k}}
# \label{eq_1}
# \end{equation}
#
#
#
# Naturally, we could describe all the kinematic variables in the fixed reference frame. But in circular motions, is convenient to define a basis with a vector in the direction of the position vector $\bf\vec{r}$. So, the vector $\bf\hat{e_R}$ is defined as:
#
#
# \begin{equation}
# {\bf\hat{e_R}} = \frac{\bf\vec{r}}{\Vert{\bf\vec{r} }\Vert}
# \label{eq_2}
# \end{equation}
#
#
# The second vector of the basis can be obtained by the cross multiplication between $\bf\hat{k}$ and $\bf\hat{e_R}$:
#
#
#
# \begin{equation}
# {\bf\hat{e_\theta}} = {\bf\hat{k}} \times {\bf\hat{e_R}}
# \label{eq_3}
# \end{equation}
#
#
#
# The third vector of the basis is the conventional ${\bf\hat{k}}$ vector.
#
#
# This basis can be used also for non-circular movements. For a 3D movement, the versor ${\bf\hat{e_R}}$ is obtained by removing the projection of the vector ${\bf\vec{r}}$ onto the versor ${\bf\hat{k}}$:
#
#
#
# \begin{equation}
# {\bf\hat{e_R}} = \frac{\bf\vec{r} - ({\bf\vec{r}.{\bf\hat{k}}){\bf\hat{k}}}}{\Vert\bf\vec{r} - ({\bf\vec{r}.{\bf\hat{k}}){\bf\hat{k}}\Vert}}
# \label{eq_4}
# \end{equation}
#
#
#
#
# Note that if the movement is on a plane, the expression above is equal to ${\bf\hat{e_R}} = \frac{\bf\vec{r}}{\Vert{\bf\vec{r} }\Vert}$ since the projection of $\bf\vec{r}$ on the $\bf\hat{k}$ versor is zero.
# ## Time-derivative of the versors ${\bf\hat{e_R}}$ and ${\bf\hat{e_\theta}}$
#
# To obtain the expressions of the velocity and acceleration vectors, it is necessary to obtain the expressions of the time-derivative of the vectors ${\bf\hat{e_R}}$ and ${\bf\hat{e_\theta}}$.
#
# This can be done by noting that:
#
#
#
# \begin{align}
# {\bf\hat{e_R}} &= \cos(\theta){\bf\hat{i}} + \sin(\theta){\bf\hat{j}}\\
# {\bf\hat{e_\theta}} &= -\sin(\theta){\bf\hat{i}} + \cos(\theta){\bf\hat{j}}
# \label{eq_5}
# \end{align}
#
#
# Deriving ${\bf\hat{e_R}}$ we obtain:
#
#
#
# \begin{equation}
# \frac{d{\bf\hat{e_R}}}{dt} = -\sin(\theta)\dot\theta{\bf\hat{i}} + \cos(\theta)\dot\theta{\bf\hat{j}} = \dot{\theta}{\bf\hat{e_\theta}}
# \label{eq_6}
# \end{equation}
#
#
# Similarly, we obtain the time-derivative of ${\bf\hat{e_\theta}}$:
#
#
#
# \begin{equation}
# \frac{d{\bf\hat{e_\theta}}}{dt} = -\cos(\theta)\dot\theta{\bf\hat{i}} - \sin(\theta)\dot\theta{\bf\hat{j}} = -\dot{\theta}{\bf\hat{e_R}}
# \label{eq_7}
# \end{equation}
#
#
# ## Position, velocity and acceleration
# ### Position
#
# The position vector $\bf\vec{r}$, from the definition of $\bf\hat{e_R}$, is:
#
#
#
# \begin{equation}
# {\bf\vec{r}} = R{\bf\hat{e_R}} + z{\bf\hat{k}}
# \label{eq_8}
# \end{equation}
#
#
# where $R = \Vert\bf\vec{r} - ({\bf\vec{r}.{\bf\hat{k}}){\bf\hat{k}}\Vert}$.
# ### Velocity
#
# The velocity vector $\bf\vec{v}$ is obtained by deriving the vector $\bf\vec{r}$:
#
#
#
# \begin{equation}
# {\bf\vec{v}} = \frac{d(R{\bf\hat{e_R}})}{dt} + \dot{z}{\bf\hat{k}} = \dot{R}{\bf\hat{e_R}}+R\frac{d\bf\hat{e_R}}{dt}=\dot{R}{\bf\hat{e_R}}+R\dot{\theta}{\bf\hat{e_\theta}}+ \dot{z}{\bf\hat{k}}
# \label{eq_9}
# \end{equation}
#
# ### Acceleration
#
# The acceleration vector $\bf\vec{a}$ is obtained by deriving the velocity vector:
#
#
#
# \begin{align}
# {\bf\vec{a}} =& \frac{d(\dot{R}{\bf\hat{e_R}}+R\dot{\theta}{\bf\hat{e_\theta}}+\dot{z}{\bf\hat{k}})}{dt}=\\\nonumber
# =&\ddot{R}{\bf\hat{e_R}}+\dot{R}\frac{d\bf\hat{e_R}}{dt} + \dot{R}\dot{\theta}{\bf\hat{e_\theta}} + R\ddot{\theta}{\bf\hat{e_\theta}} + R\dot{\theta}\frac{d{\bf\hat{e_\theta}}}{dt} + \ddot{z}{\bf\hat{k}}=\\\nonumber
# =&\ddot{R}{\bf\hat{e_R}}+\dot{R}\dot{\theta}{\bf\hat{e_\theta}} + \dot{R}\dot{\theta}{\bf\hat{e_\theta}} + R\ddot{\theta}{\bf\hat{e_\theta}} - R\dot{\theta}^2{\bf\hat{e_R}}+ \ddot{z}{\bf\hat{k}} =\\
# =&\ddot{R}{\bf\hat{e_R}}+2\dot{R}\dot{\theta}{\bf\hat{e_\theta}}+ R\ddot{\theta}{\bf\hat{e_\theta}} - {R}\dot{\theta}^2{\bf\hat{e_R}}+ \ddot{z}{\bf\hat{k}} =\\\nonumber
# =&(\ddot{R}-R\dot{\theta}^2){\bf\hat{e_R}}+(2\dot{R}\dot{\theta} + R\ddot{\theta}){\bf\hat{e_\theta}}+ \ddot{z}{\bf\hat{k}}\nonumber
# \label{eq_10}
# \end{align}
#
#
# + The term $\ddot{R}$ is an acceleration in the radial direction.
#
# + The term $R\ddot{\theta}$ is an angular acceleration.
#
# + The term $\ddot{z}$ is an acceleration in the $\bf\hat{k}$ direction.
#
# + The term $-R\dot{\theta}^2$ is the well known centripetal acceleration.
#
# + The term $2\dot{R}\dot{\theta}$ is known as Coriolis acceleration. This term may be difficult to understand intuitively. It appears when there is displacement in the radial and angular directions at the same time.
# ## Important to note
#
# The reader must bear in mind that the use of a different basis to represent the position, velocity or acceleration vectors is only a different representation of the same vector. For example, for the acceleration vector:
#
#
#
# \begin{equation}
# {\bf\vec{a}} = \ddot{x}{\bf\hat{i}}+ \ddot{y}{\bf\hat{j}} + \ddot{z}{\bf\hat{k}}=(\ddot{R}-R\dot{\theta}^2){\bf\hat{e_R}}+(2\dot{R}\dot{\theta} + R\ddot{\theta}){\bf\hat{e_\theta}} + \ddot{z}{\bf\hat{k}}=\dot{\Vert\bf\vec{v}\Vert}{\bf\hat{e}_t}+{\Vert\bf\vec{v}\Vert}^2\Vert{\bf\vec{C}} \Vert{\bf\hat{e}_n}
# \label{eq_11}
# \end{equation}
#
#
# In which the last equality is the acceleration vector represented in the path-coordinate of the particle (see http://nbviewer.jupyter.org/github/BMClab/bmc/blob/master/notebooks/Time-varying%20frames.ipynb).
# ## Example
#
# Consider a particle following the spiral path described below:
#
#
#
# \begin{equation}
# {\bf\vec{r}}(t) = (2\sqrt{t}\cos(t)){\bf\hat{i}}+ (2\sqrt{t}\sin(t)){\bf\hat{j}}
# \label{eq_12}
# \end{equation}
#
# ### Solving numerically
# In[40]:
t = np.linspace(0.01,10,30).reshape(-1,1) #create a time vector and reshapes it to a column vector
R = 2*np.sqrt(t)
theta = t
rx = R*np.cos(t)
ry = R*np.sin(t)
r = np.hstack((rx, ry)) # creates the position vector by stacking rx and ry horizontally
# In[41]:
e_r = r/np.linalg.norm(r, axis=1, keepdims=True) # defines e_r vector
e_theta = np.cross([0,0,1],e_r)[:,0:-1] # defines e_theta vector
# In[42]:
dt = t[1] #defines delta_t
Rdot = np.diff(R, axis=0)/dt #find the R derivative
thetaDot = np.diff(theta, axis=0)/dt #find the angle derivative
v = Rdot*e_r[0:-1,:] +R[0:-1]*thetaDot*e_theta[0:-1,:] # find the linear velocity.
# In[43]:
Rddot = np.diff(Rdot, axis=0)/dt
thetaddot = np.diff(thetaDot, axis=0)/dt
# In[44]:
a = ((Rddot - R[1:-1]*thetaDot[0:-1]**2)*e_r[1:-1,:]
+ (2*Rdot[0:-1]*thetaDot[0:-1] + Rdot[0:-1]*thetaddot)*e_theta[1:-1,:])
# The versors $\bf\hat{e_R}$, $\bf\hat{e_\theta}$ are plotted below at some points of the path.
# In[45]:
from matplotlib.patches import FancyArrowPatch
get_ipython().run_line_magic('matplotlib', 'inline')
plt.rcParams['figure.figsize'] = (8,8)
fig = plt.figure()
ax = fig.add_axes([0,0,1,1])
plt.plot(r[:,0],r[:,1],'.')
for i in np.arange(len(t)-2):
vec1 = FancyArrowPatch(r[i,:],r[i,:]+e_r[i,:],mutation_scale=30,color='r', label='e_r')
vec2 = FancyArrowPatch(r[i,:],r[i,:]+e_theta[i,:],mutation_scale=30,color='g', label='e_theta')
ax.add_artist(vec1)
ax.add_artist(vec2)
plt.xlim((-10,10))
plt.ylim((-10,10))
plt.grid()
plt.legend([vec1, vec2],[r'$\vec{e_r}$', r'$\vec{e_{\theta}}$'])
plt.show()
# The velocity and acceleleration of the particle are plotted below at some points of the path.
# In[46]:
from matplotlib.patches import FancyArrowPatch
get_ipython().run_line_magic('matplotlib', 'inline')
plt.rcParams['figure.figsize'] = (8,8)
fig = plt.figure()
plt.plot(r[:,0],r[:,1],'.')
ax = fig.add_axes([0,0,1,1])
for i in np.arange(len(t)-2):
vec1 = FancyArrowPatch(r[i,:],r[i,:]+v[i,:],mutation_scale=10,color='r')
vec2 = FancyArrowPatch(r[i,:],r[i,:]+a[i,:],mutation_scale=10,color='g')
ax.add_artist(vec1)
ax.add_artist(vec2)
plt.xlim((-10,10))
plt.ylim((-10,10))
plt.grid()
plt.legend([vec1, vec2],[r'$\vec{v}$', r'$\vec{a}$'])
plt.show()
# ### Solved symbolically
# In[47]:
O = sym.vector.CoordSys3D(' ')
t = sym.symbols('t')
# In[48]:
r = 2*sym.sqrt(t)*sym.cos(t)*O.i+2*sym.sqrt(t)*sym.sin(t)*O.j
r
# In[49]:
plot_parametric(r.dot(O.i),r.dot(O.j),(t,0,10))
# In[50]:
e_r = r - r.dot(O.k)*O.k
e_r = e_r/sym.sqrt(e_r.dot(O.i)**2+e_r.dot(O.j)**2+e_r.dot(O.k)**2)
# In[51]:
e_r
# In[52]:
e_theta = O.k.cross(e_r)
e_theta
# In[53]:
from matplotlib.patches import FancyArrowPatch
plt.rcParams['figure.figsize'] = (8,8)
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1])
ax.axis("on")
time = np.linspace(0,10,30)
for instant in time:
vt = FancyArrowPatch([float(r.dot(O.i).subs(t,instant)),float(r.dot(O.j).subs(t,instant))],
[float(r.dot(O.i).subs(t,instant))+float(e_r.dot(O.i).subs(t,instant)),
float(r.dot(O.j).subs(t, instant))+float(e_r.dot(O.j).subs(t,instant))],
mutation_scale=20,
arrowstyle="->",color="r",label='${{e_r}}$')
vn = FancyArrowPatch([float(r.dot(O.i).subs(t, instant)),float(r.dot(O.j).subs(t,instant))],
[float(r.dot(O.i).subs(t, instant))+float(e_theta.dot(O.i).subs(t, instant)),
float(r.dot(O.j).subs(t, instant))+float(e_theta.dot(O.j).subs(t, instant))],
mutation_scale=20,
arrowstyle="->",color="g",label='${{e_{theta}}}$')
ax.add_artist(vn)
ax.add_artist(vt)
plt.xlim((-10,10))
plt.ylim((-10,10))
plt.legend(handles=[vt,vn],fontsize=20)
plt.grid()
plt.show()
# In[54]:
R = 2*sym.sqrt(t)
# In[55]:
Rdot = sym.diff(R,t)
Rddot = sym.diff(Rdot,t)
Rddot
# In[56]:
v = Rdot*e_r + R*e_theta
# In[57]:
v
# In[58]:
a = (Rddot - R)*e_r + (2*Rdot*1+0)*e_theta
aCor = 2*Rdot*1*e_theta
aCor
# In[59]:
a
# In[60]:
from matplotlib.patches import FancyArrowPatch
plt.rcParams['figure.figsize'] = (8,8)
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1])
ax.axis("on")
time = np.linspace(0.1,10,30)
for instant in time:
vt = FancyArrowPatch([float(r.dot(O.i).subs(t,instant)),float(r.dot(O.j).subs(t,instant))],
[float(r.dot(O.i).subs(t,instant))+float(v.dot(O.i).subs(t,instant)),
float(r.dot(O.j).subs(t, instant))+float(v.dot(O.j).subs(t,instant))],
mutation_scale=20,
arrowstyle="->",color="r",label='${{v}}$')
vn = FancyArrowPatch([float(r.dot(O.i).subs(t, instant)),float(r.dot(O.j).subs(t,instant))],
[float(r.dot(O.i).subs(t, instant))+float(a.dot(O.i).subs(t, instant)),
float(r.dot(O.j).subs(t, instant))+float(a.dot(O.j).subs(t, instant))],
mutation_scale=20,
arrowstyle="->",color="g",label='${{a}}$')
vc = FancyArrowPatch([float(r.dot(O.i).subs(t, instant)),float(r.dot(O.j).subs(t,instant))],
[float(r.dot(O.i).subs(t, instant))+float(aCor.dot(O.i).subs(t, instant)),
float(r.dot(O.j).subs(t, instant))+float(aCor.dot(O.j).subs(t, instant))],
mutation_scale=20,
arrowstyle="->",color="b",label='${{a_{Cor}}}$')
ax.add_artist(vn)
ax.add_artist(vt)
ax.add_artist(vc)
plt.xlim((-10,10))
plt.ylim((-10,10))
plt.legend(handles=[vt,vn,vc],fontsize=20)
plt.grid()
plt.show()
# ## Further reading
#
# - Read pages 916-931 of the 18th chapter of the [Ruina and Rudra's book] (http://ruina.tam.cornell.edu/Book/index.html) about polar coordinates.
# ## Video lectures on the Internet
#
# - Khan Academy:
# - [Polar coordinates 1](https://www.youtube.com/watch?v=B5dOy4m6I1E)
# - [Polar coordinates 2](https://www.youtube.com/watch?v=1pQXJuWbz9w)
# - [Polar coordinates 3](https://www.youtube.com/watch?v=roSG9V3zApM)
# ## Problems
#
# 1. Find the polar basis (using the computer) for a projectile motion of a particle following the parametric equation below:
#
#
# \begin{equation}
# \vec{r}(t) = (10t-51){\bf{\hat{i}}} + \left(-\frac{9,81}{2}t^2+50t+100\right){\bf{\hat{j}}}
# \label{eq_13}
# \end{equation}
#
#
# 2. Problems from 15.1.1 to 15.1.14 from Ruina and Rudra's book,
# 3. Problems from 18.1.1 to 18.1.8 and 18.1.10 from Ruina and Rudra's book.
# ## Reference
#
# - Ruina A, Rudra P (2019) [Introduction to Statics and Dynamics](http://ruina.tam.cornell.edu/Book/index.html). Oxford University Press.