Suggested background reading:
- THE FEYNMAN LECTURES ON PHYSICS, volume 2, Section 14-5.
- When Gravity Balances the Lorentz Force, www.physics-by-dixon.com
In this article we compute the orbit of Mercury as it travels around the Sun. Two cases are examined: (1) with the Sun not spinning (Newton 2 and his force law of Universal Gravitation), and (2) with the spinning Sun engendering a gravitomagnetic field in addition to the omnipresent gravitational field.
The center of the sun is assumed to be fixed at the origin of a rectangular coordinate system in inertial space. When the sun spins CCW, its angular velocity points in the positive z-direction.
In Gravitomagnetic Theory the spinning sun is theorized to engender a dipolar gravitomagnetic field that (a) points in the positive z-direction at internal points of the Sun’s equatorial plane, and (b) points in the negative z-direction at points in the external equatorial plane. The orbit of Mercury is assumed to lie in the xy plane and, from the perspective of the positive z axis, also to be CCW. In Case 1 Mercury is theorized to be acted upon solely by a gravitational force. In Case 2 Mercury is theorized to be acted upon by a gravitational force and by a much weaker gravitomagnetic force.
Analytically the orbit is a perfect ellipse when only gravity acts (Case 1). However, in Case 2, with the small gravitomagnetic force also acting, the result should be an open ellipse.
It is expected that the computed result may not be precisely what is observed owing to certain realities. First the spinning Sun is modeled as a spinning square line mass, but in actuality it is a spinning, oblate sphere of mass with a complicated internal density and with zones that rotate at different rates. Secondly, the Sun’s spin equatorial plane does not actually coincide with the plane of Mercury’s orbit. (Mercury’s orbital plane is tilted slightly, relative to the Sun’s spin equatorial plane.) Finally, there may be small numeric computer errors that cannot be ignored, since they are comparable in size to the very faint gravitomagnetic affects.
Feynman (Reference 2) models a square electric current loop as a square line charge q, with sides a, and line charge density λ=q/4a. He assumes that this line charge circulates CCW in the xy-plane with a current, I=q/τ, where τ is the time for q to make one trip around the loop. I.e. I=qω/2π. The magnetic moment for this model is μ=Ia2, or μ=qωa2/2π. For a point r>a in the loop’s plane he finds that Bz points in the negative z-direction and has the value Bz=-μ/4πεοc2r3.
We model the spinning Sun as a square line mass of side R, and mass M, where (a) M is the mass of the Sun, and (b) R is the radius of the Sun. The line mass density is λ=M/4*R. We assume that this line mass circulates CCW with a mass current of I=λv, where v=2ω*R/π is the line mass speed around the loop.
At points in the equatorial plane, and at distances r>R, our model indicates that O=-μG/c2r3 = -MωR2G/4πc2r3. Using this value for Oz, the following Python program yields a precession of .435 arcseconds per year … a close approximation to the presently measured precession of .43 arcseconds
Following is the listing for the Python program used to compute the orbits.
import math
pi=math.pi
c=2.99792458e8 #Speed of light
N=4000000
#Constants for Sun.
M=1.9891e30 #Sun Mass
R=6.955e8 #Sun Radius
Tau=2*pi/(28*24*60**2) #Sun rotation period
Omega=2*pi*Tau #Sun rotation rate
#Constants for Mercury.
m=.33011e24 #Mass of Mercury
Peri=46001200e3 #Perihelion distance
PeriSpeed=58.98e3 #Mercury speed at perihelion
tau=87.969*86400 #Mercury orbital period
#Program constants
DT=tau/N #Time increment between orbital step computations
G=6.67408e-11 #Gravitational constant.
Obase=-M*R**2*G*Omega/(2*pi*c**2) #Constant
#Lists for gravity only.
Dis1=[] #Mercury distance from Sun
Dis1x=[]
Dis1y=[]
V1=[] #Mercury speed
V1x=[]
V1y=[]
A1=[] #Mercury acceleration
A1x=[]
A1y=[]
#Gravity plus gravitomagnetism lists.
Dis2=[]
Dis2x=[]
Dis2y=[]
V2=[]
V2x=[]
V2y=[]
A2=[]
A2x=[]
A2y=[]
#Initialize gravity only lists.
Dis1.append(Peri)
Dis1x.append(-Peri)
Dis1y.append(0.)
V1.append(PeriSpeed)
V1x.append(0.)
V1y.append(-PeriSpeed)
F1=G*M*m/Dis1[0]**2
F1x=F1
F1y=0.
A1.append(F1/m)
A1x.append(F1x/m)
A1y.append(F1y/m)
#Update gravity only lists.
for index in range(1,N+2):
Dis1x.append(Dis1x[index-1]+V1x[index-1]*DT)
Dis1y.append(Dis1y[index-1]+V1y[index-1]*DT)
Dis1.append(math.sqrt(Dis1x[index]**2+Dis1y[index]**2))
V1x.append(V1x[index-1]+A1x[index-1]*DT)
V1y.append(V1y[index-1]+A1y[index-1]*DT)
V1.append(V1x[index]**2+V1y[index]**2)
F1=G*M*m/Dis1[index]**2
F1x=F1*(-Dis1x[index])/Dis1[index]
F1y=F1*(-Dis1y[index])/Dis1[index]
A1.append(F1/m)
A1x.append(F1x/m)
A1y.append(F1y/m)
#initialize gravity plus lists.
Dis2.append(Peri)
Dis2x.append(-Peri)
Dis2y.append(0.)
V2.append(PeriSpeed)
V2x.append(0.)
V2y.append(-PeriSpeed)
F2=G*M*m/Dis2[0]**2
F2x=F2*(-Dis2x[0]/Dis2[0])
F2y=F2*(-Dis2y[0]/Dis2[0])
O=Obase/(Dis2[0]**3)
Oz=O
GMFx=m*V2y[0]*Oz
GMFy=m*(-V2x[0])*Oz
F2x=F2x+GMFx
F2y=F2y+GMFy
A2x.append(F2x/m)
A2y.append(F2y/m)
#Update grav-gravitomagnetism lists.
for index in range(1,N+2):
Dis2x.append(Dis2x[index-1]+V2x[index-1]*DT) #>0>>0>>>0
Dis2y.append(Dis2y[index-1]+V2y[index-1]*DT)
Dis2.append(math.sqrt(Dis2x[index]**2+Dis2y[index]**2))
V2x.append(V2x[index-1]+A2x[index-1]*DT)
V2y.append(V2y[index-1]+A2y[index-1]*DT)
V2.append(math.sqrt(V2x[index]**2+V2y[index]**2))
F2=G*M*m/Dis2[index]**2 #Gravitational force
F2x=F2*(-Dis2x[index]/Dis2[index])
F2y=F2*(-Dis2y[index]/Dis2[index])
O=Obase/(Dis2[index]**3)
Oz=O
GMFx=m*V2y[index]*Oz #Gravitomagnetic force
GMFy=m*(-V2x[index])*Oz
GMF=math.sqrt(GMFx**2+GMFy**2)
F1x=F1x+GMFx #Total force when Sun spins
F2y=F2y+GMFy
A2x.append(F2x/m) #Mercury acceleration
A2y.append(F2y/m)
A2.append(F2/m)
#Compute precession ... angle between grav and gravplus distances.
Precess=math.acos(Dis2x[N]/Dis1x[N])
Precess=Precess*360/(2*pi)
Precess=Precess*60**2
print("Precess=",Precess," arcsec")
print("End Of Program")