Libration

Presentation and Results

The Moon exhibits an oscillatory motion around its principal axis, with an amplitude of about 7 degrees, called longitudinal libration.
In the traditional SM (Spinning Moon) conception, this phenomenon is explained kinematically, through the combination of the eccentricity of the lunar orbit and its supposed rotation on its own axis.

In the NSM (Non‑Spinning Moon) conception, where this rotation does not exist, libration can be explained as the conjunction of two dynamical effects:

the well‑known tidal effect ¹, presented for example in https://fr.wikipedia.org/wiki/Rotation_synchrone
a previously unstudied effect, which we call the tumble effect², caused by a slight mass imbalance in favor of the hemisphere facing the Earth.

¹ In NSM the tidal effect is real, although synchronous rotation is not.
² The tidal and tumbling effects both result in a pendulum-like movement.

The program whose source code is provided below can be run in two modes: a main mode called “libration”, and an auxiliary mode called “rotation”.

Thanks to the Matplotlib library, the image of the Moon is added to the plot at regular angular intervals, and the final image at the end of the orbit is produced at the end of the program.

In “libration” mode

The program simulates the pendulum effect in an idealized case where the tidal effect is assumed to be zero and where the mass imbalance ubl is concentrated at a single point, the HDS (High Density Spot), located at point N (Nearest, origin of selenographic coordinates).
The initial pendulum angle is set to 7 degrees, and ubl is then varied at each run until the best match with reality is obtained, namely one full pendulum swing per orbit.
This result is reached for ubl = 5 per thousand, where the following figure is obtained:

In “rotation” mode

All libration effects are cancelled, which makes it possible to simulate the continuous rotation(s) of the Moon as seen in the left‑hand or right‑hand figures on the site https://science.nasa.gov/moon/tidal-locking/ (these figures appear almost identically on https://en.wikipedia.org/wiki/Tidal_locking).

left‑hand figure (“real” Moon)

right‑hand figure (“false” Moon)

The two figures make it possible to confirm dynamically the result already obtained through geometry, namely the validity of the NSM conception and the unreality of the SM conception. When the program is run with zero angular velocity omega, it generates the “real” Moon on the left, and when it is run with a specific value of omega, it generates the “false” Moon on the right.

Code

			
# coding = utf - 8
"""
                           ********************************************************
                           * SIMULATION OF THE LONGITUDINAL LIBRATION OF THE MOON *
						   *				   February 28, 2022 				  *	
                           ********************************************************
                           						Gilbert Vidal
                            Minor updates: on March 20, 2023 and Decembre 10, 2024
                           
No entries.
Change any parameter directly in its declaration.
If you do so, don't forget to register before running.
A vector begins with an upper case letter.
The unit vector of an axis begins with the upper case U.
The vector from B to C is labelled BC. But, as E is the origin, the vector EC is labelled C.
A scalar begins with a lower case letter.
The coordinates of the point P are those of the vector EP (because E is the origin of the plan).
The distance between the points B and C is labelled lBC (l stands for length).
The angles are expressed  1/ in radians for calculations
                          2/ in degrees for inputs/outputs (with suffix _d).
"""
import numpy as np
from math import cos, sin, radians, sqrt, pi
import matplotlib.pyplot as plt
import matplotlib.lines as lines
import matplotlib.path as path
import matplotlib.patches as patches
#1 - choose mode libration or rotation
mode = 0											#0: rotation mode		1: libration mode
#2 - choose figure left or right (no incidence if mode == 1)
rotation = "left"									#"left": related to the left animated figure
													#"right": related to the right animated figure
omega = 0 											#angular velocity of the Moon around
													#its centre of mass (rd/s)
if rotation == "right" and mode == 0:
	omega = -2*pi/(28*24*3600)
#VECTORIAL FUNCTIONS *********************************************************************************
def sum(V0,V1):                                     #sum of vectors V0 and V1
    return V0[0] + V1[0], V0[1] + V1[1]
def dif(V0,V1):                                     #sum of vectors V0 and V1
    return V0[0] - V1[0], V0[1] - V1[1]
def mod(V):                                         #module of the vector V
    return sqrt(V[1]**2 + V[0]**2)                         
def kV(k, V0):                                      #product of vector V0 by a scalar k
    return k*V0[0], k*V0[1]
def Vecprod(V0, V1):                                #scalar of the vectorial product of vectors V0 and V1
    return V0[0]*V1[1] - V0[1]*V1[0]
def unit(V0, V1):                                   #unitary Vector of an axis bearing (parallel to)
    V = dif(V1,V0)                                     #the vector P0P1 (= (V1 - V0))
    return kV(1/mod(V), V) 
#UNIVERSAL CONSTANTS  *********************************************************************
r = 1.737e6                         #radius of the Moon (m)               
mass = 7.34e22                      #mass of the Moon (kg)
mass_earth = 5.97e24                #mass of the Earth (kg)
g = 6.673e-11                       #gravitational constant  N·(m/kg)
v_zero = 970                        #velocity of the Moon at start (m/s)
lEC_zero = 4.06e8                   #Distance (length) at start between the mass centre C
                                    #    of the Moon and the Earth M (at apogee).
                                    #As regards the distance between Earth and Moon,
                                    #    EC and EM are supposed equal.
                                    
                                    #https://en.wikipedia.org/wiki/Lunar_distance_(astronomy)
                                    #from 356,500 km at the perigee to 406,700 km at apogee, 
                                    #https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
                                        #Perigee            3.63e8
                                        #Apogee              4.06e8  **************
                                        #Max. orbital velocity     1082 
                                        #Min. orbital velocity (km/s)    0.970  **********
                                        #revolution period (days)    27.3217
                                
#CONSTANTS *********************************************************************************
ubl = + .005                        #relative unbalance due to the HDS (high density spot)
if mode == 0: ubl = 0
mass_s = ubl*mass                   #mass of the HDS
da = .01                            #angular orbital increment of a sample (rd)
b_zero_d = 7						#pendulum angle of the Moon at start IN DEGREES
if mode == 0: b_zero_d = 0
b_zero = b_zero_d*pi/180            #pendulum angle of the Moon at start (in rd of course)
                  
m = mass_s*r/(mass)                 #abciss on the axis bd of the geographic centre of the moon 
                                    #bd is the relative axis joining the bright and dark poles.
                                    #bd is oriented from bright pole to dark pole.
                                    #Its origin is C, center of mass of the Moon.    
s = - r                             #abciss on the bd axis of the HDS
MoI = mass*r**2*(1 + ubl)           #moment of inertia (relative to C) of the Moon
#INITIALIZATION OF THE VARIABLES   ********************************************************
a = 0                               #orbital angle (rd)
b = b_zero                          #pendulum angle (rd)
if mode == 0: b = 0
db = 0							    #angular increment (for display purposes)
C = [lEC_zero, 0]                   #vector EC (coordinates of the centre of mass of the Moon (m))
V = [0, v_zero]                     #coordinates of the vector velocity of the centre of mass (m/s)
t = 0                               #time
dt = da*lEC_zero/v_zero             #duration of the orbital increment of a sample.
                                    #dt is not constant. This is only its initialization value.     
                                    #dt is NOT the increment of the main loop.
                                    #dt is calculated at each loop. The true increment is angle da.
nb_samp_total = int(2*pi/da)        #number of samples of the MAIN LOOP
nb_plots = 12                       #number of plots of the MAIN LOOP
nb_samp_between_plots = int(nb_samp_total/nb_plots)          #number of samples between 2 plots
num = 0								#number of the current plot (for display purposes)	
#FOR VISUALIZATION NEEDS ******************************************************************
kb = 5                              #ampli factor applied on the pendulum angle for visu needs
if mode == 0: kb = 1		
kr = 30                             #ampli factor applied on the radius for visu needs
r_v = kr*r							#radius of the visualized Moon
rspot_v	= r_v/5						#radius of the visualized HDS (High Density Spot)
									#(the real HDS has no radius. It is infinitely small)
s_v = -r_v + rspot_v				#abcissis on the bd axix of the centre of the visualized HDS
#GRAPHICS  (except the plotting of the Moon itself)*****************************************
fig = plt.figure(figsize=(7, 6))	#total size: 7 by 6 inches
ax = fig.add_subplot(1,1,1)
    
xy = 450e6
xmin =  - xy  ;  xmax = xy  ;  ymin =  - xy  ;  ymax = xy
s1="{:.1e}".format(omega)
titre1 = 'ubl=' + str(ubl) + '     omega=' + s1 + ' rd/s' + '     b_zero=' + str(b_zero_d) + ' deg' + '\n'
ax.set(xlim = (xmin, xmax), ylim = (ymin, ymax), title = titre1)
plt.plot([xmin/2, xmax/2], [0,0], color='yellow')  ;  plt.plot([0,0], [ymin/2, ymax/2], color='yellow')
def annot(txt, xy, xytext):
	ax.annotate(txt, xy, xytext, arrowprops=dict(facecolor='black', shrink=0.01, width=1))
def annot1(txt, xy, size):#**kwargs):
	ax.annotate(txt, xy, size=size)
if mode == 1:
	annot('E (Earth)', (0, 0), (-2e8, -1.5e8))
	annot1('The radius of the moon is enlarged ' + str(kr) + ' times', (-2.5e8, 1e8), size=10)
	annot1('The libration angle is enlarged ' + str(kb) + ' times', (-2.5e8, .5e8), size=10)
	annot('Bright-Dark axis (bd)', (3.2e8, -.6e8), (.1e8, -1e8))
	annot('high density spot (HDS)', (3.6e8, -.4e8), (.1e8, -1.65e8))
if mode == 0:
	if rotation == "left":
		annot1('omega = 0', (-1.5e8, 1e8), size=20)		
		annot1('this Moon', (-1e8, 0), size=15)
		annot1('is NOT spinning', (-1.6e8, -1e8), size=15)
	if rotation == "right":
		annot1('omega = -2.6e-6', (-2e8, 1e8), size=20)		
		annot1('this "Moon"', (-1e8, 0), size=15)
		annot1('is really spinning', (-1.5e8, -1e8), size=15)
		annot1('(clockwise)', (-1e8, -2e8), size=15)
		
#MAIN LOOP  The basic increment is da, it is not dt*****************************************
#           ********************************************************************************
 
for i_sample in range(nb_samp_total):
    
	def mapping():
	
		global lEC, lEB, Ue, Ufs, Ubd, M, Ubd_v, B, CB, a, b
		
		if  mode == 0 and rotation == "right":
			pass
			b = -a		# This line is a deliberate cheating to display the ideal expected result.
						# Comment it out (replace 'b' with '#b') if you want to highlight
						# the real figure.
		
		Ue = (cos(a), sin(a))                     #unitary vector on the axis EM   
		Ubd = (cos(a + b), sin(a + b))            #unitary vector on the axis bd
		Ubd_v = (cos(a + kb*b), sin(a + kb*b))      #unitary vector on the axis bd "magnified" (displayed)
		E = (0, 0)							#centre of the Earth
		M = C                              	#centre of the Moon,
											#it is assumed to coincide with the center of mass
		CB = kV(s, Ubd)                     	#vector CB 
		B = sum(C, CB)                     	#vector EB: coordinates of the point B
		Ufs = unit(B, E)				   	#unitary vector of the attractive force due to the HDS.	
		lEC = mod(C)                       	#distance between the Earth and the centre of mass of the Moon
		lEB = mod(B)                       	#distance between the Earth and the HDS (or B)	
    
	   
	mapping()
 
    #TRAJECTORY OF THE CENTRE OF MASS OF THE MOON
    
	f = g*mass*mass_earth/lEC**2 
	F = kV(-f, Ue)
	V = sum(V, kV(dt/mass, F))
	C = sum(C, kV(dt, V))
    #PENDULUM MOVEMENT
	fs = g*mass_s*mass_earth/lEB**2     #module of the force exercised by the Earth on the HDS
	Fs = kV(fs, Ufs)                    #vector of the force exercised by the Earth on the HDS
	moment = Vecprod(kV(s, Ubd), Fs)     #scalar of the moment due to the HDS
	
    #OUTPUTS
	if (i_sample % nb_samp_between_plots  == 0) and (i_sample < (nb_samp_total -10)):
		
		MB_v = kV(s_v, Ubd_v)						#visualized vector MB 
		B_v = sum(M, MB_v)                      #coordinates of the centre of the visualized HDS
		moon = patches.Circle(M, r_v, alpha=1, facecolor="yellow", edgecolor="black")
		spot = patches.Circle(B_v, rspot_v, alpha=1, color="black")
		centre = patches.Circle(M, r_v/10, alpha=1, color="black")
		me1 = sum(M, kV(-2*r_v, Ue))  ;  me2 = sum(M, kV(1.5*r_v, Ue))
		bd1 = sum(M, kV(-2*r_v, Ubd_v))  ;  bd2 = sum(M, kV(1.5*r_v, Ubd_v))
		axis_me1 = (me1, me2)
		axis_bd1 = (bd1,bd2)
		axis_me = patches.Polygon(axis_me1, color="blue")
		axis_bd = patches.Polygon(axis_bd1, color="red")
		ax.add_artist(moon)
		if mode == 1: ax.add_artist(spot)
		ax.add_artist(centre)
		ax.add_artist(axis_bd)
		if mode == 1: ax.add_artist(axis_me)
		
	mapping()
    
   
    #INCREMENTATIONS
    
	v_trajectory = sqrt(mod(V)**2)
	Proj = Vecprod(V, Ufs)				#projection of vector V on the perpendicular to Ufs
	dt = da*lEC/v_trajectory
	#dt = da*lEC/Proj
	omega = omega + moment*dt/MoI
	db = omega*dt
	a = a + da
	b = b + db
	t = t + dt
	
 
#FINAL
days = t / (24*3600)
print ('\n', i_sample, '   ', 'v0=', v_zero, '    ', 'total time=', '%.2f' % days,'days')
plt.show()

		

Presentation and Results

In “libration” mode

In “rotation” mode

Code

Share this:

Leave a comment Cancel reply