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