Joe, if I'm following you correctly you're saying a planet of ~3.7 Earth masses at 197.7 AU would fit the Uranus perturbation data. A planet of that mass at that distance would be cold and very dim and would require no masking nebula. It could be the "green dot". Wouldn't it be more parsimonious to make that assumption?

This lower mass (my 3.7 became 2.4 Earth masses upon improved measurement of the graph, and revision of my simple dynamical theory above) applies to a planet at ~ 52AU, not 197.7AU. My idea now is, that the phase incoherence, in the signal remaining after white noise subtraction (in Standish's graphs, Figs. 5a & 6a, of Uranus' RA) is inconsistent with Newtonian perturbation by a planet (see previous post; this idea was appended to it, after "nemesis" posted his question).

The 2.4 Earth masses might become zero if Standish's space probe data were perfected. The ephemeris based on Earth-based planetary mass measurements (the basis of Fig. 5a) gives 7.45 Earth masses: almost exactly half the mass of Uranus itself. It's as if Uranus were interacting with itself, through the intermediation of some structure or barrier that orbits at ~ 52AU.

Though my periodograms don't confirm the evidence for Barbarossa that I saw grossly in Standish's Fig. 5a, Standish's analysis doesn't disprove Barbarossa's existence. If Standish's ephemeris is right, his predecessors' ephemerides were wrong; if Standish's predecessors were fallible, then Standish is, too. Also, the ~ 52AU, ~7.45 (~14.6/2!) Earthmass *incoherent* signal in Standish's Fig. 5a, which remains after detrending and de-noising, but is mostly or completely removed by use of space probe planetary mass measurements as in Fig. 6a, indicates new physics.

With new physics, the only way to prove anything is to look, not theorize. My $100 reward offer, previously detailed on this messageboard, still stands. The reward is merely for getting someone really sufficiently equipped, to look definitively; not for looking oneself. Previously on this thread I've listed many reasons for looking. If a dense body isn't there, then likely something just as interesting is.

Regarding *Neptune's* ephemeris, Standish's Fig. 8a remarks that the ephemeris of USNO Publications, Vol. XII, agrees with Standish's, for Neptune's RA, only between 1830 & 1950. Before 1830 & again after 1950, for whatever reasons, the USNO Vol. XII ephemeris gives +20"/century greater RA for Neptune than Standish's does. Pursuant to my recent posts, Barbarossa should cause *Uranus'* time derivative of RA to deviate +/- 0.73"/(84/2/4)*pi/2*100 = +/- 10.92"/century. Tidal acceleration goes as r, time before reversal goes as r^1.5, and subtended angle goes as 1/r. So, Barbarossa should cause *Neptune's* time derivative of RA to deviate +/- 10.92 * (30.07/19.18)^1.5 = +/- 21.44"/century. This is the same deviation of Neptune's time derivative of RA, that by whatever complicated calculation, is seen in USNO Vol. XII vs. Standish.

Sorry, I misunderstood... but wouldn't it be worthwhile to look for a planet of ~4 Earth masses at ~52 AU? It would still be very dim and could have been overlooked. And it shouldn't rule out Barbarossa.

Sorry, I misunderstood... but wouldn't it be worthwhile to look for a planet of ~4 Earth masses at ~52 AU? It would still be very dim and could have been overlooked. And it shouldn't rule out Barbarossa.

Thanks for emphasizing this. Barbarossa would be the dimmer, because (198/52)^4 = 250 = 6 magnitudes, not enough to compensate for the size difference between a Uranus (i.e., 5-15 Earthmass planet, even if of Earthlike density) and a Jupiter (i.e., quantum mechanical estimate of brown dwarf size). Unfortunately, the phase incoherence of the residual signal I find in Standish's data, leaves the phase indeterminate; furthermore it argues against a planetary cause at all.

The USNO Vol. XII (1929, according to the USNO library webpage) ephemeris for Neptune's RA (according to Standish, Fig. 8a) is consistent with Barbarossa in amplitude, period and phase. The USNO mission is to produce accurate ephemerides by any and all lawful means. Before radio, a cruiser might lose all hands, if the navigator lacked accurate information on the only 3rd mag. star he could identify on a partly cloudy night. Tomorrow, after electronic countermeasures, and the electromagnetic side effects of superweapons such as nuclear bombs, the navigator may well revert to the sextant. I doubt he'd sight Neptune, but if the equivalent of a thousand taxpayer working lives are wasted because the Navy's rocket to Neptune missed, the loss is, in a way, comparable to that of a cruiser.

Fig. 8a (Neptune) resembles the error in fitting 1.5 cycles (trough-peak-trough-peak) of an 80-yr sinusoid, to a cubic curve. For the 120 yr between 1830 & 1950 (Neptune was discovered in 1846, but Neptune's ephemeris can be extrapolated backward) the fit is good, but before and after that, the ephemeris for Neptune's RA deviates roughly as a cubic curve from a sinusoid. The USNO might have tried to adjust for some small unexplained error, by zeroing the discrepant RA, with a cubic curve, at 1930 (approx. date of the catalog), 1850 (approx. beginning of Neptune data) & 1890 (midpoint)(the three roots of the sinusoid). This attempt would begin to fail markedly a quarter cycle away, before 1830 and after 1950, as Fig. 8a shows.

James DeMeo had to fly to Cleveland to rummage through the Physics Dept. at Case Western Reserve Univ. (after tactfully getting permission) to rediscover the only extant copy of Dayton Miller's data, in a cardboard box in a dusty storeroom. Everyone involved with USNO Pubs. Vol. XII is deceased. If notes exist detailing all their corrections, finding them might be like rediscovering something in the Vatican Library. The USNO mission was and is, to produce accurate ephemerides. They understood that there are systematic observation errors, calculation errors, undiscovered planets, unknown forces and unimagined "new physics". A cubic curve might be used to correct a small, incomprehensible discrepancy.

The slope of the discrepancy (after inclusion of such a cubic correction term) in Fig. 8a, at 1830 and 1950 (a sinusoid trough & peak, resp.) would approximate the maximum sinusoid slope due to Barbarossa, i.e., +21.44"/century (see my earlier post).

Not only are the period and amplitude of the unexplained sinusoid (the 1929 USNO approximation of which by a cubic, results in Standish's Fig. 8a) consistent with Barbarossa; so is the phase. One of the times at which the sinusoid should have most negative slope, is c. (1830+1950)/2 = 1890. In 1890.0, I find Neptune's heliocentric ecliptic longitude was 63.5. Barbarossa's effective heliocentric ecliptic longitude was: 173.327 (for 1987.082) - (1987.082-1890.0)/2780*360 (for orbital motion) - 175/4/2780*360 (for retardation of speed change, by a Neptune-minus-Barbarossa quarter cycle vis a vis the tidal acceleration due to Barbarossa) = 155.1. Recalling the 2nd harmonic nature of the tidal force, Neptune's RA should have most positive acceleration at 155.1-45, most negative acceleration at 155.1-135, and zero acceleration, i.e., most negative speed, at 155.1-90=65.1, agreeing perfectly, considering the rough dates (1.6deg < 1yr), with the actual 63.5.

There are many reasons why Standish might find an ephemeris which excludes Barbarossa. However, the 1929 USNO Vol. XII ephemeris, supports Barbarossa's existence, indirectly but accurately confirming the amplitude, period and phase of Barbarossa's predicted tidal effect on Neptune.

(sent three minutes ago via JPL website email form)

Attn.: Dr. P. C. Liewer, Jet Propulsion Laboratory (pls. forward)

To: Dr. E. Myles Standish (if active) or his successor (otherwise)

There is an interesting comment that has been made regarding Dr. E. Myles Standish's ephemeris. Apparently this comment only has been published on Dr. Tom Van Flandern's online messageboard (attached to Dr. Van Flandern's website, www.metaresearch.org) by Joe Keller, date today, May 15, 2008.

Sincerely, Joseph C. Keller B.A., cumlaude, Harvard, Mathematics, 1977

Update May 31, 2008: no response of any kind so far, from anyone.

Maran et al, Planetary & Space Science 45:1037-1043, 1997, solving differential equations of motion numerically using large amounts of computer time, found the nearest distance at which a Neptune-mass Planet X, could orbit without destabilizing (causing large increases in major axis within 10^6 yr) known Trans-Neptunian Objects (TNOs). For a circular orbit of i=0 or 5deg inclination, the Planet X distance for which all 14 studied TNOs become stable, is ~62.5AU according to Maran, Table 6, p. 1043. For i=10deg, it's ~58.75AU, and for i=45deg, ~54AU. Interpolating to Barbarossa's i=12deg10' (curiously, the same inclination as "Vulcan") gives 57.125AU.

Maran studied real TNOs with nonzero eccentricity & inclination. To compare Maran's Planet X to Barbarossa, in its effect on a TNO that ranks in the 93rd percentile (i.e., 1 in 14) for vulnerability to destabilization, let's consider a hypothetical TNO with zero eccentricity & inclination, but semimajor axis 45AU instead of the usual 42-43AU for a classical TNO (or 39.44AU for a plutino). As in posts above, the tidal effect is estimated by multiplying the tangential tidal force at quadrature, by the sine of the angle subtended at the sun at quadrature:

M/R^2*(sec^2(theta)-cos(theta))*cos(theta) where sin(theta)=r/R

If this effect is same for Barbarossa as for Maran's hypothetical Planet X (i.e., just enough to destabilize a 93rd percentile, 1 in 14, vulnerable TNO) then Barbarossa must have 3253 Earth masses, agreeing with the 3430 Earth masses I've attributed to Barbarossa by assuming orbital precession resonances in the outer solar system.

Summarizing, Maran's 1997 computation-heavy study of 14 TNOs, found the distance at which a Neptune-mass Planet X in circular orbit, at Barbarossa's inclination, would cause the stability limits of the Edgeworth-Kuiper belt, to be what they are. Using an elementary geometric estimate, I found that a Barbarossa-mass Planet X, in circular orbit at Barbarossa's distance and inclination, also would produce the observed limits of the Edgeworth-Kuiper belt.

(sent three minutes ago to email address listed on univ. website)

Attn.: Prof. J. P. Emerson, Queen Mary/University of London (please forward)

To: Dr. M. D. Maran

Yesterday I discovered your article, "Limitations on the Existence of a Tenth Planet". Your heavy computation found one point on the mass-radius curve constraining any possible Planet X. The actual Planet X, for which there is much evidence, lies elsewhere on this curve.

I've put this evidence on the messageboard of Dr. Tom Van Flandern's website, www.metaresearch.org, under my name, Joe Keller. Today I added a discussion of your article.

Sincerely, Joseph C. Keller, M. D. B. A., Harvard, cumlaude, Mathematics, 1977

Update May 31, 2008: no response of any kind so far, from anyone.

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.

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.

"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.8).)

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:

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.

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 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

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

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

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.

Posted (in essence) here May 15, 2008, by Joe Keller:

... Fig. 8a (Neptune) [Standish, Astronomical Journal, 1993] resembles the error in fitting 1.5 cycles (trough-peak-trough-peak) of an 80-yr sinusoid, to a cubic curve. For the 120 yr between 1830 & 1950 (Neptune was discovered in 1846, but Neptune's ephemeris can be extrapolated backward) the fit is good, but before and after that, the ephemeris for Neptune's RA deviates roughly as a cubic curve from a sinusoid. The USNO might have tried to adjust for some small unexplained error, by zeroing the discrepant RA, with a cubic curve, at 1930 (approx. date of the catalog), 1850 (approx. beginning of Neptune data) & 1890 (midpoint)(the three roots of the sinusoid). This attempt would begin to fail markedly a quarter cycle away, before 1830 and after 1950, as Fig. 8a shows. ...

Today's comment:

The Vol. XII USNO Neptune ephemeris, as plotted by Standish, has a nonsinusoidal error. Real orbits are sums of sinusoidal functions (in a way, the epicycles of Ptolemy have evolved into Fourier terms). The Vol. XII error, on the other hand, resembles that due to use of a cubic polynomial term, and therefore must be the work of man. Polynomials often are used to enhance the accuracy of slowly varying constants such as precession, so it's likely that the (basically) cubic polynomial error in USNO Vol. XII, is due to the limited range of accuracy of a polynomial used to describe empirically, 1.5 cycles of an unexplained sinusoidal residual. This is appropriate despite its limitations, because the USNO's primary mission is accuracy by any means.

Newcomb's 1866 ephemeris, applied to 1866-1889 USNO Neptune data, apparently uses an extra second order harmonic to describe Barbarossa's tidal effect. The small residual arises from the difference between the period, 0.5*164.8=82.4 yr, of Newcomb's second order harmonic term, and the actual period of Barbarossa's tidal effect, 0.5/(1/164.8-1/2780)=87.59 yr. Yesterday and today I composed and ran the following program:

REM program name: USNONEWC.NEP REM lang. QBasic REM This program finds the (linear & quadratic) time & opposition dependence REM of USNO 1866-1889 Neptune residuals vs. Newcomb's 1866 ephemeris.

PRINT : PRINT pi = 4 * ATN(1): pi180 = pi / 180: n = 0: nn = 0 DIM sunra(12): PRINT "Read sun RA" REM sunra is RA of sun in radians, on 15th of mo. in 1880, per Nautical Alm. FOR i = 1 TO 12 READ a, b REM datacheck IF a < 0 OR a > 23 OR b < 0 OR b > 59 THEN GOSUB 3000 sunra(i) = pi180 * (a * 15 + b / 4) NEXT DIM nep(100, 4): DIM x(100): DIM y(100): PRINT "Read Neptune data" REM nep's components are: time from 1866.0 in yrs; no. observations for mo.; REM residual vs. Newcomb's ephemeris in " (ave. for mo.); REM RA of Nept. in radians, at 1st observation of mo.; REM RA of sun in rad, 15th of mo., per 1880 Naut. Alm. FOR i = 1 TO 100 READ a, b, c, d, e IF a = -1 GOTO 100 REM datacheck IF a < 0 OR a > 23 OR b < 1 OR b > 12 OR d < -200 OR d > 400 THEN GOSUB 3000 IF c < 1 OR c > 15 OR e < 0 OR e > 240 THEN GOSUB 3000 n = n + 1: nn = nn + c nep(i, 1) = a + (b - .5) / 12 nep(i, 2) = c nep(i, 3) = d / 1000 * 15 nep(i, 4) = e / 4 * pi180 + pi - sunra(b)

REM adjusts for independent effect of distance from opposition REM nep(i, 3) = nep(i, 3) + .285 * nep(i, 4) ^ 2

NEXT

REM find linear & quad effect, of distance from opposition 100 PRINT "n = "; n; "; nn = "; nn FOR i = 1 TO n x(i) = nep(i, 4): y(i) = nep(i, 3) NEXT PRINT : PRINT "Effect of distance from opposition" GOSUB 1000

REM find lin & quad effect of time 200 FOR i = 1 TO n x(i) = nep(i, 1): y(i) = nep(i, 3) NEXT PRINT : PRINT "Effect of time" GOSUB 1000

REM find dependence of distance from opposition, on year of observation 300 FOR i = 1 TO n x(i) = nep(i, 1) y(i) = nep(i, 4) REM for dependence of sq. of dist. fr. opp., on yr, use line below: REM y(i) = y(i)^ 2 NEXT PRINT : PRINT "Dependence of dist. from opposition, on time" GOSUB 1000

END

REM finds corr. coeff. & best fit slope, vs. x & vs. x^2 1000 PRINT "Corr. coeff. & slope vs. x: "; GOSUB 1100 FOR i = 1 TO n y(i) = y(i) - slope * (x(i) - mx): x(i) = (x(i) - mx) ^ 2 NEXT PRINT "Detrended corr. coeff. & slope vs. (x-mx)^2: "; GOSUB 1100 RETURN

REM calc corr coeff & slope 1100 sx = 0: sy = 0: u = 0: v = 0: w = 0: m = 0 FOR i = 1 TO n wt = nep(i, 2): sx = sx + x(i) * wt: sy = sy + y(i) * wt: m = m + wt NEXT mx = sx / m: my = sy / m FOR i = 1 TO n wt = nep(i, 2): w = w + (y(i) - my) * (x(i) - mx) * wt u = u + (x(i) - mx) ^ 2 * wt: v = v + (y(i) - my) ^ 2 * wt NEXT cc = w / SQR(u * v): slope = cc * SQR(v / u) PRINT "corr. coeff "; cc; " slope "; slope; PRINT " mean x "; mx; " mean y "; my; PRINT " Fisher's z = "; SQR(n - 3) * LOG((1 + cc) / (1 - cc)) RETURN

END

REM sun RA at Greenwich mean noon on 15th of month per 1880 Nautical Alm. REM hr & min; nearest minute; zero signifies month not used 2000 DATA 19,47,0,0,0,0,0,0,0,0,0,0,7,41,9,41,11,35,13,24,15,25,17,34

REM from House of Reps. Misc. Documents v. 101, pp. B151-157; REM USNO Neptune positions 1866-1889; REM see main program for more explanation 2101 DATA 0,7,2,85,50 DATA 0,8,12,116,49 DATA 0,9,10,91,47 DATA 0,10,13,88,45 DATA 0,11,10,56,42 DATA 0,12,10,109,40 DATA 1,1,1,90,40 DATA 1,7,2,60,58 DATA 1,8,6,16.5,58 2110 DATA 1,9,10,22,56 DATA 1,10,14,36,53 DATA 1,11,4,67.5,49 DATA 1,12,7,19,48 DATA 2,8,4,11,66 DATA 2,9,7,70,64 DATA 2,10,9,64,62 DATA 2,11,12,47.5,59 DATA 2,12,6,65,57 DATA 3,1,2,80,56 2120 DATA 5,9,3,40,91 DATA 5,10,2,25,87 DATA 5,11,1,-30,83 DATA 5,12,2,-15,82 DATA 6,9,1,80,98 DATA 6,11,8,-101,93 DATA 7,10,1,20,105 DATA 7,11,8,54,101 DATA 7,12,1,-10,98 DATA 8,10,6,118,113 2130 DATA 8,11,7,84,111 DATA 8,12,9,37,108 DATA 9,1,2,115,107 DATA 9,9,5,158,125 DATA 9,10,5,90,123 DATA 9,11,11,125,120 DATA 9,12,8,182.5,117 DATA 10,9,2,95,134 DATA 10,10,8,124,131 DATA 10,11,4,172.5,128.5 2140 DATA 10,12,7,134,125 DATA 11,10,4,115,139 DATA 11,11,7,149,137 DATA 11,12,4,137.5,134 DATA 12,9,1,80,151 DATA 12,11,4,152.5,146 DATA 13,9,2,165,160 DATA 13,10,2,190,157 DATA 13,11,5,174,154 DATA 13,12,1,200,151 2150 DATA 14,1,2,170,150 DATA 14,11,1,260,162 DATA 14,12,8,287.5,161 DATA 15,10,3,177,175 DATA 15,11,4,267.5,172 DATA 15,12,8,190,170 DATA 16,1,1,200,167 DATA 16,11,1,240,181 DATA 17,9,1,260,195 DATA 17,11,5,260,192 2160 DATA 17,12,7,181,189 DATA 18,11,6,165,200 DATA 18,12,4,172.5,198 DATA 19,1,2,230,195 DATA 19,11,3,190,209 DATA 19,12,7,221,207 DATA 20,1,1,300,204 DATA 20,11,3,117,218 DATA 20,12,6,122,216 DATA 21,1,1,230,214 2170 DATA 21,11,7,221,229 DATA 21,12,6,248,226 DATA 22,1,5,248,222 DATA 22,12,4,280,235 DATA 23,1,3,280,232

DATA -1,0,0,0,0

END

REM datacheck 3000 PRINT "! incongruous data value at "; i RETURN

END

Trying to read Newcomb's mind, I think he would have fitted the data with, essentially, an empirical series of sinusoids, all harmonics of the presumed Neptune period. Surely Newcomb would have admitted an ad hoc second order harmonic term, because of the possibility of yet another remaining distant outer planet, as LeVerrier already in 1846 publicly had suspected and as Todd in 1877 would predict and seek. Because of resonances, small undiscovered planets nearer than Neptune likely would give higher harmonics too. Newcomb's series wouldn't exactly equal the theoretical Newtonian perturbed elliptical orbit. I think Newcomb would have published the orbital elements of the ellipse which most resembled his series of sinusoids, while using the empirical series itself to calculate positions.

Above, I've estimated that Barbarossa's effect on Neptune's longitude would have period 87.59 yr (so, it isn't exactly a second harmonic in period), would cross the time axis downward in 1890, and would have (semi)amplitude 0.2144"/87.5/(2*pi)=5.623". Conveniently, 1795 (Lalande's unwitting but usefully accurate Neptune position; see Rawlins, Astronomical Journal 75:856+, 1970, "The Great Unexplained Residual in the Orbit of Neptune") and 1853 (near the midpoint of the 1846 to, at latest, 1866, post-discovery data available for Newcomb's 1866 ephemeris) lie on the same level of this actual sinusoid. A constant term, and the amplitude and phase of the exact second harmonic sinusoidal correction term (implicitly) used, are three parameters which can be fitted to Lalande's point, to the 1853 point (i.e., approximate midpoint of the post-discovery data) and to the 1853 time derivative of the actual (not exactly second harmonic) sinusoid.

I did so. This is about the best fit Newcomb could achieve with the harmonics he likely used. Third and higher harmonics help little because Newcomb's base and prediction epochs are 1/8 cycle apart; for the third harmonic this is 3/8 cycle, so the values of its derivatives at the base epoch bear little relation to their values at the prediction epoch.

Using the above program, I found the actual residuals. The best-fit (using correlation coefficients) quadratic polynomial describing the residuals, is

Season-minus-opposition (i.e., signed angular distance of the observation, from opposition) also correlated with the residuals, but only through its correlation with time. On the other hand, the squared distance from opposition, did correlate, at -1.9 sigma, with the residuals, independently of time. The best-fit quadratic polynomial including this (season-minus-opposition)^2 effect, is

The r.m.s. distance from opposition was 0.619 radian. Using this r.m.s. value to get corrected, smoothed observations, I find for residuals from Newcomb's presumed best-fit exact second harmonic (discussed above)

1866.0 +0.724" 1877.5 +1.927" 1889.0 +3.976"

(The residuals according to the fit which ignored the distance of the observations from opposition, hardly differed from these.)

For comparison, the predicted residuals (assuming, as above, Barbarossa's tidal effect, and Newcomb's harmonic term best fitting it) are

1866.0 +0.552" 1877.5 +1.627" 1889.0 +2.400"

Integrating by Simpson's rule, we find that this theory, of a difference between Newcomb's exact second harmonic, and the actual slightly longer second harmonic due to Barbarossa, explains 69% of the magnitude of the area under the residual curve, which is to say, 90% of the residual variance (excluding short-term variations).

Similarly, Eckert, Brouwer & Clemence fit Neptune's orbit. Their method was to find the closest fitting Newtonian orbit. They avoided empirical correction terms. Their purpose was to determine the discrepancy implied by Lalande's 1795 postion. They were adjusting the Newtonian orbit, in effect adding a first-order harmonic only (the higher order harmonics would be small due to Neptune's small eccentricity). I found that their procedure, given Barbarossa's tidal force, should imply exactly the Lalande residual that it did. I'll write about that tomorrow.

Barbarossa's tidal effect on Neptune's longitude, spans 5/4 cycle of a sine wave from 1846 to 1956 (zeros at 1846, 1890, 1934; peaks 1868 & 1956; trough 1912). Eckert et al published in 1951, but for a simple approximation, let's call that 1956. Eckert's purpose wasn't accurate prediction; his purpose was to find discrepancy, from that Neptune orbit which would be implied by known solar system bodies. So, Eckert mainly was adjusting position, period, eccentricity and perihelion, for best fit to observed longitude. With Neptune's small eccentricity, this amounts to adding some small constant, plus a small first order harmonic, to the function theta(time).

Such a fit would use the bottom half of a first order sine wave, intersecting the second order curve of Barbarossa's perturbation, in 1868, 1912, & 1956, and tangent to it in 1912. The maximum discrepancy would be only half the (semi)amplitude of Barbarossa's perturbation of Neptune. A small constant term could halve this again, to only 5.6"/4 = 1.4" max. The fit would worsen < 1868, but these data are sparser and less accurate anyway; in Parks Library at ISU, I found no USNO Neptune observations earlier than 1866, in any of the House or Senate documents.

Barbarossa's perturbation of Neptune has period ~ 88 yr, but really, Neptune's period is only 164.8, not 88*2. For an approximation, let's fix the 1912 point of the first order fitting sinusoid. The fitting sinusoid has not only about twice the period but also about twice the amplitude of the actual Barbarossa perturbation. Thus for the hypothetical case of no Barbarossa at all, the predicted 1795 (Lalande) residual is cos(360*(1912-1795)/164.8) * 2 * 5.62" - 1.405"(from the small constant term giving least max residual 1868-1956) - 5.62" (because the midline of the first order sinusoid is at the peak of the Barbarossa sinusoid) = -9.82". For the actual Barbarossa case, 1795's longitude is perturbed sin(360*(1890-1795)/87.59) * 5.62" = +2.85".

So in this approximation Eckert should actually find a longitude residual of

-9.82"+2.85" = -6.97"

and if the main approximation error is corrected, by raising the first order sinusoid 0.2144"/yr * (1956-1951) * 0.5 (because I've already split the difference to minimize the max error), we get

-6.97" - 0.54" = -7.51"

for 1795. Eckert did find (Rawlins, Astronomical Journal 75:856+, 1970, Table II; citing Eckert, Brouwer & Clemence, Astron. Papers Am. Ephemeris 12, 1951) -6.4", assuming a negligible mass (e.g., the currently accepted 0.002 Earth mass) for Pluto. Rawlins thought Eckert's figure should be adjusted to -8.6". So the effect of Barbarossa, after partial removal by best-fit orbital adjustment, equals the Lalande residual, as estimated by Eckert et al and by Rawlins.

(This is the first email letter I've sent about this subject, that got a response from a professional astronomer besides Dr. Van Flandern. The responder is a department chair, too.)

June 15, 2008

Dear Prof. *******,

I'm a Harvard math major who has compiled much evidence regarding "Planet X", which I've discovered and named Barbarossa. I've written a computer program which analyzes Neptune positions for tidal influence; the result is that the direction and tide due to Barbarossa are roughly what Lowell thought they were. I've analyzed residuals from Eckert's (1951), the USNO Vol. XII (1929) and Newcomb's (1866) ephemerides, and found that all of them seem to be due to one and the same Barbarossa, minus whatever ad hoc fitting corrections were used in the ephemeris.

The direction and tidal strength of Barbarossa are about the same as estimated by Todd, Lowell, or Harrington. I also calculate that Standish's correction to Neptune's mass is far too small to make such residuals go away, so, I don't know why Standish doesn't find them. Maybe Standish made the graduate student do it over and over until it came out zero.

Maran's numerical experiment on the stability of Trans-Neptunian Objects, also supports a Barbarossa of this tidal strength, as do simple outer solar system orbital precession resonances I discovered myself. Articles purporting to disprove "Planet X" (Barbarossa), have big holes in their arguments, and often explicitly acknowledge that they do.

The unexplained correlation of the multipoles of the "cosmic" microwave background, with the ecliptic, caused me, motivated by a novel physical theory of the CMB, to search the USNO-B catalog online, for excess inconsistent magnitudes near the (+) CMB dipole. These might be due to automated misidentification of Barbarossa and its moons, and stars, with one another. Later I realized that there is much evidence of a moving (i.e., intra-solar system) nebula in this direction. Also IRAS might have detected Barbarossa's bow shock.

I found disappearing dots on sky survey plates yielding an accurate and believable orbit. Amateur astronomers Joan Genebriera (Spain; 16", Tenerife), Steve Riley (U.S.A.; 8", S. California) and Robert Turner (England; automated 14", Tenerife) aiming at my coordinates, prospectively captured electronic images, presumably of Barbarossa or its moons, consistent with this ephemeris. Most of these images aren't starlike, but there might be novel reasons for that.

Working backwards from my own, through Tombaugh's, now through Adams or LeVerrier's methods, I'd be ready for Herschel's method of simply looking through a telescope, except that I can't afford the time or money (many say it's easy but few do it). The Barbarossa system, Barbarossa together with its moons (which include one identified and two inferred giant planets in their own right), seems to have 11 Jupiter masses. Barbarossa itself, with ~ 8 Jupiter masses, is near the theoretical boundary between cold brown dwarf and giant planet. At equilibrium temperature (at est. 197.7 AU) it could escape IRAS detection. The apparent magnitude is only +18 or +19; this is consistent with somewhat smaller than quantum-theoretical size, albedo in the low theoretical range (like the alkali metal "black smoker" brown dwarf) and/or an obscuring nebula.

Maybe your observatory would like to look for Barbarossa. Coordinates and other information have been posted by me on the messageboard of Dr. Tom Van Flandern, at www.metaresearch.org.

For more than a year, I've contacted hundreds of persons about this. The response of messageboards, except for Dr. Van Flandern's, has been either total silence, or, usually, one form or another of expulsion. The response of professional astronomers (except for Dr. Van Flandern, who does not agree with all my ideas, but who has suggested books and articles, and has formulated some hypothetical questions) has been silence. I've written (mail, not email) government officials and Congressmen; the only response there, has been a postcard saying, basically, I'd like to help you son but you don't vote in my district.

Sincerely, Joseph C. Keller, M. D. B. A., Harvard, cumlaude, 1977

Above, I've calculated the orbit of Barbarossa assuming its mass is negligible compared to the sun's. The actual 0.0103::1 mass ratio alters the calculated Barbarossa-sun distance from 197.7 to 198.4 A.U. However, this alteration affects Barbarossa's coordinates by only ~ 2" or less.

When an outer planet perturbs an inner one (both circular orbits) the amplitude of the 2nd harmonic term of tangential tidal force, is approximately proportional to the inner radius, r, for fixed R and M. So the 2nd harmonic amplitude determines r*M/R^3, and M is known if the outer radius, R, is, or vice versa. (Based on 2nd harmonic amplitude alone, Barbarossa could be 198.4/5=39.7 AU from the sun with mass 3430/5^3=27 Earth masses, but such a planet would be only a magnitude dimmer than Neptune.)

For r/R = 0.1, the 1st harmonic amplitude (of tidal force) is only 0.012x the 2nd; for r/R = 0.5, it's 0.066x. The halved frequency of the first harmonic helps the amplitude of displacement (force twice integrated) fourfold, so for r/R = 0.5, the 1st harmonic displacement is 4*0.066=0.264x the 2nd.

Suppose Lowell, Harrington, and Maran, each in their own way, essentially determined the first harmonic of perturbing force. For the first harmonic, the formula r*M/R^3=const., is drastically inaccurate. Above, I used instead M/R^2*f(r/R), where alpha=arcsin(r/R), and f is either

(sec^2-cos)(alpha) or (sec-cos^2)(alpha).

The former expression gives the tangential tidal force at quadrature; the latter expression multiplies this by the sine of the angle at the sun (at quadrature as seen from the perturbed planet). Either of these expressions will maintain approximate proportionality to the actual 1st harmonic amplitude as r varies. For the range r/R = 0.1 to 0.5, the error is never more than 10%. The Maclaurin series in r/R of these expressions are, resp., (3/2)*(r/R)^2+(9/8)*(r/R)^4+..., and (3/2)*(r/R)^2+(3/8)*(r/R)^4+... .

The Maclaurin series in r/R, for the actual convolution for the 1st harmonic amplitude, using the series of 1st derivatives of Legendre polynomials, is found to be proportional to (3/2)*(r/R)^2+0.94*(r/R)^4+... . As expected from these series, exact computation shows that the former trigonometric expression above, becomes (for r/R increasing from 0.1 to 0.5) 4% too big, and the latter 10% too small. Above, I've shown that either expression implies that the 1st harmonic amplitudes given to Uranus by Lowell's or Harrington's, and to TNO's by Maran's estimates, are about the same as those given by Barbarossa.

The longitudes of Lowell's & Harrington's estimates, and also of my variational estimate, agree with each other, but are greater than that of Barbarossa. This might be because these methods (certainly my variational method) include all unknown gravitational and possibly pseudogravitational forces, not only Barbarossa.

On the other hand, Eckert et al (with or without Rawlins' refinement) effectively restricted attention to the difference between the 2nd harmonic due to Barbarossa, and a best-fitting 1st harmonic approximating it. Eckert's and Rawlins' values constitute a +/- 1" range which, considering the slopes of both the sinusoids in 1795, confirms my position for Barbarossa, +/- 2.5deg longitude.

The electron in a hydrogen atom has kinetic energy equal to 0.5*alpha^2 times its rest mass (i.e., rest energy), where alpha = 1/137.036, is the fine-structure constant. This constant is regarded as a "fundamental" physical constant. That is, its value is caused by some unknown, likely statistical, mechanism, which might recur for gravitational as for electrical fields.

Barbarossa's mean orbital potential energy happens to be -0.5*alpha^2, and its mean orbital kinetic energy 0.25*alpha^2, times the difference, between the potential energy at the surface of Barbarossa and the potential energy at the surface of the sun. Barbarossa itself apparently has somewhat less than 1/100 solar mass and theoretically about 1/10 solar diameter (or maybe less, which could account for its dimness). (The factor of 2 between this gravitational case, and the hydrogen atom case, might be from the gyromagnetic ratio g=2, or some other reason.) If the -0.5*alpha^2 (mean potential energy) formula is exact, then Barbarossa's semimajor axis should be

2*136.036^2*1392000km/2/149598000km /(1 - 0.00845*7/9/0.1) = 184.3 AU

using my best estimates of Barbarossa's individual mass (0.00845*7/9 solar masses) and diam (0.1 solar diameters).

Known distant (arbitrarily, > 100 AU, thus beyond the, somewhat disputed, "brown dwarf desert") brown dwarf companions roughly confirm the foregoing. Gleaned from the internet:

HD3651B, 480AU, 20-40 Jup mass, primary star is type K0V AB Pictoris, 275 AU, 13.5 Jmass, primary K2V GQ Lupi, 103 AU, 1-42 Jmass, primary K7eV CD-33 7795, 100 AU, 20 Jmass, primary M1

These detected brown dwarf companions are those that are unusually massive and hot. My theory above predicts that if Barbarossa's mass were 0.04 solar masses (42 Jmass) its distance would be 287 AU.

If the Barbarossa system contained only Barbarossa and Frey (not the Freya and Lowell I've hypothesized, to explain Frey's precession and Barbarossa+Frey's ecliptic longitude discrepancy, resp.) then my best estimate of Barbarossa's individual mass in solar units would be 0.00845*0.8771=0.0074115=1/134.9 ~ 1/137. Maybe such is the usual mass for a star's main distant hyperjovian companion.

Recently I estimated the mass of the Barbarossa system as 0.0103 solar masses, based on hypothetical resonances of orbital precessions due to Barbarossa in the outer solar system. Earlier (see above letter to Lowell family, Aug. 30, 2007) I got practically the same estimate, 0.0104 solar masses, based on hypothetical period equality, of the Edgeworth Belt's typical orbital precession due to Barbarossa and the Edgeworth Belt's typical precession due to the known solar system (these estimates assumed Barbarossa's orbit is circular with radius 197.7 AU). My theory of the CMB temperature and dipole (also mentioned in my letter to the Lowell family) gives (recalculated June 25, 2008) 0.008763 solar masses for the Barbarossa system (carried to precision O(m^3), neglecting the known planets, and assuming circular orbit with radius 198.4 AU because of the reduced-mass correction).

The precession effects are functions of mass*mean(orbital radius^(-3)). On the other hand, because Barbarossa's position seems to be within a degree of the (+) CMB dipole, the dipole is a function of instantaneous orbital radius. The orbital energy is a function of the semimajor axis. (All my radii neglect the mass of the known planets, which alters Barbarossa's radius by less than 1 part in 2000, for a given period.)

Last year, I emailed Lowell Observatory, telling them that Barbarossa's position coincided, in heliocentric longitude, with the Jupiter/Saturn heliocentric conjunction of 1981 (a "triple conjunction" as seen from Earth), if Jupiter & Saturn were corrected to mean circular orbits. I had taken Barbarossa's orbit to be circular. Alternatively, Barbarossa might have been near its apse, so that no mean circular correction to its longitude was needed; if so, then Barbarossa's apse is near 172.6 ecliptic longitude.

Explanation: r is 1981 radius, a is semimajor axis, e is eccentricity, m is mass. (1) says that Barbarossa's aphelion happened to occur at the time of the heliocentric Jupiter/Saturn conjunction (~1981.3), nearly the midpoint (~1980.8) of the 1954-2007 range of my Barbarossa identifications in online sky surveys and prospective amateur photos. This is justified above, as the reason no eccentricity correction (mean circular correction) needed to be applied to Barbarossa's orbital longitude, to equate it with the Jupiter/Saturn conjunction point. If perihelion is assumed instead of aphelion, the equations become inconsistent. (2) says that Barbarossa's average angular speed for 1954-2007 is that of a circular orbit at 198.4 AU (using the standard "reduced mass" equation). (3) comes from the polar equation of an ellipse, and Kepler's equal-areas-in-equal-times law. It's accurate to O(e^6). (4) adjusts the mass determined from precession resonances, to the true effective distance. (5) adjusts the mass determined from my theory of the CMB dipole, to the true effective distance.

This determines the semimajor axis as 184.5 AU, in perfect agreement with the quasi-atomic theory above. The instantaneous (1981 aphelion) radius is 194.8 AU. The mass (of the Barbarossa system) becomes 0.00845 solar masses. The eccentricity is 0.05295. By comparison, the eccentricities of Jupiter & Saturn are 0.048775 & 0.055723, resp. The arithmetic mean of the eccentricities of Jupiter & Saturn is 0.052249; the mean weighted by mass*(major axis)^2 is 0.052279; the mean weighted by mass*aphelion^2 is 0.052303.

To first order, Jupiter's and Saturn's directed eccentricities differ (subtracting vectorially using the law of cosines, then multiplying by 2.4833/(2.4833-1))(J:S value from Franklin & Soper, AJ 125:2678+, 2003) so as to amount to an equivalent eccentricity of 0.11009 with perihelion 138.9, for the cyclical advancement and retardation of the J:S resonance point.

The "Great Inequality" of Jupiter & Saturn, is known only to the range 900+/-60 yr (Vladimir Ladma, www.sweb.cz). This corresponds to a mean period for the advancing Jupiter/Saturn resonance, of 2700+/-180 yr, i.e., a major axis for Barbarossa (using standard "reduced mass" correction) of 185.8 to 203.1 AU. The usually given Jupiter & Saturn periods, 11.86 & 29.46 yr, give a period of 2758 yr and a major axis (using "reduced mass") of 197.3 AU.

So, Barbarossa's major axis in my model (184.5 AU), corresponds to a period at the lower range of accepted values for the "Great Inequality" (the reason for the uncertainty, is that the nearness of the J/S frequency ratio to 2.5, causes mathematical divergence and chaos). My model gives an eccentricity for Barbarossa, which is suspiciously near the mass*radius^2 weighted mean of the eccentricities of Jupiter & Saturn. The eccentricity and apse of Barbarossa differ grossly from those required to follow the J/S conjunction point, but "in the mean" (i.e., replacing all ellipses with approximating circular orbits) Barbarossa does follow (one of) the J/S conjunction points, about as nearly as possible considering Barbarossa's inclination.

The asymmetric (apparent) radial velocities on outwardly progressing shells of galaxies, approach equivalence to the CMB dipole, when the shells reach radius ~ 100 Mpc. This suggests (but doesn't prove) that the CMB dipole isn't a solar system phenomenon. Another objection to my solar system theory of the dipole, is that for higher spherical harmonics of the CMB anisotropy, prediction greatly exceeds observation (the coefficient of the 2nd spherical harmonic term is -36.81% = -1/e of the 1st, and of the 3rd spherical harmonic -43% of the 2nd).

However, today I found yet another prop for my theory that there is something very important about the distance, 52.6 AU from the sun. The gravitational acceleration due to the sun at 52.6 AU, is the geometric mean, of the acceleration due to the sun at Earth's orbit, and Hubble's constant times the speed of light (that is, 80.0 km/s/Mpc * c, or 7.7*10^(-8) cm/s^2). This is a few percent greater than most estimates of Hubble's constant, and a few percent less than most estimates of the Pioneer anomalous acceleration (which I and, in their own way, mainstream relativists before me have argued, is essentially also a manifestation of Hubble's constant).

Hi Joe, over in cosmicsurfers thread I give a link to the papers by Consoli, where he re examines all of the aether drift experiments. There's definitely something there, worthy of investigation. I do think though that its unfortunate that this is seen as some sort of slight breeze. It looks as though there's a slight anisotropy in light speed for the vacuum, about 3.45E-9 metres but a larger figure for air, about 72 metres. An entrained aether bubble round any mass object, drops the speed of the aether wind to very low values, depending on the refractive index of the material through which the light of the interferometer passes.

So rather than thinking of it as a breeze, perhaps we should think of it as a very peculiar windshield glass. Using my speed of gravity, where c^2 / b^2 = h (where b is the speed of gravity.) I had to find a scaling factor to fit the measured aether drift, it turned out that the fine structure constant was the best fit.

I'm still trying to make sense of it. The implication is that a diamagnetic material lets say, will have a different "awareness" of its motion than a paramagnetic material. Hmm, pretty bizarre.

Hi Joe and Stoat, In measuring h with a "moving-coil watt balance," electric current generates egual balanced upward motion of magnetic pressures against the downward motion of gravity. What if bulk charge at 4D scale levels much like a stretched rubber band is the graviton negative charged bec from bare electrons that emit virtual photon/positron/electron pairs, while the return wave of antigraviton bulk 4d scale charge is paired with forward motion of positive space. This looped current is conducted across 4D space instantaneously between the positive and negative regions forming the energy carrier wave between bare electrons and bare positrons. It is a chicken or egg scenario which came first? The fact that electric current generates opposing gravito-magnetic forces in the "moving-coil watt balance," in measuring h it is obvious that the two forces are entwined. John