Brief Review of Undergraduate Classical Mechanics (unfinished)

By of Ruits - 9월 26, 2020



1. d'Alembert's principle $\to$ Lagrangian mechanics

Let's consider a dynamical system with N particles & k constraints. $\delta x^i \: \left( i=1,2,\cdots,3N \right) $ is a virtual displacement, and let $q^i \left( i=1,2,\cdots,3N-k \right) $ be a set of generalized coordinates. Virtual displacement can be represented by
$$\delta x^i = \frac{\partial x^i}{\partial q^j} \delta q^j$$
There is for sure no time-dependence since virtual displacement is just an instant phenomena. d'Alembert's principle indicates that virtual work done by constraint force $R_i$ in conjugation with $\delta x^i$ is always zero :
$$ R_i \delta x^i = 0$$
Alternatively, we can express the same principle using external force $F_i^{(a)}$ like below :
$$ (F_i^{(a)}-\dot{p}_i) \delta x^i =0$$

In case with k holonomic constraints, work done by external force is
$$\delta W = F_i^{(a)} \delta x^i = \left( F_i^{(a)} \frac{\partial x^i}{\partial q^j} \right) \delta q^j$$
Here, inside the parenthesis is called 'generalized force' $Q_j$. Thus $\delta W$ is,
$$\delta W = Q_j \delta q^j$$
Then now consider the work done by $\dot{p}_i$.
$$\begin{align} \dot{p}_i \delta x^i & = \sum_{i=1}^{3N} m_i \ddot{x}_i \delta x^i = \sum_{i=1}^{3N} m_i \ddot{x}_i \frac{\partial x^i}{\partial q^j} \delta q^j \\ & = \sum_{i=1}^{3N} m_i \left( \frac{d}{dt} \left( \dot{x_i} \frac{\partial x^i}{\partial q^j} \right) - \dot{x}_i \frac{d}{dt} \frac{\partial x^i}{\partial q^j} \right) \end{align}$$
Now the total derivative of $x_i(q^j)$ is
$$\dot{x}_i = \frac{\partial x_i}{\partial q^j} \dot{q}^j + \frac{\partial x_i}{\partial t}$$
$x^i$ and $q^i$ are in linear relationship with each other, so
$$\frac{\partial \dot{x}^i}{\partial \dot{q}^j}=\frac{\partial x^i}{\partial q^j}$$
Inserting this, we obtain
$$\sum_{i=1}^{3N} m_i  \frac{d}{dt} \left( \dot{x_i} \frac{\partial x^i}{\partial q^j} \right) = \sum_{i=1}^{3N} m_i  \frac{d}{dt} \left( \dot{x_i} \frac{\partial \dot{x}^i}{\partial \dot{q}^j} \right) = \frac{d}{dt} \frac{\partial}{\partial \dot{q}^j} \left( \frac{1}{2} \sum_{i=1}^{3N} m_i {\dot{x}_i}^2 \right) = \frac{d}{dt} \frac{\partial T}{\partial \dot{q}^j}$$ 
And for the next term,
$$\frac{d}{dt} \frac{\partial x^i}{\partial q^j}= \frac{\partial}{\partial q^j} \frac{dx^i}{dt}$$
I won't prove it here, but it is proved simply.
$$\sum_{i=1}^{3N} m_i \dot{x}^i \frac{d}{dt} \frac{\partial x^i}{\partial q^j} = \sum_{i=1}^{3N} m_i \dot{x}^i \frac{\partial}{\partial q^j} \frac{d x^i}{dt} = \frac{\partial T}{\partial q^j}$$
Using this, d'Alembert's principle can be simplified like below :
$$ (F_i^{(a)}-\dot{p}_i) \delta x^i =0 \:\:\:\:\: \Rightarrow \:\:\:\:\: \left( \frac{d}{dt} \frac{\partial T}{\partial \dot{q}^j} -  \frac{\partial T}{\partial q^j} - Q_j \right) \delta q^j = 0 $$
This expression we desired to obtain is called Lagrange equation, where $Q_j$ corresponds to the $q^j$ directional constraint force. If the external force is conservative, then $Q_j$ can be expressed by
$$Q_j \equiv F_i^{(a)} \frac{\partial x^i}{\partial q^j} = -\frac{\partial U}{\partial x^i} \frac{\partial x^i}{\partial q^j} = -\frac{\partial U}{\partial q^j}$$
Thus we get well known form of Lagrange equation ($L=T-U$).
$$\frac{d}{dt} \frac{\partial L}{\partial \dot{q}^j}-\frac{\partial L}{\partial q^j}=0$$





2. Action extremization $\to$ Lagrangian mechanics

Which path $x(t)$ minimize the action $S[x(t)]$? Use variational calculus to answer this question.
$$S[x(t)]=\int_{t_0}^{t_1} L(x^i(t), \dot{x}^i(t)) dt $$
Here $i=1,2,3,\cdots,3N$, The reason for not including higher order terms more than 1st is that physical state is completely determined by coordinates and speeds. Suppose there also are k constraints.
$$f_j(x^i)=0$$
Also $j=1,2,3, \cdots, k$. Simply varying $S[x(t)]$, we get
$$\delta S =0 \:\:\:\:\: \Rightarrow \:\:\:\:\: \int dt \left( \frac{d}{dt} \frac{\partial L}{\partial \dot{x}^i}-\frac{\partial L}{\partial x^i} \delta x^i \right)= 0$$
Since all $\delta x^i$s are not independent each other, we need to vary the condition of constraints.
$$\delta f_j =0 \:\:\:\:\: \Rightarrow \:\:\:\:\: \frac{\partial f_j}{\partial x^i} \delta x^j = 0$$
Each term for index i  is the functional of only $x^i(t)$, so we can introduce some Lagrange multiplier $\lambda(t)$ to solve this equation.
$$\frac{d}{dt} \frac{\partial L}{\partial \dot{x}^i}-\frac{\partial L}{\partial x^i} = \lambda_j(t) \frac{\partial f_j}{\partial \dot{x}^i}$$
We want to know k unknown functions $\lambda_j(t)$ and 3N unknown functions $x^i(t)$. Also we have k unknown equations $f_j(x^1, x^2, \cdots, x^{3N})=0$ and 3N Lagrange equations. So we can find all unknown functions with appropriate initial conditions. In addition, the right side of the above equations means the $x^i$ directional constraint force.





3. Noether's Theorem

Sometimes it is the case that $\delta S$ vanishes for certain limited variations of the path without imposing any condition at all. When this happens, we say that $S$ has a symmetry. A symmetry of an action functional $S[x]$ is an infinitesimal transformation of the path, $x^i(t) \rightarrow x^i+\epsilon^i(x^j)$ that leaves the action invariant :
$$S[x^i+\epsilon^i(x^j(t),t)]=S[x^i(t)]$$
or, up to constant :
$$S[x^i+\epsilon^i(x^j(t),t)]=S[x^i(t)]+C$$
which means action is transformed up to total derivative of some function. $$L(x^j)=L(x^j+\epsilon^j)+\frac{d}{dt} G(x^j,t) $$
in any case, physics remains unchanged.

(1) For some symmetry transformation, regardless of the choice of path, suppose 
$$\delta L \equiv \frac{dG(x^j,t)}{dt}$$

(2) Also suppose path is transformed like below :
$$x^i \rightarrow x^i+\epsilon^i(x^j(t))$$

(3) Then Noether current (more natural concept; charge) is conserved.
$$\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{x}^i} \epsilon^i(x^j) - G(x^j,t) \right)_{Classical\:path} = 0$$

[PROOF] The action has a symmetry. so identically $\delta S \equiv G \mid ^{t_2}_{t_1}$. No equation of motion has been used. Since path transformed $x^i \rightarrow x^i+\epsilon^i(x^j(t))$, $\delta S$ can be represented also like
$$\begin{align} \delta S & \equiv \int_{t_1}^{t_2} dt \left( \frac{\partial L}{\partial x^i} \epsilon^i + \frac{\partial L}{\partial \dot{x}^i} \dot{\epsilon}^i \right) \\ & = \left. \frac{\partial L}{\partial \dot{x}^i} \epsilon^i \right|^{t_2}_{t_1} + \int_{t_1}^{t_2} dt \left( \frac{\partial L}{\partial x^i} - \frac{d}{dt} \frac{\partial L}{\partial \dot{x}^i} \right) \end{align}$$
This equation satisfies for all paths. Regardless of path selected, we here only imposed symmetry condition. Let's now select our path clasically. For the classical path, Lagrange equation holds. Thus 2nd term vanishes. Then we get
$$ \left. \left( \frac{\partial L}{\partial \dot{x}^i} \epsilon^i - G \right) \right|_{t_1}^{t_2} = 0$$
for any two times $t_1, t_2$.





4. How to construct Lagrangian?

Let's next think about why our Lagrangian is the form of $T-U$ in Classical Mechanics.

(1) Consider a freely moving particle in an inertial frame. There is no preferred point in the space (nor the time), thus there is no dependence on $\vec{x}$ or $t$. In this case we can obtain $L=L(\dot{x}^2)$. It follows 
$$\frac{\partial L}{\partial x^i}=0 \:\:\:\:\: \Rightarrow \:\:\:\:\: \frac{d}{dt} \frac{\partial L}{\partial \dot{x}^i}=0$$
Of course, this implies $\vec{v}$ is a constant.

(2) Let's do a Galilei's transformation. Take two inertial frame $K$ and $K'$, moving one respect to other with a small velocity $\vec{w}$. Then our physics must be same, thus Lagrangian must be up to total time derivative of some function.
$$L'({v'}^2) = L(v^2+2\vec{w} \cdot \vec{v}+w^2) = L(v^2)+2\frac{\partial L}{\partial {v^2}} \vec{w} \cdot \vec{v}$$
The last term is a total derivative of something.
$$\frac{\partial L}{\partial {v^2}} \vec{w} \cdot \vec{v} = \frac{d}{dt} g(\vec{x}, t)$$
From the above, we can notice that $\partial L / \partial {v^2}$ is just a constant. So we can conclude $L \propto v^2$.

(3) Now consider the case with a potential. The exact differential of $L$ is
$$dL = \frac{\partial L}{\partial q} dq + \frac{\partial L}{\partial \dot{q}} d\dot{q}$$
Suppose $\frac{\partial^2 L}{\partial q \partial \dot{q}} = 0$. (i.e. $\frac{\partial L}{\partial \dot{q}}=\frac{\partial L}{\partial \dot{q}} (\dot{q})$) Then the first term should be a function of only $q$. ($\frac{\partial L}{\partial q} = \frac{\partial L}{\partial q} (q)$) From (2), we know $\frac{\partial L}{\partial \dot{q}}=2a \dot{q}$. Take $a=m/2$, the half mass of the particle. Then we know $\frac{\partial L}{\partial q} = F(q)$.
$$dL=p(\dot{q}) d\dot{q} + F(q) dq$$
From the fact that $F=dU/dq$, we can construct $L(q,\dot{q})=T(\dot{q})-U(q)$.

(4) Then what if $\frac{\partial^2 L}{\partial q \partial \dot{q}} \neq 0$? Lets consider an example of particle in a scalar field ($\phi(\vec{x},t)$) and a vector field ($\vec{A} (\vec{x},t)$). You can also think this situation as an analogue of electromagentics.

[ Example : Constructing EM Lagrangian ]
We already know in this case Lorentz force is
$$m\ddot{\vec{x}} = q\vec{E} + \frac{q}{c} \vec{v} \times \vec{B}$$
Expanding $i$ th component,
$$\begin{align} m\ddot{x}^i &= qE_i + \frac{q}{c} \epsilon_{ijk} \dot{x}^j B^k \\ & = -q\left( \partial_i \phi -\frac{1}{c}\partial_t A_i \right) + \frac{q}{c} \epsilon_{ijk} \dot{x}^j \epsilon^{klm} \partial_l A_m \\ &=-q \partial_i \phi - \frac{q}{c} \frac{d}{dt} A_i + \frac{q}{c} \dot{x}^j \partial_j A_i +\frac{q}{c} \left( \dot{x}^j \partial_i A_j - \dot{x}^j \partial_j A_i \right) \end{align}$$
Then rearranging this, we obtain
$$\frac{d}{dt} \left( m\dot{x}_i+\frac{q}{c} A_i \right) = -q \partial_i \phi + \frac{q}{c} \dot{x}^j \partial_i A_j$$
Suppose RHS is $\partial L / \partial x^i$ and LHS is $d/dt \left( \partial L / \partial \dot{x}^i \right)$. Lets check exactness. (if not exact, we should multiply some factor)
$$\frac{\partial^2 L}{\partial \dot{x}^j \partial x^i}=\frac{q}{c} \partial_i A_j = \frac{\partial^2 L}{ \partial x^i \partial \dot{x}^j}$$
Thus in this non-conservative potential, we can find Lagrangian from pre-known Lerentz force law. (In fact, vector potential is just a differential one-form, magnetic field is two-form, and electric field is one-form. Think about field strength tensor $F_{\mu \nu}$. Not vectors strictly.)
$$L(\vec{x},\dot{\vec{x}})=\frac{1}{2} m \dot{x}^2 + \frac{q}{c} \dot{\vec{x}} \cdot \vec{A} - q \phi$$





5. Central potential

In this section, we now discuss about a potential which depends only on radial distance $r$.

(1) Reducing degrees of freedom : In principle, motion in 3D spherical coordinates, not any constraint is determined. So if we should use 3d spherical coordinates, one simple solution of Euler-Lagrange equation gives us one constraint that azimuthal angle is always constant ($\theta = \pi/2$). Also equivalently, use angular momentum conservation to reduce one degree of freedom in motion. You can easily prove $d\vec{L}/dt=0$. This means 3 degrees of freedom decreased to 2 (planar motion). So we now we can choose 2 generalized coordinate $r$ and $\theta$. Then our Lagrangian is the form of
$$L=\frac{1}{2}m \left( \dot{r}^2+r^2\dot{\theta}^2 \right) - U(r)$$
The angle $\theta$ is, at a glance, cyclic coordinate whose corresponding canonical momentum is the angular momentum of the system
$$p_{\theta}=\frac{\partial L}{\partial \theta}=mr^2\dot{\theta}$$
It was our first constant of motion we just found. let $p_{\theta}=l$. Also this means constant areal speed ($dA=\frac{1}{2}r^2 d\theta$). We need to reduce the other one degree of freedom to determine the trajectory. Let's find another constant of motion. Using Noether's theorem about time translation symmetry ($t\to t+\epsilon$). 

[ Example : Find Noether charge for the time translation symmetry ]
Regardless of any path we select, there exist $G(x^i,t)$ such that
$$\delta L = \frac{dG}{dt}$$
We can easily find $G=\epsilon L$ since $\delta L=L(t+\epsilon)-L(t)=\epsilon \frac{dL}{dt}$.
Path is transformed like below :
$$x^i(t) \to x^i(t+\epsilon) = x^i(t)+\dot{x}^i(t) \epsilon$$
Thus conserved current is
$$h(x,\dot{x})=\frac{\partial L}{\partial x^i} \dot{x}^i -L$$
$h$ is called 'Energy function'. Be aware that this is different from Hamiltionian.

We found another conserved quantity called 'Energy'. Now we can reduce one more degree of freedom.

(2) Finding the trajectory equation : We have 2 constants of motion.
$$E=\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2)+U(r) \:\:\:\:\:\: l=mr^2\dot{\theta}$$
Solving this differential equation is not very hard. The simple equation is
$$\frac{dr}{dt}=\sqrt{ \frac{2}{m} \left( E-U(r)\right) -\frac{l^2}{m^2r^2}}$$
We want to know $\theta$ vs $r$ relation. Using the fact that $dt=\frac{mr^2}{l}d\theta$,
$$\theta -\theta_0 = \int_{r_0}^{r} \frac{l dr / r^2}{\sqrt{ 2m\left( E-U(r)\right) -l^2/r^2 }}$$
If the variable of integration is changed to $u=1/r$,
$$\theta-\theta_0= \int_{u_0}^{u} \frac{du}{\sqrt{ \frac{2mE}{l^2} - \frac{2mU(r)}{l^2}-u^2 }}$$
In case of $U(r)=ar^{n+1}=au^{-(n+1)}$, this is integrable in terms of simple functions only in certain cases. The result can be expressed in terms of trigonometric functions when $n=1,-2,-3$ and can be expressed in terms of elliptic functions when $n=5,3,0,-4,-5,-7$. The other exponents can be expressed in terms of hypergeometric function (Kummer function). Further info about Kummer function : 
(3) Conditions for closed orbit (Bertrand's theorem) : The only central forces that result in closed orbits for all bound particles are the inverse-square force  ($U \sim 1/r$) and Hooke's force ($U \sim r^2$). Simple proof can be found in Goldstein's Classical Mechanics.

(4) Kepler problem (n=-2) : First consider an 2-body problem. Lagrangian is
$$L=\frac{1}{2}m_1 \dot{\vec{r}}_1^2 + \frac{1}{2}m_2 \dot{\vec{r}}_2^2 - U(\left| \vec{r}_1 - \vec{r}_2 \right|)$$
Changing the variable :
$$(\vec{r}_1 ,\vec{r}_2) \rightarrow \left( \vec{r}=\vec{r}_1-\vec{r}_2, \vec{R}=\frac{m_1 \vec{r}_1 + m_2 \vec{r}_2}{m_1+m_2} \right)$$
and then our Lagrangian is reduced to freely moving center of mass + relative central potential problem ($\mu$ : reduced mass) :
$$L=\frac{1}{2}(m_1+m_2)\dot{\vec{R}}^2 + \frac{1}{2} \mu \dot{\vec{r}}^2 -U(r)$$
Then now only thing deserve to deal with is relative motion. Consider $U(r)=-\alpha/r$. The indefinite integral is the form
$$\int \frac{dx}{\sqrt{a+bx+cx^2}}=\frac{1}{\sqrt{-c}} \cos^{-1} \left( -\frac{b+2cx}{\sqrt{b^2-4ac}}\right)$$
WLOG, set $\theta_0=0$ when $r=r_0$. Then the solution we can obtain is
$$r=\frac{p}{1+e\cos \theta}$$
where $p=\frac{l^2}{m\alpha}$ is called 'semi-latus rectum'. Geometric meaning is like below :

Fig 1. Semi-latus rectum

and $e=\sqrt{1+2El^2/m\alpha^2}$ is an eccentricity. I'll skip here the concept of geometry of conic section. Some key relations are written below. In orbital dynamics, product of $G$ and $M$ is called $\mu=GM=\alpha/m$. And specific angular momentum per mass is called $h=l/m$.

[ Some useful relations in Keplerian orbit ]
  a. $p=a(1-e^2)$ : def of semi-latus rectum
  b. $T=2\pi \sqrt{a^3/\mu}$ : period
  c. $b=a\sqrt{1-e^2}$ : semi-minor axis vs semi-major axis
  d. $r_p=a(1-e)$, $r_a=a(1+e)$ : distance to periapsis / apoapsis
  e. $h^2=\mu p$ : angular momentum per mass vs semi-latus rectum
  f. $E=-\mu /2a$ : energy vs semi-major axis
Using e & f relations, you can transform constants of motion set $(l, E) \leftrightarrow (a, e)$ whichever you want. $a$ & $e$ is purely geometrical and $l$ & $E$ is purely dynamical constants.

Briefly speaking, $r_p$ and $r_a$ can be obtain by effective potential graph.

Fig 2. Effective potential

Period can be obtained by areal speed.
$$\frac{dS}{dt}=\frac{1}{2}r^2 \dot{\theta} = \frac{l}{2m} \:\:\:\:\: \rightarrow \:\:\:\:\: \int_{0}^{\pi a b} dS=\frac{l}{2m} P$$
Another minor concept is LRL vector. $\vec{A}=\vec{p} \times \vec{L}$, which is also a constant. It always directs towards periapsis. This conserved quantity became an important key for Pauli to borrow the idea of orbital mechanics to solve Schrodinger equation of hydrogen atom in 1926. If you're interested in it, read below :

I'll introduce some more concepts in orbital dynamics. First one is eccentric anomaly $E$. $E$ is mainly used to parameterize the equation of ellipse $x^2/a^2 + y^2/b^2=1$. $x=a\cos E$ and $y=b\sin E$.
Fig 3. Eccentric anomaly

Eccentric anomaly $E$ and true anomaly $\theta$ have a following simple relationship.
$$\sqrt{1-e}\tan \frac{\theta}{2}=\sqrt{1+e} \tan \frac{E}{2}$$
and trajectory equation would be
$$r=\frac{a(1-e^2)}{1+e\cos \theta}\:\:\:\:\:\Leftrightarrow\:\:\:\:\:r=a(1-e\cos E)$$
Finally, after some algebra, $r$ vs $t$ relation can be re-written in $E$ vs $t$ relation
$$t=\sqrt{\frac{m}{2}} \int_{r_p}^{r} \frac{dr}{\sqrt{\frac{\alpha}{r}-\frac{l^2}{2mr^2}+E}} \:\:\:\:\:\Leftrightarrow \:\:\:\:\: t=\sqrt{\frac{ma^3}{\alpha}} \int_0^E (1-e\cos E) dE$$
(the left E is Energy and right one is eccentric anomaly. Sorry for confusing notation.)
If we integrate 0 to $2\pi$, can obtain period again. We want to find the angle which is linearly increasing about time. Let's define $M=E-e\sin E$. This is called 'mean anomaly', which means literally increasing uniformly.
$$M=E-e\sin E \:\:\:\:\rightarrow \:\:\:\: \frac{dM}{dt}=\sqrt{\frac{\alpha}{ma^3}} $$
Then using this angle, we can find true anomaly $\theta$ when time is $t$. This problem is called 'Kepler problem'.

(5) $1/r^2$ potential (n=-3) : Suppose radial force is given $F(r)=-k/r^2 - q/r^3\:\:(k,q>0)$. From Euler-Lagrange equations, we can obtain $r-\theta$ 2nd order ODE.
$$\frac{d^2}{d\theta^2} \left( \frac{1}{r} \right) +\frac{1}{r}=-\frac{mr^2}{l^2}F(r)$$
changing variable $u=1/r$ reduces to
$$\frac{d^2u}{d\theta^2}+\left( 1-\frac{mq}{l^2} \right)u=\frac{mk}{l^2}$$
If $q<l^2/m$, the solution is just a strange conic section :
$$\frac{1}{r}=\frac{mk}{l^2-mq} \left[ 1+a \cos\left( \sqrt{1-\frac{mq}{l^2}}(\theta-\theta_0) \right) \right]$$
If $q=l^2/m$, the solution is infinitely converging spiral :
$$\frac{1}{r}=\frac{mk}{l^2} \theta^2 + A \theta + B$$
If $q>l^2/m$, the solution is hyperbolic function :
$$\frac{1}{r}=\frac{mk}{l^2-mq} \left[ 1+A\exp\left( \sqrt{\frac{mq}{l^2}-1}(\theta-\theta_0) \right) + B \exp\left( -\sqrt{\frac{mq}{l^2}-1}(\theta-\theta_0) \right) \right]$$

(6) Relativistic Effect : Read the post I wrote before. Directly corresponds to the test of general relativity.





6. Non-inertial frame

If our system is observed in non-inertial reference frame, Newtonian mechanics doesn't seem to be valid anymore. This section is a bridge to connect lots of engineering physics (satellite attitude control, rigid body dynamics, etc.). Let's briefly review. First, we adopt our basis non-inertial. So time derivative of the basis is now not a constant. Linearly accelerating frame is too trivial to deal with, so we only consider rotating frame here.
$$\frac{d\vec{r}}{dt}=\frac{d}{dt}(x^i \vec{e}_i) = \frac{dx^i}{dt} \vec{e}_i + x^i \frac{d\vec{e}_i}{dt}$$
The last term is time derivative of basis vector. How can we describe? Suppose an rotating angular velocity to be $\vec{w}$. Thinking about it simply drawing a picture. The time derivative of our basis should be
$$\frac{d\vec{e}_i}{dt} = \vec{w} \times \vec{e}_i$$
I will write derivative of only components with prime(') notation.
$$\frac{d'\vec{r}}{dt'} \equiv \frac{dx^i}{dt} \vec{e}_i$$
Then full description can be written as
$$\frac{d\vec{r}}{dt}=\frac{d'\vec{r}}{dt'}+\vec{w} \times \vec{r}$$
Displacement is a (mathematical) vector, thus vector itself doesn't change at all. Then the second derivative is just a little more messy. (here I assumed $vec{w}$ is constant.)
$$\begin{align} \frac{d}{dt}\frac{d\vec{r}}{dt} & =\frac{d'}{dt'} \left( \frac{d'\vec{r}}{dt'}+\vec{w} \times \vec{r} \right) + \vec{w} \times \left( \frac{d'\vec{r}}{dt'}+\vec{w} \times \vec{r} \right) \\ & = \frac{{d'}^2 \vec{r}}{d{t'}^2} + 2\vec{w} \times \frac{d' \vec{r}}{dt'} + \vec{w} \times \vec{w} \times \vec{r} \end{align}$$
To see an artificial force-like effect, rearranging all these terms with non-inertial acceleration to locate LHS.
$$\vec{a}' = \frac{\vec{F}}{m} - \underbrace{ 2\vec{w} \times \frac{d' \vec{r}}{dt'} }_{\text{Coriolis term}} - \underbrace{ \vec{w} \times \vec{w} \times \vec{r} }_{\text{centrifugal term}}$$
Here $\vec{a}'$ is the (fake ; not a fake in fact) acceleration in our rotating frame. Here's an easy example of free-fall in the Earth.

[ Example : Free-fall problem revisited ]
Consider a situation as shown below. Suppose $|\vec{v}| \gg wR$ and $r'/R \ll 1$. If we drop a massive particle at the height of $h$ when time $t=0$, where will this particle be when it hits the ground ($z'=0$)?

Fig 4. Free fall example

Our expression of acceleration is :
$$\vec{a}' = \frac{\vec{F}}{m} -  2\vec{w} \times \frac{d' }{dt'} \left( \vec{R} + \vec{r}' \right)  -  \vec{w} \times \vec{w} \times \left( \vec{R}+\vec{r}' \right)$$
where $\vec{F}=-mg \hat{z}'$, $\vec{w} \times \vec{w} \times (\vec{R}+\vec{r}') \simeq -w^2 R \cos\lambda \left( -\sin \lambda \hat{x}' - \cos \lambda \hat{z}' \right)$. Here is the other assumption $v \gg wR$, so we neglect centrifugal term since
$$\left| \vec{w} \times \vec{w} \times \vec{R} \right| \ll \left| \vec{w} \times \vec{v}' \right| $$
So our problem is to solve 
$$\vec{a}'=-g\hat{z}'-2\vec{w} \times \vec{v}'$$
Compared to the first gravitational term, Coriolis term is sufficiently small. So we can treat it as a perturbation. Our zero-th order ($w^0$ order) solution is
$$\vec{r}'_{(0)}=\left( 0,0,h-\frac{1}{2}gt^2 \right) $$
Then the next order ($w^1$ order) can be obtained to solve the equation shown below
$$\vec{a}'_{(1)}=-2 \vec{w} \times \vec{v}'_{(0)} \:\:\:\:\: \Rightarrow \:\:\:\:\: \vec{r}'_{(1)}=\frac{1}{3}wgt^3 \cos \lambda \hat{y}'$$
Also higher order solution can be obtained by solving
$$\vec{a}'_{(n+1)}=-2 \vec{w} \times \vec{v}'_{(n)}$$
Thus our full solution is
$$\vec{r}'=\vec{r}'_{(0)} (w^0) + \vec{r}'_{(1)} (w^1) + \cdots = \left( h-\frac{1}{2}gt^2 \right) \hat{z}' + \left( \frac{1}{3} wgt^3 \cos \lambda \right) \hat{y}' + \mathcal{O}(w^2)$$





7. Rigid body

In the rigid body system, all particles consistituting rigid body satisfy $\vec{v}=\vec{w} \times \vec{r}$. It is an example of holomonic constraints. In arbitrary dimension, angular momentum is defined as the form of 2-vector :
$$\vec{\vec{L}}=\vec{x} \wedge \vec{p} \:\:\:\:\: \Leftrightarrow \:\:\:\:\: L^{ij}=x^i p^j$$
And in 3 spatial dimension, volume form(Levi-civita) contracted with 2-vector yields an 1-form to be the traditional vector notation. (volume duality)
$$\underline{L} = \underline{\underline{\underline{\epsilon}}} \cdot \vec{\vec{L}} \:\:\:\:\: \Leftrightarrow \:\:\:\:\: L_i = \epsilon_{ijk} L^{jk}$$
Then we can calculate angular momentum by hand.
$$\begin{align} L_i & = \epsilon_{ijk}L^{jk} = \epsilon_{ijk} x^j p^k = m \epsilon_{ijk} x^j \sqrt{g} \varepsilon^{klm} w_l x_m \\ &= m (\delta_i^l \delta_j^m - \delta_i^m \delta_j^l) x^j x_m w_l =m (w_i x^j x_j - x_i x^j w_j)  \\ &= m(x^2 \delta_{ij} -x_i x_j) w^j \equiv I_{ij} w^j \end{align}$$

(1) Thus we can define the inertia tensor by :
$$\begin{align} I_{ij} & = \sum_n m_{(n)} \left( x_{(n)}^2 \delta_{(n)ij}- x_{(n)i} x_{(n)j} \right) & \:\:\:\:\: \text{(discrete system)} \\ & =\int_V \rho (\vec{x}) \left( x^2 \delta_{ij} - x_i x_j \right) dV & \:\:\:\:\: \text{(continuous system)} \end{align}$$

In general, it is very hard to calculate 9 components of inertial tensor. So we need to change our frame to make it easier. Mathematically it always can be diagonalized orthogonally. Thus we can find orthogonal set of basis which makes our inertia tensor much more easier.

(2) Solving the eigenvalue equation
$$\mathbf{I} \cdot \vec{w} = I \vec{w}$$
we can take $\vec{e}_i=\vec{w}_i/|\vec{w_i}|$ as an orthonormal basis. then
$$\vec{w}=(w^1, w^2, w^3) \:\:\:\:\:\:\:\: \mathbf{I}=\begin{pmatrix} I_1 & 0 & 0 \\ 0& I_2 & 0 \\ 0 & 0 & I_3 \end{pmatrix} \:\:\:\:\:\:\:\: \vec{L}=(I_1w^1, I_2 w^2, I_3 w^3)$$

(3) There is another key tool to make problems easier is to take the center of mass to be the origin of our body frame. Then total kinetic wrt inertial frame would be separated into 2 parts.
$$T_{inertial} = \frac{1}{2}MV^2 + \frac{1}{2} \sum_{(n)} m_{(n)} (\vec{w} \times \vec{r}_{(n)})^2 = T_{trans}+T_{rot}$$

(4) Parallel-axis theorem states
$$ I_{ij}'=I_{ij ; \text{CM}}+M(d^2 \delta_{ij} -d_i d_j)$$

If center of body frame is not CM, we can transform inertia tensor using this theorem. Let's now show the inertia tensor is really a tensor.

[ Notes : Inertia tensor is really a tensor ]
First consider an angular momentum transformation rule. Note that angular momentum is not a vector, here we contracted with volume form, so it is 1-form.
$$L'_i=\epsilon'_{ijk} L'^{jk}=\sqrt{g'} \varepsilon_{ijk}L'^{jk}=\frac{\partial x^p}{\partial x'^i} \epsilon_{pqr} L^{qr} =\frac{\partial x^p}{\partial x'^i} L_p$$
And angular velocity is a vector.
$$w'^j=\frac{\partial x'^j}{\partial x^q} w^q$$
Since $L_i=I_{ij} w^j$, the inertia tensor transformed like rank-2 covariant tensor.
$$I'_{ij}=\frac{\partial x^p}{\partial x'^i} \frac{\partial x^q}{\partial x'^j} I_{pq}$$

(5) Newtonian rigid body mechanics : Now we're ready to find the equation of motion of rigid body. First, our reference frame is fixed body frame. So we need to consider time-derivative of the basis vector. In our non-inertial frame,
$$\vec{N}=\frac{d\vec{L}}{dt}=\frac{d' \vec{L}}{dt'}+\vec{w} \times \vec{L}$$
Suppose inertia tensor doesn't change with time.
$$N_i = I_{ij}\dot{w}^j-\epsilon_{ijk} w^j I^{k}_{\:\:l} w^l$$
If we take our basis to be diagonalize the inertia tensor, we get
$$I_x \dot{w}_x -w_y w_z (I_y-I_z)=N_x \:\:\:\:\:(x,y,z \:\text{cyclic})$$


(6) Euler angles : Since $\vec{w}$ is not coordinates, but a time derivative of some angular coordinates. Our goal is to extend to Lagrangian mechanics, so we need to define some proper angular coordinates wrt inertial coordinates : people call this 'Euler angles'. There are also many other angular systems that indicates the principal axis of the rigid body, but 'Euler angles' are mainly used many applied fields. See the figure below.

Fig 5. Euler angles

We can find each time derivative of $\theta$, $\phi$, $\psi$ is directional angular velocity.
$$\begin{align} \dot{\vec{\phi}} & = ( \dot{\phi} \sin \theta \sin \psi ) \hat{x}' - ( \dot{\phi} \sin \theta \cos \psi ) \hat{y}' + ( \dot{\phi} \cos \theta ) \hat{z}' \\ \dot{\vec{\theta}} & = ( \dot{\theta} \cos \psi ) \hat{x}' + ( \dot{\theta} \sin \psi ) \hat{y}' \\ \dot{\vec{\psi}} &= \dot{\psi} \hat{z}' \end{align}$$
And $\vec{w}=\dot{\vec{\phi}}+\dot{\vec{\theta}} + \dot{\vec{\psi}}$. Sorting by components, we obtain
$$\begin{align}w_1 & = \dot{\phi} \sin \theta \sin \psi + \dot{\theta} \cos \psi \\ w_2 & = \dot{\phi} \sin \theta \cos \psi - \dot{\theta} \sin \psi \\ w_3 &= \dot{\phi} \cos \theta + \dot{\psi} \end{align}$$

(7) Lagrangian prescription : We just specified 3 independent angular coordinates. We saw above kinetic term can be rewritten $T=T_{trans,CM}+T_{rot}$ when our origin of frame is CM. If not, in general, full contribution of kinetic is $T=T_{rot}$ wrt our non-moving arbitrary origin. Rotational kinetic is also can be written with Euler angles.
$$T_{rot}=\frac{1}{2}(I_1 w_1^2 + I_2 w_2^2 + I_3 w_3^2) = \frac{1}{2} I_1 \left( \dot{\phi}^2 \sin^2 \theta + \dot{\theta}^2 \right) + \frac{1}{2} I_3 \left( \dot{\phi} \cos \theta + \dot{\psi} \right)^2$$
and if our total Lagrangian is the form of $U(\phi, \theta, \psi)$, $\phi$ and $\psi$ are cyclic. Thus we get 2 independent constants of motion.
$$\frac{\partial L}{\partial \dot{\phi}}=p_{\phi}, \;\;\;\;\; \frac{\partial L}{\partial \dot{\psi}}=p_{\psi}$$
And as we saw before, it also has time translation symmetry. So we get final constant of motion 'Energy'. In principle, 3 constants of motion fully describe the dynamics of our system. One of the 3 constants are replaced instead of $\theta$-directional equation of motion. Lets look at an example of motion of a top.

[ Example : Motion of a symmetric top ]
Consider a rigid top whose moment of inertia satisfies $I_1=I_2\neq I_3$. See the figure 6.


Fig 6. Symmetric Top & Effective potential

Then our Lagriangian is
$$L(\phi,\theta,\psi) = \frac{1}{2} I_1 \left( \dot{\phi}^2 \sin^2 \theta + \dot{\theta}^2 \right) + \frac{1}{2} I_3 \left( \dot{\phi} \cos \theta + \dot{\psi} \right)^2 - Mgl\cos \theta $$
and 3 constants of motion we get are
$$ \begin{align} \text{(1)} &\;\;\; \frac{\partial L}{\partial \dot{\phi}}=\left( I_1 \sin^2 \theta + I_3 \cos^2 \theta \right) \dot{\phi} + I_3 \dot{\psi} \cos \theta = p_{\phi} \\ \text{(2)} &\;\;\; \frac{\partial L}{\partial \dot{\psi}} = I_3 \left( \dot{\psi} + \dot{\phi} \cos \theta \right)  =p_{\psi} \\ \text{(3)} &\;\;\;  \frac{1}{2} I_1 \left( \dot{\phi}^2 \sin^2 \theta + \dot{\theta}^2 \right) + \frac{1}{2} I_3 \left( \dot{\phi} \cos \theta + \dot{\psi} \right)^2 + Mgl\cos \theta = E \end{align}$$
The last one is replacement of $\theta$-directional Euler-Lagrange equation. We can rewrite (1) & (2) and inserting into (3) to leave only $\theta$ dependence.
$$ \begin{align} \text{(1)} &\;\;\; \dot{\phi}=\frac{p_{\phi}-p_{\psi} \cos \theta}{I_1 \sin^2 \theta} \\ \text{(2)} & \;\;\; \dot{\psi} = \frac{p_{\psi}}{I_3}-\frac{(p_{\phi}-p_{\psi} \cos \theta) \cos \theta}{I_1 \sin^2 \theta} \\ \text{(3)} &\;\;\; E=\frac{1}{2}I_1 \dot{\theta}^2 + \frac{(p_{\theta}-p_{\phi} \cos \theta)^2}{2 I_1 \sin^2 \theta} + \frac{p_{\psi}^2}{2I_3} + Mgl \cos \theta \end{align}$$
Then we can deform (3) to solve the exact solution
$$t(\theta)=\int \frac{d\theta}{\sqrt{2(E-U_{eff}(\theta))/I_1}}, \;\;\;\; U_{eff}(\theta)=\frac{(p_{\theta}-p_{\phi} \cos \theta)^2}{2 I_1 \sin^2 \theta} + \frac{p_{\psi}^2}{2I_3} + Mgl \cos \theta $$
Rather than finding a solution, let's guess the motion through the effective potential.
If Energy minimized ($\theta=\theta_m$), $\theta_m$ satisfies
$$\left.\frac{dU_{eff}}{d\theta} \right|_{\theta=\theta_m}=\frac{-B^2 \cos \theta_m + p_{\psi}B\sin^2 \theta_m}{I_1 \sin^3 \theta_m} - Mgl \sin \theta_m = 0, \:\:\: B=p_{\phi}-p_{psi}\cos \theta_m$$
Solving this equation, 
$$B=\frac{p_{\psi}\sin^2 \theta_m}{2\cos\theta_m} \left( 1 \pm \sqrt{1-\frac{4MglI_1 \cos\theta_m}{p_{\psi}^2}} \right)$$
This should be real. Thus $\psi$-directional angular momentum needed to rotate is at least
$$p_{\psi}^2 \geq 4Mgl I_1 \cos \theta_m$$
If Energy is $E=E'$, $\theta$ has two turning points ($\theta_1, \theta_2$). In this case $\hat{z}'$ rotate about $\hat{z}$ axis with precession. If $\dot{\phi}>0$ for all time, it has precession like 7(a). If $\dot{\phi}=0$ at some moment, it has precession like 7(b). If $\dot{\phi}<0$ for some time interval, it contains retrograde motion 7(c).
Fig 7. (a) anterograde precession / (b) peak / (c) retrograde precession





8. Hamiltonian mechanics

If you want to use dynamical degrees of freedom $(p,q)$ instead of $(q,\dot{q})$, we now should conduct a Legendre transformation.
$$\begin{align} dL(q^i, \dot{q}^i) &=\frac{\partial L}{\partial q^i}dq^i + \frac{\partial L}{\partial \dot{q}^i} d\dot{q}^i + \frac{\partial L}{\partial t} dt \\ &= \dot{p}_i dq^i + p_i d\dot{q}^i + \frac{\partial L}{\partial t} dt \\ &=d(p_i \dot{q}^i)-\dot{q}^i dp_i + \dot{p}_i dq^i + \frac{\partial L}{\partial t} dt  \\ d(p_i \dot{q}^i -L) & =\dot{q}^i dp_i - \dot{p}_i dq^i - \frac{\partial L}{\partial t} dt\end{align}$$
We can define another quantity $H=pq-L$ parameterized with $p$ and $q$, which is called 'Hamiltonian'. It is different from Energy function $h(q,\dot{q}$. Fundamentally its parameters cover $q$ and $\dot{q}$. Think about thermodynamic potentials.
$$dU(S,V,N)=TdS-PdV+\mu dN$$
$$dF(T,V,N)=-SdT-PdV+\mu dN$$
It is also an example of Legendre transformation. Anyway, If we define Hamiltonian like this, we can get Hamilton's equation of motion.
$$\dot{q}^i = \frac{\partial H}{\partial p_i},\;\;\;\;\;\; \dot{p}_i=-\frac{\partial H}{\partial q^i}$$
This is the 2 sets of 1st order differential equation while Lagrangian mechanics yields a set of 2nd order differential equation. Let's consider total time-derivative of the Hamiltonian.
$$\frac{dH}{dt}=\frac{\partial H}{\partial p_i} dp_i + \frac{\partial H}{\partial q^i} dq^i + \frac{\partial H}{\partial t}=\dot{q}^i \dot{p}_i - \dot{p}_i \dot{q}^i + \frac{\partial H}{\partial t}=\frac{\partial H}{\partial t}$$
This is a key concept of Hamiltonian mechanics. If Hamiltonian doesn't depend on time explicitly, where we call it has a time-translational symmetry, $H$ is conserved quantity. We saw above that Noether current(charge) about time-translation symmetry was Energy function.

We can also construct Hamiltonian mechanics from action extremization.
$$S=\int L dt =\int dt \left(p_i\dot{q}^i-H \right)$$
Varying the action, where our boundary variation is always 0, we get
$$\delta S = \int dt \left( \dot{q}^i \delta p_i + p_i \delta \dot{q}^i - \frac{\partial H}{\partial p_i} \delta p_i - \frac{\partial H}{\partial q^i} \delta q^i \right) = 0$$
and after integration by parts,
$$\int dt \left( \left( \dot{q}^i - \frac{\partial H}{\partial p_i} \right) \delta p_i - \left( \dot{p}_i + \frac{\partial H}{\partial q_i} \right) \delta q^i \right) =0$$
from which we can derive Hamilton's equations of motion also.





9. Liouville's theorem

Let's consider a vector like below :
$$\vec{X} = (q^1, q^2, \cdots, q^N, p_1, p_2, \cdots, p_N) \equiv (q^i, p_i)$$
The future of this vector is fully determined with 2N initial conditions. It means that in the hyperspace of $(q^i, p_i)$ its streamline is fully determined. The time derivative is
$$\vec{V}=(\dot{q}^i,\dot{p}_i)=\left(\frac{\partial H}{\partial p_i},-\frac{\partial H}{\partial q^i} \right)$$
This forms a velocity field. Then the divergence of this velocity field is
$$\nabla \cdot \vec{V}= \frac{\partial^2 H}{\partial q^i \partial p_i} - \frac{\partial^2 H}{\partial p_i \partial q^i}=0$$
This means dynamical stream on hyperspace is a kind of incompressible fluid. We call this hyperspace $(q^i, p_i)$ 'phase space'. A phase space distribution $\rho(q^i,p_i)$ forms the infinitesimal probability $\rho d^n q d^n p$ that our system will be found in the infinitesimal volume $d^n q d^n p$. In other word, incompressibility condition is
$$\frac{d \rho}{dt}=\frac{\partial \rho}{\partial q^i} \dot{q}^i + \frac{\partial \rho}{\partial p_i}\dot{p}_i + \frac{\partial \rho}{\partial t}=0$$
Thus we conclude that
Incompressibility in phase space $\Leftrightarrow$ Probability conserved
Also $\rho$ obeys continuity equation in phase space :
$$\frac{\partial \rho}{\partial t} + \partial_{q^i} (\rho \dot{q}^i) + \partial_{p_i} (\rho \dot{p}_i)=0$$
This is the result of the incompressibility in the phase space. Thus $\vec{J}=(\rho,\rho \dot{q}^i, \rho \dot{p}_i)$ is a conserved current.





10. Canonical transformation

Let's consider an coordinate transformation which makes Hamilton's equation of motion invariant. In the Lagrangian mechanics we studied physics is invariant when the action is transformed up to constant. In this case, we call this transformation is 'canonical' or 'canonically transformed'. Assume a canonical transformation is defined
$$Q^i(t)=Q^i(q^j,t)$$
Alternatively, inverse transform also can be defined.
$$q^i(t)=q^i(Q^j,t)$$
Then the action is transformed up to constant.
$$S[Q^i(t)]=s[q^i(t)]+Const$$
It means Lagriangian is up to total time derivative of some function $G(q^i,p_i,t)$. Here I used capital notation is for transformed coordinate, and small letter is for the original coordinate.
$$p_i \dot{q}^i - h(q^i, p_i, t) = P_i Q^i - H(Q^i, P_i, t) + \frac{dG}{dt}$$
where $G$ is described using only 2 parameters of $p$, $P$, $q$, $Q$. For example, if we take $q$ and $Q$ as 2 parameters to describe generating function, then $G=G(q^i, Q^i, t)$. Multiplying above expression by $dt$ on both sides, we get
$$p_i dq^i - P_i dQ^i + (H-h)dt = dG \;\;\;\;\;\;\cdots \text{(1)}$$

(case 1) $G=G_1(q^i, Q^i, t)$
Since $dG=(\partial G/ \partial q^i) dq^i + (\partial G/ \partial Q^i) dQ^i + (\partial G/ \partial t) dt$, If $G(q^i, Q^i, t)$ is given, corresponding canonical transformation can be obtained by followings :
$$p_i=\frac{\partial G_1}{\partial q^i} \;\;\;\;\;\; P_i=-\frac{\partial G_1}{\partial Q^i} \;\;\;\;\;\; H=h+\frac{\partial G}{\partial t}$$
we can find the appropriate coordinate transformation rule $(q^i, p_i) \leftrightarrow (Q^i, P_i)$ and new Hamiltonian $H$.

(case 2) $G=G_2(q^i, P_i, t)$
We can rewrite (1) as shown below :
$$d(G+P_i Q^i)=p_i dq^i + Q^i dP_i +(H-h)dt$$
If we define $G'(q^i, P_i, t)=G+P_iQ^i$ and it is given properly, then the (2N+1) equations should be
$$p_i=\frac{\partial G'}{\partial q^i} \;\;\;\;\;\; Q_i=\frac{\partial G'}{\partial P_i} \;\;\;\;\;\; H=h+\frac{\partial G'}{\partial t}$$

(case 3) $G=G_3(p_i, Q^i, t)$
We can rewrite (1) as shown below :
$$d(G-q^i p_i)=-q^i dp_i - P_i dQ^i + (H-h)dt$$
If we define $G'(p_i, Q^i, t)=G-q^i p_i$ and it is given properly, then the (2N+1) equations should be
$$q^i = -\frac{\partial G'}{\partial p_i} \;\;\;\;\;\; P_i=-\frac{\partial G'}{\partial Q^i} \;\;\;\;\;\; H=h+\frac{\partial G'}{\partial t}$$

(case 4) $G=G_4(p_i, P_i, t)$
We can rewrite (1) as shown below :
$$d(G-q^i p_i + Q^i P_i)=-q^i dp_i + Q^i dP_i + (H-h)dt$$
If we define $G'(p_i, P_i, t)=G-q^i p_i + Q^i P_i$ and it is given properly, then the (2N+1) equations should be
$$q^i = -\frac{\partial G'}{\partial p_i} \;\;\;\;\;\; Q^i=\frac{\partial G'}{\partial P_i} \;\;\;\;\;\; H=h+\frac{\partial G'}{\partial t}$$

You can feel it is useless, but in some complex system, it plays sometimes a role of a cheat. Look at this example.

[ Example : SHO ]
Consider a Hamiltonian of the simple harmonic oscillator.
$$h(q,p,t)=\frac{p^2}{2m}+\frac{mw^2q^2}{2}$$
Generating function is given by
$$G(q,Q,t)=\frac{mwq^2}{2} \cot Q$$
We want to find canonical transform $(q,p) \to (Q,P)$ with a generating function $G(q,Q,t)$. For our transformation to be canonical, it must satisfy
$$p=\frac{\partial G}{\partial q}=mwq\cot Q \;\;\;\;\;\; P=-\frac{\partial G}{\partial Q}=\frac{mwq^2}{2}\csc^2 Q$$
Equivalently,
$$q=\sqrt{\frac{2P}{mw}}\sin Q \;\;\;\;\;\; p=\sqrt{2mwP} \cos Q$$
Then our new Hamiltonian is
$$H(Q,P,t)=wP$$
In this coordinates, $Q$ is a cyclic coordinate whose momentum $P$ is a constant of motion. $H$ has no explicit time dependence, so $H$ itself is a constant of motion. In conclusion, constant $P$ indicates conserved energy $P=E/w$.
$$\dot{Q}=\frac{\partial H}{\partial P} = w \;\;\;\;\;\Rightarrow\;\;\;\;\; Q=wt+\delta$$
Here we get the other constant of motion $\delta$ which means a initial phase. The original solution of this system can be recovered using the above.
$$q=\sqrt{\frac{2P}{mw}}\sin Q=\sqrt{\frac{2E}{mw^2}}\sin (wt+\delta)$$



추가예정 : Hamilton-Jacobi / Poisson bracket / Dirac bracket
예상 업로드일 : 11월 중

  • Share:

You Might Also Like

0 개의 댓글