Proof that Gravitational Constant is not constant

They were using this formula
PGRb = -(192 * Pi * (G ^ (5 / 3)) / (5 * (c ^ 5))) * ((Pb / (2 * Pi)) ^ -(5 / 3)) * (1 + ((73 / 24) * e ^ 2) + ((37 / 96) * e ^ 4)) * ((1 - (e ^ 2)) ^ -(7 / 2)) * (M1 * M2 * ((M1 + M2)) ^ -(1 / 3))

Which is the integral (is this the right word?) of the change in orbital period per second, based on the period, and the eccentricity of the binary.
If I get the link to this formula could someone check the formula was derived correctly?
http://iopscience.iop.org/0004-637X/722/2/1030/pdf/0004-637X_722_2_1030.pdf
TIMING MEASUREMENTS OF THE RELATIVISTIC BINARY PULSAR PSR B1913+16

Formula 3 is on the 3rd page.

 

Log in or Sign up to hide all adverts.

From that article 3.3. "Gravitational Radiation Damping - According to general relativity a binary star system should radiate energy in the form of gravitational waves. Peters &
Matthews (1963) showed that the resulting rate of change in orbital period should be ....
PGRb = -(192 * Pi * (G ^ (5 / 3)) / (5 * (c ^ 5))) * ((Pb / (2 * Pi)) ^ -(5 / 3)) * (1 + ((73 / 24) * e ^ 2) + ((37 / 96) * e ^ 4)) * ((1 - (e ^ 2)) ^ -(7 / 2)) * (M1 * M2 * ((M1 + M2)) ^ -(1 / 3))"

Now can anyone source that "Peters & Matthews (1963) paper? "Gravitational Radiation from Point Masses in a Keplerian Orbit" Please.

 

Log in or Sign up to hide all adverts.

Strangely enough the above formula gave a perfect answer when G was set on 6.66110E-11 when the result of formula 3 was used in formula 6 giving a result of 0.99999 rather than 0.99720 when the standard value of G applied 6.67384E-11. Is this another indicator that "G" in the binary is lower than in our solar system?

 

Log in or Sign up to hide all adverts.

Now going to put this lower value of G back into the orbital period formula and see what the value of "r" (the sum of the semi major axis) becomes.
1.9480635068E+09 calculated Sum of Semi Major Axis in meters using GB 1.9491384000E+09 Previous calculated Sum of Semi Major Axis in meters using GB.
That is a difference of 1.0748932000E+06 meters. Could Weisberg and Taylor be that far out in their estimates?

Using these new G and r figures back into the dr/dT formula give resultant orbital decay of 0.297 meters/year. So it has not helped to explain the difference in expected orbital period change.

Is it possible the negative potential energy is being paid back at a rate whereby the G declines sufficient per year that the orbital period shortens by 76.5 microseconds but it only appears to move 0.297 m/y but in actual fact moves more like 3.54 m/y?
So there is a difference in distance (3.26 m/y) to account for, and from the Newtonian gravitational force formula we can estimate the amount of gravitational potential energy it loses over that distance and assuming it also pays back to "G" let's see if that is a possible solution.

 

This was like taking out a loan to repay debt. One debt is partially paid off but you are getting closer and closer to the fiscal cliff. Sounds a bit like the USA but the interest has to be paid yearly and it seems like a good idea to have quantitative easing but will it work? Can the binary stars dance their way apart before the big collapse occurs?
http://en.wikipedia.org/wiki/United_States_fiscal_cliff

 

Getting excited by the figures coming out of the equations. I won't say too much at this stage, but it looks like the change in Gravitational Potential energy due to the declining Gravitational Constant of the Binary "GB" and "GB2" is coming out nearly equal to the value calculated for the Gravitational Radiation Energy produced for the same time period.
So that means Gravitational Radiation Energy could be the same amount of energy as is paid back to the negative Gravitational Potential energy (due to the declining G force), so if this is true, and I'm just trying to understand what might be happening, is that the Gravitational Radiation Energy is not lost as such but is balanced internally within the binary.

I could still be wrong on this, for the output to the macro is a bit messy, but this could be what it means.

 

Excitement back to minus 1 again. False alarm.

 

The silver lining for you is you'll probably never know how stupid this is. 24 posts to yourself in a row. They're ought to be a law.

 

There ought to be a law against counting posts. Complaining about it should be be even worse.

 

Strange thing was that going back through the macro checking on the results "GB" was smaller than "GB2". Now that has my theory rattling for I was expecting the Gravitational force to be declining and paying back the negative gravitational potential energy but it initially appears to be the other way around. Is there another mistake in my math or is this real?

This is beginning to get interesting again for that was the expected result if the gravity interaction moved at the speed of light, the binary stars are always moving into a region of higher G and would therefore be accelerating. I'll have to refresh my memory of how that works again. (For there to be an acceleration you need a force and a force creates an action and a reaction so you can't get acceleration out of thin air! Something else would need to decelerate.)

 

I Found the error and now I have a solution, that unfortunately I can't express as a formula but as concepts and this is backed up with an Excel macros. Now the Excel sheet only easily works with numbers to 15 significant digits so I have gone to the limit of the Excel program refining the Gravitational Constant (now a variable) the sum of the semi major axes, and the combined amount the stars fall in toward each other on a yearly basis, that would satisfy their orbital period change and the amount of gravitational energy produced.

The nearest I can refine the movement (orbital decay) was -5.19156873220307 m/year (another attempt -5.191568732211110E+00 m/year fully macro'ed)
, but even this had -1.957805145262480E+25 Joules of energy unaccounted for.
But this wasn't bad considering the binary was producing -2.0682747571E+31 Joules/year Gravitational radiant energy.
4.1365534299E+31 Joules/year are repaid to the negative potential energy debt to lower the Gravitational force between them.
The total negative potential energy debt between the stars was -2.7077872656E+41 joules so only a minute fraction is being repaid yearly.

I'll tidy up the macro so you can all run it and see the results for yourself soon.

 

This macro will set up a special sheet to run on so it is entirely safe. Once you have run it once change theses 4 lines if you need to see it run again or delete the special sheet.
Change from this:

Sheets.Add
ActiveSheet.Name = "Finding_Delta_r"
Range("A1").Select
' Sheets("Finding_Delta_r").Select
Change to this (by adding 3 comment blocks and removing 1 comment block):

'Sheets.Add
' ActiveSheet.Name = "Finding_Delta_r"
'Range("A1").Select
Sheets("Finding_Delta_r").Select

It scans for the limit of motion required to keep all the measured values correct.

It hasn't got all the variables declared but it will still run OK. Takes about 10 seconds run-time on my machine.


Sub Finding_Delta_r()
'
' Finding_Delta_r Macro
' Macro recorded 14/11/2012 by Robittybob1
'
Dim thrs As Double
Dim T As Double
Dim Tc As Double
Dim r As Double
Dim drInt As Double 'from dr_dt_GB Macro
Dim drGRN As Double 'from EQ3_Power Macro
Dim e As Double '(eccentricity -0.617131)
Dim K2 As Double
Dim drAvg As Double 'dr Averaged
Dim dr As Double
Dim dr_dt As Double
Dim GB As Double
Dim c As Double
Dim M1 As Double
Dim M2 As Double
Dim mu As Double '(std gravitational parameter (mu = G*(M1+M2))
Dim a As Double 'semi major axis 1,950,100,000 m
'T = 27906.97959 '(7.75*60*60)based on measured data 27906.979587552
Twt = 27906.979587552


Sheets.Add
ActiveSheet.Name = "Finding_Delta_r"
Range("A1").Select
' Sheets("Finding_Delta_r").Select


Columns("A:C").Select
Range("A13").Activate
Selection.NumberFormat = "0.000000000000000E+00"
Selection.ColumnWidth = 41.86
Range("C32").Select
Selection.Interior.ColorIndex = 4
Range("A19").Select
Selection.Interior.ColorIndex = 4
Range("B19").Select
Selection.Interior.ColorIndex = 4
Range("B1").Select
With Selection.Interior
.ColorIndex = 4
.Pattern = xlSolid
End With
Range("A13").Select
ActiveCell.FormulaR1C1 = "=1949138400"

Range("B13").Select
ActiveCell.FormulaR1C1 = _
"=""Calculated Sum of Semi Major Axis in meters Weisberg and Taylor"""
Range("A17").Select
ActiveCell.FormulaR1C1 = "=6.67213237560131E-11"
Range("B18").Select
ActiveCell.FormulaR1C1 = _
"=""calculated GB2 from orbital measured period and measured A binary"""
Range("A19").Select
ActiveCell.FormulaR1C1 = "=R[-2]C-R[-1]C"
Range("B19").Select
ActiveCell.FormulaR1C1 = "=""GB2 is higher (-) / lower (+) than GB"""
Range("A26").Select
ActiveCell.FormulaR1C1 = "=R[4]C-R[5]C"
Range("C30").Select
ActiveCell.FormulaR1C1 = "=RC[-2]-R[1]C[-2]"
Range("C32").Select
ActiveCell.FormulaR1C1 = "=R[21]C[-2]-RC[-2]"
Range("C33").Select
ActiveCell.FormulaR1C1 = "=R[20]C[-2]-R[-1]C[-2]"
Range("C32").Select
Selection.Copy
Range("C25").Select
Selection.PasteSpecial Paste:=xlPasteValues, Operation:=xlNone, SkipBlanks _
:=False, Transpose:=False
Range("D27").Select
ActiveWindow.SmallScroll Down:=12
Range("A51").Select
Application.CutCopyMode = False
ActiveCell.FormulaR1C1 = "=R[-2]C-R[-1]C"
Range("B51").Select
ActiveCell.FormulaR1C1 = "=""differnece in mechanical energy"""
Range("A52").Select
ActiveCell.FormulaR1C1 = "=R[-20]C/2"
Range("B52").Select
ActiveCell.FormulaR1C1 = _
"=""half the potential change goes to kinetic energy"""
Range("A53").Select
ActiveCell.FormulaR1C1 = "=R[-2]C-R[-16]C"
Range("B53").Select
ActiveCell.FormulaR1C1 = _
"=""Mechanical energy to adjust after GRN accounted for"""
Range("B51").Select
ActiveCell.FormulaR1C1 = "=""difference in mechanical energy"""
Range("B1").Select


Q = 10
O = 10

For Counter30 = 1 To 16
XYZ = 0

Q = Q / 10
O = O / 10
QO = Q - 2 * O
If Counter30 > 1 Then
QO = QO + PQ
End If

'L00P

'Do Until Range("C32") < 0
Do Until XYZ > 0
'Range("B1") = Range("B1") + 0.000000000001
QO = QO - Q
Range("B1") = QO
Range("C52") = QO
'drGRN = Range("B1")
drGRN = QO

'drGRN = Range("B1")
'G = 6.67384 * 10 ^ -11 'as measured on Earth
GB = 6.6721323756E-11 'calculated GB from orbital measured period and measured A binary
r = 1949138400 ' sum of semi major axes
c = 299792458
M1 = 1.4398 * 1.9891 * 10 ^ 30
M2 = 1.3886 * 1.9891 * 10 ^ 30
ab = 1949138400 'abinary = is the sum of the semi-major axes
Pi = 3.14159265358979
e = 0.617131
mu = G * (M1 + M2)
Range("A13") = ab
Range("B13") = "Calculated Sum of Semi Major Axis in meters Weisberg and Taylor"


'solve for Gb

T = 27906.979587552 '"Measured T (in secs Weisberg and Taylor)"
T2 = 27906.9795110896 'Measured shorter period after 1 year, T (in secs Weisberg and Taylor)"
GB = ab ^ 3 * (2 * Pi / T) ^ 2 / (M1 + M2)

Range("A17") = GB
Range("B17") = " calculated GB from orbital measured period and measured A binary"
'ab2 = (ab + drGRN)'because drGRN is a negative number add it.
ab2 = (ab + drGRN)
GB2 = (ab + drGRN) ^ 3 * (2 * Pi / T2) ^ 2 / (M1 + M2)
Range("A18") = GB2
Range("B18") = " forced reduction in GB2 from orbital measured period and measured A binary"

'Potential energy difference

U = -GB * M1 * M2 / ab
Range("A30") = U
Range("B30") = "U"

Ub = -GB2 * M1 * M2 / (ab + drGRN) '(ab - drAvg)
Range("A31") = Ub
Range("B31") = "Ub"

Ubn = -GB * M1 * M2 / (ab + drGRN)
Range("A29") = Ubn
Range("B29") = "Ubn"

Range("A32") = (U - Ub)
Range("B32") = "(U- Ub) . some GPE goes to KE, find KE"
Range("A33") = (Ub - Ubn)
Range("B33") = "(Ub - Ubn). GPE if GB does not change"
'Range("A34") = Range("A33") - Range("A32")
'Range("A34") = (Ub - Ubn) - (U- Ub)

'Range("B34") = "A33-A32 change in GPE"

n = 35

'Gravitational energy changes

k_dEdt = (32 / 5) * (GB ^ 4 / c ^ 5) * ((M1 * M2) ^ 2) * (M1 + M2) 'k_dEdt is not the complete formula /r^5 is still part of it.
Power1 = -k_dEdt / r ^ 5 'Power is in Watts
Range("A" & n) = Power1 * 31557600
Range("B" & n) = "GRN Power1 * 31557600"
n = n + 1

k_dEdt = (32 / 5) * (GB2 ^ 4 / c ^ 5) * ((M1 * M2) ^ 2) * (M1 + M2) 'k_dEdt is not the complete formula /r^5 is still part of it.
Power2 = -k_dEdt / ((r + drGRN) ^ 5) ' more power radiated the shorter the r value, this is looking at 3.5 m closer.
Range("A" & n) = Power2 * 31557600
Range("B" & n) = "GRN Power2* 31557600 drGRN m closer"
n = n + 1

Range("A" & n) = (Power1 / 2 + Power2 / 2) * 31557600
Range("B" & n) = "(Power1/2 + Power2/2) * 31557600 =Average energy"
n = n + 2

F1 = GB * M1 * M2 / r ^ 2
F2 = GB2 * M1 * M2 / (r + drGRN) ^ 2
Energy = (F1 / 2 + F2 / 2) * drGRN 'average force times distance
Power3 = 0.5 * Energy / 31557600

Range("A" & n) = Power1
Range("B" & n) = "GRN Power1"
n = n + 1

Range("A" & n) = Power2
Range("B" & n) = "GRN Power2"
n = n + 1

Range("A" & n) = Power1 - Power2
Range("B" & n) = "GRN Power1 - GRN Power2"
n = n + 1

Range("A" & n) = Power3
Range("B" & n) = "GPE Power3 Joules/second"
n = n + 1

Power3 = 0.5 * Energy
Range("A" & n) = Power3
Range("B" & n) = "Power3 Joules/Year"
n = n + 1

Power3 = 0.5 * Energy / 31557600
Range("A" & n) = Power3 / (Power1 / 2 + Power2 / 2)
Range("B" & n) = "Power3/(Power1/2+Power2/2)" 'Proportion accounted for
n = n + 1

n = n + 4


'Emech = -Gb*M1*m2/(2*ab)

Emech = -GB * M1 * M2 / (2 * ab)
Range("A" & n) = Emech
Range("B" & n) = "Emech"
n = n + 1

'Emech2 = -Gb2*M1*m2/(2*ab2)

Emech2 = -GB2 * M1 * M2 / (2 * ab2)
Range("A" & n) = Emech2
Range("B" & n) = "Emech2"
n = n + 1

Z = Z + 1
If Z = 1000 Then
GoTo 99
End If
' xyz = (Emech -Emech2)- ((Power1/2+Power2/2)*31557600)- (U - Ub)
XYZ = (Emech - Emech2) - ((Power1 / 2 + Power2 / 2) * 31557600) - (U - Ub)
Range("C32") = XYZ
Range("C32") = XYZ
Range("C30") = XYZ 'previous xyz
Loop
PQ = QO + Q

Next
99 Range("C52") = QO
'Range("A2:E34").Select
'Selection.Insert Shift:=xlDown
'
'
End Sub

 

The crunch line of the above macro is "XYZ = (Emech - Emech2) - ((Power1 / 2 + Power2 / 2) * 31557600) - (U - Ub)"

Which is basically saying you need to account for all the energy changes and you are looking for the number "XYZ" that is closest to zero that one can find, but as I said already Excel is limited to 15 significant numbers. But you can see if you were able to run it on a program that would handle 45 digits you could find an exact "GB2" value where all the measured parameters and the above energy equation are satisfied as in "(Emech - Emech2) - ((Power1 / 2 + Power2 / 2) * 31557600) - (U - Ub) = 0"

Has anyone attempted to run the macro? Are you having trouble? If so let me know so I can tell you what has to be done to get it to work for you.

Here is how to start: Open up a new Excel workbook. Go to Tools > macros >record new macro and then start, but then stop recording immediately.
Go to macros > view macros and open "Macro1" Paste the following macro over top of the Macro1 template. There might be a better way but that is how I would do it.

 

Now that it can be shown that the gravitational attraction between two masses can be altered and especially altered when the masses are radiating gravitational energy, we are given a clue as to how gravity works. As a simple guess it is like magnetism in that it can be more (stronger) or less (weaker) without changing the mass of the material.
When the gravitational attraction between two masses declines part of the gravitational potential energy debt incurred has to be repaid, so you get a breakdown in the dr/dt formula, which only accounts for the amount of movement due to the loss of momentum from the gravitational energy component, but that loss also reduces G so a huge amount of extra movement is required to fund that. The dE/dt, and dp/dt formulas holds true, and they correctly account for the period change (so Weisberg and Taylor were OK on that).
But with the macro a solution was recovered that altered both "r" and "dr" and "G" (called GB2 in the macro).

It will take a bit of thought to even attempt to explain what could be happening at the subatomic/quantum level.
The only thing that hasn't held is the postulate that G is constant throughout all time and space.

 

An interesting extension of this project would be to continue to follow what will happen as the Hulse Taylor Binary stars continue their orbital decay. We have the verified Formula 3 from the Weisberg and Taylor paper that very accurately follows the orbital period change. So each year we could recalculate the orbital period for year (n) and the new Gb(n) and see how long the collapse will take. Note: As the Gb(n) weakens maybe the final collapse could be slowed, but then when they do combine into the black-hole we could have the Gravitational constant for this black-hole (GBH) . Maybe we can show that time is not zero at the event horizon of a black hole! Every time I read that time is stopped by a black-hole it never rings true with me, so one way or the other we will find out.
I could streamline the above macro and make it loop till the sum of the semi major axes equals zero. I would also need to introduce formula 3 to keep track of the orbital period changes.
Is there anyone out there who is able to help?

We'd need it to loop real quickly for the computer would have to run day and night for 9.5 years if the calculation took 1 sec/y if the period to collapse was the expected 300 million years.
We could make it a little less accurate and look at increments of 1,000 years. and get an answer in 3.5 days. So it is not impossible.

 

I streamlined and stripped the macro and got it running at 450 Hz. Running out of possible improvements but it still might be able to be tweaked up to 1000 Hz. The output value changes depending on whether I do the maths prior to entering the variables into the macro, maybe it will be more accurate if I put all the fractions in as unsolved values rather than as a calculated value. For example one equation has a component "(32 / 5)" but this may give a different result than if I was to put in 6.4E+00, (purely a bad example). This is a bit of a surprise to me.

PS: I'll have to settle on 575 Hz. 6 days run-time!
Still have to introduce formula 3 into the macro so that will slow it down again, but I have feeling once we start looking at the movement from year to year, it will be progressively larger each year so we won't have to scan from the 0 (zero) each time. So that will be a major saving of time, also instead of going to 15 decimal places well could limit it to 8 for example.

 

Still trying to get the introduction of formula 3 meaningful, but progress is being made. It will have several large measurement errors in it but even if I can get it to give me a sensible estimate maybe we can figure out a way of dealing with the errors.
Thanks for reading the thread everyone.

 

For those who work with macros try this beauty. Was Einstein wrong? Maybe time doesn't change but gravity does instead? Press escape "ESC" if it goes on and on. I haven't speed tested it completely but it looked quite useful. It is set to run the " first million years" of a possible 300 million year to full orbital decay.
Didn't take too long, but it locked up first time so I've turned screen updating off " Application.ScreenUpdating = False"
And slowed it down at one point in case it was rushing it too fast by the line ActiveSheet.calculate

Sub Finding_Delta_r_stripped_Plus3()
'
' Finding_Delta_r Macro
' Macro recorded 14/11/2012 by Robittybob1
'

'Dim r As Double
Dim drGRN As Double 'from EQ3_Power Macro
Dim GB As Double
Dim c As Double
Dim c5 As Double '6.4/c^5
Dim M1 As Double
Dim M2 As Double
Dim MX As Double
Dim MX2 As Double 'MX^2

Dim MpM As Double
Dim Emech As Double
Dim Emech2 As Double
Dim Power1 As Double
Dim Power2 As Double
Dim Pb As Double
Dim PGRb As Double
Dim e As Double
Dim U As Double
Dim Ub As Double
'Dim Ubn As Double
Dim Gb2 As Double
Dim ab2 As Double
Dim T As Double
Dim T2 As Double
'Dim G As Double
Dim Pi As Double 'Pi2
Dim QO As Double
Dim Q As Double
Dim Counter30 As Double
Dim CounterX As Double
Dim xyz As Double
Dim PQ As Double
Dim ab As Double

Dim k_dEdt As Double
Sheets("Finding_Delta_r").Select
Application.ScreenUpdating = False


GB = 6.6721323756E-11 'calculated GB from orbital measured period and measured A binary
c = 299792458
M1 = 1.4398 * 1.9891 * 10 ^ 30
M2 = 1.3886 * 1.9891 * 10 ^ 30
MX = M1 * M2
MpM = M1 + M2
MX2 = MX ^ 2
c5 = 6.4 / c ^ 5
Pb = 27906.97959 '(roughly 7.75*60*60)based on measured data 27906.979587552
ab = 1949138400 'abinary = is the sum of the semi-major axes
Pi = 3.14159265358979
Pi2 = 2 * Pi
e = 0.617131
T = 27906.979587552 '"Measured T (in secs Weisberg and Taylor)"
PQ = 0

For CounterX = 1 To 1000000
Q = 10
PGRb = -(192 * Pi * (GB ^ (5 / 3)) / (5 * (c ^ 5))) * ((T / (Pi2)) ^ -(5 / 3)) * (1 + ((73 / 24) * e ^ 2) + ((37 / 96) * e ^ 4)) * ((1 - (e ^ 2)) ^ -(7 / 2)) * (MX * ((MpM)) ^ -(1 / 3)) * 31557600
ActiveSheet.calculate
T2 = T + PGRb

For Counter30 = 1 To 8 '16
xyz = 0
Q = Q / 10
QO = PQ + Q


Do Until xyz > 0
QO = QO - Q
drGRN = QO
drGRN = QO

'solve for Gb
GB = ab ^ 3 * (Pi2 / T) ^ 2 / (MpM)
ab2 = (ab + drGRN)
Gb2 = (ab2) ^ 3 * (Pi2 / T2) ^ 2 / (MpM)

'Potential energy difference

U = -GB * MX / ab

Ub = -Gb2 * MX / (ab2)

'Gravitational energy changes

k_dEdt = c5 * (GB ^ 4) * ((MX2)) * (MpM)
Power1 = -k_dEdt / ab ^ 5 'Power is in Watts confusing ab with r

k_dEdt = c5 * (Gb2 ^ 4) * ((MX2)) * (MpM)
Power2 = -k_dEdt / ((ab + drGRN) ^ 5)



F1H = (GB * MX / ab ^ 2) / 2
F2H = (Gb2 * MX / (ab + drGRN) ^ 2) / 2

'Emech = -Gb*M1*m2/(2*ab)

Emech = -GB * MX / (2 * ab)

Emech2 = -Gb2 * MX / (2 * ab2)

xyz = (Emech - Emech2) - ((Power1 / 2 + Power2 / 2) * 31557600) - (U - Ub)
Loop

PQ = QO + Q
Next

T = T2
GB = Gb2
ab = ab + QO

99 Range("C8") = QO
Next
Application.ScreenUpdating = True

End Sub

 

After the first million years it was still moving at roughly the same rate so an estimate of years to collapse just based on the first million years, is 378,541,083 years

 

Now going for the big one ! Computer will have to run for a day. It is pretty swish now with only a fraction of the lines to begin with.

Sub Finding_Delta_r_stripped_Plus3V2()
'
' Finding_Delta_r Macro
' Macro recorded 14/11/2012 by Robittybob1
'

'Dim r As Double
Dim drGRN As Double 'from EQ3_Power Macro
Dim GB As Double
Dim c As Double
Dim c5 As Double '6.4/c^5
Dim M1 As Double
Dim M2 As Double
Dim MX As Double
Dim MX2 As Double 'MX^2

Dim MpM As Double
Dim Emech As Double
Dim Emech2 As Double
Dim Power1 As Double
Dim Power2 As Double
Dim Pb As Double
Dim PGRb As Double
Dim e As Double
Dim U As Double
Dim Ub As Double
'Dim Ubn As Double
Dim Gb2 As Double
Dim ab2 As Double
Dim T As Double
Dim T2 As Double
'Dim G As Double
Dim Pi As Double 'Pi2
Dim Pi2 As Double
Dim QO As Double
Dim Q As Double
Dim Counter30 As Double
Dim CounterX As Double
Dim xyz As Double
Dim PQ As Double
Dim ab As Double

Dim k_dEdt As Double
Sheets("Finding_Delta_r").Select
Application.ScreenUpdating = False


GB = 6.6721323756E-11 'calculated GB from orbital measured period and measured A binary
c = 299792458
M1 = 1.4398 * 1.9891 * 10 ^ 30
M2 = 1.3886 * 1.9891 * 10 ^ 30
MX = M1 * M2
MpM = M1 + M2
MX2 = MX ^ 2
c5 = 6.4 / c ^ 5
Pb = 27906.97959 '(roughly 7.75*60*60)based on measured data 27906.979587552
ab = 1949138400 'abinary = is the sum of the semi-major axes
Pi = 3.14159265358979
Pi2 = 2 * Pi
e = 0.617131
T = 27906.979587552 '"Measured T (in secs Weisberg and Taylor)"
PQ = 0

k_dEdt = c5 * (GB ^ 4) * ((MX2)) * (MpM)
Power1 = -k_dEdt / ab ^ 5 'Power is in Watts

U = -GB * MX / ab

Emech = -GB * MX / (2 * ab)

'For CounterX = 1 To 1000000
Do Until ab <= 0
Q = 10
PGRb = -(192 * Pi * (GB ^ (5 / 3)) / (5 * (c ^ 5))) * ((T / (Pi2)) ^ -(5 / 3)) * (1 + ((73 / 24) * e ^ 2) + ((37 / 96) * e ^ 4)) * ((1 - (e ^ 2)) ^ -(7 / 2)) * (MX * ((MpM)) ^ -(1 / 3)) * 31557600
ActiveSheet.Calculate
T2 = T + PGRb

For Counter30 = 1 To 10 '16
xyz = 0
Q = Q / 10
QO = PQ


Do Until xyz > 0
QO = QO - Q
'QO = PQ - Q
drGRN = QO
'drGRN = QO

'solve for Gb
'GB = ab ^ 3 * (Pi2 / T) ^ 2 / (MpM)
ab2 = (ab + drGRN)
Gb2 = (ab2) ^ 3 * (Pi2 / T2) ^ 2 / (MpM)

'Potential energy difference

' U = -GB * MX / ab

Ub = -Gb2 * MX / (ab2)

'Gravitational energy changes

' k_dEdt = c5 * (GB ^ 4) * ((MX2)) * (MpM)
' Power1 = -k_dEdt / ab ^ 5 'Power is in Watts

k_dEdt = c5 * (Gb2 ^ 4) * ((MX2)) * (MpM)
Power2 = -k_dEdt / ((ab + drGRN) ^ 5)



' F1H = (GB * MX / ab ^ 2) / 2
' F2H = (Gb2 * MX / (ab + drGRN) ^ 2) / 2

'Emech = -Gb*M1*m2/(2*ab)

'Emech = -GB * MX / (2 * ab)

Emech2 = -Gb2 * MX / (2 * ab2)

xyz = (Emech - Emech2) - ((Power1 / 2 + Power2 / 2) * 31557600) - (U - Ub)

Loop

PQ = QO + Q

' Here here


Next
Power1 = Power2
U = Ub
T = T2
GB = Gb2
Emech = Emech2
'ab = ab + QO
ab = ab2

99 Range("C8") = QO
' Next
Loop
Application.ScreenUpdating = True
Range("A1") = CounterX
Range("A2") = Power2
Range("A3") = Emech
Range("A4") = Gb2
ActiveWorkbook.Save

End Sub