Requiem for Relativity

Harrington Corroborated Lowell

RS Harrington, Astronomical Journal 96:1476+, 1988, made a Planet X prediction comparable to Lowell's. A post above, has details of my interpretation of Lowell's final prediction, with a reference to Lowell's final description of his work. Harrington's trial orbits ranged in semimajor axis from 30 to 80 AU, vs. Lowell's two extreme apse cases whose semimajor axis averaged 44AU. Harrington's trial perturbations didn't improve Neptune's fit much, so, like Lowell, he relied on Uranus; but Harrington's Uranus data spanned 1833-1988, vs. 1715-1914 for Lowell's.

In Figs. 1 & 2, p. 1477, Harrington's densest clusters of best fit Planet X positions, are for 1988 at RA 16h50m Decl -30; and for 1930 at RA 14h30m Decl -26. My approx. slide rule extrapolation of these along a constant speed great circle to 1914, is RA 208 Decl -23. Lowell's final prediction amounted to ecliptic longitude 202, ecl. latitude small and indeterminate; if ecliptic latitude zero, this would be RA 200 Decl -9.

Pulsar Timing Confirms Dayton Miller's Ether Drift Axis


Data. I included all the millisecond pulsars (MSPs) from Taylor's 1995 catalog (online, VizieR) and from the current (2008) online ATNF pulsar catalog, that met these criteria:

1. P < 30ms.
2. Pdot/P < 8*10^(-18) (cgs units) (This removed six MSPs which formed a second modal peak in Pdot/P; three of these were obvious statistical outliers vis a vis Pdot.)
3. Proper Motion given; and, not "0" in Taylor's catalog (to avoid difficult or conflicting measurements).
4. Not aligned with any globular cluster (according to my check vs. Bica's 2006 globular cluster catalog on VizieR).

All MSPs had distances given. All Pdot/P were positive. For the five acceptable MSPs in Taylor's catalog I used Taylor's values. Seventeen more came from the ATNF catalog, for N=22 total.


Methods. Except for the farthest pulsars, theoretical acceleration due to rotation of the Milky Way is small, and I neglected it in the main study. I did use the Shklovsky correction term (to correct acceleration for transverse motion)(see Taylor, ApJ 411:674, 1993; Zakamska & Tremaine, AJ 130:1939, 2005). I found that Pdot/P, thus corrected for transverse motion, correlates with the cosine of the angle formed between the line of sight to the pulsar, and the ether drift vector identified by Miller in his 1933 Reviews of Modern Physics article. That is, looking toward the direction which Miller thought our solar system to be moving through the ether, one observes greater Pdot/P, as if pulsars in that direction are accelerating away from us. This suggests that our solar system is braking relative to the ether.

Experimenting with various linear combinations of the raw Pdot/P term and the Shklovsky term, on various subsamples, discovered a better correlation, when the Shklovsky term is multiplied by about -1/2 (but the difference in correlation, vs. that for -1/3, is tiny), and then the entire corrected Pdot/P is multiplied by -1 for Pdot/P < 2.43*10^(-18) (s/s)/s (N=10) and multiplied by +1 for Pdot/P > 2.43*10^(-18) /s. I call this overall +/- 1 factor, the "reversal" option: small positive Pdot/P really is negative, when corrected for the Hubble inflation of the universe; upstream in the ether, big Pdot/P looks bigger and small (really negative, when corrected for inflation) Pdot/P looks smaller (really, when corrected, bigger negative). For this small sample, all coefficients between 2.2 & 2.6, inclusive, are equivalent, but I write 2.43, because 2.43*10^(-18)(cm/s)/cm, i.e., /s, equals 75 (km/s)/Mpc, a recent determination of the Hubble constant from a meta-analysis of all the most accurate determinations.

There is a theory for the choice, - 1/3; see below. I call this use of a -1/3 factor in the Shklovsky term, the "bent ether" option: the apparent transverse velocity might really be due to a lightpath along a fiber of ether which usually is progressively unbending; something is added to Pdot/P to compensate for this pseudo-acceleration toward the observer.


Results. Correlation coefficient, normalized Fisher z, & direction:

usual "etherless" Shklovsky correction: cc=0.397, N=22, z=1.83, p=0.034 (one-tailed), (l,b)=(301,- (0.1 radian search grid)

"Reversal + bent ether, -1/2 factor" option: cc=0.482, N=22, z=2.29, p=0.011 (one-tailed), (l,b)=(305.4,-25.6) (0.004 radian local search grid)

Dayton Miller's ether drift: (l,b)=(282.0,-35.2)


Discussion. The "reversal + bent ether" option gives considerably better correlation, and also its best axis is closer to Miller's ether drift (25deg away). In the "bent ether" option, the pathlength increment isn't plus theta^2 / 2, as in textbook optics; it's minus theta^2 / 6, as for a slightly curved line that unbends. Two of the three Taylor catalog MSPs disqualified for excessive Pdot/P (those with the 2nd & 3rd biggest Pdot/P of the three, but not with obvious statistical outlier Pdot, and also among the nearest MSPs to Earth) have precisely equal Pdot/P, 1.11497*10^(-17) in cgs units, when the (-1/3 or - 1/2) factor by which the transverse velocity correction term is multiplied in my "bent ether" option, is changed to very nearly (-1/5) (precisely, -0.1988). Thus the MSPs, seem to correspond to the cases of a slowly straightening quadratic or cubic planar curve, whose end tangents lie on right circular cones whose altitudes are the line between the ends, and whose vertices are the endpoints, of the line of sight. As one cone rotates with respect to the other, the curve's tangents at its endpoints, are coplanar in two cases: one approximable by a quadratic, and one by a cubic, Maclaurin polynomial.


Preliminary finding. I was led to the foregoing, by a preliminary finding based on an exact solution for Taylor's five included MSPs only. One of these had Pdot/P = 3.3*10^(-18); the other four all had Pdot/P < 2.5*10^(-18). Adjustment of the factor by which the Shklovsky term is multiplied, and fitting to the best axis, amounts to adjusting three parameters, which often will suffice to fit five points to perfect correlation. A factor of exactly -1/3, gave the perfect fit, a constant plus cosine. The factor became exactly -1/3, when 20% of the Milky Way rotational acceleration correction term was included. This confirms that visible matter and only visible matter (99% stars) produces the real Milky Way gravity. As I've discussed on this messageboard, and also in my article submitted Feb. 2002 to Aircraft Engineering & Aerospace Technology, most of Oort's law and the like, is ether effect, not really due to motion.

For these five, the unique perfectly fitting axis of most *negative* (least positive) corrected Pdot/P, is galactic coords (l,b)=(276.3,-22.6). Though the sign of the effect on Pdot/P is reversed, this is only ~30deg from the axis found for the larger sample above. Dayton Miller's final ether drift coords (from a webpage of Dr. James deMeo, the geographer) were RA 4h54m Decl -70deg33', which assuming B1900.0 coords., are by Martindale's online coordinate transformation utility, (l,b)=(282.0,-35.2); this pulsar axis and Miller's axis differ only 13.5deg. (WW Campbell's approximate compromise solar motion antapex was RA 6h Decl -30, which is, again assuming B1900.0 coords, (l,b)=(236.2,-22..)

The maximum of this perfectly fitted cosine function, differs only 1% from the common Pdot/P value for the two abovementioned Taylor catalog pulsars which had been excluded because of big Pdot/P. (These were found to have equal Pdot/P when adjusted using a -1/5 factor, i.e., unbending cubic curve ether fiber, for the Shklovsky term.) One of these two excluded "cubic fiber" pulsars lies near the Miller axis and one near its antipode, but both exhibit big Pdot/P as if they lay at whichever end of the axis gives them the biggest Pdot/P.


Determination of Barbarossa acceleration. Barbarossa's 1983.2 position was (l,b)=(265.9,+48.4); the (preliminary, perfect fit to five Taylor pulsars) pulsar and Barbarossa axes differ 71.5deg. The (preliminary fit to five pulsars) pulsar axis roughly interpolates the Miller & Barbarossa axes along a great circle. Barbarossa's gravity would cause a cosine-dependent solar acceleration of (semi)amplitude a/c=0.5216 * 10^(-17) in cgs units. In those units, the five pulsars' cosine-dependent apparent acceleration has (semi)amplitude 1.775 * 10^(-17). If the pulsar axis resulted from the arithmetic sum of Barbarossa gravitational and Miller ether drift effects, the axes wouldn't add vectorially. The spherical harmonic addition formula (see, inter alia, Jahnke & Emde) would be needed. However, the underlying physics are largely unknown. There might be an underlying physical process which manifests accurately as vectorial addition:

1.775 - 0.5216 * cos71.5 = 1.609
0.5216 * sin71.5 / 1.609 = tan17

and 17 is close to the actual 13.5deg separation of the (preliminary) pulsar & Miller axes. So, this estimate of the pulsar axis, found by fitting exactly to a small number of pulsars (which I discovered the night of May 19, 2008), corroborates the existence of both the Miller ether drift, and of Barbarossa.

You'll get an "A" on your term paper. I guarantee it.

The foregoing post would make a good term paper for a class like Statistics 101. I'll help you with it, if you'll tell the teacher that I helped you. If you don't get an "A", I'll talk to the teacher myself.

Taylor's 1995 pulsar catalog (online, "VizieR" website) has only about a page of millisecond (P < 0.025s) pulsars. Only a few remain after those lacking Pdot (the zero entry means no data), lacking Proper Motion (click on the VizieR line number to see), or aligned with globular clusters (check against the globular cluster catalog) are discarded.

Your paper can assess the effect on statistical significance, of discarding 3 of the 8 pulsars, of adjusting one parameter (the coefficient of the transverse motion "Shklovsky correction term"; Cowling, MNRAS 204:1237+, 1983), and of choosing the axis in space that gives the best correlation (8-3-1-2 = 2 d.o.f., and there's always one line through 2 points; but this isn't proof that the correlation isn't significant). The Monte Carlo method could be used to assess statistical significance: make up fictitious pulsars with data randomly chosen in the same range as the real pulsar data. Follow the same procedure as above, but fail to get the same level of significance. Repeat, say, 10 times. I haven't tried this myself, and would like to see what happens.

*******

The pulsar distances in Taylor's catalog are determined, at least usually, from signal dispersion with, at least to a considerable extent, the assumption of constant interstellar "free electron density" (Taylor, ApJ 411:674+, 1993). This unique method available for determining pulsar distance, might result in unique accuracy. The "free electron density" might really be constant or as good as constant, if it is involved with the phenomenon which Miller, Galaev, etc., have called the ether.

*********

I'll be busy for awhile rechecking Standish's claim about Neptune, using original USNO data I'm copying by hand from bound volumes. Would someone with access to such a library, email (through this messageboard) or mail me (Joe Keller, POB 9122, Ames, IA 50014 USA) copies of four pages from the 1905 American Ephemeris? (It's missing from ISU.)

Needed:

Feb. 1905: "THE SUN'S" "AT GREENWICH MEAN [not apparent] NOON" table with RA, Decl, etc., and probably on the other side of the sheet, "THE SUN'S" "AT GREENWICH MEAN [not apparent] NOON" table with Logarithm of Radius Vector, etc.

Dec. 1905: ".

Update May 31, 2008: no responses re ephemeris request; 5/29, I requested the volume by interlibrary loan.

Plan for Even More Rebuttal of Standish's Claim

If I knew all the radii from the sun to Neptune, then the tabulated data would tell me all about Neptune's position, and I could check what net unknown acceleration Neptune underwent after known tidal accelerations (i.e., acceleration of Neptune minus acceleration of the sun) are subtracted. I'll consider a segment of Neptune's orbit for which the tangential tidal acceleration due to Barbarossa should be especially big.

I'll determine Neptune's radii from Hamilton's principle (numerically extremizing the integral of Lagrangian*dt with respect to variation of all the Neptune-sun radii), assuming that there is no unknown acceleration. If there really isn't any unknown acceleration, then this result also will extremize the integral of L*dt with respect to variation of all Neptune's celestial coordinates (which are known) (after extremization with respect to the unknown radii has been performed). The actual unknown acceleration will be that which causes the integral of L*dt to be extremized with respect to Neptune's celestial coordinates, as soon as it is extremized with respect to Neptune's radii.

In detail, I'll vary the radii at the endpoints of the orbit, over a generous range compatible with the ephemeris. For each pair of endpoint radii, I'll start with ephemeris radii at all the points, then vary them by adding trial half or whole Fourier sine terms (1/2, 1, 3/2, 2, 5/2, etc., sine waves in the interval) of 1/200,000 radius amplitude (analogous to 1 arcsec). For a given sine term, two such trials, plus the original, will allow a quadratic estimate of the coefficient needed for extremization. I can go through all the sines until the period is smaller than the gap between the two nearest points, then repeat the cycle until there is enough convergence. (I'll approximate the orbit with a cubic spline, to compute each action integral.)

This extremization effectively incorporates the scatter in the data. That is, an error in the coordinate data, will cause the extremal action integral to be less small, but the extremal curve remains, in a sense, the best approximation to the actual orbit.

Then, for the same pair of endpoint radii, I'll do the foregoing three more times, once assuming each theoretical component of Barbarossa's tidal influence (at the midpoint of the orbit segment). Each of the four results then can be tested by altering each coordinate point in turn by 1" in RA and in Decl, and evaluating each new action integral (again interpolating with cubic splines). Since Barbarossa's actual influence has been more or less omitted, there is no extremization with respect to the coordinates, and the measured changes in the action integrals are linear functions of the unknown vector influence.

So, the components of Barbarossa's influence can be estimated by minimizing, the sum of squares of the changes in the action integrals (extremized w.r.t. the intermediate radii) when each data coordinate is varied in turn. For convenience, the projection on the orbital plane, of Barbarossa's influence, can be plotted as a vector field w.r.t. the assumed endpoint radii. Also the endpoint radii found from various ephemerides can be marked on this vector field, to get error bars for the result.

Summarizing: the influence of Barbarossa is that which causes the orbit of Neptune, extremized w.r.t. radii according to Hamilton's principle assuming perfect accuracy of the data coordinates, to be also most nearly extremized w.r.t. the data coordinates.

Update May 31, 2008: Following essentially the above plan, I'm more than half done writing the program, 21 monitor screens-full of Basic code so far. When I'm done, I'll post the entire program to this messageboard.

Here's the program I mentioned May 22 & 31 ("notepad" is useful for opening these files):


REM program name: USNO.BAS or USNONEPT.UNE
REM This program in QBasic, uses Hamilton's principle to search for
REM tidal forces acting on Neptune, based on season's widest from opposition
REM USNO 2d ser. vol. IX pt. 1 positions 1903-1911; n=14.
REM run time per i5,j5 step = 12min on Intel Pentium4 1.6 GHz (9 steps = 108min);
REM step = 3.6hr on IBM 486 (1994).

REM 1454 dbl precision vars in arrays, 65 dbl prec solo; 15 integer; 1534 tot
REM initialize consts; everything double precision
100 PRINT : PRINT : pi# = ATN(1) * 4: pi180# = pi# / 180: twelveinv# = 1 / 12
REM "clight" in AU/hr; other consts in Earthmass-AU-day units
REM clight & AU fr. Columbia Enc.
clight# = 3600 * 2.99792458# / 1496.0497#
REM grav fr. 2006 CODATA; masses fr. TP Snow 1983
grav# = 6.67428 * 5.974 * 10 ^ (-8 + 27 - 3 * 13 + 2 * (1 + 3)) * (2.4 * 3.6) ^ 2 / 1.4960497# ^ 3
REM include Merc,Venus,Mars,asteroids & centaurs with sunmass
n = 14: DIM xyz(6, n, 3) AS DOUBLE: DIM xcom(n, 3) AS DOUBLE
DIM ne(6, n, 4) AS DOUBLE: DIM kin(n) AS DOUBLE: DIM m(6) AS DOUBLE
DIM deriv(n, 4) AS DOUBLE: DIM delta(n, 4) AS DOUBLE
DIM qdel(n) AS DOUBLE: DIM qdelinv(n) AS DOUBLE: DIM denom(n) AS DOUBLE
DIM sq(n) AS DOUBLE: DIM diffsq(n) AS DOUBLE: DIM avedel(n) AS DOUBLE
REM more dbl prec var arrays def'd @ lines 210, 220, 240
m(1) = 17.2: m(2) = 332943 + .056 + .851 + .107 + .02
REM include 4,5 moons with Jup,Sat resp.; Luna with Earth
m(4) = 318.1 + (8.89 + 4.85 + 14.9 + 10.7) / 7.35 / 81
m(5) = 95.2 + (13.5 + .076 + .105 + .244 + .188) / 7.35 / 81
m(6) = 14.6: m(3) = 1 + 1 / 81: mm0# = m(2) + m(3) + m(4) + m(5) + m(6)
REM mm012# adjusts kinetic energy, to finite mass of sun et al
mm012# = 1 / (1 + m(1) / mm0#)
mbarbarossa# = 3430: tidebarbarossa# = grav# * mbarbarossa# / mm012# / 197.7 ^ 3
incrth# = 1 / 10 ^ .75 * pi180# / 3600: incrr# = incrth# * 30.07: icounter = 0

REM get times, RAs & Decls of Neptune from data bloc 2000
200 DIM celest(6, n, 3) AS DOUBLE: DIM time(n) AS DOUBLE
REM 1st entry in celest denotes resp. Nept,Sun,Earth,J,Sat,U
REM 3rd denotes resp. RA in rad, Decl, radius in AU
REM Nov. 16, 1858, noon = JD 2,400,000; time(i)=JD-2,400,000
juliannoon1900# = 46 + 365 * 41 + 10
FOR i = 1 TO n
READ a#, b#, c#, d#, e#, f#, g#, h#, ee#, ff#, gg#, hh#
time(i) = a# + juliannoon1900#
equcorr# = -.038 - (.066 - .038) * INT((d# + .1) / 1907)
celest(1, i, 1) = (e# * 15 + f# / 4 + (g# + h# + equcorr#) / 240) * pi180#
i1 = 1: IF ee# < 0 THEN i1 = -1
celest(1, i, 2) = (ee# + i1 * ff# / 60 + (i1 * gg# + hh#) / 3600) * pi180#
NEXT
gap142# = time(14) - time(2)
gap62# = time(6) - time(2): gap146# = gap142# - gap62#
gap102# = time(10) - time(2): gap1410# = gap142# - gap102#

REM get RAs & Decls, radii, derivatives of sun from data bloc 2100
210 DIM rasundot(n) AS DOUBLE: DIM declsundot(n) AS DOUBLE: DIM radiussundot(n) AS DOUBLE
FOR i = 1 TO n
READ a#, b#, c#, d#, e#, f#, r#, rasundot(i), declsundot(i), radiussundot(i)
celest(2, i, 1) = (a# * 15 + b# / 4 + c# / 240) * pi180#
i1 = 1: IF d# < 0 THEN i1 = -1
REM need to write, e.g., S 0deg 17', as -0.00000001, 17, etc.
celest(2, i, 2) = (d# + i1 * e# / 60 + i1 * f# / 3600) * pi180#
celest(2, i, 3) = 10 ^ (r# - 10)
NEXT

REM Neptune-sun dists
220 DIM nsd(n) AS DOUBLE: DIM nsd0(n) AS DOUBLE
REM roughly interpolated Naut. Alm. Neptune-sun dists, all Feb or early March
nsd0(2) = 10 ^ 1.475745#: nsd0(6) = 10 ^ 1.4760451#
nsd0(10) = 10 ^ 1.4764035#: nsd0(14) = 10 ^ 1.47661265#

REM Jup, Sat, Ur positions
GOSUB 3260
REM done reading data


REM heart of program
240 DIM acin(5, 5) AS DOUBLE: DIM actint(n, 2, 7) AS DOUBLE
DIM co(6, 7) AS DOUBLE: DIM hess(6) AS DOUBLE: DIM lambda(3) AS DOUBLE

FOR i5 = -1 TO 1

nsd(2) = nsd0(2) + i5 * incrr#

FOR j5 = -1 TO 1

nsd(14) = nsd0(14) + j5 * incrr#
FOR i = 1 TO 6
FOR j = 1 TO 7
co(i, j) = 0
NEXT: NEXT
PRINT i5, j5
FOR i6 = 0 TO 6
PRINT i6; " ";
REM i6>0 signifies entry in Hessian for tidal potential
FOR i7 = 1 TO n
FOR j7 = 1 TO 2
REM i7 & j7 signify node and component of perturbation of observed coords
actint(i7, j7, i6) = 0
FOR i1 = 1 TO -2 STEP -3
celest(1, i7, j7) = celest(1, i7, j7) + incrth# * i1
actionint0# = 10 ^ 12
FOR i9 = -2 TO 2
nsd(6) = nsd0(6) + (i9 + (gap146# * i5 + gap62# * j5) / gap142#) * incrr#
FOR j9 = -2 TO 2
nsd(10) = nsd0(10) + (j9 + (gap1410# * i5 + gap102# * j5) / gap142#) * incrr#
REM Lagrange cubic interpol for nsd
GOSUB 1220
REM Nept-Earth dists
GOSUB 1240
REM 1-time aberration corr (sun coords, time fns, only)
IF icounter = 0 THEN GOSUB 3250
IF icounter = 0 THEN GOSUB 3600
icounter = 1
REM geocentric rect coords
GOSUB 1300
REM c.o.m. (excluding Nept) rect coords
GOSUB 1320
REM Neptunian ecliptic coords
GOSUB 1340
k1 = 1: k3 = 3
REM quadratic numerical differentiation
GOSUB 1500
REM calculate action integral, 1st-order Euler-Maclaurin summation
GOSUB 1550
acin(i9 + 3, j9 + 3) = actionint#
NEXT j9: NEXT i9
REM est actionintegral min quadratically
GOSUB 1600
actint(i7, j7, i6) = actint(i7, j7, i6) + actionint0# * (i1 + .25 * (ABS(i1) - i1))
NEXT i1
IF i6 > 0 THEN actint(i7, j7, i6) = actint(i7, j7, i6) - actint(i7, j7, 0)
celest(1, i7, j7) = celest(1, i7, j7) + incrth#
NEXT j7: NEXT i7
NEXT i6
FOR i7 = 1 TO n
FOR j7 = 1 TO 2
FOR i6 = 1 TO 6
co(i6, 7) = co(i6, 7) - actint(i7, j7, 0) * actint(i7, j7, i6)
FOR j = 1 TO 6
co(i6, j) = co(i6, j) + actint(i7, j7, j) * actint(i7, j7, i6)
NEXT j: NEXT i6: NEXT: NEXT
REM solve 6x6 linear syst
GOSUB 1800
REM diagonalize 3x3 Hessian matrix
GOSUB 1850
NEXT j5: NEXT i5

990 END

REM Lagrange cubic to find rest of n=14 from 2d,6th,10th,14th
1220 FOR i = 1 TO n
IF i = 2 OR i = 6 OR i = 10 OR i = 14 GOTO 1228
nsd(i) = 0
FOR j = 2 TO 14 STEP 4
u# = 1
FOR k = 2 TO 14 STEP 4
IF k = j GOTO 1224
u# = u# * (time(i) - time(k)) / (time(j) - time(k))
1224 NEXT k
nsd(i) = nsd(i) + u# * nsd(j)
NEXT j
1228 NEXT i
RETURN

REM law of cosines
1230 bb# = celest(2, i, 2): dd# = celest(2, i, 1): r2# = celest(2, i, 3)
cs# = SIN(aa#) * SIN(bb#) + COS(aa#) * COS(bb#) * COS(cc# - dd#)
REM r1^2 -r2*cs*r1+r2^2-r3^2=0; use quadratic eqn
r1# = .5 * (r2# * cs# + SQR((r2# * cs#) ^ 2 - 4 * (r2# ^ 2 - r3# ^ 2)))
RETURN

REM Neptune-Earth dists
1240 FOR i = 1 TO n
aa# = celest(1, i, 2): cc# = celest(1, i, 1): r3# = nsd(i)
GOSUB 1230
celest(1, i, 3) = r1#
NEXT
RETURN

REM convert to geocentric rect. coords.
1300 FOR j = 1 TO 6
IF j = 3 GOTO 1308
FOR i = 1 TO n
cs# = COS(celest(j, i, 2))
xyz(j, i, 3) = SIN(celest(j, i, 2)) * celest(j, i, 3)
xyz(j, i, 1) = cs# * COS(celest(j, i, 1)) * celest(j, i, 3)
xyz(j, i, 2) = cs# * SIN(celest(j, i, 1)) * celest(j, i, 3)
NEXT i
1308 NEXT j
RETURN

REM convert to c.o.m. (excluding Nept) rect. coords.
1320 FOR i = 1 TO n
xcom(i, 1) = 0: xcom(i, 2) = 0: xcom(i, 3) = 0
FOR j = 2 TO 6
IF j = 3 GOTO 1326
FOR k = 1 TO 3
xcom(i, k) = xcom(i, k) + m(j) * xyz(j, i, k)
NEXT k
1326 NEXT j
FOR k = 1 TO 3
xcom(i, k) = xcom(i, k) / mm0#
NEXT k
NEXT i
FOR i = 1 TO n
xyz(3, i, 1) = 0: xyz(3, i, 2) = 0: xyz(3, i, 3) = 0
FOR j = 1 TO 6
FOR k = 1 TO 3
xyz(j, i, k) = xyz(j, i, k) - xcom(i, k)
NEXT: NEXT: NEXT
RETURN

REM convert to Neptunian ecliptic coords; 1st time pt. = 0 long.
REM 1st & nth def. plane; c.o.m. (excluding Nept) is ctr
REM 4th dimension of "ne" is Lagrangian, not time
1340 a1# = xyz(1, 1, 1): a2# = xyz(1, 1, 2): a3# = xyz(1, 1, 3)
b1# = xyz(1, n, 1): b2# = xyz(1, n, 2): b3# = xyz(1, n, 3)
maginv# = 1 / SQR(a1# ^ 2 + a2# ^ 2 + a3# ^ 2)
a1# = a1# * maginv#: a2# = a2# * maginv#: a3# = a3# * maginv#
GOSUB 1350
b1# = c1#: b2# = c2#: b3# = c3#
GOSUB 1350
c1# = -c1#: c2# = -c2#: c3# = -c3#
FOR j = 1 TO 6
FOR i = 1 TO n
x# = xyz(j, i, 1) * a1# + xyz(j, i, 2) * a2# + xyz(j, i, 3) * a3#
y# = xyz(j, i, 1) * c1# + xyz(j, i, 2) * c2# + xyz(j, i, 3) * c3#
z# = xyz(j, i, 1) * b1# + xyz(j, i, 2) * b2# + xyz(j, i, 3) * b3#
nee# = ATN(y# / x#) + .5 * pi# * (x# - ABS(x#)) / x#
IF nee# < 0 THEN nee# = nee# + pi#
ne(j, i, 1) = nee#
ne(j, i, 2) = ATN(z# / SQR(x# ^ 2 + y# ^ 2))
ne(j, i, 3) = SQR(x# ^ 2 + y# ^ 2 + z# ^ 2)
NEXT: NEXT
RETURN

REM normalized vector cross product
1350 c1# = a2# * b3# - b2# * a3#
c2# = a3# * b1# - b3# * a1#
c3# = a1# * b2# - b1# * a2#
maginv# = 1 / SQR(c1# ^ 2 + c2# ^ 2 + c3# ^ 2)
c1# = c1# * maginv#: c2# = c2# * maginv#: c3# = c3# * maginv#
RETURN

REM quadratic numerical differentiation
1500 FOR i = 2 TO n
FOR k = k1 TO k3
delta(i, k) = ne(1, i, k) - ne(1, i - 1, k)
NEXT: NEXT
FOR i = 2 TO n - 1
FOR k = k1 TO k3
deriv(i, k) = (delta(i, k) * qdel(i) + delta(i + 1, k) * qdelinv(i)) * denom(i)
NEXT: NEXT
FOR k = k1 TO k3
deriv(1, k) = 2 * delta(2, k) * denom(1) - deriv(2, k)
deriv(n, k) = 2 * delta(n, k) * denom(n) - deriv(n - 1, k)
NEXT
RETURN

REM calc action integral
REM calc. kin=v^2/2="T"
1550 FOR i = 1 TO n
sph# = deriv(i, 2) ^ 2 + (COS(ne(1, i, 2)) * deriv(i, 1)) ^ 2
kin(i) = mm012# * .5 * (deriv(i, 3) ^ 2 + ne(1, i, 3) ^ 2 * sph#): NEXT
REM calc pot#="U" & action=T-U
1554 FOR i = 1 TO n: pot# = 0: FOR j = 2 TO 6: r2# = 0: FOR k = 1 TO 3
del# = xyz(j, i, k) - xyz(1, i, k): r2# = r2# + del# ^ 2
NEXT k: pot# = pot# - m(j) / SQR(r2#): NEXT j
IF i6 = 0 GOTO 1558
j = 2
IF i6 < 4 THEN LET j = 1
IF i6 = 6 THEN LET j = 3
k = 2
IF i6 = 1 THEN LET k = 1
IF (i6 = 3 OR i6 > 4) THEN LET k = 3
k1 = 2
IF j = k THEN LET k1 = 1
pot# = pot# - tidebarbarossa# * xyz(1, i, j) * xyz(1, i, k) * k1
1558 ne(1, i, 4) = kin(i) - pot#: NEXT i
REM 1st order Euler-Maclaurin integration
s# = 0: ss# = 0: k1 = 4: k3 = 4
GOSUB 1500
REM 0th order step = trapezoidal rule
FOR i = 1 TO n: s# = s# + ne(1, i, 4) * avedel(i): NEXT
REM 1st order step
FOR i = 1 TO n: ss# = ss# - deriv(i, 4) * diffsq(i): NEXT
actionint# = s# - twelveinv# * ss#
RETURN

REM est actionintegral min quadratically
1600 FOR i = 1 TO 5
FOR j = 2 TO 4
b# = acin(i, j): a# = acin(i, j - 1): c# = acin(i, j + 1)
GOSUB 1650
NEXT: NEXT
FOR j = 1 TO 5
FOR i = 2 TO 4
b# = acin(i, j): a# = acin(i - 1, j): c# = acin(i + 1, j)
GOSUB 1650
NEXT: NEXT
RETURN

REM Newton ("Bessel") quad interp
1650 d# = c# - 2 * b# + a#
IF d# < 0 GOTO 1658
actionint# = b# - .125 * (a# - c#) ^ 2 / d#
IF actionint# < actionint0# THEN LET actionint0# = actionint#
1658 RETURN

REM solve 6x6 linear syst for 3x3 Hessian matrix
REM row elim
1800 FOR i = 1 TO 5
FOR j = i + 1 TO 6
a# = co(j, i) / co(i, i)
FOR k = i + 1 TO 7
co(j, k) = co(j, k) - co(i, k) * a#
NEXT: NEXT: NEXT
REM back subst
FOR i = 6 TO 1 STEP -1
FOR j = 6 TO i + 1 STEP -1
co(i, 7) = co(i, 7) - co(i, j) * hess(j)
NEXT j
hess(i) = co(i, 7) / co(i, i)
NEXT i
RETURN

REM diagonalize 3x3 Hessian matrix
REM characteristic poly
1850 a1# = hess(1): b1# = hess(2): c1# = hess(3)
b2# = hess(4): c2# = hess(5): c3# = hess(6): a2# = b1#: a3# = c1#: b3# = c2#
c# = -(a1# * b2# * c3# + a2# * b3# * c1# + a3# * b1# * c2# - a3# * b2# * c1# - a1# * b3# * c2# - a2# * b1# * c3#)
a# = -(a1# + b2# + c3#)
b# = a1# * b2# + a1# * c3# + b2# * c3# - a3# * c1# - b3# * c2# - a2# * b1#
REM Cardan-Tartaglia, Dummit & Foote pp. 611-614
p# = (3 * b# - a# ^ 2) / 3: q# = (2 * a# ^ 3 - 9 * a# * b# + 27 * c#) / 27
d# = -4 * p# ^ 3 - 27 * q# ^ 2: th# = ATN(-SQR(d# / 27) / q#)
r# = -13.5 * q# / COS(th#): th# = th# / 3: r# = (ABS(r#)) ^ (1 / 3) * r# / ABS(r#)
r2# = r# - 3 * p# / r#
lambda(1) = r2# * COS(th#) / 3: lambda(3) = r2# * COS(th# + 2 * pi# / 3) / 3
lambda(2) = -lambda(1) - lambda(3)
FOR i = 1 TO 3
lambda(i) = lambda(i) - a# / 3
REM PRINT lambda(i),
NEXT
PRINT
REM find char. vectors
FOR i = 1 TO 3
a# = lambda(i)
GOSUB 1880
NEXT
RETURN

REM Cramer's rule for char. vectors
1880 PRINT a#; ": 1, ";
d# = (b2# - a#) * (c3# - a#) - b3# * c2#
f# = (-a2# * (c3# - a#) + a3# * c2#) / d#
g# = ((b2# - a#) * (-a3#) + b3# * a2#) / d#
PRINT f#; ", "; g#
RETURN

1990 END

2000 REM data bloc 2000
REM Neptune positions; days from noon, Jan1,1900
REM Washington Mean Time; ".0" = "noon"; mm/dd/yyyy
REM RA h,m,s, observer corr; Decl deg,',", observer corr
REM B1900.0 coords; p. A CLXIX
REM observer corr fr. pp. A CLXIV,CLXV; unlisted observers get ave. corr
REM equinox corr to RA, fr. p. A CLIV: 1903-1906, -0.038s; 1907-1911, -0.066s
2001 DATA 1382.7,10,14.7,1903,6,25,41.95,-.041,22,14,44.4,-.56
2002 DATA 1518.3,2,27.3,1904,6,13,52.01,.014,22,21,15.9,.06
2003 DATA 1816.5,12,21.5,1904,6,30,7.56,.014,22,14,19.5,.06
2004 DATA 1875.4,2,18.4,1905,6,24,2.55,.014,22,19,40.5,.06
2005 DATA 2171.6,12,11.6,1905,6,41,13.97,-.0185,22,8,22.9,-.14
2006 DATA 2252.3,3,2.3,1906,6,33,13.57,-.0185,22,17,18.3,-.14
2007 REM 2256.3,3,6.3,1906,6,33,5.23,-.0185,22,17,34.3,-.14
REM 2007 is alternative to 2006
2008 DATA 2873.6,11,13.6,1907,7,3,25.29,.021,21,48,15.3,.01
2009 DATA 2978.4,2,26.4,1908,6,52,57.42,-.061,22,4,0.0,-.12
2010 DATA 3263.6,12,7.6,1908,7,11,10.07,-.061,21,39,25.0,-.12
2011 DATA 3319.4,2,1.4,1909,7,4,46.15,.021,21,50,23.6,.01
2012 DATA 3637.6,12,16.6,1909,7,20,2.24,.021,21,27,15.6,.01
2013 DATA 3685.4,2,2.4,1910,7,14,29.53,-.061,21,38,5.7,-.12
2014 DATA 4001.6,12,15.6,1910,7,29,54.23,.021,21,11,11.1,.01
2015 DATA 4052.4,2,4.4,1911,7,24,4.02,-.041,21,23,54.9,-.56

2100 REM data bloc 2100
REM sun RA & Decl apparent for GMT instead of WMT; aberration will compensate
REM sun dist instantaneous for GMT
REM linear interp from Nautical Almanacs
REM 1905 interpolated 365.25 from '04&'06 (missing from ISU lib.)
REM moon approx full at both 1905 dates; so, RA OK, radius corr made
REM hourly rates, from nearest whole 1906 dates ~ equivalent mod 365.25;
REM Feb&Dec hourly rates direct from 1906 Naut. Alm.; others extrapolated
2101 DATA 13,14,49.11,-8,4,44.5,9.9987079,9.291,-55.70,-48.71
2102 DATA 22,37,51.60,-8,38,66.4,9.9958609,9.425,56.21,43.1
2103 DATA 17,58,50.86,-23,26,52.25,9.9927975,11.104,-.29,-11.9
2104 DATA 22,7,29.42,-11,33,24.07,9.99508278,9.632,53.15,40.5
2105 DATA 17,15,20.17,-23,3,3.26,9.99315232,11.020,-11.94,-18.8
2106 DATA 22,51,2.25,-7,19,33.9,9.9962700,9.355,57.10,44.1
2107 REM DATA 23,5,55.95,-5,47,16.5,9.9967047,9.263,58.26,45.3
REM 2107 goes with 2007
2108 DATA 15,12,29.87,-17,54,3.4,9.9953083,10.081,-43.995,-48.95
2109 DATA 22,34,32.50,-8,58,49.7,9.9957800,9.449,55.89,42.8
2110 DATA 16,57,32.80,-22,40,31.7,9.9933602,10.946,-16.44,-21.9
2111 DATA 20,59,34.11,-17,5,2.2,9.9936809,10.174,42.98,26.5
2112 DATA 17,36,9.59,-23,20,22.6,9.9929852,11.082,-6.17,-15.7
2113 DATA 21,2,37.93,-16,52,6.7,9.9937518,10.174,42.98,26.5
2114 DATA 17,30,39.49,-23,16,53.7,9.9930158,11.073,-7.34,-16.4
2115 DATA 21,9,46.36,-16,21,9.7,9.9938729,10.104,44.43,28.2

2600 REM data bloc 2600
REM Jup RA, Decl, &Jup-sun dist from Naut. Alm. with interpolation
REM in my head; accuracy ~0.1 deg heliocentric ecliptic long.;
REM 1905 (lines 4&5) (Naut. Alm. missing) interpolated from neighboring yrs
REM with ~1 mo. different dates --> accuracy ~1 deg helio. ecl. long.

2601 DATA 23,5,15,-7,29,0,0.69609
2602 DATA 23,59,45,-1,13,0,0.69498
2603 DATA 1,17,2,6,42,9,0.69579
2604 DATA 1,54,0,9,0,0,0.6995
2605 DATA 4,0,0,7,0,0,0.7
2606 DATA 3,47,51,19,20,0,0.70412
2607 REM DATA 3,49,56,19,28,0,0.70423
REM 2607 goes with 2007
2608 DATA 9,2.1,0,7,23.4,0,0.722038
2609 DATA 8,31.1,0,19,46.0,0,0.725235
2610 DATA 11,1,41,7,24,0,0.73170
2611 DATA 10,58.5,0,7,59,0,0.73270
2612 DATA 12,45,0,-3,29,0,0.73628
2613 DATA 12,55,44,-4,23,0,0.73650
2614 DATA 14,19.5,0,-12,45,0,0.73577
2615 DATA 14,46,0,-14,46,0,0.73530

2700 REM data bloc 2700
REM same as bloc 2600 but for Sat.

2701 DATA 20,20,42.8,-20,13.4,0,0.9979512
2702 DATA 21,9,52,-17,8.8,0,0.99678
2703 DATA 21,21,30,-16,39,12.3,.99384
2704 DATA 22,7,0,-14,0,0,0.993
2705 DATA 22,0,0,-11,0,0,0.99
2706 DATA 22,33,58.7,-10,42,0,0.988779
2707 REM DATA 22,35,48.8,-10,32,0,0.988729
REM 2707 goes with 2007
2708 DATA 23,29,53.5,-5,48,47,0.980465
2709 DATA 23,53.6,0,-2,58,0,0.979005
2710 DATA 0,16,21.83,-.000000001,57.8,0,0.9750043
2711 DATA 0,26.2,0,0,19.5,0,0.97420
2712 DATA 1,4,23,4,2.25,0,0.96983
2713 DATA 1,10.75,0,4,45,0,0.969189
2714 DATA 1,55.25,0,9,1,0,0.96517
2715 DATA 1,58,0,9,30.6,0,0.964565

2800 REM data bloc 2800
REM same as bloc 2600 but for Ur.

2801 DATA 17,27,3,-23,23,0,1.2842568
2802 DATA 19,57,18,-23,37,10.45,1.2848027
2803 DATA 18,0,22,-23,39,30,1.28598
2804 DATA 19,16,0,-23,34,0,1.286
2805 DATA 18,20,0,-23,55,0,1.29
2806 DATA 18,33,51.7,-23,30,0,1.287696
2807 REM DATA 18,34,26.1,-23,29.5,0,1.287712
REM 2807 goes with 2007
2808 DATA 18,43,44,-23,25,40.7,1.2901014
2809 DATA 19,8.2,0,-22,54,0,1.2905008
2810 DATA 19,6.25,0,-22,59.0,0,1.2915801
2811 DATA 19,20,0,-22,35,0,1.291782
2812 DATA 19,25,0,22,27.3,0,1.292948
2813 DATA 19,37,0,-22,1.5,0,1.29312
2814 DATA 19,42,0,-21,53.3,0,1.29423
2815 DATA 19,54,0,-21,22,0,1.29440

2990 END

REM aberration corr; USNO long. 77.067W
3250 FOR i = 1 TO n
lag# = celest(1, i, 3) / clight# - 77.067 / 15
celest(2, i, 3) = celest(2, i, 3) / 10 ^ (lag# * radiussundot(i) / 10 ^ 7)
time(i) = time(i) - lag# / 24
r2# = celest(2, i, 3)
celest(2, i, 1) = celest(2, i, 1) - (lag# - r2# / clight#) * rasundot(i) / 240 * pi180#
celest(2, i, 2) = celest(2, i, 2) - (lag# - r2# / clight#) * declsundot(i) / 3600 * pi180#
NEXT
RETURN

REM JSU positions
3260 FOR j = 4 TO 6
3270 FOR i = 1 TO n
3280 READ a#, b#, c#, d#, e#, f#, r#
i1 = 1: IF d# < 0 THEN i1 = -1
cc# = (a# * 15 + b# / 4 + c# / 240) * pi180#
aa# = (d# + i1 * e# / 60 + i1 * f# / 3600) * pi180#
celest(j, i, 1) = cc#: celest(j, i, 2) = aa#
r3# = 10 ^ r#
GOSUB 1230
celest(j, i, 3) = r1#
NEXT: NEXT
RETURN

REM time fns
3600 FOR i = 2 TO n
delta(i, 4) = time(i) - time(i - 1): sq(i) = delta(i, 4) ^ 2
NEXT
FOR i = 2 TO n - 1
qdel(i) = delta(i + 1, 4) / delta(i, 4): qdelinv(i) = 1 / qdel(i)
denom(i) = 1 / (delta(i, 4) + delta(i + 1, 4))
diffsq(i) = sq(i + 1) - sq(i): avedel(i) = .5 * (delta(i, 4) + delta(i + 1, 4))
NEXT
denom(1) = 1 / delta(2, 4): denom(n) = 1 / delta(n, 4)
diffsq(1) = sq(2): diffsq(n) = -sq(n)
avedel(1) = .5 * delta(2, 4): avedel(n) = .5 * delta(n, 4)
gap62# = time(6) - time(2): gap146# = time(14) - time(6)
gap1410# = time(14) - time(10): gap102# = time(10) - time(2)
gap142# = time(14) - time(2)
RETURN

3990 END


Update July 3, 2008: Here is the more accurate data bloc incorporating the 1905 Nautical Almanac as well.

2000 REM data bloc 2000
REM Neptune positions; days from noon, Jan1,1900
REM Washington Mean Time; ".0" = "noon"; mm/dd/yyyy
REM RA h,m,s, observer corr; Decl deg,',", observer corr
REM B1900.0 coords; p. A CLXIX
REM observer corr fr. pp. A CLXIV,CLXV; unlisted observers get ave. corr
REM equinox corr to RA, fr. p. A CLIV: 1903-1906, -0.038s; 1907-1911, -0.066s
2001 DATA 1382.7,10,14.7,1903,6,25,41.95,-.041,22,14,44.4,-.56
2002 DATA 1518.3,2,27.3,1904,6,13,52.01,.014,22,21,15.9,.06
2003 DATA 1816.5,12,21.5,1904,6,30,7.56,.014,22,14,19.5,.06
2004 DATA 1875.4,2,18.4,1905,6,24,2.55,.014,22,19,40.5,.06
2005 DATA 2171.6,12,11.6,1905,6,41,13.97,-.0185,22,8,22.9,-.14
2006 DATA 2252.3,3,2.3,1906,6,33,13.57,-.0185,22,17,18.3,-.14
2007 REM 2256.3,3,6.3,1906,6,33,5.23,-.0185,22,17,34.3,-.14
REM 2007 is alternative to 2006
2008 DATA 2873.6,11,13.6,1907,7,3,25.29,.021,21,48,15.3,.01
2009 DATA 2978.4,2,26.4,1908,6,52,57.42,-.061,22,4,0.0,-.12
2010 DATA 3263.6,12,7.6,1908,7,11,10.07,-.061,21,39,25.0,-.12
2011 DATA 3319.4,2,1.4,1909,7,4,46.15,.021,21,50,23.6,.01
2012 DATA 3637.6,12,16.6,1909,7,20,2.24,.021,21,27,15.6,.01
2013 DATA 3685.4,2,2.4,1910,7,14,29.53,-.061,21,38,5.7,-.12
2014 DATA 4001.6,12,15.6,1910,7,29,54.23,.021,21,11,11.1,.01
2015 DATA 4052.4,2,4.4,1911,7,24,4.02,-.041,21,23,54.9,-.56

2100 REM data bloc 2100
REM sun RA & Decl apparent for GMT instead of WMT; aberration will compensate
REM sun dist instantaneous for GMT
REM linear interp from Nautical Almanacs
REM hourly rates, from nearest whole 1906 dates ~ equivalent mod 365.25;
REM Feb&Dec hourly rates direct from 1906 Naut. Alm. (but 1905 exact);
REM others extrapolated from 1906 alm.

2101 DATA 13,14,49.11,-8,4,44.5,9.9987079,9.291,-55.70,-48.71
2102 DATA 22,37,51.60,-8,38,66.4,9.9958609,9.425,56.21,43.1
2103 DATA 17,58,50.86,-23,26,52.25,9.9927975,11.104,-.29,-11.9
2104 DATA 22,6,47.92,-11,37,12.18,9.99505484,9.636,52.96,39.2
2105 DATA 17,13,56.13,-23,1,33.54,9.99317072,11.011,-12.11,-19.6
2106 DATA 22,51,2.25,-7,19,33.9,9.9962700,9.355,57.10,44.1
2107 REM DATA 23,5,55.95,-5,47,16.5,9.9967047,9.263,58.26,45.3
REM 2107 goes with 2007
2108 DATA 15,12,29.87,-17,54,3.4,9.9953083,10.081,-43.995,-48.95
2109 DATA 22,34,32.50,-8,58,49.7,9.9957800,9.449,55.89,42.8
2110 DATA 16,57,32.80,-22,40,31.7,9.9933602,10.946,-16.44,-21.9
2111 DATA 20,59,34.11,-17,5,2.2,9.9936809,10.174,42.98,26.5
2112 DATA 17,36,9.59,-23,20,22.6,9.9929852,11.082,-6.17,-15.7
2113 DATA 21,2,37.93,-16,52,6.7,9.9937518,10.174,42.98,26.5
2114 DATA 17,30,39.49,-23,16,53.7,9.9930158,11.073,-7.34,-16.4
2115 DATA 21,9,46.36,-16,21,9.7,9.9938729,10.104,44.43,28.2

2600 REM data bloc 2600
REM Jup RA, Decl, &Jup-sun dist from Naut. Alm. with interpolation
REM in my head (but 1905 exact);
REM accuracy ~0.1 deg heliocentric ecliptic long.;

2601 DATA 23,5,15,-7,29,0,0.69609
2602 DATA 23,59,45,-1,13,0,0.69498
2603 DATA 1,17,2,6,42,9,0.69579
2604 DATA 1,40,48.04,9,19.17,0,0.6964711
2605 DATA 3,48,25.00,19,1.98,0,0.7020857
2606 DATA 3,47,51,19,20,0,0.70412
2607 REM DATA 3,49,56,19,28,0,0.70423
REM 2607 goes with 2007
2608 DATA 9,2.1,0,7,23.4,0,0.722038
2609 DATA 8,31.1,0,19,46.0,0,0.725235
2610 DATA 11,1,41,7,24,0,0.73170
2611 DATA 10,58.5,0,7,59,0,0.73270
2612 DATA 12,45,0,-3,29,0,0.73628
2613 DATA 12,55,44,-4,23,0,0.73650
2614 DATA 14,19.5,0,-12,45,0,0.73577
2615 DATA 14,46,0,-14,46,0,0.73530

2700 REM data bloc 2700
REM same as bloc 2600 but for Sat.

2701 DATA 20,20,42.8,-20,13.4,0,0.9979512
2702 DATA 21,9,52,-17,8.8,0,0.99678
2703 DATA 21,21,30,-16,39,12.3,.99384
2704 DATA 21,47,36.91,-14,3.308,0,0.99320164
2705 DATA 22,1,37.76,-13,44.916,0,0.98977142
2706 DATA 22,33,58.7,-10,42,0,0.988779
2707 REM DATA 22,35,48.8,-10,32,0,0.988729
REM 2707 goes with 2007
2708 DATA 23,29,53.5,-5,48,47,0.980465
2709 DATA 23,53.6,0,-2,58,0,0.979005
2710 DATA 0,16,21.83,-.000000001,57.8,0,0.9750043
2711 DATA 0,26.2,0,0,19.5,0,0.97420
2712 DATA 1,4,23,4,2.25,0,0.96983
2713 DATA 1,10.75,0,4,45,0,0.969189
2714 DATA 1,55.25,0,9,1,0,0.96517
2715 DATA 1,58,0,9,30.6,0,0.964565

2800 REM data bloc 2800
REM same as bloc 2600 but for Ur.

2801 DATA 17,27,3,-23,23,0,1.2842568
2802 DATA 19,57,18,-23,37,10.45,1.2848027
2803 DATA 18,0,22,-23,39,30,1.28598
2804 DATA 18,14,5.90,-23,37,56.8,1.2862177
2805 DATA 18,15,16.35,-23,39,45.0,1.28738096
2806 DATA 18,33,51.7,-23,30,0,1.287696
2807 REM DATA 18,34,26.1,-23,29.5,0,1.287712
REM 2807 goes with 2007
2808 DATA 18,43,44,-23,25,40.7,1.2901014
2809 DATA 19,8.2,0,-22,54,0,1.2905008
2810 DATA 19,6.25,0,-22,59.0,0,1.2915801
2811 DATA 19,20,0,-22,35,0,1.291782
2812 DATA 19,25,0,22,27.3,0,1.292948
2813 DATA 19,37,0,-22,1.5,0,1.29312
2814 DATA 19,42,0,-21,53.3,0,1.29423
2815 DATA 19,54,0,-21,22,0,1.29440

2990 END

1903-1911 USNO Neptune Observations Support Barbarossa

Standish reported that Uranus positions increased in accuracy in about 1911. Lowell used Uranus data from until about 1914. Thus these 1903-1911 USNO Neptune data, which I've found in the Iowa State Univ. Library, have the advantage of being no more accurate than Lowell's Uranus data. So, if Barbarossa is implied by these data, then Lowell might well have correctly deduced Barbarossa (or rather, the roughly gravitationally equivalent, smaller but nearer, Planet X) from his Uranus data.

These data are from original reports, not summaries possibly degraded by omitting details. The data are homogeneous in their method of collection and presentation.

So far, I've used only the first and last USNO observations of each opposition season; the greatest Earth parallax gives the most accurate triangulation of position. I used March 2.3, 1906, instead of March 6.3, because the other last-of-season oppositions were in Feb. Other things being equal, I planned to favor observations by observers whose personal corrections were given, but as I recall, this didn't become an issue. There were no observations from the 1906-1907 season, so late 1903 - early 1911 encompassed (8-1)*2=14 first- & last-of-season observations.

The data table had small footnotes telling the user to make observer and equinox corrections, which I looked up in the lengthy introduction. I wonder if these were on the magnetic tape from which Standish says he got his USNO Neptune, etc., data.

The 1905 ephemeris is missing from the ISU library, so my planetary positions (other than Earth, which easily could be interpolated, & Neptune) for that year are especially rough. This book is still on interlibrary loan order; when it comes, I can improve my 1905 data.

My program assumes that the center of mass of the known solar system w.r.t. Barbarossa, is unaffected by the variations in Neptune's trajectory. This makes the potential term tidal, i.e., quadratic in the sun-to-Neptune (more precisely, rest-of-known-solar-system-c.o.m.-to-Neptune) distance, and described by a Hessian quadratic form. The six independent entries of the Hessian matrix are scaled so that the perfect Barbarossa effect (in the limit where the ratio of Neptune's to Barbarossa's orbital radius, is zero) is a "1" and two "-1/2"s, on the diagonal (of the eventually diagonalized Hessian), and all other entries "0". These entries were optimized by a 6x6 system of linear equations, to minimize the sum of squares of derivatives (more precisely, central first differences) of the minimum Hamilton's principle action integral, w.r.t. alterations of the 14x2=28 observed celestial coordinates. (This minimum action integral is minimized w.r.t. variations of the Neptune-sun radius at two midpoints, with cubic or at least quadratic interpolation; in Neptunian ecliptic coordinates so second time derivatives, of coordinates, are small.)

The implied average rounding error of the USNO Neptune Declinations is 0.025", of the RAs 0.0375", and of the Nautical Almanac Neptune-sun distances, 0.25/10^7 Briggs log unit which is to the radius as 0.012" is to a radian, i.e., equivalent to 0.012". (The error due to Triton's displacement of Neptune is at most 0.01".)

I tried varying two midcourse radii (for action integral minimization) and observed coordinates (for finding which extraneous tidal force minimized the constraint they imposed) in step sizes ranging from 0.1" to 2", or the equivalent (and twice that, for the midcourse radii). The 1.2" and 2" step sizes, seemingly implied a much bigger extraneous tidal force than did the smaller 0.1", 0.1778" = 1"/10^0.75, 0.3162"=1"/10^0.5, 0.5623"=1"/10^0.25, 1.0", 1.05", 1.1", 1.15", or 1.17" step sizes.

Standish (p. 2001, 1st col.) says Uranus observations before 1911 have 1.2" "scatter"; presumably the scatter of Neptune observations is the same. This is exactly the step size, at and above which, my variational method diverges to unbelievably big characteristic values.

The step size which minimized the magnitude of the matrix trace (ideally +1-1/2-1/2=0), was 1.05", which gave characterisitic values +0.4918, +0.0076, and -0.3821. This and all smaller step sizes, yielded roughly the same characteristic vector, for the largest-magnitude (which also was the most positive) characteristic value. This vector, with a rough correction for Neptune's displacement from the sun, corresponds to RA 210.9 Decl -36.6 (modulo 180). With 1.1" and larger step sizes, the characteristic vector began to fluctuate much more.

[Update June 30, 2008: I got the missing 1905 Nautical Almanac on interlibrary loan from BYU via the Bertha Bartlett library in Story City. The improved data are appended to the USNO.BAS program above. Rerunning the above USNO.BAS program after correcting the sun, Jupiter, Saturn & Uranus positions for 1905, I find, for the stepsizes 0.1", 1"/10^0.75, 1"/10^0.5, 1"/10^0.25, 1.0", 1.05", 1.1", 1.15", & 1.2" (a subset of those above) that again the magnitude of the trace is minimized by the 1.05" stepsize. Again, the magnitude of the characteristic values for the 1.2" stepsize, is much bigger than for smaller stepsizes; and somewhat as before, the direction of the characteristic vector corresponding to Barbarossa's position, fluctuates more for stepsize 1.15" and above. Stepsize 1.05" gave characteristic values +0.5632, +0.0256, and -0.5572. The characteristic direction corresponding to the large (+) characteristic value, is (with the same rough correction for Neptune's orbital radius) RA 213.5 Decl -40.8 (modulo 180).]

[Update July 2: The result oscillates rapidly as a function of stepsize. Stepsize 1.050001" gives trace +0.9, bigger than the usual trace range +0.3 to +0.5; 1.051" gives very big characteristic values; 1.049" gives a small trace like 1.050", but the trace doesn't interpolate, even between 1.049" & 1.050". Generally however the characteristic vector corresponding to the large (+) characteristic value, almost always lies in the same octant (modulo 180), and the characteristic vector obtained with stepsize 1.04", RA 236.2 Decl -34.6, differs only about 17 deg (mod 180) from that obtained with stepsize 1.05". When variation of the Neptune radii is omitted altogether, 3 of 3 stepsizes tried, give very big characteristic values.]

[Update July 3: With the 1.05" stepsize, the other 8 endpoint radii choices in the program, i.e., ~(-1, 0, or +1)/200000 of the contemporary ephemeris Neptune radius, failed to give a trace any closer to zero. Despite double precision, this program suffers greatly from rounding error and numerical instability, but it is nevertheless remarkable that the stepsize runs giving tidal force Hessians in best agreement with Laplace's equation (i.e., zero trace) also point toward the direction given by Lowell and Harrington.]

Harrington's best 1988 and 1930 positions (see my post above, discussing Harrington's agreement with Lowell), extrapolated linearly to 1910.5 (the midpoint of Harrington's Uranus data, on which he relied) give RA 206, Decl -25. The midpoint of my 1903-1911 data is actually about 1907.5, so +3.0 years of Barbarossa's presumed orbital motion would change my [updated] best direction to RA 213.8, Decl -41.0 for 1910.5. Thus Harrington's Planet X direction (relying on Uranus data) essentially differs only 17 deg (p=0.044, for modulo 180) from the direction I determine variationally from 14 USNO Neptune observations 1903-1911.