- Thank you received: 0
Requiem for Relativity
15 years 9 months ago #15775
by Maurol
Replied by Maurol on topic Reply from Mauro Lacy
<blockquote id="quote"><font size="2" face="Verdana, Arial, Helvetica" id="quote">quote:<hr height="1" noshade id="quote"><i>Originally posted by nemesis</i>
<br />Thanks, Maurol. You describe the tropical year as a convenience - it sounds more like a fudge factor to me. I don't see the problem with the solstices changing calendar dates, especially over such a long time frame.
<hr height="1" noshade id="quote"></blockquote id="quote"></font id="quote">
I suppose you mean equinoxes. Well, if I recall correctly, one of the Popes with a funny name starting with G didn't think the same. There must be an historical reason why the calendar is aligned with the equinoxes; you can surely google about it.
Oh, one "detail" that I got wrong: 13000 years in the future, the winter solstice will be around 18-20 of January. It was the other way around: shorter calendars will need more days to indicate the same moment.
<br />Thanks, Maurol. You describe the tropical year as a convenience - it sounds more like a fudge factor to me. I don't see the problem with the solstices changing calendar dates, especially over such a long time frame.
<hr height="1" noshade id="quote"></blockquote id="quote"></font id="quote">
I suppose you mean equinoxes. Well, if I recall correctly, one of the Popes with a funny name starting with G didn't think the same. There must be an historical reason why the calendar is aligned with the equinoxes; you can surely google about it.
Oh, one "detail" that I got wrong: 13000 years in the future, the winter solstice will be around 18-20 of January. It was the other way around: shorter calendars will need more days to indicate the same moment.
Please Log in or Create an account to join the conversation.
- Joe Keller
- Offline
- Platinum Member
Less
More
- Thank you received: 0
15 years 9 months ago #15776
by Joe Keller
Replied by Joe Keller on topic Reply from
Synchronized Precession
Above, I note that the presumed Barbarossa/Frey objects show a binary orbit agreeing with Kepler's 2nd law to 0.125%. What about the changing aspect, of the binary orbit, due to the solar orbit? I used a kind of average ellipse which might compensate to first order, but the second order term is as large as (0.0965 radian)^2 * 0.5 = 0.47%. Maybe the phase was favorable. Anyway, there is another possible explanation: synchronized precession.
Let's compare the Sun/Earth/Luna system to the Sun/Barbarossa/Frey system. The precession period of the latter, would be 19yr
* 200^3 for weaker Sun tide
* (1/400)^2 for greater binary radius given same angular momentum
* sqrt(3000) for more massive central object
* sqrt(400) for the increase in angular momentum with binary radius
= 1,000,000 yr.
Let's replace the sun with a second, Saturn-mass moon, Freya, 2 AU from Barbarossa. Freya could precess Frey's orbit in 1,000,000 * (1/0.0003) * (2/200)^3 = 3000 yr. If this precession came to be synchronized with Barbarossa's solar orbit, the binary orbit's aspect would not change.
Above, I note that the presumed Barbarossa/Frey objects show a binary orbit agreeing with Kepler's 2nd law to 0.125%. What about the changing aspect, of the binary orbit, due to the solar orbit? I used a kind of average ellipse which might compensate to first order, but the second order term is as large as (0.0965 radian)^2 * 0.5 = 0.47%. Maybe the phase was favorable. Anyway, there is another possible explanation: synchronized precession.
Let's compare the Sun/Earth/Luna system to the Sun/Barbarossa/Frey system. The precession period of the latter, would be 19yr
* 200^3 for weaker Sun tide
* (1/400)^2 for greater binary radius given same angular momentum
* sqrt(3000) for more massive central object
* sqrt(400) for the increase in angular momentum with binary radius
= 1,000,000 yr.
Let's replace the sun with a second, Saturn-mass moon, Freya, 2 AU from Barbarossa. Freya could precess Frey's orbit in 1,000,000 * (1/0.0003) * (2/200)^3 = 3000 yr. If this precession came to be synchronized with Barbarossa's solar orbit, the binary orbit's aspect would not change.
Please Log in or Create an account to join the conversation.
15 years 9 months ago #23400
by Maurol
Replied by Maurol on topic Reply from Mauro Lacy
<blockquote id="quote"><font size="2" face="Verdana, Arial, Helvetica" id="quote">quote:<hr height="1" noshade id="quote"><i>Originally posted by Maurol</i>
Oh, one "detail" that I got wrong: 13000 years in the future, the winter solstice will be around 18-20 of January. It was the other way around: shorter calendars will need more days to indicate the same moment.
<hr height="1" noshade id="quote"></blockquote id="quote"></font id="quote">
And if we consider leap seconds, in 13000 years the calendar will fall behind between 8 and 12 days and fraction. We cannot know exactly in advance, because you cannot know with an advance bigger than 6 months if and when a leap second will be needed
That certainly doesn't look as a very good way to keep time, don't you think?
So, considering leap seconds, the winter solstice will fall between January 6-10 in the year 13000.
But don't worry, leap seconds are here not to stay
Oh, one "detail" that I got wrong: 13000 years in the future, the winter solstice will be around 18-20 of January. It was the other way around: shorter calendars will need more days to indicate the same moment.
<hr height="1" noshade id="quote"></blockquote id="quote"></font id="quote">
And if we consider leap seconds, in 13000 years the calendar will fall behind between 8 and 12 days and fraction. We cannot know exactly in advance, because you cannot know with an advance bigger than 6 months if and when a leap second will be needed
That certainly doesn't look as a very good way to keep time, don't you think?
So, considering leap seconds, the winter solstice will fall between January 6-10 in the year 13000.
But don't worry, leap seconds are here not to stay
Please Log in or Create an account to join the conversation.
- Joe Keller
- Offline
- Platinum Member
Less
More
- Thank you received: 0
15 years 9 months ago #23401
by Joe Keller
Replied by Joe Keller on topic Reply from
REM program to fit Barbarossa/Frey orbits to 3 data points, predict others
REM initialize constants
PRINT : PRINT
REM GOTO 9000
pi# = 4 * ATN(1): pi180# = pi# / 180: pi2# = 2 * pi#: crit = 0
kk# = 0 - (pi2# / 365.25#) ^ 2 * 332447 / 332001: corr# = 1 / 1.01#
REM kk# is solar accel. const. in AU/day^2, including known planets
REM dimension variables; double precision throughout
REM positive integration weights for 1954, 1986, 1987:
DIM w(3) AS DOUBLE
w(1) = 1 / 4: w(2) = 3 / 8: w(3) = 3 / 8
DIM t(6) AS DOUBLE: DIM ra(6) AS DOUBLE: DIM decl(6) AS DOUBLE
DIM x(6) AS DOUBLE: DIM y(6) AS DOUBLE: DIM z(6) AS DOUBLE
DIM r(6) AS DOUBLE
DIM xe(6) AS DOUBLE: DIM ye(6) AS DOUBLE: DIM ze(6) AS DOUBLE
DIM re(6) AS DOUBLE
DIM rab(6) AS DOUBLE: DIM declb(6) AS DOUBLE: DIM raf(6) AS DOUBLE
DIM declf(6) AS DOUBLE
DIM nax(6) AS DOUBLE: DIM nay(6) AS DOUBLE: DIM naz(6) AS DOUBLE
DIM x0(6) AS DOUBLE: DIM y0(6) AS DOUBLE: DIM z0(6) AS DOUBLE
DIM theta(6, 6) AS DOUBLE
REM get times
GOSUB 1000
dt1# = 1 / (t(2) - t(1)): dt2# = 1 / (t(3) - t(2)): dt3# = 2 / (t(3) - t(1))
REM get coordinates
GOSUB 1100
REM get barycentric Earth positions
GOSUB 1200
REM PRINT "Done initializing; now searching..."
REM best mass ratio and radius
REM for circular orbit using A&A2, B&B3, C11&C, is q0#=.8747#; r00#=209.2#
REM and rp#=-.19 is best change from rp#=0
REM best using work and general orbit is q0#=.8770#; r00#=208.3#; rp0#=-25.9#
rel0# = 10 ^ 9
nn = 1
FOR i = 1 TO nn
q# = .877# + 0 * .0001# * i: p# = 1 - q#
FOR r0# = 208.3# TO 208.3# STEP .1
FOR rp# = -25.9 TO -25.9 STEP .1
rt# = rp# * pi2# / 3024 / 365
REM 3024 is an earlier approximation of the actual period; this isn't critical
REM getting c.o.m. coords. by weighted ave. of Barbarossa & Frey coords.,
REM causes < 0.5" error
FOR j = 1 TO 5
ra(j) = rab(j) * q# + raf(j) * p#
decl(j) = declb(j) * q# + declf(j) * p#
zz# = SIN(decl(j)): cs# = COS(decl(j))
xx# = cs# * COS(ra(j)): yy# = cs# * SIN(ra(j))
csthe# = (xx# * xe(j) + yy# * ye(j) + zz# * ze(j)) / re(j)
r1# = r0# + rt# * (t(j) - t(1))
REM cosine rule: rr# ^ 2 + 2 * rr# * re(j) * csthe# + re(j) ^ 2 - r1#^2 = 0
REM use quadratic equation for rr#
b# = 2 * re(j) * csthe#: c# = re(j) ^ 2 - r1# ^ 2
rr# = (0 - b# + SQR(b# ^ 2 - 4 * c#)) * .5#
x(j) = xe(j) + rr# * xx#
y(j) = ye(j) + rr# * yy#
z(j) = ze(j) + rr# * zz#: r(j) = r1#
NEXT j
REM one-time light time correction
IF crit = 0 THEN GOSUB 1050
crit = 1
REM Minimize difference between Newtonian accel.
REM and accel. of trial curve.
REM accel. of trial curve
GOSUB 2000
REM Newtonian accel.
GOSUB 2100
REM work done by difference
GOSUB 2200
IF rel# < rel0# THEN GOSUB 3000
NEXT rp#: NEXT r0#
REM IF i = 1 THEN PRINT i; "/"; nn; "done...";
REM IF i > 1 THEN PRINT i; "/"; nn; "...";
NEXT i
PRINT : PRINT "unexplained work, AU/day^2*AU:"; rel0#
PRINT "best mass ratio, initial radius, & rp# parameter:"
PRINT q0#, r00#, rp0#
PRINT "t(3) radius, midrange of t(1) & t(3)";
rt0# = rp0# * pi2# / 3024 / 365
rmid# = r00# + .5# * rt0# * (t(3) - t(1))
PRINT r00# + 2 * (rmid# - r00#), rmid#
PRINT "theoretical circular period"
PRINT rmid# ^ 1.5# * SQR(332000 / 332447)
FOR i = 2 TO 3
FOR j = 1 TO i - 1
cs# = (x0(i) * x0(j) + y0(i) * y0(j) + z0(i) * z0(j)) / r0(i) / r0(j)
theta(i, j) = ATN(SQR(1 - cs# ^ 2) / cs#)
PRINT "sector "; j; i; ": ";
PRINT theta(i, j)
PRINT "radians/day"; " ";
om# = theta(i, j) / (t(i) - t(j))
PRINT om#
IF i = 3 AND j = 1 THEN GOSUB 3100
NEXT j: NEXT i
PRINT theta(2, 1) + theta(3, 2) - theta(3, 1)
REM GOTO 990
PRINT "first order elliptical extrapolation"
xx# = x0(3) / r0(3): yy# = y0(3) / r0(3): zz# = z0(3) / r0(3)
dt# = t(3) - t(1)
x0# = x0(1) / r0(1): y0# = y0(1) / r0(1): z0# = z0(1) / r0(1)
cs# = xx# * x0# + yy# * y0# + zz# * z0#
xx# = xx# - cs# * x0#: yy# = yy# - cs# * y0#: zz# = zz# - cs# * z0#
rr# = SQR(xx# ^ 2 + yy# ^ 2 + zz# ^ 2)
xx# = xx# / rr#: yy# = yy# / rr#: zz# = zz# / rr#
th0# = ATN(SQR(1 - cs# ^ 2) / cs#)
FOR i = 4 TO 6
th# = th0# * (t(i) - t(1)) / dt#
th# = th# / (1 + rt0# * (t(i) - .5# * (t(3) + t(1))) / rmid#)
cs# = COS(th#): sn# = SIN(th#)
r# = r00# + rt0# * (t(i) - t(1))
xpred# = r# * (cs# * x0# + sn# * xx#)
ypred# = r# * (cs# * y0# + sn# * yy#)
zpred# = r# * (cs# * z0# + sn# * zz#)
REM PRINT "barycentric c.o.m. X, Y, Z for i-th time": PRINT
REM PRINT xpred#, ypred#, zpred#
REM PRINT "geocentric X, Y, Z"
xpred# = xpred# - xe(i): ypred# = ypred# - ye(i): zpred# = zpred# - ze(i)
REM PRINT xpred#, ypred#, zpred#
REM PRINT "J2000 celestial coords."
decl# = ATN(zpred# / SQR(xpred# ^ 2 + ypred# ^ 2)) / pi180#
ra# = ATN(ypred# / xpred#) / pi180# + 180
PRINT i; "-th point:"; ra#, decl#
PRINT "RA"; INT(ra# / 15); "h"; 4 * (ra# - INT(ra# / 15) * 15); "m"
PRINT "Decl"; -INT(-decl#); "deg"; 60 * (-decl# - INT(-decl#)); "arcmin"
NEXT i
990 END
REM Get Julian date - 2400000.
1000 t(1) = 34798.5# + 8 / 24 + 5 / 1440
t(2) = 46504.5# + 14 / 24 + 18 / 1440
REM Change time from B plate to C plate.
t(3) = t(2) + 365 - 28 - 31 + 16 + 3 / 24 - 16 / 1440
tbarbarossa3# = 54184.5# + 42 / 1440: tfrey3# = 54191.5# + 7 / 24 + 39 / 1440
REM Interpolate t(4) using est. of mass ratio.
t(4) = .877# * tbarbarossa3# + .123# * tfrey3#
REM Guess that 12/22/08 epoch, t(5), is approx. time on meridian;
REM Barbarossa was near stationarity anyway.
t(5) = 54822.5# + 12 / 24 + 47 / 1440
t(6) = 54851.5# + 9 / 24 + 48.5# / 1440
RETURN
REM light time correction
1050 FOR j = 1 TO 6
t(j) = t(j) - rr# * 1.496 * 10 ^ 13 / (2.99793 * 10 ^ 10)
NEXT j
RETURN
REM get J2000 geocentric Barbarossa & Frey coordinates
REM "A" plate (1954); objects A2 & A
1100 raf(1) = pi180# * (11 * 15 + 2 / 4 + 25.16# / 240)
declf(1) = 0 - pi180# * (5 + 56 / 60 + 11.3# / 3600)
rab(1) = pi180# * (11 * 15 + 3 / 4 + 12.4# / 240)
declb(1) = 0 - pi180# * (5 + 58 / 60 + 9 / 3600)
REM "B" plate (1986); objects B3 & B
raf(2) = pi180# * (11 * 15 + 16 / 4 + 51.55# / 240)
declf(2) = 0 - pi180# * (7 + 49 / 60 + 41.1# / 3600)
rab(2) = pi180# * (11 * 15 + 16 / 4 + 56.07# / 240)
declb(2) = 0 - pi180# * (7 + 55 / 60 + 14.3# / 3600)
REM "C" plate (1987); objects C & C11
raf(3) = pi180# * (11 * 15 + 18 / 4 + 3.18# / 240)
declf(3) = 0 - pi180# * (7 + 58 / 60 + 46.1 / 3600)
rab(3) = pi180# * (11 * 15 + 18 / 4 + .41 / 240)
declb(3) = 0 - pi180# * (8 + 1 / 60 + 57.7 / 3600)
REM Genebriera March 25, 2007
rab(4) = pi180# * (11 * 15 + 26 / 4 + 22.2# / 240)
declb(4) = 0 - pi180# * (9 + 4 / 60 + 59 / 3600)
REM Riley April 1, 2007
raf(4) = pi180# * (11 * 15 + 26 / 4 + 25 / 240)
declf(4) = 0 - pi180# * (8 + 57 / 60 + 26 / 3600)
REM alternate coords. for Riley's "Frey" (see messageboard May 7, 2007)
REM raf(4) = pi180# * (11 * 15 + 26 / 4 + 24.6 / 240)
REM declf(4) = 0 - pi180# * (8 + 57 / 60 + 48.5 / 3600)
REM Genebriera's "Frey" (see messageboard May 15, 2007)
REM raf(4) = pi180# * (11 * 15 + 26 / 4 + 31.8 / 240)
REM declf(4) = 0 - pi180# * (9 + 0 / 60 + 11 / 3600)
REM Dec. 22, 2008 U. of Iowa photo
rab(5) = pi180# * (11 * 15 + 28 / 4 + 22.08# / 240)
declb(5) = 0 - pi180# * (9 + 16 / 60 + 6.4# / 3600)
raf(5) = pi180# * (11 * 15 + 29 / 4 + 4.66# / 240)
declf(5) = 0 - pi180# * (9 + 7 / 60 + 2.3# / 3600)
RETURN
REM Get J2000 barycentric Earth XYZ position (barycenter excludes Barbarossa).
REM In this program, "barycentric" = sun + traditional planets.
REM all but 1954, from Astronomical Almanac
1200 A# = 14 / 24 + 18 / 1440: i = 2
GOSUB 1220
A# = 17 / 24 + 2 / 1440: i = 3
GOSUB 1220
REM as in subroutine 1000, weighted time for 2007
A# = 42 / 1440 * .877# + (7 + 7 / 24 + 39 / 1440) * (1 - .877#): i = 4
GOSUB 1220
REM as in subroutine 1000, guess that 12/2008 epoch was meridian time
A# = 12 / 24 + 47 / 1440: i = 5
GOSUB 1220
A# = 9 / 24 + 48.5# / 1440: i = 6
GOSUB 1220
GOSUB 1250
RETURN
1220 b# = 1 - A#
READ xe1#, xe2#, ye1#, ye2#, ze1#, ze2#
xe(i) = b# * xe1# + A# * xe2#
ye(i) = b# * ye1# + A# * ye2#
ze(i) = b# * ze1# + A# * ze2#
re(i) = SQR(xe(i) ^ 2 + ye(i) ^ 2 + ze(i) ^ 2)
RETURN
REM find 1954 barycentric Earth position
1250 msun# = 332000: mjup# = 318.1#: msat# = 95.2#: mur# = 14.6#: mnep# = 17.2
me# = 1: mtot# = msun# + mjup# + msat# + mur# + mnep# + me#
sx# = 0: sy# = 0: sz# = 0
REM interpolate sun coordinates
A# = 8 / 24 + 5 / 1440: b# = 1 - A#
raa# = 22 * 15 + 30 / 4 + 33.54 / 240: rab# = 22 * 15 + 34 / 4 + 20.69 / 240
ra# = (raa# * b# + rab# * A#) * pi180#
decla# = 9 + 22 / 60 + 8.3 / 3600: declb# = 8 + 59 / 60 + 52 / 3600
decl# = 0 - (decla# * b# + declb# * A#) * pi180#
ra# = .9899423000000001#: rb# = .9901818999999999#
r# = ra# * b# + rb# * A#: w# = msun#
GOSUB 1270
REM planets interpolated to nearest RA second, or Decl arcminute
REM Jupiter
ra# = pi180# * (5 * 15 + 2 / 4 + 41 / 240)
decl# = pi180# * (22 + 29 / 60 + 51 / 3600)
r# = 4.846#: w# = mjup#
GOSUB 1270
REM Saturn
ra# = pi180# * (14 * 15 + 30 / 4 + 57 / 240)
decl = 0 - pi180# * (12 + 9 / 60 + 22 / 3600)
r# = 9.308#: w# = msat#
GOSUB 1270
REM Uranus
ra# = pi180# * (7 * 15 + 24 / 4 + 19 / 240)
decl# = pi180# * (22 + 31 / 60 + 46 / 3600)
r# = 18.001#: w# = mur#
GOSUB 1270
REM Neptune
ra# = pi180# * (13 * 15 + 38 / 4 + 22 / 240)
decl# = 0 - pi180# * (8 + 22 / 60 + 9 / 3600)
r# = 29.662#: w# = mnep#
GOSUB 1270
REM convert to 1950 ecliptic coords.
theta# = 23.45# * pi180#
GOSUB 1290
REM convert to 2000 ecliptic coords.
theta# = pi180# * 50.2619# / 3600 * 50
cs# = COS(theta#): sn# = SIN(theta#): sx0# = sx#
sx# = cs# * sx# - sn# * sy#
sy# = cs# * sy# + sn# * sx0#
REM convert to 2000 celestial coords.
theta# = 0 - 23.45# * pi180#
GOSUB 1290
xe(1) = 0 - sx# / mtot#
ye(1) = 0 - sy# / mtot#
ze(1) = 0 - sz# / mtot#
re(1) = SQR(xe(1) ^ 2 + ye(1) ^ 2 + ze(1) ^ 2)
RETURN
1270 z# = r# * SIN(decl#): cs# = COS(decl#)
x# = r# * cs# * COS(ra#): y# = r# * cs# * SIN(ra#)
sx# = sx# + w# * x#
sy# = sy# + w# * y#
sz# = sz# + w# * z#
RETURN
1290 cs# = COS(theta#): sn# = SIN(theta#): sz0# = sz#
sz# = cs# * sz# - sn# * sy#
sy# = cs# * sy# + sn# * sz0#
RETURN
REM find the accel. of the parabolas fitting 1st thru 3rd points
2000 vx1# = (x(2) - x(1)) * dt1#
vy1# = (y(2) - y(1)) * dt1#
vz1# = (z(2) - z(1)) * dt1#
vx2# = (x(3) - x(2)) * dt2#
vy2# = (y(3) - y(2)) * dt2#
vz2# = (z(3) - z(2)) * dt2#
REM correct accel. to fixed frame
ax# = (vx2# - vx1#) * dt3# * corr#
ay# = (vy2# - vy1#) * dt3# * corr#
az# = (vz2# - vz1#) * dt3# * corr#
RETURN
REM find Newtonian accel. at 1st thru 3rd points
2100 FOR j = 1 TO 3
dr# = kk# / r(j) ^ 3
nax(j) = x(j) * dr#
nay(j) = y(j) * dr#
naz(j) = z(j) * dr#
NEXT j
RETURN
REM weighted ave. sunward accel. & difference from observed
2200 avenax# = 0: avenay# = 0: avenaz# = 0
FOR j = 1 TO 3
avenax# = avenax# + w(j) * nax(j)
avenay# = avenay# + w(j) * nay(j)
avenaz# = avenaz# + w(j) * naz(j)
NEXT j
rel# = SQR((avenax# - ax#) ^ 2 + (avenay# - ay#) ^ 2 + (avenaz# - az#) ^ 2)
REM PRINT "subroutine 2200 --> Newtonian vs. apparent accel."
REM PRINT SQR(avenax# ^ 2 + avenay# ^ 2 + avenaz# ^ 2)
REM PRINT SQR(ax# ^ 2 + ay# ^ 2 + az# ^ 2)
REM work done by unexplained force
rel# = (avenax# - ax#) * (x(3) - x(1))
rel# = rel# + (avenay# - ay#) * (y(3) - y(1))
rel# = rel# + (avenaz# - az#) * (z(3) - z(1))
rel# = ABS(rel#)
RETURN
REM update best
3000 rel0# = rel#: i0 = i: q0# = q#: r00# = r0#: rp0# = rp#
FOR m = 1 TO 5
x0(m) = x(m): y0(m) = y(m): z0(m) = z(m): r0(m) = r(m)
NEXT m
RETURN
3100 PRINT "apparent circular period"
PRINT pi2# / om# / 365.25
RETURN
REM '86
DATA -.992473,-.99431,.097062,.081274,.042041,.035195
REM less accurate Earth XYZ for '87 "C" plate (see below)
REM DATA -.66451,-.66451,.6939,.6939,.30093,.30093
DATA -.647098,-.660316,.689355,.678874,.298895,.29435
REM '07
DATA -.992917,-.991887,-.057633,-.073429,-.025095,-.031943
REM '08
DATA -.008637,-.026129,.906672,.906369,.39304,.392909
REM '09
DATA -.493191,-.508316,.786488,.778431,.340941,.337449
REM data from April 2008 program
REM (times of plates or photos, Sun positions,
REM coords. of identified disappearing dots)
REM plate A
DATA 1954.154,.9900229, 22,34, 27.6,-8,-59,-7
REM object A2
DATA 11,2,25.47,-5,-56,-12
REM object A
DATA 11,3,12.4,-5,-58,-9
REM plate B
DATA 1986.201,.9945905,23,40,30.9,-2,-10,-49
REM object B3
DATA 11,16,51.55,-7,-49,-41.1
REM object B
DATA 11,16,56.07,-7,-55,-13.2
REM plate C
DATA 1987.08215,.9852006,20,53,6.9,-17,31,13
REM object C11 (fewer than 11 found; named for easy recall)
DATA 11,18,.41,-8,-1,-57.7
REM object C
DATA 11,18,3.18,-7,-58,-46.1
REM plate D (optical infrared)
REM DATA 1997.16711,.9914845,22,57,21.1,-6,-39,-54
REM DATA 11,22,32.9,-8,-26,-56
REM DATA 11,22,16.77,-8,-29,-32.6
REM DATA 11,22,1.33,-8,-34,-36.9
REM initialize constants
PRINT : PRINT
REM GOTO 9000
pi# = 4 * ATN(1): pi180# = pi# / 180: pi2# = 2 * pi#: crit = 0
kk# = 0 - (pi2# / 365.25#) ^ 2 * 332447 / 332001: corr# = 1 / 1.01#
REM kk# is solar accel. const. in AU/day^2, including known planets
REM dimension variables; double precision throughout
REM positive integration weights for 1954, 1986, 1987:
DIM w(3) AS DOUBLE
w(1) = 1 / 4: w(2) = 3 / 8: w(3) = 3 / 8
DIM t(6) AS DOUBLE: DIM ra(6) AS DOUBLE: DIM decl(6) AS DOUBLE
DIM x(6) AS DOUBLE: DIM y(6) AS DOUBLE: DIM z(6) AS DOUBLE
DIM r(6) AS DOUBLE
DIM xe(6) AS DOUBLE: DIM ye(6) AS DOUBLE: DIM ze(6) AS DOUBLE
DIM re(6) AS DOUBLE
DIM rab(6) AS DOUBLE: DIM declb(6) AS DOUBLE: DIM raf(6) AS DOUBLE
DIM declf(6) AS DOUBLE
DIM nax(6) AS DOUBLE: DIM nay(6) AS DOUBLE: DIM naz(6) AS DOUBLE
DIM x0(6) AS DOUBLE: DIM y0(6) AS DOUBLE: DIM z0(6) AS DOUBLE
DIM theta(6, 6) AS DOUBLE
REM get times
GOSUB 1000
dt1# = 1 / (t(2) - t(1)): dt2# = 1 / (t(3) - t(2)): dt3# = 2 / (t(3) - t(1))
REM get coordinates
GOSUB 1100
REM get barycentric Earth positions
GOSUB 1200
REM PRINT "Done initializing; now searching..."
REM best mass ratio and radius
REM for circular orbit using A&A2, B&B3, C11&C, is q0#=.8747#; r00#=209.2#
REM and rp#=-.19 is best change from rp#=0
REM best using work and general orbit is q0#=.8770#; r00#=208.3#; rp0#=-25.9#
rel0# = 10 ^ 9
nn = 1
FOR i = 1 TO nn
q# = .877# + 0 * .0001# * i: p# = 1 - q#
FOR r0# = 208.3# TO 208.3# STEP .1
FOR rp# = -25.9 TO -25.9 STEP .1
rt# = rp# * pi2# / 3024 / 365
REM 3024 is an earlier approximation of the actual period; this isn't critical
REM getting c.o.m. coords. by weighted ave. of Barbarossa & Frey coords.,
REM causes < 0.5" error
FOR j = 1 TO 5
ra(j) = rab(j) * q# + raf(j) * p#
decl(j) = declb(j) * q# + declf(j) * p#
zz# = SIN(decl(j)): cs# = COS(decl(j))
xx# = cs# * COS(ra(j)): yy# = cs# * SIN(ra(j))
csthe# = (xx# * xe(j) + yy# * ye(j) + zz# * ze(j)) / re(j)
r1# = r0# + rt# * (t(j) - t(1))
REM cosine rule: rr# ^ 2 + 2 * rr# * re(j) * csthe# + re(j) ^ 2 - r1#^2 = 0
REM use quadratic equation for rr#
b# = 2 * re(j) * csthe#: c# = re(j) ^ 2 - r1# ^ 2
rr# = (0 - b# + SQR(b# ^ 2 - 4 * c#)) * .5#
x(j) = xe(j) + rr# * xx#
y(j) = ye(j) + rr# * yy#
z(j) = ze(j) + rr# * zz#: r(j) = r1#
NEXT j
REM one-time light time correction
IF crit = 0 THEN GOSUB 1050
crit = 1
REM Minimize difference between Newtonian accel.
REM and accel. of trial curve.
REM accel. of trial curve
GOSUB 2000
REM Newtonian accel.
GOSUB 2100
REM work done by difference
GOSUB 2200
IF rel# < rel0# THEN GOSUB 3000
NEXT rp#: NEXT r0#
REM IF i = 1 THEN PRINT i; "/"; nn; "done...";
REM IF i > 1 THEN PRINT i; "/"; nn; "...";
NEXT i
PRINT : PRINT "unexplained work, AU/day^2*AU:"; rel0#
PRINT "best mass ratio, initial radius, & rp# parameter:"
PRINT q0#, r00#, rp0#
PRINT "t(3) radius, midrange of t(1) & t(3)";
rt0# = rp0# * pi2# / 3024 / 365
rmid# = r00# + .5# * rt0# * (t(3) - t(1))
PRINT r00# + 2 * (rmid# - r00#), rmid#
PRINT "theoretical circular period"
PRINT rmid# ^ 1.5# * SQR(332000 / 332447)
FOR i = 2 TO 3
FOR j = 1 TO i - 1
cs# = (x0(i) * x0(j) + y0(i) * y0(j) + z0(i) * z0(j)) / r0(i) / r0(j)
theta(i, j) = ATN(SQR(1 - cs# ^ 2) / cs#)
PRINT "sector "; j; i; ": ";
PRINT theta(i, j)
PRINT "radians/day"; " ";
om# = theta(i, j) / (t(i) - t(j))
PRINT om#
IF i = 3 AND j = 1 THEN GOSUB 3100
NEXT j: NEXT i
PRINT theta(2, 1) + theta(3, 2) - theta(3, 1)
REM GOTO 990
PRINT "first order elliptical extrapolation"
xx# = x0(3) / r0(3): yy# = y0(3) / r0(3): zz# = z0(3) / r0(3)
dt# = t(3) - t(1)
x0# = x0(1) / r0(1): y0# = y0(1) / r0(1): z0# = z0(1) / r0(1)
cs# = xx# * x0# + yy# * y0# + zz# * z0#
xx# = xx# - cs# * x0#: yy# = yy# - cs# * y0#: zz# = zz# - cs# * z0#
rr# = SQR(xx# ^ 2 + yy# ^ 2 + zz# ^ 2)
xx# = xx# / rr#: yy# = yy# / rr#: zz# = zz# / rr#
th0# = ATN(SQR(1 - cs# ^ 2) / cs#)
FOR i = 4 TO 6
th# = th0# * (t(i) - t(1)) / dt#
th# = th# / (1 + rt0# * (t(i) - .5# * (t(3) + t(1))) / rmid#)
cs# = COS(th#): sn# = SIN(th#)
r# = r00# + rt0# * (t(i) - t(1))
xpred# = r# * (cs# * x0# + sn# * xx#)
ypred# = r# * (cs# * y0# + sn# * yy#)
zpred# = r# * (cs# * z0# + sn# * zz#)
REM PRINT "barycentric c.o.m. X, Y, Z for i-th time": PRINT
REM PRINT xpred#, ypred#, zpred#
REM PRINT "geocentric X, Y, Z"
xpred# = xpred# - xe(i): ypred# = ypred# - ye(i): zpred# = zpred# - ze(i)
REM PRINT xpred#, ypred#, zpred#
REM PRINT "J2000 celestial coords."
decl# = ATN(zpred# / SQR(xpred# ^ 2 + ypred# ^ 2)) / pi180#
ra# = ATN(ypred# / xpred#) / pi180# + 180
PRINT i; "-th point:"; ra#, decl#
PRINT "RA"; INT(ra# / 15); "h"; 4 * (ra# - INT(ra# / 15) * 15); "m"
PRINT "Decl"; -INT(-decl#); "deg"; 60 * (-decl# - INT(-decl#)); "arcmin"
NEXT i
990 END
REM Get Julian date - 2400000.
1000 t(1) = 34798.5# + 8 / 24 + 5 / 1440
t(2) = 46504.5# + 14 / 24 + 18 / 1440
REM Change time from B plate to C plate.
t(3) = t(2) + 365 - 28 - 31 + 16 + 3 / 24 - 16 / 1440
tbarbarossa3# = 54184.5# + 42 / 1440: tfrey3# = 54191.5# + 7 / 24 + 39 / 1440
REM Interpolate t(4) using est. of mass ratio.
t(4) = .877# * tbarbarossa3# + .123# * tfrey3#
REM Guess that 12/22/08 epoch, t(5), is approx. time on meridian;
REM Barbarossa was near stationarity anyway.
t(5) = 54822.5# + 12 / 24 + 47 / 1440
t(6) = 54851.5# + 9 / 24 + 48.5# / 1440
RETURN
REM light time correction
1050 FOR j = 1 TO 6
t(j) = t(j) - rr# * 1.496 * 10 ^ 13 / (2.99793 * 10 ^ 10)
NEXT j
RETURN
REM get J2000 geocentric Barbarossa & Frey coordinates
REM "A" plate (1954); objects A2 & A
1100 raf(1) = pi180# * (11 * 15 + 2 / 4 + 25.16# / 240)
declf(1) = 0 - pi180# * (5 + 56 / 60 + 11.3# / 3600)
rab(1) = pi180# * (11 * 15 + 3 / 4 + 12.4# / 240)
declb(1) = 0 - pi180# * (5 + 58 / 60 + 9 / 3600)
REM "B" plate (1986); objects B3 & B
raf(2) = pi180# * (11 * 15 + 16 / 4 + 51.55# / 240)
declf(2) = 0 - pi180# * (7 + 49 / 60 + 41.1# / 3600)
rab(2) = pi180# * (11 * 15 + 16 / 4 + 56.07# / 240)
declb(2) = 0 - pi180# * (7 + 55 / 60 + 14.3# / 3600)
REM "C" plate (1987); objects C & C11
raf(3) = pi180# * (11 * 15 + 18 / 4 + 3.18# / 240)
declf(3) = 0 - pi180# * (7 + 58 / 60 + 46.1 / 3600)
rab(3) = pi180# * (11 * 15 + 18 / 4 + .41 / 240)
declb(3) = 0 - pi180# * (8 + 1 / 60 + 57.7 / 3600)
REM Genebriera March 25, 2007
rab(4) = pi180# * (11 * 15 + 26 / 4 + 22.2# / 240)
declb(4) = 0 - pi180# * (9 + 4 / 60 + 59 / 3600)
REM Riley April 1, 2007
raf(4) = pi180# * (11 * 15 + 26 / 4 + 25 / 240)
declf(4) = 0 - pi180# * (8 + 57 / 60 + 26 / 3600)
REM alternate coords. for Riley's "Frey" (see messageboard May 7, 2007)
REM raf(4) = pi180# * (11 * 15 + 26 / 4 + 24.6 / 240)
REM declf(4) = 0 - pi180# * (8 + 57 / 60 + 48.5 / 3600)
REM Genebriera's "Frey" (see messageboard May 15, 2007)
REM raf(4) = pi180# * (11 * 15 + 26 / 4 + 31.8 / 240)
REM declf(4) = 0 - pi180# * (9 + 0 / 60 + 11 / 3600)
REM Dec. 22, 2008 U. of Iowa photo
rab(5) = pi180# * (11 * 15 + 28 / 4 + 22.08# / 240)
declb(5) = 0 - pi180# * (9 + 16 / 60 + 6.4# / 3600)
raf(5) = pi180# * (11 * 15 + 29 / 4 + 4.66# / 240)
declf(5) = 0 - pi180# * (9 + 7 / 60 + 2.3# / 3600)
RETURN
REM Get J2000 barycentric Earth XYZ position (barycenter excludes Barbarossa).
REM In this program, "barycentric" = sun + traditional planets.
REM all but 1954, from Astronomical Almanac
1200 A# = 14 / 24 + 18 / 1440: i = 2
GOSUB 1220
A# = 17 / 24 + 2 / 1440: i = 3
GOSUB 1220
REM as in subroutine 1000, weighted time for 2007
A# = 42 / 1440 * .877# + (7 + 7 / 24 + 39 / 1440) * (1 - .877#): i = 4
GOSUB 1220
REM as in subroutine 1000, guess that 12/2008 epoch was meridian time
A# = 12 / 24 + 47 / 1440: i = 5
GOSUB 1220
A# = 9 / 24 + 48.5# / 1440: i = 6
GOSUB 1220
GOSUB 1250
RETURN
1220 b# = 1 - A#
READ xe1#, xe2#, ye1#, ye2#, ze1#, ze2#
xe(i) = b# * xe1# + A# * xe2#
ye(i) = b# * ye1# + A# * ye2#
ze(i) = b# * ze1# + A# * ze2#
re(i) = SQR(xe(i) ^ 2 + ye(i) ^ 2 + ze(i) ^ 2)
RETURN
REM find 1954 barycentric Earth position
1250 msun# = 332000: mjup# = 318.1#: msat# = 95.2#: mur# = 14.6#: mnep# = 17.2
me# = 1: mtot# = msun# + mjup# + msat# + mur# + mnep# + me#
sx# = 0: sy# = 0: sz# = 0
REM interpolate sun coordinates
A# = 8 / 24 + 5 / 1440: b# = 1 - A#
raa# = 22 * 15 + 30 / 4 + 33.54 / 240: rab# = 22 * 15 + 34 / 4 + 20.69 / 240
ra# = (raa# * b# + rab# * A#) * pi180#
decla# = 9 + 22 / 60 + 8.3 / 3600: declb# = 8 + 59 / 60 + 52 / 3600
decl# = 0 - (decla# * b# + declb# * A#) * pi180#
ra# = .9899423000000001#: rb# = .9901818999999999#
r# = ra# * b# + rb# * A#: w# = msun#
GOSUB 1270
REM planets interpolated to nearest RA second, or Decl arcminute
REM Jupiter
ra# = pi180# * (5 * 15 + 2 / 4 + 41 / 240)
decl# = pi180# * (22 + 29 / 60 + 51 / 3600)
r# = 4.846#: w# = mjup#
GOSUB 1270
REM Saturn
ra# = pi180# * (14 * 15 + 30 / 4 + 57 / 240)
decl = 0 - pi180# * (12 + 9 / 60 + 22 / 3600)
r# = 9.308#: w# = msat#
GOSUB 1270
REM Uranus
ra# = pi180# * (7 * 15 + 24 / 4 + 19 / 240)
decl# = pi180# * (22 + 31 / 60 + 46 / 3600)
r# = 18.001#: w# = mur#
GOSUB 1270
REM Neptune
ra# = pi180# * (13 * 15 + 38 / 4 + 22 / 240)
decl# = 0 - pi180# * (8 + 22 / 60 + 9 / 3600)
r# = 29.662#: w# = mnep#
GOSUB 1270
REM convert to 1950 ecliptic coords.
theta# = 23.45# * pi180#
GOSUB 1290
REM convert to 2000 ecliptic coords.
theta# = pi180# * 50.2619# / 3600 * 50
cs# = COS(theta#): sn# = SIN(theta#): sx0# = sx#
sx# = cs# * sx# - sn# * sy#
sy# = cs# * sy# + sn# * sx0#
REM convert to 2000 celestial coords.
theta# = 0 - 23.45# * pi180#
GOSUB 1290
xe(1) = 0 - sx# / mtot#
ye(1) = 0 - sy# / mtot#
ze(1) = 0 - sz# / mtot#
re(1) = SQR(xe(1) ^ 2 + ye(1) ^ 2 + ze(1) ^ 2)
RETURN
1270 z# = r# * SIN(decl#): cs# = COS(decl#)
x# = r# * cs# * COS(ra#): y# = r# * cs# * SIN(ra#)
sx# = sx# + w# * x#
sy# = sy# + w# * y#
sz# = sz# + w# * z#
RETURN
1290 cs# = COS(theta#): sn# = SIN(theta#): sz0# = sz#
sz# = cs# * sz# - sn# * sy#
sy# = cs# * sy# + sn# * sz0#
RETURN
REM find the accel. of the parabolas fitting 1st thru 3rd points
2000 vx1# = (x(2) - x(1)) * dt1#
vy1# = (y(2) - y(1)) * dt1#
vz1# = (z(2) - z(1)) * dt1#
vx2# = (x(3) - x(2)) * dt2#
vy2# = (y(3) - y(2)) * dt2#
vz2# = (z(3) - z(2)) * dt2#
REM correct accel. to fixed frame
ax# = (vx2# - vx1#) * dt3# * corr#
ay# = (vy2# - vy1#) * dt3# * corr#
az# = (vz2# - vz1#) * dt3# * corr#
RETURN
REM find Newtonian accel. at 1st thru 3rd points
2100 FOR j = 1 TO 3
dr# = kk# / r(j) ^ 3
nax(j) = x(j) * dr#
nay(j) = y(j) * dr#
naz(j) = z(j) * dr#
NEXT j
RETURN
REM weighted ave. sunward accel. & difference from observed
2200 avenax# = 0: avenay# = 0: avenaz# = 0
FOR j = 1 TO 3
avenax# = avenax# + w(j) * nax(j)
avenay# = avenay# + w(j) * nay(j)
avenaz# = avenaz# + w(j) * naz(j)
NEXT j
rel# = SQR((avenax# - ax#) ^ 2 + (avenay# - ay#) ^ 2 + (avenaz# - az#) ^ 2)
REM PRINT "subroutine 2200 --> Newtonian vs. apparent accel."
REM PRINT SQR(avenax# ^ 2 + avenay# ^ 2 + avenaz# ^ 2)
REM PRINT SQR(ax# ^ 2 + ay# ^ 2 + az# ^ 2)
REM work done by unexplained force
rel# = (avenax# - ax#) * (x(3) - x(1))
rel# = rel# + (avenay# - ay#) * (y(3) - y(1))
rel# = rel# + (avenaz# - az#) * (z(3) - z(1))
rel# = ABS(rel#)
RETURN
REM update best
3000 rel0# = rel#: i0 = i: q0# = q#: r00# = r0#: rp0# = rp#
FOR m = 1 TO 5
x0(m) = x(m): y0(m) = y(m): z0(m) = z(m): r0(m) = r(m)
NEXT m
RETURN
3100 PRINT "apparent circular period"
PRINT pi2# / om# / 365.25
RETURN
REM '86
DATA -.992473,-.99431,.097062,.081274,.042041,.035195
REM less accurate Earth XYZ for '87 "C" plate (see below)
REM DATA -.66451,-.66451,.6939,.6939,.30093,.30093
DATA -.647098,-.660316,.689355,.678874,.298895,.29435
REM '07
DATA -.992917,-.991887,-.057633,-.073429,-.025095,-.031943
REM '08
DATA -.008637,-.026129,.906672,.906369,.39304,.392909
REM '09
DATA -.493191,-.508316,.786488,.778431,.340941,.337449
REM data from April 2008 program
REM (times of plates or photos, Sun positions,
REM coords. of identified disappearing dots)
REM plate A
DATA 1954.154,.9900229, 22,34, 27.6,-8,-59,-7
REM object A2
DATA 11,2,25.47,-5,-56,-12
REM object A
DATA 11,3,12.4,-5,-58,-9
REM plate B
DATA 1986.201,.9945905,23,40,30.9,-2,-10,-49
REM object B3
DATA 11,16,51.55,-7,-49,-41.1
REM object B
DATA 11,16,56.07,-7,-55,-13.2
REM plate C
DATA 1987.08215,.9852006,20,53,6.9,-17,31,13
REM object C11 (fewer than 11 found; named for easy recall)
DATA 11,18,.41,-8,-1,-57.7
REM object C
DATA 11,18,3.18,-7,-58,-46.1
REM plate D (optical infrared)
REM DATA 1997.16711,.9914845,22,57,21.1,-6,-39,-54
REM DATA 11,22,32.9,-8,-26,-56
REM DATA 11,22,16.77,-8,-29,-32.6
REM DATA 11,22,1.33,-8,-34,-36.9
Please Log in or Create an account to join the conversation.
- Joe Keller
- Offline
- Platinum Member
Less
More
- Thank you received: 0
15 years 9 months ago #15777
by Joe Keller
Replied by Joe Keller on topic Reply from
The previous post is my new computer program for demonstrating the progress of the Barbarossa/Frey c.o.m. from 1954 through 1986 to 1987. It also allows predictions to be made for the geocentric c.o.m. position at any future time, if modified with appropriate Julian date and Earth XYZ coords.
Mauro's computer instructions worked: after right-clicking on the top bar, I had to left-click edit->selectall; then left-click on the screen; then highlight the text on the screen to be copied; then left-click edit->copy. A few of the lines might become wrapped, but otherwise it looks OK. As it appears in the messageboard now, none of the lines are wrapped.
Mauro's computer instructions worked: after right-clicking on the top bar, I had to left-click edit->selectall; then left-click on the screen; then highlight the text on the screen to be copied; then left-click edit->copy. A few of the lines might become wrapped, but otherwise it looks OK. As it appears in the messageboard now, none of the lines are wrapped.
Please Log in or Create an account to join the conversation.
15 years 9 months ago #15778
by nemesis
Replied by nemesis on topic Reply from
[/quote]
I suppose you mean equinoxes. Well, if I recall correctly, one of the Popes with a funny name starting with G didn't think the same. There must be an historical reason why the calendar is aligned with the equinoxes; you can surely google about it.
[/quote]
Is there any particular reason the equinoxes are preferred to the solstices? It's like the lunar cycle - a lunar month is traditionally reckoned from new moon to new moon but could just as well go from full moon to full moon. That makes more sense to me because it's easier to see when the moon is full.
I suppose you mean equinoxes. Well, if I recall correctly, one of the Popes with a funny name starting with G didn't think the same. There must be an historical reason why the calendar is aligned with the equinoxes; you can surely google about it.
[/quote]
Is there any particular reason the equinoxes are preferred to the solstices? It's like the lunar cycle - a lunar month is traditionally reckoned from new moon to new moon but could just as well go from full moon to full moon. That makes more sense to me because it's easier to see when the moon is full.
Please Log in or Create an account to join the conversation.
Time to create page: 0.758 seconds