Celestial Navigation using last years Almanac

[jdc]

physical model versus efficient computation?

...All of the methods can be made less accurate by reducing the number of periodic terms used, or by reducing the number of significant digits used at each calculation stage.

Van Flandern and Pulkkinen's process used about 45 periodic terms for the Sun (Earth) position to give accuracy of about 0.02°. They use this method because it is the same method used for the Moon and other planets which required a (fairly) similar number of terms. I'm not aware of any update done to their work, or corrected figures....

L0 = 280.46646 + 36000.76983T + 0.0003032T^2
M = 357.52911 + 35999.05029T - 0.0001537T^2
e = 0.016708634 - 0.000042037T - 0.0000001267T^2
C = (1.914602 - 0.004817T - 0.000014T^2)sin(M) + (0.019993 - 0.000101T)sin(2M) + 0.000289sin(3M)
R = 1.000001018(1-e^2)/(1+e.cos(M+C))
O = L0 + C
L = O - 0.00569 - 0.00478sin(125.04 - 1934.136T)
p = 23°26'21".448 - 46".8150T - 0".00029T^2 + 0".001813T^3
...

You've clearly gone into this in much greater detail than I, but I _think_ there's a difference between what I was describing, ie modeling the motions of the planets and moon using Newton's (or Keppler's) laws and the method you describe above, which looks to be using polynomial or harmonic series to make an arbitrarily close approximation, but isn't a physical model per-se.

The latter seems to me to have a few disadvantages: it's is only useful if the coefficients have been derived first and it's no good at getting an understanding of the magnitude or source of perturbations - at least they are not explicit. However I admit it has one huge advantage: it's much more computationally efficient provided you know the coefficients. Do you know how Van Flandern and Pulkkinen's derived them?

I have to admit my interest is a kind of pig-headed attempt at self sufficiency: my 'project' is to see what I could predict for the year to come based solely on observations I could, in principle at least, make by myself with only a watch and a telescope or sextant.

 

[danielefua]

The latter seems to me to have a few disadvantages: it's is only useful if the coefficients have been derived first and it's no good at getting an understanding of the magnitude or source of perturbations - at least they are not explicit. However I admit it has one huge advantage: it's much more computationally efficient provided you know the coefficients. Do you know how Van Flandern and Pulkkinen's derived them?

I thought I mentioned it: Van Flandern and Pulkkinnen found the constants to put in the series by what is called a "best fit" procedure - a branch of numerical math particularly suitable for large computers. They start by assuming "wise" guesses for the constants based on what is suggested by the theory and iteratively find new guesses that tend to minimize the "distance" between data obtained by actual observations and the same data simulated by the series expansions themselves. Much easier to say than to do, believe me...

Daniel

 

[AntarcticPilot]

I thought I mentioned it: Van Flandern and Pulkkinnen found the constants to put in the series by what is called a "best fit" procedure - a branch of numerical math particularly suitable for large computers. They start by assuming "wise" guesses for the constants based on what is suggested by the theory and iteratively find new guesses that tend to minimize the "distance" between data obtained by actual observations and the same data simulated by the series expansions themselves. Much easier to say than to do, believe me...

Daniel

Um. How did they avoid local minima in the iterative process? Did they use simulated annealing or something analogous?

 

[Sybarite]

Um. How did they avoid local minima in the iterative process? Did they use simulated annealing or something analogous?

The simple way to do it: crossing oceans without charts, almanacs, instruments, not even a watch in a 20' boat they built themselves!

http://www.youtube.com/watch?v=shXAUu_WFmg

http://www.sansboussole.com/

 

[AntarcticPilot]

The simple way to do it: crossing oceans without charts, almanacs, instruments, not even a watch in a 20' boat they built themselves!

http://www.youtube.com/watch?v=shXAUu_WFmg

http://www.sansboussole.com/

As all mariners did before the (practically simultaneous) invention of the Chronometer and publication of practical Lunar tables.

Really, all that is necesary is to understand the limits of whatever methods you are using, and navigate accordingly. So, if I know I know my latitude but am not sure of my longitude (the usual situation prior to the mid 18th century), I plan my voyage to I can run down the Easting (or westing) for the last part of my voyage. Or I take a departure from a location known to be on the same longitude as my destination and sail dies north or south.

 

[John Barry]

You've clearly gone into this in much greater detail than I, but I _think_ there's a difference between what I was describing, ie modeling the motions of the planets and moon using Newton's (or Keppler's) laws and the method you describe above, which looks to be using polynomial or harmonic series to make an arbitrarily close approximation, but isn't a physical model per-se.

The latter seems to me to have a few disadvantages: it's is only useful if the coefficients have been derived first and it's no good at getting an understanding of the magnitude or source of perturbations - at least they are not explicit. However I admit it has one huge advantage: it's much more computationally efficient provided you know the coefficients. Do you know how Van Flandern and Pulkkinen's derived them?

I have to admit my interest is a kind of pig-headed attempt at self sufficiency: my 'project' is to see what I could predict for the year to come based solely on observations I could, in principle at least, make by myself with only a watch and a telescope or sextant.

The method I gave was based on the mean orbit, not the Van Flandern and Pulkkinen work. Their work produced a series of periodic elements (VSOP87 works the same way). Of course, they all use the mean orbit as a base.

You can use Kepler's equation, but I think that the solution I gave is simpler.

Kepler's equation E = M +e.sin(E) is trancendental, and needs to be solved by iteration.

For 12:00 (ET) of 1/6/13 T is 0.13415471, T^2 = 0.01799749. T^3 = .00241445

L0 (mean longitude of the Sun) = 280.46646 + 4829.67284 + 0.00000 = 5110.13930 = 70.13930
M (mean anomaly of the Sun) = 357.52911 + 4829.44215 + 0.00000 = 5186.97126 = 146.97126
e (eccentricity or the orbit) = 0.016708634 - 0.000005639 - 0.000000002 = 0.016702993

The T^2 functions are not part of a harmonic series, they're an allowance for precession. For dates before 2100 you can ignore them.

E = M + e.sin(E) iterated is stable (three decimals) after three cycles. Then calculate v (true anomaly) from tan(v/2) = sqrt((1+e)/(1-e))tan(E/2). Then C (equation of centre) is v - M.

However, as the eccentricity is small (<0.02) you can calculate C directly from the formula I gave. This gives C= 1.025245 for the data above, the Kepler equation gives 1.025242. These are the same within the limits of the precision we're using.

The other parts of the calculation are:
R is the Sun's radius vector
O is the Sun's true longitude (true geometric longitude referred to the mean equinox of the date)
L is the Sun's apparent longitude (referred to the true equinox of the date, corrected for nutation and aberation)
p is the obliquity to the ecliptic.

(for accuracy in GHA Aries, you should also include a correction for nutation in right ascension which I omitted, currently about +0.003°, reducing to -0.003 over the next 4 years.)

Best of luck with your watch, telescope and sextant.

 

[danielefua]

Um. How did they avoid local minima in the iterative process? Did they use simulated annealing or something analogous?

The formula by Pulkkinen and Van Flandern rely mostly on the monumental work of the famous Canadian astronomer Simon Newcomb dated back in 1898; I am not sure he cared very much about those arguments. You may find a good review paper here: http://iau-comm4.jpl.nasa.gov/XSChap8.pdf.

Daniel

 

[jdc]

The method I gave was based on the mean orbit, not the Van Flandern and Pulkkinen work. Their work produced a series of periodic elements (VSOP87 works the same way). Of course, they all use the mean orbit as a base.

You can use Kepler's equation, but I think that the solution I gave is simpler.

Kepler's equation E = M +e.sin(E) is trancendental, and needs to be solved by iteration.

For 12:00 (ET) of 1/6/13 T is 0.13415471, T^2 = 0.01799749. T^3 = .00241445

L0 (mean longitude of the Sun) = 280.46646 + 4829.67284 + 0.00000 = 5110.13930 = 70.13930
M (mean anomaly of the Sun) = 357.52911 + 4829.44215 + 0.00000 = 5186.97126 = 146.97126
e (eccentricity or the orbit) = 0.016708634 - 0.000005639 - 0.000000002 = 0.016702993

The T^2 functions are not part of a harmonic series, they're an allowance for precession. For dates before 2100 you can ignore them.

E = M + e.sin(E) iterated is stable (three decimals) after three cycles. Then calculate v (true anomaly) from tan(v/2) = sqrt((1+e)/(1-e))tan(E/2). Then C (equation of centre) is v - M.

However, as the eccentricity is small (<0.02) you can calculate C directly from the formula I gave. This gives C= 1.025245 for the data above, the Kepler equation gives 1.025242. These are the same within the limits of the precision we're using.

The other parts of the calculation are:
R is the Sun's radius vector
O is the Sun's true longitude (true geometric longitude referred to the mean equinox of the date)
L is the Sun's apparent longitude (referred to the true equinox of the date, corrected for nutation and aberation)
p is the obliquity to the ecliptic.

(for accuracy in GHA Aries, you should also include a correction for nutation in right ascension which I omitted, currently about +0.003°, reducing to -0.003 over the next 4 years.)

Best of luck with your watch, telescope and sextant.

Many, many thanks for a comprehensive and intelligent reply: it was exactly what I was asking.

In various versions of my program I used both iteration and graphical methods to invert E = M + e.sin(E) to get M, which would have been a pain on a calculator but no problem in a high level language like matlab (which is what I used).

As for my "watch, telescope and sextant" project, it's fairly successful I think, or at least I do actually use it. I split the computations into two distinct steps:

1. Run physical model, calculating and printing a table of values of GHA, DEC, GHA_Aries and earth-sun distance (in AU) for 12:0UTC for every day of the year
2. Import this table into an html + JavaScript program, which then only has to do quadratic interpolation for times other than 12:00. I added sight reduction and altitude corrections.

I think this approach of splitting the computation of ephemerides from subsequent brute-force interpolation, although using more storage than the polynomial method, is more extensible to data for other bodies; for instance I added data for Venus and Moon, just importing the tables without my having generated them (cheating I admit). It's not an elegant piece of coding but, as I said, it does seem to work.

PS: To Daniel, I'm not at all sure that LMS methods really require 'large computers', at least not by today's standards; I doubt that anything more powerful than a pdp11 was used when the original work was carried out, and astronomical functions are well behaved so 'steepest descent' would have worked ok I'd have thought.

 

[AntarcticPilot]

PS: To Daniel, I'm not at all sure that LMS methods really require 'large computers', at least not by today's standards; I doubt that anything more powerful than a pdp11 was used when the original work was carried out, and astronomical functions are well behaved so 'steepest descent' would have worked ok I'd have thought.

Thanks - it was I who raised the issue of local minima in the solutions, so thanks for elaborating on that.

I'm sure you know, but just to elaborate a bit for others, the numerical techniques for solving this kind of problem have been around for a very long time; there has even been a standard library of computer functions (the Numerical Algoriths Group (NAG) library) since the 1980s, or perhaps even before (I was using it in the early 80s). And I have a massively useful book on my bookshelf called "Numerical Recipes", which dates from 1986, and includes all the techniques required to solve this kind of problem. The rather under-powered netbook I'm writing this on is probably more powerful than the majority of computers available in 1986!