Monday, February 29, 2016

The 132 AD Autumnal Equinox Assuming Ptolemy's Latitude for Alexandria


  On p. 100 of Ptolemy, The Geography, Ptolemy lists the latitude of Alexandria a 31°. The smallest step in the table is 5 minutes of arc or 1/12th of a degree which corresponds to a round off error of 1/24th of a degree. If one assumes a spherical earth the latitude is the altitude of the pole at a particular location so it can be obtained by observation. I recomputed the 132 AD autumnal equinox using Ptolemy's latitude and got about 4 pm on  Sep 23. Since the astronomical day starts at noon this would have been Sep 24 according to that calendar.


A 1/24th degree round off error might account for Ptolemy's value of 2 pm for the time of the equinox.


correction (Mar 1): the value for the day of year in the 1st image should read doyAE = 267.681. Forgot to change the value for the month in the function.

Determining When An Equinox Occurs


  With the Vernal Equinox is less than a month away it is a good time to use the HORIZONS data to show how one can find the time of an equinox. This example starts with hourly azimuth values for the latter half of September 132. The time is represented by the day of the month plus a fraction which includes the hour of the day in UT and the difference in longitude of Alexandria from the prime meridian. A cubic spline fit is used to obtain functions which will calculate the azimuth and altitude of the Sun for an arbitrary time. The time of midday, t_md, is defined as when the Sun reaches a azimuth of 180° and is found by linear interpolation. This is equivalent to finding the time when the Sun crosses the meridian. The functions previously found are used to check the azimuth and calculate the altitude, alt_md, of the sun at noon.


The positions of the altitude at noon appear as a set of points along a line decreasing by about the same amount each day as shown by the red line in the figure below. The blue line represents the altitude of the equator which is assumed to be 90 degrees minus the latitude of Alexandria. The coefficients of the line passing through the altitudes for midday are found with the line function which does a least squares fit. The equinox occurs when the Sun crosses the equator which is easily computed using the coefficients of the line.


If the value for the altitude is 3 minutes of arc too high the calculated time for the equinox will occur about 2 hours earlier.


One can see how sensitive the results are to errors in measurement.

Sunday, February 28, 2016

Comparison of Ptolemy's Observations With Those of Horizons


  One can compare Ptolemy's observations of the solsticial and equinoctial points with data from JPL's HORIZONS online ephemeris. To avoid comparing apples and oranges we need to use methods similar to those of Ptolemy.  One has to set HORIZONS to compute azimuth and elevation using the standard atmosphere to get the apparent positions of the Sun and then determine the time in Alexandria when the Sun crosses the local meridian on days near equinoxes and solstices. One can then use linear fits to determine the local time of an equinox or solstice. The computed times are sensitive to the precision used to measure angles and an error of three minutes of arc can shift the time for an event by a couple of hours.


Although the time of the autumnal equinox for 132 AD differs slightly from stated time of Sep 25, 2pm on p. 168 of Toomer, Ptolemy's Almagest, the results compare fairly well with the intervals for the passage through a quadrant considering the accuracy with which Ptolemy worked.

The difference between the altitudes of the midday Sun at the solstices can be used to estimate the obliquity of the ecliptic which appears to be about 23°40'.


Ptolemy gives a value for the difference in the Sun's altitude at the solstices of 11/83rds of a complete circle which corresponds to the value shown above. It's difficult to determine exact times since the Sun altitude is relatively stationary at the solstices and one needs to estimate the altitude of the equator to get the times for the equinoxes. The errors for the intervals are twice those of the estimates of the time of an individual event. Ptolemy used Hipparchus' times of passage for the quadrants and method to deduce the anomaly. They may not have been correct for the year 132 AD.

Saturday, February 20, 2016

JPL HORIZONS AD 132 Phenomena (Corrected)


When fatigue starts setting in your error rate goes up and that may account for some of the calendar creep. One of my settings for HORIZONS was J2000.0 which doesn't use the reference frame for the dates calculated. There is also a switch for the Julian date and decimal RA and DEC make it easier to check the time of the equinoxes and solstices. I also used apparent solar time and got a better fit to Ptolemy's data.

     Date__(UT)__HR:MN   JDUT     R.A._(a-appar)_DEC.

      0132-Mar-19 00:00 1769348.5  357.63697  -1.03602
      0132-Mar-20 00:00 1769349.5  358.53070  -0.64442
      0132-Mar-21 00:00 1769350.5  359.42378  -0.25287
      0132-Mar-22 00:00 1769351.5    0.31632   0.13854
      0132-Mar-23 00:00 1769352.5    1.20843   0.52971

      0132-Jun-21 00:00 1769442.5   87.33879  23.65461
      0132-Jun-22 00:00 1769443.5   88.37944  23.66893
      0132-Jun-23 00:00 1769444.5   89.42051  23.67628
      0132-Jun-24 00:00 1769445.5   90.46193  23.67666
      0132-Jun-25 00:00 1769446.5   91.50360  23.67006
      0132-Jun-26 00:00 1769447.5   92.54542  23.65648

      0132-Sep-22 00:00 1769535.5  178.07833   0.84227
      0132-Sep-23 00:00 1769536.5  178.99124   0.44213
      0132-Sep-24 00:00 1769537.5  179.90462   0.04163
      0132-Sep-25 00:00 1769538.5  180.81856  -0.35911
      0132-Sep-26 00:00 1769539.5  181.73312  -0.75999
      0132-Sep-27 00:00 1769540.5  182.64839  -1.16091

      0132-Dec-19 00:00 1769623.5  266.83065 -23.64465
      0132-Dec-20 00:00 1769624.5  267.94346 -23.66328
      0132-Dec-21 00:00 1769625.5  269.05640 -23.67395
      0132-Dec-22 00:00 1769626.5  270.16930 -23.67667
      0132-Dec-23 00:00 1769627.5  271.28202 -23.67144
      0132-Dec-24 00:00 1769628.5  272.39439 -23.65827

Using cubic splines gives the following estimates for the equinoxes. The solstices were found by estimating the zero for the differences.

      Mar 21 15:30   Vernal Equinox
      Jun 23   01:18   Summer Solstice
      Sep 24  02:30   Autumnal Equinox
      Dec 21  08:13   Winter Solstice

The astronomical day starts at noon so the morning hours would be on one day and the afternoon hours on the next. If the dates indicated are ordinary Julian calendar dates then Sept 25 12:00 would be Julian date 1769539.0 as previously determined. There is still seems to be a little bit of discrepancy between HORIZONS' and Ptolemy's dates and times.

It might help if we could develop a "personal equation" for Ptolemy as an observer. It might help to characterize the Alexandrians of his time and the members of the Museum.

edit (Feb 28): corrected the time of the summer solstice to 01:18 which is what my saved data indicates.

Thursday, February 18, 2016

JPL HORIZONS 132 Autumnal Equinox


  I checked JPL's HORIZONS online ephemeris for the date and time of the 132 AD autumnal equinox and received the following results from which the time of the equinox can be determined.

      Date__(UT)__HR:MN     R.A._(ICRF/J2000.0)_DEC
      0132-Aug-28 19:00     11 59 46.47 +00 01 17.9
      0132-Aug-28 20:00     11 59 55.46 +00 00 18.7
      0132-Aug-28 21:00     12 00 04.46 -00 00 40.5
      0132-Aug-28 22:00     12 00 13.45 -00 01 39.7

      18.7/59.2*60=19.0

      interpolated equinox date & time
      132 Aug 28 20:19 UT

This is the descending node since the declination is decreasing as it passes through 0° DEC as expected for the autumnal equinox. On p. 168 of Toomer, Ptolemy's Almagest Ptolemy states,

      "Now we have established that, among the first of the equinoxes 
      observed by us, one of the most accurately determined was the 
      autumnal equinox which occurred in the seventeenth year of 
      Hadrian, on Athyr [III] 7 in the egyptian calendar [132 Sept. 25], 
      about 2 equinoctial hours after noon."

JPL's equinox is about a month prior to Ptolemy's. It appears that somehow we got the date wrong. What could be the explanation? On p. 9 Toomer says that Ptolemy used the 365 day Egyptian calendar. The Julian calendar has on average 365.25 days per year so there were about 132/4 = 33 leap days which may account for some of the difference. Each leap day would shift the Egyptian calendar one day forward relative to the Julian calendar. The dates given in HORIZONS are Julian dates and adding 33 days would take us to Sept. 30. The equinox is still a few days off and Ptolemy's datum is probably just an estimate. If Ptolemy didn't take into account the fact that the diffraction of the atmosphere affects apparent positions his pole would be shifted slightly north as well as his equator if he assumed it was 90° south of the pole. Since the Sun descends in declination at the autumnal equinox his estimate of the date and time would be a little early. Ptolemy was concerned about diffraction later in life so he may have become aware of the problem after the Almagest was written. An error in Ptolemy's anomaly could also affect the times of the equinoxes that of perigee.

Supplemental (Feb 18): While researching calendar creep I came across this Amagest Ephemeris Calculator which gives Toomer's Julian calendar date for Hathor 7. The discrepancy in the date of the autumnal equinox is a puzzle.

Using Ptolemy's Anomaly Table to Find the Time of Perigee


  As indicated in the last blog Ptolemy's anomaly table can be used to the time it takes to go from one longitude to another. To find the time interval between the autumnal equinox and perigee one can use linear interpolation. Ptolemy's anomaly table is found on page 167 of Toomer, Ptolemy's Almagest. The first two columns are mean motion, M below, and the third column is the anomaly, α below, in degrees. The apparent longitude is the sum of the sum M + α  and one looks for two longitudes in the table that straddle some arbitrary longitude, λ, then do the interpolation.


One gets essentially the same result found previously. Finding the time from perigee is a special case but the time between any two longitudes will be the difference of the individual times to/from perigee.

Wednesday, February 17, 2016

Did Ptolemy Have a Time Interval Table or Formula?


  The time interval formula is all that is needed to compute the times of other events if one has a known position and time for the Sun like Ptolemy's autumnal equinox.


Did Ptolemy a table of time differences from the time of perigee that allowed him to compute the time for passage between two longitudes? Ptolemy's anomaly table may be such a device. Ptolemy gave the angular position of perigee so one can compute the angular distance from perigee for a given longitude. Using the anomaly table in reverse one can determine the corresponding mean motion and the time interval from perigee. Doing this for two longitudes allows one to determine the time it would take to go from one to the other. So if the date and time of one longitude is known one can compute the date and time of another.

Supplemental (Feb 18): The Julian dates assume UT which is Greenwich Mean Time but the latitude of Alexandria is about 30° E so it is 1/12th of a day ahead of UT and this means that one needs to subtract 0.083 from the times above to correct them.

Ptolemy's Ecliptic Longitude Formula & Event Julian Dates


  The longitude formula assumed in the last blog gives incorrect ecliptic longitudes and we need to be more careful about the Julian dates of the equinoxes, solstices and perigee. The times for the quadrant passages were correct since they involved 2nd differences.


Since the time of the autumnal equinox is known, the Julian dates of the vernal equinox and solstices can be computed. The epoch is chosen to be the beginning of the year 132 and is given for comparison.


The anomaly is measured from the time of perigee, tp, and this longitude is given by Ptolemy. Since the time of the autumnal equinox is also given we can determine the date and time of perigee.


We now have the correct longitude formula and can compute the true anomaly for a set of times of the year.


This data can then be used to fit an ellipse to Ptolemy's anomaly. As one can see there is a systematic error of about half a degree for the residuals of the linear formula involving the true anomaly and its rate of change but this is below the 1 minute of arc precision that Ptolemy appears to have been working at.


The values of σ and τ allow us to compute the Keplerian elements as before.


Tuesday, February 16, 2016

Keplerian Elements From Ptolemy's Anomaly for the Sun


  Ptolemy gave a table of the Sun's anomaly that was accurate to about 1 arc minute in his Almagest along with other data from which we can make an estimate of the Sun's "orbit."
We start with a figure of the equinoxes and solstices and the positions of perigee and apogee in the plane of the ecliptic. In Toomer, Ptolemy's Almagest, p. 153, Ptolemy states that the apogee was 24.5° in advance of the summer solstice which would place it at Gem 5.5° in his twelve 30° divisions of the ecliptic measured from the vernal equinox. Similarly the perigee is at Sag 5.5°.


The positions of perigee and apogee were determined by the nonuniform rate of motion of the Sun being slowest at apogee and fastest at perigee. He also gives Hipparchus' values for  the time it took to go from the vernal equinox to the summer solstice (94½ days) and the time to go from the summer solstice to the autumnal equinox (92½ days). He says that his observations gave similar values and from a fit of his anomaly we see that this is so.


The time it takes to traverse a quadrant is greatest for the quadrant containing the apogee and least in the quadrant containing perigee. The value for λ0 was chosen so that the autumnal equinox occurs at 132 Sep 25 2pm. Ptolemy's anomaly only has sin terms in it and is accurate to about 1 arc minute as the residuals for the fit show.


We can use Ptolemy's ecliptic longitude to compute an "orbit" for the Sun in the same way that we did using the longitude obtained with the MICA data. The epoch chosen was the beginning of the year, 132 Jan 1, noon, and we can use the the longitude to estimate the angle and time of perigee, then calculate the true anomaly, ν, at intervals of one day.


Using the linear relation between the cosine of ν and the square root of its rate of change we can estimate σ and τ which in turn allows us to estimate the eccentricity, e, and the other Keplerian elements. The residuals for fit of the formula are a few minutes of arc.


One can find Ptolemy's value for obliquity of the ecliptic on p. 63 of Toomer's book. The summary of the estimated Keplerian elements is as follows.


Most of the changes in orientation of the elliptical orbit are due to the rotation of the equatorial reference frame in the plane of the ecliptic. One would expect some of the difference in the value for the eccentricity e to be due to Ptolemy's longitude and the use of an eccentric circle but some of it may be due to perturbations also. The astronomical unit, AU, used here is based on the elliptical orbit for Ptolemy's data but a change in the mean motion would suggest a change in the mean radius since the areal velocity is fairly constant.

edit (Feb 17): Corrected the date and time for the autumnal equinox but the formula assumed for the longitude above proved to be incorrect. See the next blog for more details.

Tuesday, February 9, 2016

What Does Linear Least Squares Do?


  The equations for the linear least squares coefficients are derived from the assumption that the sum of the square of the errors is a minimum. This may not be the actual case as the following example shows. The result is slightly biased with the fit line (dotted) being slightly oblique to the true line (solid).


A calculation shows that for the fit the 1st moment of the errors, δ'', is zero while this is untrue for the assumed errors, δ. This result is quite general and can be proven assuming the equations for the fit.


One can also show the unexpected result that the average of the errors is also zero starting from the first of the fit equations.


Both conclusions are due to the original assumption that the sum of the square of the errors is a minimum and may be in conflict with the true values. For random errors the conclusions are approximately true and are not in conflict with the use of linear least squares for an estimate of the coefficients.

Residuals for the Solar Longitude Fit


  The residuals for the longitude fit in the last blog have arc second accuracy. Here are the final correction for the longitude fit and resulting residuals.



The periodic sin function present appears to be associated with the Moon orbital motion and can be removed leaving a long term residual.


This may be due to planetary perturbations and/or errors in the fit. The lunar synodic period works quite well with time measured in days. The fit of Sun's motion with the lunar disturbance of the Earth's position added appears to have about 10 arc second accuracy.

Tuesday, February 2, 2016

Finding the Mean Motion by Averaging


  One does not have to use least squares to find the mean motion of the Sun along the ecliptic. The 3 years of data for the ecliptic longitude obtained from MICA give a fairly good estimate of the mean motion, n. One can interpolate to find the time it takes to move an integral multiple of 360° then divide to estimate the mean motion.


Using this estimate of the mean motion it's a simple task to find the nonuniform motion.


The longer the interval of time the better the estimate of the mean motion as can be seen in the following plot. Notice that the difference between the successive averages and the original estimate returns to zero at the end of the years.


Judging by Ptolemy's value of the mean motion found in Toomer, Ptolemy's Almagest, p. 140, it hasn't changed much in 2000 years.


Ptolemy's ecliptic longitude was measured from the Vernal Equinox. To find the mean motion in a fixed reference frame like J2000.0 one has to subtract the precession rate of the Equatorial reference frame. A similar numerical value is used to define the astronomical unit (AU) which is the radius of a circular orbit that has this value for its mean motion. The AU is a trifle smaller than the Earth mean distance from the Sun.

Monday, February 1, 2016

A Problem With Least Squares Fits and a Fix


  One sometimes encounters a problem when using least squares to find the best fit of a formula to a set of data. A simple example is fitting a line to the Sun's ecliptic longitude as a function of time.


Notice that the nonuniform motion, δ, slopes slightly downward when one would expect just a sinusoidal function. The problem is that least squares warped the error raising it on the left and lowering it on the right to reduce the variance. If one could fit the smaller nonuniform motion first one might avoid the problem but one has to fit the linear portion first to estimate the mean motion, n. What one really wants to do is separate the uniform motion from the nonuniform motion. This can be done by successive residual fits that alternates between finding uniform and nonuniform motion in the residuals. The function correct_α does 10 iterations of this here.


The nonuniform motion is level after the correction and the variance, V, is slightly larger. Here are the programs that were used to do the corrections.


The first function, Δα, fits the nonuniform motion in a set of residuals then uses an new residual to extract some uniform motion from it. The second function, correct_α, repeats this process nitr times.