Earth's Gravitational Binding Energy
Dave Typinski, January 2009

Earth's Gravitational Binding Energy

(* Universal gravitational constant, NIST-CODATA, (m^3)·(kg^(-1))·(s^(-2)) *)
G=6.67428×(10^(-11));

(* Speed of light, NIST-CODATA, m·(s^(-1))*)
clight=299792458;

(* Earth mass, NASA, kg *)
(Subscript[Mass, accepted])=5.9736×(10^24) ;

(* Average Earth density, kg·(m^(-3)) *)
(Subscript[ρ, avg])=5513.94 ;

(* Earth's equatorial radius in meters- WGS 84 reference ellipsoid via Wikipedia *)
(Subscript[r, equatorial])=6.378137×(10^6); 

(* inverse flattening (1/f) - WGS 84 reference ellipsoid via Wikipedia *)
f=(1/298.257223563);

(Subscript[r, polar])=(Subscript[r, equatorial]) (1-f);

(Subscript[r, volumetric])=((((Subscript[r, equatorial])^2) (Subscript[r, polar]))^(1/3));

(* Quadratic fit to Preliminary Reference Earth Model (PREM), Dziewonski & Anderson, 1981 *)
(* Note that the PREM uses (Subscript[r, volumetric]) as the radius of the Earth *)

(* Quadratic density function constants *)
a={-2.177275×(10^(-10)),
-2.412312×(10^(-10)),
0,
-3.092201×(10^(-11)),
0,
0,
0,
0,
0,
0,
0};

b={1.911026×(10^(-8)),
1.397596×(10^(-4)),
-5.000714×(10^(-4)),
-2.444117×(10^(-4)),
-2.328571×(10^(-4)),
-1.260300×(10^(-3)),
-5.970556×(10^(-4)),
1.086873×(10^(-4)),
0,
0,
0};

c={1.308848×(10^4),
1.234623×(10^4),
7.306695×(10^3),
6.782253×(10^3),
5.319659×(10^3),
1.124903×(10^4),
7.108267×(10^3),
2.690961×(10^3),
2.900×(10^3),
2.600×(10^3),
1.020×(10^3)};

h={1.2215×(10^6),
3.4800×(10^6),
3.6300×(10^6),
5.7010×(10^6),
5.7710×(10^6),
5.9710×(10^6),
6.1510×(10^6),
6.3466×(10^6),
6.3560×(10^6),
6.3680×(10^6),
(Subscript[r, volumetric])};
h[[0]]=0;

(* Gravitational potential energy at r using average density *)
(Subscript[U, gavg])[R_]:=(16/15)G (π^2) ((Subscript[ρ, avg])^2) (R^5);
(Subscript[E, bind-avg])[R_]:=-(Subscript[U, gavg])[R];

(* Gravitational potential energy at r using manual integration *)
Ui=Table[0,{i,1,11}];
Do[
Ui[[i]]=-(π^2) G ((16/27) a[[i]]^2 (h[[i-1]]^9-h[[i]]^9)+(4/3) a[[i]] b[[i]] (h[[i-1]]^8-h[[i]]^8)+(16/21) (2 a[[i]] c[[i]]+b[[i]]^2) (h[[i-1]]^7-h[[i]]^7)+(16/9) b[[i]] c[[i]] (h[[i-1]]^6-h[[i]]^6)+(16/15) c[[i]]^2 (h[[i-1]]^5-h[[i]]^5));
,{i,1,11}
];

(Subscript[U, g])[R_]:=-(π^2) G ((16/27) a[[1]]^2 (h[[0]]^9-R^9)+(4/3) a[[1]] b[[1]] (h[[0]]^8-R^8)+(16/21) (2 a[[1]] c[[1]]+b[[1]]^2) (h[[0]]^7-R^7)+(16/9) b[[1]] c[[1]] (h[[0]]^6-R^6)+(16/15) c[[1]]^2 (h[[0]]^5-R^5))/;h[[0]]≤R≤ h[[1]];
(Subscript[U, g])[R_]:=Ui[[1]]-(π^2) G ((16/27) a[[2]]^2 (h[[1]]^9-R^9)+(4/3) a[[2]] b[[2]] (h[[1]]^8-R^8)+(16/21) (2 a[[2]] c[[2]]+b[[2]]^2) (h[[1]]^7-R^7)+(16/9) b[[2]] c[[2]] (h[[1]]^6-R^6)+(16/15) c[[2]]^2 (h[[1]]^5-R^5))/;h[[1]]≤R≤ h[[2]];
(Subscript[U, g])[R_]:=(∑_(i=1)^2)Ui[[i]]-(π^2) G ((16/27) a[[3]]^2 (h[[2]]^9-R^9)+(4/3) a[[3]] b[[3]] (h[[2]]^8-R^8)+(16/21) (2 a[[3]] c[[3]]+b[[3]]^2) (h[[2]]^7-R^7)+(16/9) b[[3]] c[[3]] (h[[2]]^6-R^6)+(16/15) c[[3]]^2 (h[[2]]^5-R^5))/;h[[2]]<R≤ h[[3]];
(Subscript[U, g])[R_]:=(∑_(i=1)^3)Ui[[i]]-(π^2) G ((16/27) a[[4]]^2 (h[[3]]^9-R^9)+(4/3) a[[4]] b[[4]] (h[[3]]^8-R^8)+(16/21) (2 a[[4]] c[[4]]+b[[4]]^2) (h[[3]]^7-R^7)+(16/9) b[[4]] c[[4]] (h[[3]]^6-R^6)+(16/15) c[[4]]^2 (h[[3]]^5-R^5))/;h[[3]]<R≤ h[[4]];
(Subscript[U, g])[R_]:=(∑_(i=1)^4)Ui[[i]]-(π^2) G ((16/27) a[[5]]^2 (h[[4]]^9-R^9)+(4/3) a[[5]] b[[5]] (h[[4]]^8-R^8)+(16/21) (2 a[[5]] c[[5]]+b[[5]]^2) (h[[4]]^7-R^7)+(16/9) b[[5]] c[[5]] (h[[4]]^6-R^6)+(16/15) c[[5]]^2 (h[[4]]^5-R^5))/;h[[4]]<R≤ h[[5]];
(Subscript[U, g])[R_]:=(∑_(i=1)^5)Ui[[i]]-(π^2) G ((16/27) a[[6]]^2 (h[[5]]^9-R^9)+(4/3) a[[6]] b[[6]] (h[[5]]^8-R^8)+(16/21) (2 a[[6]] c[[6]]+b[[6]]^2) (h[[5]]^7-R^7)+(16/9) b[[6]] c[[6]] (h[[5]]^6-R^6)+(16/15) c[[6]]^2 (h[[5]]^5-R^5))/;h[[5]]<R≤ h[[6]];
(Subscript[U, g])[R_]:=(∑_(i=1)^6)Ui[[i]]-(π^2) G ((16/27) a[[7]]^2 (h[[6]]^9-R^9)+(4/3) a[[7]] b[[7]] (h[[6]]^8-R^8)+(16/21) (2 a[[7]] c[[7]]+b[[7]]^2) (h[[6]]^7-R^7)+(16/9) b[[7]] c[[7]] (h[[6]]^6-R^6)+(16/15) c[[7]]^2 (h[[6]]^5-R^5))/;h[[6]]<R≤ h[[7]];
(Subscript[U, g])[R_]:=(∑_(i=1)^7)Ui[[i]]-(π^2) G ((16/27) a[[8]]^2 (h[[7]]^9-R^9)+(4/3) a[[8]] b[[8]] (h[[7]]^8-R^8)+(16/21) (2 a[[8]] c[[8]]+b[[8]]^2) (h[[7]]^7-R^7)+(16/9) b[[8]] c[[8]] (h[[7]]^6-R^6)+(16/15) c[[8]]^2 (h[[7]]^5-R^5))/;h[[7]]<R≤ h[[8]];
(Subscript[U, g])[R_]:=(∑_(i=1)^8)Ui[[i]]-(π^2) G ((16/27) a[[9]]^2 (h[[8]]^9-R^9)+(4/3) a[[9]] b[[9]] (h[[8]]^8-R^8)+(16/21) (2 a[[9]] c[[9]]+b[[9]]^2) (h[[8]]^7-R^7)+(16/9) b[[9]] c[[9]] (h[[8]]^6-R^6)+(16/15) c[[9]]^2 (h[[8]]^5-R^5))/;h[[8]]<R≤ h[[9]];
(Subscript[U, g])[R_]:=(∑_(i=1)^9)Ui[[i]]-(π^2) G ((16/27) a[[10]]^2 (h[[9]]^9-R^9)+(4/3) a[[10]] b[[10]] (h[[9]]^8-R^8)+(16/21) (2 a[[10]] c[[10]]+b[[10]]^2) (h[[9]]^7-R^7)+(16/9) b[[10]] c[[10]] (h[[9]]^6-R^6)+(16/15) c[[10]]^2 (h[[9]]^5-R^5))/;h[[9]]<R≤ h[[10]];
(Subscript[U, g])[R_]:=(∑_(i=1)^10)Ui[[i]]-(π^2) G ((16/27) a[[11]]^2 (h[[10]]^9-R^9)+(4/3) a[[11]] b[[11]] (h[[10]]^8-R^8)+(16/21) (2 a[[11]] c[[11]]+b[[11]]^2) (h[[10]]^7-R^7)+(16/9) b[[11]] c[[11]] (h[[10]]^6-R^6)+(16/15) c[[11]]^2 (h[[10]]^5-R^5))/;h[[10]]<R≤ h[[11]];

(Subscript[E, bind])[R_]:=-(Subscript[U, g])[R];

E1=(Subscript[E, bind-avg])[(Subscript[r, volumetric])];
E2=(Subscript[E, bind])[(Subscript[r, volumetric])];

M1=((-E1)/(clight^2));
M2=((-E2)/(clight^2));
F1=(M1/(Subscript[Mass, accepted]));
F2=(M2/(Subscript[Mass, accepted]));

Print["AVERAGE Density Method:"]
Print["The Earth's gravitational binding energy is ",E1," J."]
Print["The Earth's gravitational binding energy mass equivalent is ",M1," kg."]
Print["The Earth's gravitational binding energy mass fraction is ",F1]
Print[""]
Print["VARIABLE Density Method:"]
Print["The Earth's gravitational binding energy is ",E2," J."]
Print["The Earth's gravitational binding energy mass equivalent is ",M2," kg."]
Print["The Earth's gravitational binding energy mass fraction is ",F2]
Print[""]
Print["DIFFERENCE Analysis:"]
Print["The relative error in the average density calculations is ",((E1-E2)/E2) 100,"%."]
Print[""]


Plot[{-(Subscript[E, bind])[(10^3) x] (10^(-32)),-(Subscript[E, bind-avg])[(10^3) x] (10^(-32))},{x,0,(10^(-3)) (Subscript[r, volumetric])},PlotRange->{0,2.4},PlotStyle->{{RGBColor[0,0,.7],Thin},{Gray,Dashed,Thin}},Axes->None,GridLines->None,Frame->True,RotateLabel-> True,FrameStyle->{{Directive[Thin,FontFamily->Helvetica,14],Thin},{Directive[Thin,FontFamily->Helvetica,14],Thin}},FrameLabel->{Style["Radius, km",FontSize->14,FontFamily->Helvetica],Style["Binding Energy, J×(10^32)",FontSize->14,FontFamily->Helvetica],Style["-Gravitational Binding Energy of the Earth as a Function of Distance From Center",Bold,FontSize->16,FontFamily->Helvetica],""}]

The Moon's Gravitational Binding Energy

(* Universal gravitational constant, NIST-CODATA, (m^3)·(kg^(-1))·(s^(-2)) *)
G=6.67428×(10^(-11));

(* Speed of light, NIST-CODATA, m·(s^(-1))*)
clight=299792458;

(* Lunar mass, NASA, kg *)
(Subscript[Mass, accepted])=7.349×(10^22) ;

(* Lunar radii in meters *)
(Subscript[r, volumetric])=1.7371×(10^6);
(Subscript[r, core])=3.5×(10^5);

(* Average lunar density, kg·(m^(-3)) *)
(Subscript[ρ, avg])=((Subscript[Mass, accepted])/((4/3) π ((Subscript[r, volumetric])^3)));
 
(* Average lunar core density, kg·(m^(-3)) *)
(Subscript[ρ, core])=7200.;

(* Find average density of lunar mantle, kg·(m^(-3)) *)
(Subscript[ρ, mantle])=(((Subscript[Mass, accepted])-(Subscript[ρ, core]) (4/3) π ((Subscript[r, core])^3))/((4/3) π (((Subscript[r, volumetric])^3)-((Subscript[r, core])^3))));

Print["The lunar core average density is ",(Subscript[ρ, core])," kg·(m^(-3))"]
Print["The lunar mantle average density is ",(Subscript[ρ, mantle])," kg·(m^(-3))"]
Print["The lunar average density is ",(Subscript[ρ, avg])," kg·(m^(-3))"]
Print[""]


(* Quadratic density function constants *)
a={0,
0};

b={0,
0};

c={(Subscript[ρ, core]),
(Subscript[ρ, mantle])};

h={(Subscript[r, core]),
(Subscript[r, volumetric])};
h[[0]]=0;

(* Gravitational potential energy at r using average density *)
(Subscript[U, gavg])[R_]:=-(16/15)G (π^2) ((Subscript[ρ, avg])^2) (R^5);
(Subscript[E, bindavg])[R_]:=(Subscript[U, gavg])[R];

(* Gravitational potential energy at r using manual integration *)
Ui=Table[0,{i,1,2}];
Do[
Ui[[i]]=-(π^2) G ((16/27) a[[i]]^2 (h[[i-1]]^9-h[[i]]^9)+(4/3) a[[i]] b[[i]] (h[[i-1]]^8-h[[i]]^8)+(16/21) (2 a[[i]] c[[i]]+b[[i]]^2) (h[[i-1]]^7-h[[i]]^7)+(16/9) b[[i]] c[[i]] (h[[i-1]]^6-h[[i]]^6)+(16/15) c[[i]]^2 (h[[i-1]]^5-h[[i]]^5));
,{i,1,2}
];

(Subscript[U, g])[R_]:=-(π^2) G ((16/27) a[[1]]^2 (h[[0]]^9-R^9)+(4/3) a[[1]] b[[1]] (h[[0]]^8-R^8)+(16/21) (2 a[[1]] c[[1]]+b[[1]]^2) (h[[0]]^7-R^7)+(16/9) b[[1]] c[[1]] (h[[0]]^6-R^6)+(16/15) c[[1]]^2 (h[[0]]^5-R^5))/;h[[0]]≤R≤ h[[1]];
(Subscript[U, g])[R_]:=Ui[[1]]-(π^2) G ((16/27) a[[2]]^2 (h[[1]]^9-R^9)+(4/3) a[[2]] b[[2]] (h[[1]]^8-R^8)+(16/21) (2 a[[2]] c[[2]]+b[[2]]^2) (h[[1]]^7-R^7)+(16/9) b[[2]] c[[2]] (h[[1]]^6-R^6)+(16/15) c[[2]]^2 (h[[1]]^5-R^5))/;h[[1]]≤R≤ h[[2]];

(Subscript[E, bind])[R_]:=-(Subscript[U, g])[R];

E1=(Subscript[E, bindavg])[(Subscript[r, volumetric])];
E2=(Subscript[E, bind])[(Subscript[r, volumetric])];

M1=((-E1)/(clight^2));
M2=((-E2)/(clight^2));
F1=(M1/(Subscript[Mass, accepted]));
F2=(M2/(Subscript[Mass, accepted]));

Print["AVERAGE Density Method:"]
Print["The Moon's gravitational binding energy is ",E1," J."]
Print["The Moon's gravitational binding energy mass equivalent is ",M1," kg."]
Print["The Moon's gravitational binding energy mass fraction is ",F1]
Print[""]
Print["DUAL Density Method:"]
Print["The Moon's gravitational binding energy is ",E2," J."]
Print["The Moon's gravitational binding energy mass equivalent is ",M2," kg."]
Print["The Moon's gravitational binding energy mass fraction is ",F2]
Print[""]
Print["DIFFERENCE Analysis:"]
Print["The relative error in the average density calculations is ",((E1-E2)/E2) 100,"%."]
Print[""]


Plot[{-(Subscript[E, bind])[(10^3) x] (10^(-29)),-(Subscript[E, bindavg])[(10^3) x] (10^(-29))},{x,0,(10^(-3)) (Subscript[r, volumetric])},PlotRange->{0,1.4},PlotStyle->{{RGBColor[0,0,.7],Thin},{Gray,Dashed,Thin}},Axes->None,GridLines->None,Frame->True,RotateLabel-> True,FrameStyle->{{Directive[Thin,FontFamily->Helvetica,14],Thin},{Directive[Thin,FontFamily->Helvetica,14],Thin}},FrameLabel->{Style["Radius, km",FontSize->14,FontFamily->Helvetica],Style["Binding Energy, J×(10^29)",FontSize->14,FontFamily->Helvetica],Style["-Gravitational Binding Energy of the Moon as a Function of Distance From Center",Bold,FontSize->16,FontFamily->Helvetica],""}]

A Neutron Star's Gravitational Binding Energy, Average Density Method

(* Universal gravitational constant, NIST-CODATA, (m^3)·(kg^(-1))·(s^(-2)) *)
G=6.67428×(10^(-11));

(* Speed of light, NIST-CODATA, m·(s^(-1))*)
clight=299792458;

(* Solar mass, NASA, kg *)
(Subscript[Mass, solar])=1.9891×(10^30) ;
(Subscript[Mass, accepted])=1.5 (Subscript[Mass, solar]);

(* Neutron star radius in meters *)
(Subscript[r, volumetric])=6×(10^3);

(* Average lunar density, kg·(m^(-3)) *)
(Subscript[ρ, avg])=((Subscript[Mass, accepted])/((4/3) π ((Subscript[r, volumetric])^3)));
 
Print["The neutron star average density is ",(Subscript[ρ, avg])," kg·(m^(-3))"]
Print[""]

(* Gravitational potential energy at r using average density *)
(Subscript[U, gavg])[R_]:=-(16/15)G (π^2) ((Subscript[ρ, avg])^2) (R^5);
(Subscript[E, bindavg])[R_]:=(Subscript[U, gavg])[R];

E1=(Subscript[E, bindavg])[(Subscript[r, volumetric])];

M1=((-E1)/(clight^2));

F1=(M1/(Subscript[Mass, accepted]));

Print["AVERAGE Density Method:"]
Print["The Moon's gravitational binding energy is ",E1," J."]
Print["The Moon's gravitational binding energy mass equivalent is ",M1," kg."]
Print["The Moon's gravitational binding energy mass fraction is ",F1]
Print[""]

Plot[{-(Subscript[E, bind])[(10^3) x] (10^(-46)),-(Subscript[E, bindavg])[(10^3) x] (10^(-46))},{x,0,(10^(-3)) (Subscript[r, volumetric])},PlotRange->{0,6.1},PlotStyle->{{RGBColor[0,0,.7],Thin},{Gray,Dashed,Thin}},Axes->None,GridLines->None,Frame->True,RotateLabel-> True,FrameStyle->{{Directive[Thin,FontFamily->Helvetica,14],Thin},{Directive[Thin,FontFamily->Helvetica,14],Thin}},FrameLabel->{Style["Radius, km",FontSize->14,FontFamily->Helvetica],Style["Binding Energy, J×(10^46)",FontSize->14,FontFamily->Helvetica],Style["-Gravitational Binding Energy of a 1.5 Solar Mass, 6 km radius Neutron Star as a Function of Distance From Center",Bold,FontSize->16,FontFamily->Helvetica],""}]