Home

On the Spinning Sun-Engendered Mercury Orbital Precession of .43 Arc Seconds per Year

 References and suggested background reading:

 (1) Roots of Gravitomagnetic Theory. (Article on this website.)

(2) THE FEYNMAN LECTURES ON PHYSICS, volume 2, Section 14-5.

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 Sun spinning and (according to GMT) 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 GMT 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 also 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 an opposing gravitomagnetic force.

Analytically the orbit is a perfect ellipse when only gravity acts (Case 1). However, in Case 2, with the small (and opposing) 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 reality 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 effects.

Feynman(Reference 2) models a square electric current loop as a square line charge q, with sides a, and line charge density
l=q/4a. He assumes that this line charge circulates CCW in the xy-plane with a current, I=q/t, where t is the time for q to make one trip around the loop. I.e. I=qw/(2p). The magnetic moment for this model is m=Ia2, or m=qwa2/2p. For a point r>a in the loop's plane he finds that B points in the negative z direction and has the value Bz=-m/(4peoc2r3) = -qwa2/2p(4peo)c2r3.

We model the spinning Sun as (a) a square line mass M, where M is the mass of the Sun, and (b) with side a equal to R, the radius of the Sun, and (c) with a line mass density of
l=M/4R). We assume that this line mass circulates CCW with a mass current of I=lv, where v=2wR/p is the line mass speed around the loop.

At points in the equatorial plane, and at distances r>R, our model indicates that Oz=-
mG/c2r3 = -MwR2G/2pc2r3. Given the observed complexities of the real Sun, this model-based value for Oz does not produce an orbit whose precession matches the observed one of .43 arc seconds per year. Since the outer regions of the Sun are less dense than the more central ones, it might be guessed that the value of Oz in the real case must be somewhat less than what our model indicates. Indeed one can experiment with different adjustment values for Oz to produce a precession that approximates the observed .43 arc seconds. It turns out that an adjustment factor of .959 results in an Oz that closely produces the required .43 arc seconds per year.

Appendix. Python Program.

#Python 3.3.4 (v3.3.4:7ff62415e426, Feb 10 2014, 18:12:08) [MSC v.1600 32 bit (Intel)] on win32

#Type "copyright", "credits" or "license()" for more information.

#class Precession:

import math

#Constants for Sun.

MassOfSun=1.9891e30

RadiusOfSun=6.955e8

OmegaOfSun=2.*math.pi/(25.71*24*60**2) #Based on the model.

#Constants for Mercury.

MercuryMass=.33011e24

PerihelionDistance=46001200e3

PerihelionSpeed=58.98e3

OrbitalPeriod=87.969*86400 #seconds

#Program constants

N=2**20

DT=OrbitalPeriod/(N-1) #Time increment between orbital step computations.

SpeedOfLight=2.99792458e8

G=6.67408e-11 #Gravitational constant.

#First run the orbit computing code for a non-spinning Sun, and determine the computer error.

#Do first loop.

print("computing. Please wait.")

Dis=[]

DisX=[]

DisY=[]

DisX.append(PerihelionDistance)

DisY.append(0.)

Dis.append(math.sqrt(DisX[0]**2+DisY[0]**2))

g=G*MassOfSun/Dis[0]**2   #g is imaginary since MassOfSun is imaginary.

gX=g

gY=0.

FX=-gX*MercuryMass   #Gravitational force. Both gX and MercuryMass are imaginary.

FY=0.

AX=FX/MercuryMass

AY=FY/MercuryMass

VX=AX*DT/2

VY=PerihelionSpeed+(AY*DT/2)

#Set up for t=DT step.

DisX.append(DisX[0]+VX*DT)

DisY.append(DisY[0]+VY*DT)

Dis.append(math.sqrt(DisX[1]**2+DisY[1]**2))

g=G*MassOfSun/Dis[1]**2

gX=g*(-DisX[1]/Dis[1])

gY=g*(-DisY[1]/Dis[1])

FX=gX*MercuryMass   

FY=gY*MercuryMass

AX=FX/MercuryMass

AY=FY/MercuryMass

VX=VX+AX*DT

VY=VY+AY*DT

#Now do the rest of the orbits.

for index in range(2,N+1):

    DisX.append(DisX[index-1]+VX*DT)

    DisY.append(DisY[index-1]+VY*DT)

    Dis.append(math.sqrt(DisX[index]**2+DisY[index]**2))

    g=G*MassOfSun/Dis[index]**2

    gX=g*(-DisX[index]/Dis[index])

    gY=g*(-DisY[index]/Dis[index])

    FX=gX*MercuryMass

    FY=gY*MercuryMass

    AX=FX/MercuryMass

     AY=FY/MercuryMass

    VX=VX+AX*DT

    VY=VY+AY*DT

#The following instructions determine computer numeric error.

#The idea is that the orbit around a non-spinning Sun is a closed ellipse with no precession.

NumericErrorX=DisX[N]-PerihelionDistance

NumericErrorY=DisY[N]

NumericError=math.sqrt(NumericErrorX**2+NumericErrorY**2)

print("NumericError=",NumericError)   #For reference.

#Now compute the orbit around a spinning Sun.

#Given the vagaries of the Sun, a small adjustment is required to get a precession of

#approximately .43 arc secs.

OzAdjust=.959  

print("Computing. Please wait.")

Dis=[]

DisX=[]

DisY=[]

DisX.append(PerihelionDistance)

DisY.append(0.)

Dis.append(math.sqrt(DisX[0]**2+DisY[0]**2))

g=G*MassOfSun/Dis[0]**2   #g is imaginary since MassOfSun is imaginary.

gX=g

gY=0.

FX=-gX*MercuryMass   #Gravitational force. Both gX and MercuryMass are imaginary.

FY=0.

Oz=OzAdjust*MassOfSun*RadiusOfSun**2*G*OmegaOfSun/(2*math.pi*SpeedOfLight**2*Dis[0]**3)

#Oz is negative imaginary in the Sun's outer equatorial plane.

GMFX=-MercuryMass*PerihelionSpeed*Oz #Mercury Mass and Oz are both imaginary.

#GMFX is real positive.

GMFY=MercuryMass*0*Oz

TotalFX=FX+GMFX #Points toward Sun, but less than force when Sun does not spin.

TotalFY=FY+GMFY

AX=TotalFX/MercuryMass   #Newton 2. Mercury Mass is inertial and hence real.

AY=TotalFY/MercuryMass

VX=AX*DT/2

VY=PerihelionSpeed+(AY*DT/2)

#Set up for t=DT step.

DisX.append(DisX[0]+VX*DT)

DisY.append(DisY[0]+VY*DT)

Dis.append(math.sqrt(DisX[1]**2+DisY[1]**2))

g=G*MassOfSun/Dis[1]**2

gX=g*(-DisX[1]/Dis[1])

gY=g*(-DisY[1]/Dis[1])

FX=gX*MercuryMass #In these equations MercuryMass is inertial, and hence real.

FY=gY*MercuryMass

Oz=-OzAdjust*MassOfSun*RadiusOfSun**2*G*OmegaOfSun/(2*math.pi*SpeedOfLight**2*Dis[1]**3)

GMFX=-MercuryMass*VY*Oz

GMFY=MercuryMass*VX*Oz

TotalFX=FX+GMFX #Points toward Sun.

TotalFY=FY+GMFY

AX=TotalFX/MercuryMass

AY=TotalFY/MercuryMass

VX=VX+AX*DT

VY=VY+AY*DT

#Now do rest of orbits.

for index in range(2,N+1):

    DisX.append(DisX[index-1]+VX*DT)

    DisY.append(DisY[index-1]+VY*DT)

    Dis.append(math.sqrt(DisX[index]**2+DisY[index]**2))

    g=G*MassOfSun/Dis[index]**2

    gX=g*(-DisX[index]/Dis[index])

    gY=g*(-DisY[index]/Dis[index])

    FX=gX*MercuryMass

    FY=gY*MercuryMass

   Oz=-  OzAdjust*MassOfSun*RadiusOfSun**2*G*OmegaOfSun/(2*math.pi*SpeedOfLight**2*Dis[index]**3)

    GMFX=-MercuryMass*VY*Oz

    GMFY=MercuryMass*VX*Oz

    TotalFX=FX+GMFX #Points toward Sun, but less than force when Sun does not spin.

    TotalFY=FY+GMFY

    AX=TotalFX/MercuryMass

    AY=TotalFY/MercuryMass

    VX=VX+AX*DT

    VY=VY+AY*DT

CorrectDisX=DisX[N]-NumericErrorX

CorrectDisY=DisY[N]-NumericErrorY

CorrectDis=math.sqrt(CorrectDisX**2+CorrectDisY**2)

 theta=math.acos(PerihelionDistance/CorrectDis)   #In radians.

theta=theta*(360/(2*math.pi))    #Converts radians to degrees.

theta=theta*60**2   #Converts degrees to arc seconds.

theta=4.15*theta   #There are 4.15 orbits of Mercury for each Earth year.

print("Precession",theta)   #Precession per Earth year.

print("End Of Program")