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

The Schröter effect or phase anomaly of Venus

The Schröter effect was first coined by Sir Patrick Moore and was named after Johann Schröter who first noted it when he was observing Venus in the 1790's. In essence, the effect is an anomaly that is observed when determining the phase of Venus. When the planet is approaching its Greatest Elongation East (GEE) in the evening sky dichotomy is observed to take place a few days earlier than predicted whilst when the planet is approaching Greatest Elongation West (GEW) in the morning sky dichotomy is observed to take place a few days later than expected. I think it is now generally accepted that this is a genuine optical effect caused by the thick atmosphere of Venus and is not due to observer error. 

I don't think it is too hard to understand in a hand-waving way why this occurs. Most predictions about the phase of a planet or moon assume that the body is a sphere with a well-defined surface. This is largely to help make predictions about phase relatively easy to compute (see my previous blog entry). To aid with this description I include my first diagram from my earlier post:-

Venus is shown here in the gibbous phase and let us imagine that it is in the evening sky. For this to be the case Venus would have appeared from behind the Sun after a superior conjunction and then steadily increased its angular distance from the Sun as it approaches its GEE. Its phase initially would be almost full and then it would be less and less gibbous until dichotomy (half phase) is reached around GEE. In this picture we assume that Venus has a hard spherical surface and that light from the Sun can just reach the edge of the planet at the terminator (see, for example, point T).

Now imagine that Venus has an upper atmosphere that is like a thick haze. Light can penetrate this haze to a certain degree and be reflected from it, but the well-defined surface is now obscured by the haze and lower atmosphere cloud. Light destined for a point on the terminator T enters the atmosphere at a point ahead of T towards O. As it penetrates the atmosphere some of the light is scattered by the haze and so the amount of light that would have reached T is much reduced. The effect is to soften the terminator and pull its apparent edge towards O.

At the 'poles' of the planet at P and P' something else happens. Light enters the atmosphere ahead of point P and P' as it does at T. The same scattering and light reduction occurs here too due to the haze. Light that would have terminated at P or P' on a hard surface continues diffuse into the haze beyond P and P' into the dark. However, unlike at T the effect is amplified because we are looking at the atmosphere edge on - we are seeing this forward scattering over a wider depth. The net effect is that we don't see the terminator move inwards at P and P' as we did at T.

I hope you can see from my description that the overall effect is that the terminator is pushed by this scattering process to be a little ahead (sunwards) of its predicted position. This means that as we approach EGE the observed phase is always less than its predicted phase and it is this that causes dichotomy to be observed a few days early.

Now if Venus is in the morning sky then just after inferior conjunction it is a very thin crescent. On the following days its angular separation from the Sun increases until GEW is reached. As this happens the crescent phase of Venus gradually thickens until dichotomy occurs. Again because Venus has a thick atmosphere the terminator of Venus will be a little bit sunwards of where it is predicted to be. This means that the observed phase of Venus will be less than the predicted value. This will cause dichotomy to be seen a few days later than expected.

For a full description of this model of how the atmosphere affects the phase of Venus you can read the paper in the BAA journal by Anthony Mallama published in 1996. I think the idea that the observed phase is altered by scattering is compelling and it is supported by the fact that the effect is more pronounced in blue and ultraviolet light rather than yellow (see, for example, this nice observation by David Basey in 2017). It is well-known that blue and ultraviolet is more scattered than yellow light. Also, the idea that light scattering in the upper atmosphere can diffuse into the dark portion of the visible disc is supported by the observation of the extension of the cusps of the terminator into areas which are predicted to be dark when Venus is a very thin crescent (see, for example, this image by John Sussenbach in 2017).

For further discussions of the Schröter effect have a look at this article by William Sheehan in 2017.

All text and images © Duncan Hale-Sutton 2025