# Point particle on rotating curve

## Point particle on arbitrary curve


Let us assume that point with mass $m=1$ moves freely on a flat curve which spins around vertical axis with angular velocity $\omega_0$.
In general case we have a curve given by an implicit equation in a form:
$$ f(x,y) = 0 $$ 

We can treat the problem in 2d in the rotating frame of reference. However, we need to take into account additional force in non-inertial system which acts in $x$ direction:

$$F_x = m x \omega_0^2.$$






```python
load('cas_utils.sage')
```


```python
var('t')
var('g')
var('w0', latex_name='\omega_0')
xy_names = [('x','x'),('y','y')]
```


```python
to_fun, to_var = make_symbols(xy_names)
```

dAlembert principle takes the following form:


```python
dAlemb = (X.diff(t,2)-x*w0^2)*dx + (Y.diff(t,2)+g)*dy
```


```python
showmath(dAlemb.subs(to_var))
```

We will need differential form of contraints, which can also be obtained by computer algebra. 
We will use implicit equation of a circle of unit radius. In this case we can compare solutions obtained in generic procedure with one computer in generalized coordinates.


```python
f = x^2+ (y-1)^2 - 1

plt_constraints = implicit_plot(f,(x,-1,1),(y,0,2),figsize=5,color='gray')
plt_constraints
```

We will use constaints to connect variations $\delta x$ and $\delta y$ with each other, therefore we use formula for complete differential:


```python
df = diff(f,x)*dx+diff(f,y)*dy
showmath(df)
```

We solve above equation for $\delta x$ and substitute it to dAlembert equation. Since there is only one variation left, then its coeficcient must be zero. It is a first differential equation.



```python
eq1 = dAlemb.subs(df.solve(dx)).coefficient(dy).subs(to_var)
showmath(eq1)
```

The second equation can be obtained, for example by differentiation of contraints in time two times.


```python
eq2 = f.subs(to_fun).diff(t,2).subs(to_var)
showmath(eq2)
```

We solve above system for $\ddot x$ and $\ddot y$:


```python
sol = solve( [eq1,eq2], [xdd,ydd])[0]
showmath( sol[0] )
```


```python
showmath( sol[1] )
```

We are ready to set numerical values of parameters and solve system of ODEs:


```python
times = srange(0,24,0.01)

ode = [xd,yd] + [ sol_.rhs().subs({w0:4.2,g:9.81}) for sol_ in sol]

x0 = 0.6
y0 = f.subs({x:x0}).solve(y)[1].rhs().n()
dx0 = 0.0
dy0 = df.subs({x:x0,y:y0,dx:dx0}).solve(dy)[0].rhs().n()

numsol = desolve_odeint(ode,[x0,y0,dx0,dy0], times, [x,y,xd,yd])
p = line(zip(numsol[:,0],numsol[:,1]), color='red')  
p2 = line(zip(times,numsol[:,0]), figsize=5, ymin=-1, ymax=1)

(p+plt_constraints).show(title='phase space')
p2.show(title=r'$x(t)$')
```

**Experiment with Sage!**

Change parameters to see if other types of solution can be found!

### Explicit form of contraints

If contraits are in explicit for we can simplify the procedure.


```python
# f = 1/2*(x)^2 - y
f = 1-sqrt(1-x^2) - y
```


```python
df = diff(f,x)*dx+diff(f,y)*dy
```


```python
eq1 = dAlemb.subs(df.solve(dx)).coefficient(dy).subs(to_var)
showmath(eq1)
```


```python
eq2 = f.subs(to_fun).diff(t,2).subs(to_var)
showmath(eq2)
```


```python
sol = solve( [eq1,eq2], [xdd,ydd])[0]
showmath( sol[0] )
```

We see that the first equation does not depent on $y$. We can solve it numericaly.


```python
y_on_curve = y.subs(f.solve(y)[0])
```


```python
y_fun = fast_callable(y_on_curve, vars = [x])
```


```python
y_fun(numsol[:,0])
```


```python
ode = [xd,sol[0].rhs().subs({w0:4.2,g:9.81})]

times = srange(0,24,0.01)
numsol = desolve_odeint(ode,[.6,.0],times,[x,xd])
p = line(zip(numsol[:,0],y_fun(numsol[:,0])),color='red')
p2 = line(zip(times,numsol[:,0]),figsize=5,ymin=-1,ymax=1)

(p+plt_constraints).show(title='phase space')
p2.show(title=r'$x(t)$')

```

We see that it was possible to obtain the same solution. However, the implicit form allows to solve the system for  $y>1$. The explicit formulation will fail in wuch a case.

## Point particle on rotating circle


If particle moves on the rotating circle, we can use the symmetry and transform dAlembert principle to  polar coordinates in rotating frame.


```python
var('x y t')
var('w0 l g')

xy_names = ['x','y']
uv_names = [('phi','\phi')]

to_fun, to_var = make_symbols(xy_names,uv_names)
```


```python
x2u = {x:l*sin(phi),y:1-l*cos(phi)}
```


```python
showmath(x2u)
```


```python
transform_virtual_displacements(xy_names,uv_names,verbose=True)
```


```python
dAlemb = (x.subs(x2u).subs(to_fun).diff(t,2)-w0^2*x.subs(x2u))*dx_polar + \
         (y.subs(x2u).subs(to_fun).diff(t,2)+g)*dy_polar
dAlemb = dAlemb.subs(to_var)
show(dAlemb)
```


```python
sol = dAlemb.expand().coefficient(dphi).trig_simplify().solve(phidd) 
show( sol)
```

### Efffective potential

In polar coordinates we can analyze the system without solving. Namely we can observe that the generalized force depends only on $\phi$, thus:
\begin{equation}
{\ddot \phi} = \underbrace{\frac{l {\omega_0}^{2} \cos\left({\phi}\right) \sin\left({\phi}\right) - g \sin\left({\phi}\right)}{l}}_{-\displaystyle\frac{\partial U_{eff}(\phi)}{\partial\phi}}
\end{equation}

we can interpret the equation of motion as a motion in an effective potential. 

Let's extrat this potential and see how it behaves for different parameters:


```python
Ueff = -sol[0].rhs().expand().integrate(phi)
showmath(Ueff)
```


```python
Ueff = -sol[0].rhs().expand().subs({g:9.81,l:1}).integrate(phi)
showmath(Ueff)
```


```python
plot( [Ueff.subs(w0==w0_) for w0_ in [0,4,8,12]], (phi,-pi,pi),\
     legend_label=[0,4,8,12], figsize=(6,3),\
     tick_formatter=[pi,None],ticks=[pi/2,None])
```

### Numerical solutions


```python
arcsin(0.4)
```


```python
times = srange(0,24,0.01)
pars = {l:1, w0:4.2, g:9.81}

x_phi = fast_callable( x.subs(x2u).subs(pars),vars=[phi])
y_phi = fast_callable( y.subs(x2u).subs(pars),vars=[phi])

ode = [phid,sol[0].rhs().subs(pars)]
show(ode)
numsol = desolve_odeint(ode,[arcsin(0.6),0],times,[phi,phid])
p = line(zip(x_phi(numsol[:,0]),y_phi(numsol[:,0])),color='red')
p2 = line(zip(times,x_phi(numsol[:,0])),figsize=4,ymin=-pi/2,ymax=pi/2)

(p+plt_constraints).show()
p2.show()
```

## Lagrange approach

### Rotating system of coordinates

In a rotating system of coordinates, the bead moves on the circle in $(x,y)$ plane. At the same time there is a rotation with angular velocity $\omega_0$ around vetical axis. It means that a bead will have addtional velocity, perpedicular to $(x,y)$ plane, which magnitute is:
$$
v_{rot} = x\omega_0 = l {\omega_0} \sin\left({\phi}\right)
$$


```python
Ek = 1/2*sum([x_.subs(x2u).subs(to_fun).diff(t).subs(to_var)^2 for x_ in [x,y]])
Ek += 1/2*(x.subs(x2u)*w0)^2
Ek = Ek.trig_simplify()
showmath(Ek)
```


```python
Ep = g*y.subs(x2u)
showmath(Ep)
```


```python
L = Ek-Ep
showmath(L)
```


```python
EL = L.diff(phid).subs(to_fun).diff(t).subs(to_var) - L.diff(phi)
```


```python
showmath(EL)
```


### Code generation 

We can readily generate code which can be used in external programs.

En example can be found in [simulation and 3d vis](beadk3d.ipynb)




```python
from sympy import ccode
oderhs = sol[0].rhs()
ccode(oderhs._sympy_())
```

\newpage
