Real use of approximate solution to Kepler's Equation

I wanted to try out the 2nd order approximate solution to Kepler's Equation in my binary star program to see how well it compared to the numerical solution. To this end I have replaced the code in the program from the line 

X = M + Zsin M/(1 - Zcos M)

to the line

Lbl 2

(inclusive) with the following code

J = Z sin M

K = 1 - Z cos M

X = M + (-K + √(K² + 2J²))/J

Previously I had a prediction for the binary star Xi Bootis using data from my Norton's 2000.0 Star Atlas. I now want to use the same data to compare the results I had before with the results from the modified program. I get the following:-

 

Numerical                    2nd order solution

Year : Sep. : PA            Sep. : PA

2024.0 : 4.85 : 290.4    4.84 : 290.3

2034.0 : 3.90 : 270.5    3.91 : 270.7

2044.0 : 3.08 : 238.7    3.09 : 239.3

2054.0 : 2.63 : 191.3    2.63 : 191.6

2064.0 : 2.16 : 125.5    2.16 : 125.5

2074.0 : 2.76 : 53.5      2.78 : 52.9

2084.0 : 4.23 : 21.3      4.24 : 21.0

2094.0 : 5.47 : 5.2        5.47 : 5.1

2104.0 : 6.36 : 354.4    6.36 : 354.5

2114.0 : 6.92 : 346.0    6.92 : 346.0

2124.0 : 7.19 : 338.5    7.19 : 338.5

2134.0 : 7.20 : 331.2    7.20 : 331.2

2144.0 : 6.96 : 323.8    6.96 : 323.8

2154.0 : 6.50 : 315.5    6.49 : 315.5

2164.0 : 5.83 : 305.6    5.82 : 305.5

2174.0 : 4.98 : 292.7    4.98 : 292.7

I think the approximate solution does pretty well considering that the eccentricity of these stars orbit is 0.512 which is probably at the limit of what this solution can cater for. The largest error in separation is 0.02 arc seconds but mostly they are within 0.01. The largest error in PA is 0.6 degrees and this occurs in 2044 and 2074. 

If I had the situation where it wasn't practicable to seek a numerical solution then I think this approximate one would do very well.

 

All text and images © Duncan Hale-Sutton 2025

A further refinement to an approximate solution to Kepler's Equation

In my previous post I described the importance of Kepler's Equation in orbital dynamics and came up with an approximate solution which covers the cases where the eccentricity of the orbit is less than or approximately equal to 0.2. I wondered if we could do slightly better than this and come up with a solution for larger values of the eccentricity.

If we go back to equation (4) of that last post and consider a second order approximation to sin Δ and cos Δ of the form sin Δ ≈ Δ and cos Δ ≈ 1 - Δ²/2 then equation (4) becomes

Δ ≈  e (sin M (1 - Δ²/2) + Δ cos M )

Rearranging this we obtain a quadratic in Δ

½ (e sin M) Δ² + (1 - e cos M) Δ - e sin M = 0                . . . (6)

(I have assumed equality for now). If we write this in the form of

a Δ² + b Δ + c = 0

with a = ½ (e sin M), b = (1 - e cos M) and c =  - e sin M then the two solutions are 

Δ = (- b ± √(b² - 4ac))/2a

that is

Δ = (- (1 - e cos M) ± √((1 - e cos M)² + 2e² sin² M))/e sin M                . . . (7)

Note that (1 - e cos M)² + 2e² sin² M is the sum of two squares and so is always greater than or equal to zero. This means that the solutions in (7) are real. Unfortunately, the solutions are not defined when M = 0, ±π, ±2π as then sin M = 0 and they may become unstable as M tends to any of these values.

A bit of further analysis shows that only one of these two solutions is viable and that is

Δ ≈ (- (1 - e cos M) + √((1 - e cos M)² + 2e² sin² M))/e sin M                . . . (8)

I have compared this solution to the true values of Δ over all values of M and for e = 0.2 the maximum deviation is less than 0.1%. For e = 0.4 the deviation is less than 0.8% and e = 0.6 the deviation is less than 3.1%. This shows that this approximation can probably be successfully used where e ≤ 0.5.

All text and images © Duncan Hale-Sutton 2025 

Kepler's Equation

In my previous description of a program that I wrote to calculate binary star orbits I mentioned Kepler's Equation which can be stated as:-

E - e sin E = M                . . . (1)

(see section 43 on page 84 of Practical Astronomy with your Calculator by Peter Duffett-Smith). E is the eccentric anomaly, M is the mean anomaly and e is the eccentricity of the orbit. From E we can calculate the true anomaly ν. For a pair of binary stars, the mean anomaly is essentially how far round star B has gone round star A assuming that it orbited at a constant rate. If we count the number of years since the epoch of periastron (the closest approach for which data is recorded) and divide this by the orbital period then the remainder is the fraction of an orbit that has been completed since the last periastron. Multiply this by 2π and you get the mean anomaly in radians. The true anomaly is the angle taking into consideration that star B moves in an ellipse about star A, A is at a focus of this ellipse and star B does not move at a constant rate.

As you can see equation (1) is not easy to solve for E if you know M and e and that is why people resort to numerical solutions (see page 85 for Duffett-Smith's book). I wondered if there was an approximate solution which didn't rely on a numerical solution. Firstly, if we write 

Δ = E - M                        . . . (2)

Then equation (1) becomes

Δ = e sin E                     . . . (3)

If you look at page 122 of Duffett-Smith's book you can see plots Δ of versus M for various values of e. Note that the modulus of Δ is always less than e which is less than 1. This comes from equation (3); since -1 ≤ sin E ≤ 1 this implies -e ≤ e sin E ≤ e (e > 0), so from (3) -e ≤ Δ ≤ e. 

This lead me to an idea for an approximate solution to (1). Substituting for E = M + Δ in equation (3) we get 

Δ = e sin (M + Δ) = e (sin M cos Δ + cos M sin Δ)                . . . (4)

If Δ is much less than 1 then sin Δ ≈ Δ and cos Δ ≈ 1 (to first order). Thus (4) becomes

Δ ≈ e (sin M + Δ cos M)

Rearranging we get

Δ ≈ e sin M / (1 - e cos M)            . . . (5)

I have found that this solution does pretty well if e ≤ 0.2 (the maximum difference between this solution and the true estimate of Δ is less than 2.1% over all values of M).

All text and images © Duncan Hale-Sutton 2025

Observation of Mars on the 2nd February 2025

Last Sunday night we had some clear but very cold weather. I had been thinking for a while that I ought to get my Maksutov-Cassegrain out to have a look at Mars as it is so high up in the sky at the moment. Rather than spend an hour trying to set the telescope up properly I just decided to use the hand controller to keep Mars in view. This is what I had been doing with my Venus observations and it is so much quicker to get going. It is just as well because only ten minutes after I had started everything froze up (including the correcting objective lens)!

My main priority was to do a drawing of Mars rather than try to photograph it. I looked out some coloured pencils before I started and chose those which I thought would match the colours I was expecting. It was still quite a difficult thing to do. Fortunately, I worked out that I could perch on some step ladders next to the telescope so I could draw and look down the eyepiece at the same time. Holding a torch, a notebook and various pencils whilst the air temperature was below freezing is no easy task. However, I succeeded and here is the page out of my notebook:-

Mars reached opposition on January 16th and it has now moved away from us decreasing in size from about 14.6 arc seconds to 13.6. Its phase has correspondingly reduced but at 99% this was hardly noticeable. My Orion OMC 140 has a focal length of 2 metres which together with a 9mm eyepiece gives a magnification of about 222x. Even at this magnification Mars appears to be pretty small! I did try and include a x2 Barlow lens but this didn't really improve matters and it seemed that a brighter image at 222x was a better option.

The main features to be seen were the northern polar cap (which was delineated by a dark area close to its southern edge) and a large area of darker colour to the south. Not knowing exactly what I was looking at I have found a useful tool produced by Sky and Telescope - this is their Mars profiler. You can use the tool to set the date and time of your observation and to choose what telescope set-up you have. Using the date and time of 2nd February 2025 at 19:15 UT and Mirror Reversed (I was using a Maksutov-Cassegrain with a star diagonal) you get a projected map of what I was looking at. Unfortunately, it doesn't show it as a sphere so features at the edge of the map will appear foreshortened. 

What I saw does sort of match up with this projected map. The central red circle on the map is the part of Mars that points directly at Earth and this is placed over an apparently featureless red area called Tharsis which is actually a vast volcanic plateau which is home to the largest volcanoes in the solar system. With more resolution it would be possible to see white clouds associated with the tops of these mountains. The darker area to the south of Tharsis is probably Solis Lacus and shows up well on my drawing. To the left and right of this are other dark areas of Mare Erythraeum and Mare Sirenum. There is a bit of a gap between these dark areas and another dark area hinted in my drawing called Niliacus Lacus. It is to the far left of my drawing. The gap is called Chryse Planitia.

For a further explanation of the Mars profiler look at this article in Sky and Telescope. Note that the profiler also gives the position angle (PA) of the north pole and the central-meridian longitude (CM) which for my observation were 345 and 89 degrees respectively. The PA is measured from North through East and explains why the northern polar cap appears to be tilted 15 degrees westwards in my drawing. The CM is a longitude measurement on Mars and shows which longitude is pointing directly at Earth.

 All text and images © Duncan Hale-Sutton 2025

Update on the prediction for the relative positions of the binary star Xi Bootis

In my previous post I described a program to calculate the relative positions of a binary star for any given date. Using this program I can now verify the prediction that I had obtained previously using Roger Wesson's online calculator (or perhaps, more honestly, I can ensure that my code is performing the same as Roger's). Using the parameters for the Xi Bootis in that previous post Roger's utility predicted that separation and PA of Xi Boo on the 1st May 2024 (2024.333) was 4.81 arcseconds and 289.8 degrees, respectively. Running my code using the same input values I get that the PA is 289.8163749 degrees and the separation is 4.814400629 arcseconds which agrees exactly with Roger's results to the same number of decimal places as he quotes.

To make sure that my program works in a wide variety of cases I have also recomputed the predicted positions of these stars over the period 2024 to 2174 in gaps of 10 years (again refer to my earlier post) and, again, to the same number of decimal places I see no difference between the results of my code and those of Roger's. This is all good!

I have since realised that the BAA produces tables of ephemerides of visual binaries in the BAA Handbook (see page 111 for the 2025 book). I notice that the orbital elements they use come from the Sixth Catalog of Orbits of Visual Binary Stars at the US Naval Observatory. For Xi Bootis (ADS number 9413) they give the following parameters:-

Orbital Period P (years) : 152.9614

Date of periastron T : 1909.6213

Semi-major axis of orbit a (arc seconds): 4.93454

Eccentricity of orbit e : 0.51385

Inclination of orbit to plane of sky i (degrees) : 140.453

Argument of periastron ω (degrees) : 25.492

PA of ascending node Ω (degrees) :168.795

Using these parameters for Xi Bootis in my program for the 1st May 2024 (2024.333) I get that the PA is 290.9 degrees and the separation 4.97 arcseconds (as opposed to 289.8 and 4.81 found previously). This is slightly better in agreement with my observation for this date of 291 +/- 3 degrees and 5.5 +/- 0.4 arcseconds, respectively.

Finally, just to check my results against those of the BAA, using these parameters I get for 2025.0 that the PA is 289.8 degrees and the separation is 4.91 arcseconds. The BAA gets 289.7 and 4.91. For 2026.0 I get 288.2 and 4.82 and the BAA gets 288.1 and 4.82. So the BAA seems off by 0.1 in the PA each time.

All text and images © Duncan Hale-Sutton 2025

A program to calculate binary star orbits

Back in June last year I wrote a blog entry about predicting the relative positions of the binary star Xi Bootis. At that time I relied on an online calculator provided by the astronomer Roger Wesson. I thought it would be a good idea to have a means of computing these values for myself and towards that end I have now written some code for a programmable calculator (a Casio fx-4500P). This is not a particularly modern calculator (it dates from 1989/1990) but it is quite useful as it can store up to 1103 programmable steps and has 26 standard memory locations. The algorithm I have used to calculate the binary star positions comes from Practical Astronomy with your Calculator by Peter Duffet-Smith. I bought this book in 1981 but it has served me well over the years. The section (59) on binary star orbits can be found on page 130. This is my program:-

T"PERIOD"

E"PERIASTR"

Y"YEAR"

Z"E"

A"A"

L"LONG PERI"

I"I"

H"PA OF NODE"

F = (π/180)

X = (Y - E)/T

M = 2π(X - Int X)

X = M + Zsin M/(1 - Zcos M)

Lbl 1

D = X - Zsin X - M 

Abs D ≤ 1E-6 → Goto 2 ∆

W = D/(1 - Zcos X)

X = X - W

Goto 1

Lbl 2

N = 2tanˉ¹ (√((1 + Z)/(1 - Z))tan (X/2))

R = A(1 - Zcos X)

Q = sin (N + FL)cos FI

P = cos (N + FL)

X = tanˉ¹ (Q/P)

P < 0 → X = X + π: Goto 3 ∆

Q < 0 → X = X + 2π ∆

Lbl 3

O = (X/F) + H

O > 360 → O = O - 360 ∆

O < 0 → O = O + 360 ∆

O"PA" = O ▲

K"SEP" = RP/cos X ▲

Here is an explanation of some of the items in this code. The first 8 lines are the input parameters for the binary star as follows:-

 

T - the orbital period in years

E - the epoch or date of periastron

Y - the date of observation (in a decimal of a year)

Z - the eccentricity of the orbit e

A - the semi-major axis of the orbit a (in arc seconds)

L - the longitude or argument of the periastron ω (in degrees)

I - the inclination of the orbit to the plane of the sky i (in degrees)

H - the position angle (PA) of the ascending node Ω (in degrees)

The last two lines are the output values for the stars:-

O - the position angle θ of star B relative to star A (in degrees)

K - the separation ρ of star B from star A (in arc seconds)

The calculator should be set to run in radians. The symbol → is an implication sign and it forms part of a logic test. If the statement before it is true the statement after it is actioned otherwise the next line of the code is run. The symbol ∆ terminates this logic test. The symbol ▲ stops the processing at this point and displays the contents of the memory variable before the equals sign. The section of code between Lbl 1 and Lbl 2 is an iterative solution to Kepler's Equation to obtain the true anomaly (more about this later).

All text and images © Duncan Hale-Sutton 2025

Observation of Venus on the 25th January 2025

Last Saturday we had a clear patch of weather and so late afternoon, just after sunset, I thought I would have a look at Venus again to see how things have changed since half phase (dichotomy) was reached. Here is the page out of my notebook:-

The angular distance between Venus and the Sun is now decreasing but it is still pretty much in the same part of the sky as before. At the time of the observation it stood at 23 degrees above the horizon in the SSW. There was a bit of thin cloud about which glowed an orange colour in the sunset. I used the same set up as previously but due to how Venus is orientated the terminator was much more tilted than before and so I ended up just drawing it as I saw it in the eyepiece.

As you can see from my drawing the phase is now looking much more crescent in appearance. It has noticeably gone from convex to concave. This time I paid much more attention to the terminator itself and it is clear that it isn't a well defined boundary but has some fuzziness to it. I have tried to represent this in my drawing. What's more I thought I could at times detect a bit of cloud shading near the terminator - this is where there are subtle differences in the brightnesses of the clouds on Venus. Again I have very tentatively put what I saw in my drawing. One other thing was that the cusps looked somewhat brighter than the other parts of the illuminated side.

I had drawn myself some new phase diagrams on the computer to show what the predicted phase would look like from 40 to 50% in steps of 1%. With these on my phone I could compare them directly with what I saw at the telescope and I estimated that the phase was 41%. The predicted phase was a bit more than this at 42.5%, a difference of 1.5%. Here is what 42% should look like:-

All text and images © Duncan Hale-Sutton 2025