Foucault's pendulum

19 Jul 2023 - fubar - Sreekar Guddeti

scipy-Logo Invented by Leon Foucault in 1851, it demonstrates that the Earth has a rotation axis.

System

The pendulum has a bob suspended by a wire that can only apply a tensile force $\vec{\mathbf{T}}$. An additional force due to gravity $\vec{\mathbf{F}}_\mathrm{g}$ acts on the bob.

Inertial frame

In an inertial frame $F_\mathrm{A}$ fixed at the centre of earth, the pendulum is located at the surface of the earth given by the coordinates $\vec{\mathbf{R}}= R_0 \hat{ \mathbf{R}} + \Theta_0\hat{\mathbf{\Theta}} + \Phi_0 \hat{\mathbf{\Phi}}$, where $R_0$ is the radius of the Earth, $(\Theta_0, \Phi_0)$ is the polar latitude and longitude, $(\hat{\mathbf{R}}, \hat{\mathbf{\Theta}}, \hat{\mathbf{\Phi}}) $ is the spherical polar coordinates. The angular velocity of the earth is $\vec{\mathbf{\Omega}} = \Omega_0 \hat{\mathbf{R}}(\Theta_0=0)$, where $\Omega_0$ is the angular speed.

Equation of motion

In frame $F_\mathrm{A}$, therefore, there are two forces acting on the bob with mass $m$. Hence the net force $\vec{\mathbf{F}}_\mathrm{A}$ is given by

\[\vec{\mathbf{F}}_\mathrm{A} = \vec{\mathbf{T}} + \vec{\mathbf{F}}_\mathrm{g}. \tag{1a}\]

The equation of motion is given by

\[\vec{\mathbf{F}}_\mathrm{A} = m \vec{\mathbf{a}}_\mathrm{A},\tag{1b}\]

where $\vec{\mathbf{a}}_\mathrm{A}$ is the acceleration of bob.

The unit vectors are also rotating with

\[\vec{\mathbf{\Omega}} = \Omega_0 (\cos \Theta_0 \hat{\mathbf{R}} - \sin \Theta_0 \hat{\mathbf{\Theta}}) \tag{1c}\]

Orbiting and rotating frame

For a frame $F_\mathrm{B}$ located at the base of pendulum, let the coordinate axes be $(\hat{\mathbf{r}}, \hat{\mathbf{\theta}}, \hat{\mathbf{z}})$ where the base plane of the cylindrical coordinates aligns with the $R=R_0$ plane of $F_\mathrm{A}$ and the $\hat{\mathbf{r}}(\theta=0)$ unit vector aligns with the $\hat{\mathbf{\Theta}}$ unit vector. As a result, the origin of $F_\mathrm{B}$ orbits as well as rotates around $F_\mathrm{A}$ with $\Omega_0$.

If the acceleration of bob in $F_\mathrm{B}$ is $\vec{\mathbf{a}}_\mathrm{B}$, then

\[\vec{\mathbf{a}}_\mathrm{A} = \vec{\mathbf{a}}_\mathrm{B} + \vec{\mathbf{a}}_\mathrm{BA} + \vec{\mathbf{a}}_\mathrm{Coriolis} + \vec{\mathbf{a}}_\mathrm{centrifugal} \tag{2a},\]

where $\vec{\mathbf{a}}\mathrm{BA}$ is the acceleration of origin of $F\mathrm{B}$ due to orbiting and given by

\[\vec{\mathbf{a}}_\mathrm{BA} = - \Omega_0 \vec{\mathbf{R}}_\mathrm{eff} , \tag{2b}\]

where $\vec{\mathbf{R}}_\mathrm{eff}$ is the effective radius of bob from the Earth’s axis and given by

\[\vec{\mathbf{R}}_\mathrm{eff} = R_\mathrm{eff} \hat{\mathbf{R}}_\mathrm{eff}, \tag{2c}\]

where

\[\begin{array}{} R_\mathrm{eff} &= R_0 \sin \Theta_0 , \\ \hat{\mathbf{R}}_\mathrm{eff} &= \sin \Theta \hat{\mathbf{r}} + \cos \Theta \hat{\mathbf{\Theta}}, \end{array}\]

and $\vec{\mathbf{a}}_\mathrm{Coriolis}$ is the Coriolis acceleration arising due to relative velocity and given by

\[\vec{\mathbf{a}}_\mathrm{Coriolis} = 2 (\vec{\mathbf{\Omega}} \times \vec{\mathbf{v}}_\mathrm{B}), \tag{2d}\]

where $\vec{\mathbf{v}}\mathrm{B}$ is the velocity of the particle in $F\mathrm{B}$ and $\vec{\mathbf{a}}_\mathrm{centrifugal}$ is the centripetal acceleration arising due to relative radial extent and given by

\[\vec{\mathbf{a}}_\mathrm{centripetal} = \vec{\mathbf{\Omega}} \times (\vec{\mathbf{\Omega}} \times \vec{\mathbf{r}}_\mathrm{B}), \tag{2e}\]

$\vec{\mathbf{r}}\mathrm{B}$ is the position of the particle in $F\mathrm{B}$.

Reduced equation of motion

In frame $F_\mathrm{A}$, the equation of motion is given by

\[\vec{\mathbf{F}}_\mathrm{B} = m \vec{\mathbf{a}}_\mathrm{B}, \tag{3a}\]

where $\vec{\mathbf{F}}_\mathrm{B}$ is the net force. Using Equations 1b and 2a, the net force is given by

\[\vec{\mathbf{F}}_\mathrm{B} = \vec{\mathbf{F}}_\mathrm{A} + \vec{\mathbf{F}}_\mathrm{pseudo}, \tag{3b}\]

where $\vec{\mathbf{F}}_\mathrm{pseudo}$ is the additional pseudo or fictitious force experienced and given by

\[\vec{\mathbf{F}}_\mathrm{pseudo} = - m \vec{\mathbf{a}}_\mathrm{BA} - m \vec{\mathbf{a}}_\mathrm{Coriolis} - m \vec{\mathbf{a}}_\mathrm{centripetal}. \tag{3c}\]

Transformation of coordinate systems

In terms of $(\hat{\mathbf{r}}, \hat{\mathbf{\theta}}, \hat{\mathbf{z}})$, the $(\hat{\mathbf{R}}, \hat{\mathbf{\Theta}}, \hat{\mathbf{\Phi}}) $ is given by

\[\begin{array}{} \hat{\mathbf{R}} &= \hat{\mathbf{z}},\\ \hat{\mathbf{\Theta}} &= \cos \Theta_0 \hat{\mathbf{R}} - \sin \Theta_0 \hat{\mathbf{\Theta}}, \\ \hat{\mathbf{\Phi}} &= \sin \Theta_0 \hat{\mathbf{R}} + \cos \Theta_0 \hat{\mathbf{\Theta}}. \end{array} \tag{4}\]

Small angle approximation

If we assume that the length of the pendulum $L_0$ is large enough compared to the maximum radial motion of bob, then we can approximate the motion to the $z=0$ plane. Under this approximation, we can neglect the $z-$components of position, velocity and acceleration of bob in $F_B$.

Further under this approximation, the $r-$ and $\theta-$ components have an angular velocity $\vec{\mathbf{\omega}}\mathrm{B}$ and acceleration $\vec{\mathbf{\alpha}}\mathrm{B}$ only along the $z-$ component and given by

\[\begin{array} \vec{\mathbf{\omega}} &= \omega_\mathrm{B} \hat{\mathbf{z}}, \\ \vec{\mathbf{\alpha}} &= \alpha_\mathrm{B} \hat{\mathbf{z}}, \end{array} \tag{5a}\]

As a result, the position, velocity and acceleration is reduced to

\[\begin{array} \vec{\mathbf{r}}_\mathrm{B} &= r \hat{\mathbf{r}}, \\ \vec{\mathbf{v}}_\mathrm{B} &= \dot{r} \hat{\mathbf{r}} + r \omega \hat{\mathbf{\theta}}, \\ \vec{\mathbf{a}}_\mathrm{B} &= (\ddot{r} - r \omega^2 ) \hat{\mathbf{r}} + (r \alpha + 2 \dot{r}\omega) \hat{\mathbf{\theta}}. \end{array} \tag{5b}\]

It is to be noted that the subscript B in the radius $f$ and angular speed $\omega$ and angular acceleration $\alpha$ is dropped for brevity.

Degrees of freedom

Hence, under the small angle approximation, we have two degrees of freedom, the radial extent $r$ and the angular deviation $\theta$. Hence the generalized coordinates are

\[\vec{\mathbf{q}} = \{r, \theta\},\]

and the corresponding velocities are

\[\dot{\vec{\mathbf{q}}} = \{\dot{r}, \omega\}.\]

The objective is to find the equations for the rate of change of the generalized coordinates and their corresponding velocities.

Free Body diagram

Let the tension $\vec{\mathbf{T}}$ be

\[\vec{\mathbf{T}} = T_\mathrm{r} \hat{\mathbf{r}} + T_\mathrm{\theta} \hat{\mathbf{\theta}} + T_\mathrm{z} \hat{\mathbf{z}}.\]

The force due to gravity is

\[\vec{\mathbf{F}}_\mathrm{g} = -m g \hat{\mathbf{z}}.\]

Using the Equations 2 and 5b, the terms in the pseudo force are given by

\[\begin{array}{} - m \vec{\mathbf{a}}_\mathrm{BA} &= m \Omega^2 R_\mathrm{eff} \left( \cos \Theta \cos \theta \hat{\mathbf{r}} - \cos \Theta \sin \theta \hat{\mathbf{\theta}} + \sin \Theta \hat{\mathbf{z}} \right), \\ - m \vec{\mathbf{a}}_\mathrm{Coriolis} &= 2 m \Omega \left( r \omega \cos \Theta \hat{\mathbf{r}} - \dot{r} \cos \Theta \hat{\mathbf{\theta}} + \overline{\dot{r}\sin\Theta \sin\theta + r \omega \sin \Theta \cos \theta} \hat{\mathbf{z}} \right), \\ - m \vec{\mathbf{a}}_\mathrm{centripetal} &= m r \Omega^2 \left( \overline{\cos^2 \Theta + \sin^2 \Theta \sin^2 \theta} \hat{\mathbf{r}} + \sin^2\Theta \cos \theta \sin \theta \hat{\mathbf{\theta}} + \sin \Theta \cos^2 \theta \hat{\mathbf{z}} \right). \end{array}\]

Hence the resolution of the reduced equation of motion into its components is given by

\[\left( \begin{array}{} T_\mathrm{r} + m \Omega^2 R_\mathrm{eff} \cos \Theta \cos \theta + 2 m \Omega r \omega \cos \Theta + m r \Omega^2 \overline{\cos^2 \Theta + \sin^2 \Theta \sin^2 \theta} \\ T_\mathrm{\theta} - m \Omega^2 R_\mathrm{eff}\cos \Theta \sin \theta - 2 m \Omega \dot{r} \cos \Theta + m r \Omega^2 \sin^2\Theta \cos \theta \sin \theta \\ T_\mathrm{z} - m g + m \Omega^2 R_\mathrm{eff} \sin \Theta + 2 m \Omega \overline{\dot{r}\sin\Theta \sin\theta + r \omega \sin \Theta \cos \theta} + m r \Omega^2 \sin \Theta \cos^2 \theta \end{array} \right) = \left( \begin{array}{} \ddot{r} - r \omega^2 \\ r \alpha + 2 \dot{r}\omega \\ 0 \end{array} \right) \tag{6}\]

Equation of constraint

Due to the tensile nature of $\vec{\mathbf{T}}$,

\[T_\mathrm{\theta} = 0 \tag{7a}\]

and

\[T_\mathrm{z} \tan \varphi = - T_\mathrm{r}, \tag{7b}\]

where $\varphi$ is the angular deviation of the wire from the mean position and given by

\[\varphi = \frac{r}{L_0}. \tag{7c}\]

Under small angle approximation,

\[\tan \varphi \approx \varphi.\]

Analytical solution

Using expressions for $T_\mathrm{r}, T_\mathrm{z}$ from Equation 6 and subsituting in Equations. 7a, 7b we have

\[\left( \begin{array}{} \ddot{r} \\ r \alpha \end{array} \right) = \left( \begin{array}{} a_\mathrm{r} + r \omega^2 \\ a_\mathrm{\theta} - 2 \dot{r}\omega \end{array} \right)\]

where $m a_\mathrm{r}, m a_\mathrm{\theta}$ are respectively the radial and tangential components of the reduced net force and given by

\[\left( \begin{array}{} a_\mathrm{r} \\ a_\mathrm{\theta} \end{array} \right) = \left( \begin{array}{} - \frac{r}{L_0} g + R_\mathrm{eff} \Omega^2 \left( \frac{r}{L_0} \sin \Theta + \cos \Theta \cos \theta \right) + 2 \Omega \left( \frac{r}{L_0} \overline{ \dot{r} \sin \Theta \sin \theta + r \omega \sin \Theta \cos \theta } + r \omega \cos \Theta \right) + r \Omega^2 \left( \frac{r}{L_0} \sin \Theta \cos^2 \theta + \overline{ \cos^2 \Theta + \sin^2 \Theta \sin^2 \theta } \right) \\ - R_\mathrm{eff} \Omega^2 \cos \Theta \sin \theta - 2 \Omega \dot{r} \cos \Theta + r \Omega^2 \sin^2 \Theta \cos\theta \sin \theta. \end{array} \right)\]

If we denote the generalized coordinates and the corresponding velocities as generic variables

\[\{ r, \theta, \dot{r}, \omega \} = \{y_1, y_2, y_3, y_4 \},\]

then then the rate of change of these variables is given by

\[\left( \begin{array}{} \dot{y}_1 \\ \dot{y}_2 \\ \dot{y}_3 \\ \dot{y}_4 \end{array} \right) = \left( \begin{array}{} y_3 \\ y_4 \\ a_\mathrm{r} + y_1y_4^2 \\ \frac{1}{y_1} \left(a_\mathrm{\theta} - 2 y_3y_4 \right) \end{array} \right),\]

where

\[\left( \begin{array}{} a_\mathrm{r} \\ a_\mathrm{\theta} \end{array} \right) = \left( \begin{array}{} - \frac{y_1}{L_0} g + R_\mathrm{eff} \Omega^2 \left( \frac{y_1}{L_0} \sin \Theta + \cos \Theta \cos y_2 \right) + 2 \Omega \left( \frac{y_1}{L_0} \overline{ y_3 \sin \Theta \sin y_2 + y_1 y_4 \sin \Theta \cos y_2 } + y_1 y_4 \cos \Theta \right) + y_1 \Omega^2 \left( \frac{y_1}{L_0} \sin \Theta \cos^2 y_2 + \overline{ \cos^2 \Theta + \sin^2 \Theta \sin^2 y_2 } \right) \\ - R_\mathrm{eff} \Omega^2 \cos \Theta \sin y_2 - 2 \Omega y_3 \cos \Theta + y_1 \Omega^2 \sin^2 \Theta \cos y_2 \sin y_2. \end{array} \right)\]

Hence, we have a set of four first-order non-linear differential equations in four variables that can be succintly in the vector form

\[\dot{\vec{\mathbf{y}}} = \tilde{\mathrm{D}}\vec{\mathbf{y}}, \tag{8}\]

where $\vec{\mathbf{y}} = (y_1, y_2, y_3, y_4)$ is generalized position in the coordinate space and $\tilde{\mathrm{D}}$ is the differential operator. If we know the initial condition $\vec{\mathbf{y}}^0$, we can in principle solve the first order differential equation.

Reality checks

In order to test the veracity of the above analysis, let us compare the form of the equations in limiting conditions with physical behaviour.

$\Omega_0 =0$ limit

If the Earth is not rotating, then the motion of the bob should be akin to that of a simple pendulum in 2D.

\[\left( \begin{array}{} a_\mathrm{r} \\ a_\mathrm{\theta} \end{array} \right) = \left( \begin{array}{} - \frac{r}{L_0} g \\ 0 \end{array} \right).\]

Numerical solution

We can solve Equation 8 using Scipy package’s interpolate’ module solve_ivp function.

import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d
import matplotlib as mpl
import matplotlib.pyplot as plt

# magic command to convert matplotlib plots into Base64 string images

%matplotlib notebook 

# Constants

RADIUS = 6.4e6 # Radius of earth (in m)
DAY = 24 * 3600 # time period (in s)
OMEGA = 2*np.pi/DAY
LATITUDE = 61.4/180*np.pi # polar latitude (in rad)
MASS = 1 # mass of pendulum (in kg)
LENGTH = 100 # length of pendulum (in m)
GRAVITY = 9.8 # acceleration due to gravity (in ms-2)

# Default particle trajectory conditions

STATE_0 = (1, 0, 0, 0)
T_SPAN = (0, 100)
MIN_TIMESTAMPS = T_SPAN[1]*10

class Particle:
    """
    This particle is under the influence of a radial force.
    Particle's attributes are -
        - Mass
        - Initial coordinates
        - External force
        - time interval
    Particle's methods are -
        - solution
        - trajectory
    """
    def __init__(self, force, mass=MASS, state0=STATE_0, t_span=T_SPAN):
        self.mass = mass
        self.state0 = state0
        self.t_span = t_span
        self.force = force
        self.ts, self.states = self.solution()

    def motion(self, t, state):
        """
        The velocity vector for a particle under a 2D radial force.
        The parameter t is used by the ODE solver from scipy.
        """
        m = self.mass
        # generalized coordinates
        y1, y2, y3, y4 = state
        
        # force components
        f_r, f_t = self.force(state)

        v1 = y3
        v2 = y4
        v3 = f_r/m + y1*y4**2
        v4 = 1/y1*(f_t/m - 2*y3*y4)
        
        # generalized velocity
        v = (v1, v2, v3, v4)

        return v
    
    def solution(self):
        """
        ODE solver from scipy.
        """
        global MIN_TIMESTAMPS
        max_step = self.t_span[1]/MIN_TIMESTAMPS
        solver = solve_ivp(self.motion, self.t_span, self.state0, max_step=max_step)
        ts, states = solver.t, solver.y
        
        return ts, states
    
    def trajectory(self):
        """
        Trajectory of motion.
        """
        # time stamps
        ts = self.ts
        
        # Generalized coordinates.
        radii, thetas, speeds, omegas = self.states
        
        # time profiles of coordinates.
        
        fig, ax = plt.subplots()
        ax.plot(ts, radii, label='r')
        ax.plot(ts, thetas%(2*np.pi), label=r'$\theta$')
        ax.plot(ts, speeds, label='v')
        ax.plot(ts, omegas, label=r'$\omega$')
        ax.axhline(0, color ='k')
        ax.legend()
        
        # Cartesian coordinates
        
        xs = [radius*np.cos(theta) for radius, theta in zip(radii, thetas)]
        ys = [radius*np.sin(theta) for radius, theta in zip(radii, thetas)]
        
        # Interpolate the data
        f = interp1d(xs, ys, kind='cubic')
        xs_int = xs
        ys_int = f(xs_int)
        
        # Trajectory
        fig1, ax1 = plt.subplots()
        ax1.plot(xs, ys, 'o', label='trajectory')
        ax1.plot(xs_int, ys_int, '-')
        ax1.axhline(0, color ='k')
        ax1.axvline(0, color ='k')
        ax1.set_aspect('equal')
        ax1.legend()
        

Forces

Any force in general depends on the coordinates ${ r, \theta, \dot{r}, \omega }$. Let us define commonly experienced forces and test the validity of the trajectory solver of the Particle class as defined above.

Particle under Hooke’s force

To test the viability of the motion solver, let us test the motion of a particle in 2D under Hooke’s force

\[\vec{\mathbf{F}}_\mathrm{Hooke} = -r\]

The solution is a Simple Harmonic Motion (SHM) in 2D with angular speed $\omega = 1$ as shown below.

# hooke's law

def hooke(state):
    """
    The force in general depends on the coordinates (r, theta)=(y1, y2) and the velocity (v, omega)=(y3, y4).
    This is a radial force.
    """
    y1, y2, y3, y4 = state
    
    # force components in polar coordinates.
    f_r = -1*y1
    f_t = 0
    
    force = (f_r, f_t)
    
    return force

# Particle under hooke's law

state_0 = (1, 0, 0, 0)
t_span = (0, 20)
MIN_TIMESTAMPS = T_SPAN[1]*10

p = Particle(hooke, state0=state_0, t_span=t_span)

p.trajectory()

The time period of SHM is given by

\[T_\mathrm{SHM} = 2 \pi / \omega.\]

For a particle having a tangential component of velocity, we have SHM in 2D.

Tangential trajectory

For a particle having a tangential component of velocity, we have SHM in 2D. Consider a particle with initial conditions ${ r_0, \theta_0, \dot{r}_0, \omega_0 } = {1, 0, 1, 2}$ so that it also has zero component of tangential velocity, thereby, resulting in a tangential trajectory.

# tangential trajectory

state0 = (1, 0, 1, 2)
t_span = (0, 20)
MIN_TIMESTAMPS = t_span[1]*10
p_hooke1 = Particle(hooke, state0=state0, t_span=t_span)

p_hooke1.trajectory()

NOTE: The time period is still the same.

Particle under gravitation force

To check the robustness of the motion solver in polar coordinates, let us test the motion of a particle in 2D under the gravitational force

\[\vec{\mathbf{F}}_\mathrm{Newton} = - \frac{1}{r^2}\]

The solution is a conic as shown below:

def newton(state):
    """
    inverse square force that is radial in nature.
    """
    
    y1, y2, y3, y4 = state
    
    # force components in polar coordinates.
    f_r = -1/y1**2
    f_t = 0
    
    force = (f_r, f_t)
    
    return force

Escape speed

If the particle has sufficient kinetic energy to overcome the gravitational potential particle, it will escape. The minimum initial kinetic energy required for escape must be such that at infinity the total energy is zero so that

\[\frac{1}{2} m v_0^2 - \frac{1}{r_0^2} = T_\infty + U_\infty = 0,\]

where $v_0$ is the initial speed given by

\[v_0^2 = \dot{r}^2 + r^2\omega^2.\] \[\therefore \quad v_0 = \frac{1}{r_0}\sqrt{\frac{2}{m}}.\]

m = 1
r_0 = 1
r_dot = 0
v_0 = np.sqrt(2/m)/r_0 

v_0

1.4142135623730951

Therefore, the escape speed is 1.414.

Bounded particle

To be bounded, therefore, consider a particle with initial conditions ${ r_0, \theta_0, \dot{r}_0, \omega_0 } = {1, 0, 0, 1.3}$ so that it has initial speed less than escape speed, thereby, resulting in a bounded trajectory.

# tangential trajectory

state0 = (1, 0, 0, 1.3)
t_span = (0, 20)
MIN_TIMESTAMPS = t_span[1]*10
p_newton = Particle(newton, state0=state0, t_span=t_span)

p_newton.trajectory()

NOTE: The trajectory for a bounded particle is elliptical.

Unbounded particle

To be unbounded, therefore, consider a particle with initial conditions ${ r_0, \theta_0, \dot{r}_0, \omega_0 } = {1, 0, 0, 1.5}$ so that it has initial speed greater than escape speed, thereby, resulting in a unbounded trajectory.

state0 = (1, 0, 0, 1.5)
t_span = (0, 6)
MIN_TIMESTAMPS = t_span[1]*10
p_newton1 = Particle(newton, state0=state0, t_span=t_span)

p_newton1.trajectory()

NOTE: The trajectory for an unbounded particle is parabolic.

Foucault’s pendulum

The constants relevant to the motion are

radius of Earth $R$
angular speed of Earth $\Omega$
polar latitude of location of pendulum $\theta$
mass of pendulum $m$
length of pendulum $L$
acceleration due to gravity $g$

The net force acting on the Foucault pendulum is described by Equation 8. The solution is computed below.

# foucault force

def foucault(state):
    """
    Foucault force that is the net combination of gravity and pseudo force experience due to observer located at the location.
    """
    y1, y2, y3, y4 = state
    
    # force components in polar coordinates.
    global MASS
    global RADIUS
    global LATITUDE
    global OMEGA
    global GRAVITY
    global LENGTH
    
    m = MASS
    R = RADIUS
    th = LATITUDE
    w = OMEGA
    g = GRAVITY
    l = LENGTH
    
    r_eff = R * np.sin(th)
    f_cf = m*r_eff*w**2
    dy1= y1/l
    
    a_g = -1*g*dy1
    a_o = r_eff*w**2 * (dy1*np.sin(th) + np.cos(th)*np.cos(y2))
    a_c = 2*w*(dy1*(y3*np.sin(th)*np.sin(y2) + y1*y4*np.sin(th)*np.cos(y2)) + y1*y4*np.cos(th))
    a_cf = y1*w**2*(dy1*np.sin(th)*np.cos(y2)**2 + (np.cos(th)**2 + np.sin(th)**2 * np.sin(y2)**2))
    
    
    a_r = a_g + a_o + a_c + a_cf
    
    a_t = -r_eff*w**2*np.cos(th)*np.sin(y2) - 2*w*y3*np.cos(th) + y1*w**2*np.sin(th)**2*np.cos(y2)*np.sin(y2)
    
    f_r = m*a_r
    f_t = m*a_t
    
    force = (f_r, f_t)
    
    return force
    

# particle under Foucault force
state0 = (1, 0, 0, 0)
t_span = (0, 3600)
MIN_TIMESTAMPS = t_span[1]*10
pendulum = Particle(foucault, state0=state0, t_span=t_span)

# trajectory of pendulum

pendulum.trajectory()

Stable foucault pendulum

NOTE: At zero crossing, the solver becomes unstable due to the presence $1/r$ term.

Hence, it is recommended to provide a slight tangential component to the initial velocity to avoid zero crossing.

# tangential trajectory

state0 = (1, 0, 0, 0.05)
t_span = (0, 100)
MIN_TIMESTAMPS = t_span[1]*10
stable_pendulum = Particle(foucault, state0=state0, t_span=t_span)

# trajectory of pendulum

stable_pendulum.trajectory()

The approximate time period of oscillation within the so-called plane of oscillation is roughly equal to the time period of SHM motion

\[T_\mathrm{SHM} = 2 \pi \sqrt{\frac{L_0}{g}}\]

# approximate time period of oscillation

T = 2*np.pi*np.sqrt(LENGTH/GRAVITY)

20.07089923154493

Faster rotation

The effect of Earth’s rotation is not evident due to the low angular speed. Let us artificially speed up the Earth’s rotation by 10 times.

OMEGA = 2*np.pi/DAY*10
state0 = (1, 0, 0, 0.05)
t_span = (0, 360)
MIN_TIMESTAMPS = t_span[1]*10
p_at_10w = Particle(foucault, state0=state0, t_span=t_span)

p_at_10w.trajectory()

Hence, the plane of oscillation rotates in a clockwise manner.

Period of gyration

How to calculate it analytically?

Effect of latitude

As the mean position is off-centred from origin, let us look at the effect of latitude on the mean position

particle at north pole

At the north pole, the non-inertial frame $F_\mathrm{B}$ only rotates and does not orbit. Hence the only pseudo forces are due to the Coriolis and centripetal acceleration of $F_\mathrm{B}$.

… … (needs analysis)

The angular speed $\omega_\mathrm{B}$ turns out to be

\[\omega_\mathrm{B} = - \Omega\]

To demonstrate this condition numerically, let us solve the equation of motion until $\Delta T = 1/4 T_0$, where

\[T_0 = \frac{2\pi}{\Omega}\]

is one Earth Day. In this time, the plane of oscillation should rotate by $\pi/2$.

LATITUDE = 0
OMEGA = 2*np.pi/DAY
state0 = (10, 0, 0, 0.05)
t_span = (0, 6*3600)
MIN_TIMESTAMPS = 6*3600*10
p_at_0 = Particle(foucault, state0=state0, t_span=t_span)

p_at_0.trajectory()

The plane of oscillation has rotated clockwise by $\pi/2$. Hence, the plane of oscillation has angular speed of $-\Omega$.

particle at equator

At the equator, the non-inertial frame $F_\mathrm{B}$ rotates as well as orbits. Hence pseudo forces are due to the accelerating origin, Coriolis acceleration and centripetal acceleration of $F_\mathrm{B}$. However, the equation of motion is reduced to

...
... (needs analysis)

The angular speed on the average should be zero.

To demonstrate this condition numerically, let us solve the equation of motion until $\Delta T = 1/4 T_0$.

LATITUDE = np.pi/2
state0 = (10, 0, 0, 0.05)
t_span = (0, 3600)
MIN_TIMESTAMPS = 3600*10
p_at_90 = Particle(foucault, state0=state0, t_span=t_span)

p_at_90.trajectory()

particle at 61.4 $^{o}$ latitude

Consider the pendulum at the new Parliament building of India. The latitude is $61.4 ^{o}$.

Parliament building

The angular speed on the average should be $\Omega \cos\Theta$. To demonstrate this condition numerically, let us solve the equation of motion until

\[\Delta T = \frac{1}{4} \frac{T_0}{\cos \Theta}.\]

In this time, the plane of oscillation should rotate by $\pi/2$.

LATITUDE = 61.4/180*np.pi
state0 = (10, 0, 0, 0.05)
t_span = (0, DAY/np.cos(LATITUDE)/4)
MIN_TIMESTAMPS = t_span[1]*10
p_at_parliament = Particle(foucault, state0=state0, t_span=t_span)

p_at_parliament.trajectory()

Effect of $\omega_0$ on the trajectory

With no initial angular speed, in the absence of rotation of Earth, the motion is simple 2D harmonic motion. Let us check how the motion changes with a finite initial angular speed

$\omega_0 = 0$

OMEGA = 0
state0 = (1, np.pi/2, 0, 0)
t_span = (0, 100)
MIN_TIMESTAMPS = t_span[1]*10
p_at_w_0 = Particle(foucault, state0=state0, t_span=t_span)

p_at_w_0.trajectory()

/tmp/ipykernel_20147/3747781934.py:65: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig, ax = plt.subplots()

Cartesian coordinates

To avoid the instability in integration at cross-over, let us reformulate the equation of motion in the Cartesian coordinates.

class CartesianParticle:
    """
    This particle is under the influence of a 2D force.
    Particle's attributes are -
        - Mass
        - Initial coordinates
        - External force
        - time interval
    Particle's methods are -
        - solution
        - trajectory
    """
    def __init__(self, force, state0, t_span, mass=MASS):
        self.mass = mass
        self.state0 = state0
        self.t_span = t_span
        self.force = force
        self.ts, self.states = self.solution()

    def motion(self, t, state):
        """
        The velocity vector for a particle under a 2D radial force.
        The parameter t is used by the ODE solver from scipy.
        """
        m = self.mass
        # generalized coordinates
        y1, y2, y3, y4 = state
        
        # force components
        f_x, f_y = self.force(state)

        v1 = y3
        v2 = y4
        v3 = f_x/m 
        v4 = f_y/m
        
        # generalized velocity
        v = (v1, v2, v3, v4)

        return v
    
    def solution(self):
        """
        ODE solver from scipy.
        """
        global MIN_TIMESTAMPS
        max_step = self.t_span[1]/MIN_TIMESTAMPS
        solver = solve_ivp(self.motion, self.t_span, self.state0, max_step=max_step)
        ts, states = solver.t, solver.y
        
        return ts, states
    
    def trajectory(self):
        """
        Trajectory of motion.
        """
        # time stamps
        ts = self.ts
        
        # Generalized coordinates.
        xs, ys, v_xs, v_ys = self.states
        
        # time profiles of coordinates.
        
        fig, ax = plt.subplots()
        ax.plot(ts, xs, label='x')
        ax.plot(ts, ys, label='y')
        ax.plot(ts, v_xs, label=r'$v_x$')
        ax.plot(ts, v_ys, label=r'$v_y$')
        ax.axhline(0, color ='k')
        ax.legend()
        
        # Interpolate the data
        f = interp1d(xs, ys, kind='cubic')
        xs_int = xs
        ys_int = f(xs_int)
        
        # Trajectory
        fig1, ax1 = plt.subplots()
        ax1.plot(xs, ys, 'o', label='trajectory')
        ax1.plot(xs_int, ys_int, '-')
        ax1.axhline(0, color ='k')
        ax1.axvline(0, color ='k')
        ax1.legend()
        ax1.set_aspect('equal')
        

def cartesian_foucault(state):
    """
    Foucault force that is the net combination of gravity and pseudo force experience due to observer located at the location.
    """
    y1, y2, y3, y4 = state
    
    # force components in polar coordinates.
    global MASS
    global RADIUS
    global LATITUDE
    global OMEGA
    global GRAVITY
    global LENGTH
    
    m = MASS
    R = RADIUS
    th = LATITUDE
    w = OMEGA
    g = GRAVITY
    l = LENGTH
    
    s, c = np.sin(th), np.cos(th)
    r_eff = R * s
    r_sq = y1**2 + y2**2
    
    
    r_a_r = -1*r_sq/l*(g - r_eff*w**2*s - 2*w*y4*s) + r_eff*w**2*c*y1 + 2*w*(-1*y3*y2 + y4*y1)*c
    r_a_th = -1*r_eff*w**2*c*y2 - 2*w*(y3*y1 + y4*y2)*c + w**2*s**2*y1*y2
    
    a_x = (y1*r_a_r - y2*r_a_th)/r_sq
    a_y = (y2*r_a_r + y1*r_a_th)/r_sq
    
    f_x = m*a_x
    f_y = m*a_y
    
    force = (f_x, f_y)
    
    return force
    

LATITUDE = 0
state0 = (10, 0, 0, 0.01)
t_span = (0, 6*3600)
MIN_TIMESTAMPS = t_span[1]*10
cartesian_particle = CartesianParticle(cartesian_foucault, state0, t_span)

cartesian_particle.trajectory()