Programming and Technical Aspects

This is linked from Solar System Graphics, which introduces the subject and shows how to use SatGraph. This node discusses SatGraph and technical aspects for anyone interested. The code can be accessed from the various sections discussed below but the actual Visual Basic source code may be harder to obtain.


SatGraph program


The source code discussed here can be obtained from the file main.txt

frmPlanet Form



This routine sets up a fixed string array for the months and another for telescope apertures to aid the input for the user. It also reads data from two files into arrays. From the file Plan2000II.dtf are input 17 values for each of the 8 planets. The contents of this file and the corresponding array, m_dblPlanets( , ) are described below in the section “Routine btnCalculateSolar”. It also inputs from a file, PlanetPerturbII.txt, the perturbations in longitude, eccentricity, perihelion, and semimajor axis for Jupiter and Saturn discussed in MeeusI, Chapter 24. These data are stored in arrays m_dblPL( , ), m_dblPecc( , ), m_dblPH( , ), and m_dblPr( , ).

Subroutine Longitude2

This is called in btnCalculateSolar to calculate the longitude and distance from the Sun of the Earth and the planets. It is also called in Satellites_Load to calculate the longitude and distance in the plane of the satellite orbit. When used in CalculateSolar, the input values are the planet array, m_dblPlanets, the planet number and the days elapsed since the epoch. When used in Satellites_Load, the input values are the satellite array, m_dblSatellites, the satellite number, and the number of days elapsed. For the major satellites, and also for the planets when perturbations have been calculated, input values include the minor corrections for longitude, eccentricity, perihelion, and semimajor axis.

Kepler’s equation is solved iteratively, using (22.3) in MeeusI. Then the true anomaly is calculated using (25.1) in MeeusI and used to calculate the longitude. Expression (25.2) in MeeusI is used to find the radius.

As discussed in the routine CalculateSolar, the most important perturbations for Jupiter, Saturn, Uranus, and Neptune are calculated (see MeeusI, Chapter 24). These perturbations, in longitude and perihelion, are passed to Longitude2 and incorporated in the results.

Routine btnCalculateTimeElements

In this, the user either enters the Julian Date or the Universal Date and Time. In the former case, it is simple to calculate the days elapsed since the epoch 2000 by subtracting 2451545 from the Julian date. In the latter case, a subroutine DaysSinceEpoch is called to find the days since the epoch and this is then added to 2451545 to create the Julian Date. The Siderial time is also calculated.

Subroutine DaysSinceEpoch

This inputs Year, Month, Day, Hours, Minutes, Seconds and outputs Siderial Time and the days and part thereof elapsed since 2000.

Routine btnCalculateSolar

Many of the values read from the file Plan2000II.dtf and stored in the corresponding positions in the array m_dblPlanets are used by this routine. Hence, the data are described here.

They are based on epoch 2000.

Elements 0 to 10 are obtained from Table 30.a and Table 30.b of MeeusII or ESII, p704.

Element Constant
0 a (in AU)
1 e eccentricity
2 change of e per century
3 i inclination of orbit
4 change of i per century
5 L longitude
6 change of longitude per century
7 LAN Longitude of ascending node
8 change of LAN per century
9 PI (LAN + omega)
10 change of PI per century
11 diameter
12 brightness factor
13 R inclination of equator to ecliptic
14 change of R per century
15 LAN of equator on ecliptic
16 change of LAN per century

The variables used in MeeusI, Ch 24, “Planets: Principal Perturbations”, are set up.

The ecliptic longitude and orbital radius of the Earth are calculated, using subroutine Longitude2. The “approximate” longitude and orbital radius of the planet are calculated - approximate because light-time is unknown. The subroutine ConverttoPlane is used to calculate the longitude, latitude , and projected radius in the plane of the ecliptic. Subroutine GeoCoords is used to calculate longitude, latitude, projected radius, and “ radius” i.e. distance , as seen from Earth. The distance is then used to calculated the delay due to the finite speed of light, or “light time”.

It is deemed very important that the calculations for Jupiter and Saturn be precise because of their use for satellite calculations. Both MeeusI and MeeusII have chapters on these planetary “perturbations” - Chapter 24 in MeeusI and Chapter 31 in MeeusII. The latter is more up-to-date, more complex, and presumably more accurate than the former. However, SatGraph uses the former. I have done the MeeusII calculations in an Excel spreadsheet and it would not be very difficult to incorporate them into SatGraph at a later date if deemed desirable. This will be discussed in the technical section.

The data in arrays m_dblPL( , ), m_dblPecc( , ), m_dblPH( , ), and m_dblPr( , ) which were created in Form1_Load are used to calculate the minor changes in longitude, eccentricity, perihelion, and semimajor axis for Jupiter and Saturn. Refer to the section Technical Aspects for more details on how this has been coded. The corresponding minor changes for Uranus and Neptune are also calculated, but the data is built directly into the program rather than read from a file.

Longitude2 is again called, with these minor changes passed to it for Jupiter, Saturn, Uranus, or Neptune to get the more-precise values which allow for light time, and ConverttoPlane and GeoCoords again called.

After the true longitude, latitude, and radius vector have been calculated, it is necessary to add other terms ( latitude only for Saturn, p.116 MeeusI, and all three for Uranus, p.117, and for Neptune, p. 119).

Subroutine ConverttoPlane

See MeeusII, Chapter 32, p.209 and also pp295-6

Subroutine GeoCoords

This subroutine calculates the coordinates (longitude, latitude, and distance) of a body as seen from Earth from the coordinates around another body e.g. Sun or a planet. See MeeusII, Chapter 32, p209.


The source code discussed here can be obtained from the file frmLocalCoords

MeeusI, Chapter 8, gives formulae for transformation of ecliptical coordinates (longitude and latitude) to equatorial coordinates (Right Ascension (RA) and Declination), given the inclination of the equator of the Earth, in ((8.3) and (8.4).

Subroutine altaz is used to find Hour Angle, Altitude, and Azimuth.

Apparent angular diameter, phase, and visual magnitude are calculated for the planet.

Apparent angle of equator as seen from Earth and apparent angle of equator as seen from Sun (for later use) are calculated.

Subroutine altaz

Input parameters are RA, Declination, Latitude of observer on Earth, Longitude of observer, and Siderial Time. Output parameters are Hour Angle, Altitude, and Azimuth. Refer to MeeusI, Chapter 8.

Physical observations of Mars, Jupiter, and the rings of Saturn are mainly taken from ESI and MeeusII. First, we’ll discuss Mars. See the technical section for a comparison of the calculations used in the two books.


Refer to ESI, pp.334-337 and MeeusII, Chapter 41.

a0 and d0 are calculated using MeeusII, p272

It is desirable to calculate the declination of the Sun (Ds)and the declination of the Earth (De), both relative to Mars. To do this, it is first necessary to calculate Delta(Δ), Omega (Ω), and I given on the right hand side of Table 11.2 on p336 of ESI (ESI just presents these in the table, but does not produce the calculations). The calculations for De and Ds are shown on p334 and again in Example 11.7 on p337.

The Right Ascension ( Ae) of the Earth as seen from Mars is calculated using the expressions at the top of page 334 and Example 11.7.

The Longitude of Central Meridian (LCM) is calculated from the expressions at the bottom of page 336. The LCM may be useful to know what (hazy) features to expect to be visible.

If we could see much detail on Mars (unlikely for an amateur), Ds would be useful regarding shadows.


Refer to ESI, pp338-340 and MeeusII, Chapter 42.

a0 and d0 are calculated using MeeusII, p278. In order to find the LCM of the two rotational regions, W1 and W2 are setup (Step 2 of MeeusII, p278).

Ds, and De are calculated as for Mars , but using MeeusII (Steps 9,10 and 12). LCMI and LCMII are calculated. The position angle of the axis, P (see MeeusII p271) is also calculated for Jupiter. It comes from (41.4) in MeeusII or, more simply, the first equation at the top of page 334 in ESI.


Refer to ESI, pp362-366 and MeeusII, Chapter 44, p301.

Here, we calculate:

B, the Saturnicentric latitude of the Earth, referred to the plane of the rings

B’, the Saturnicentric latitude of the Sun, referred to the plane of the rings

the correction to visual magnitude due to the rings (brighter when more open).

B determines the openness of the rings and B’ the amount of sunlight on the rings. If B and B’ are of opposite sign, the rings are invisible because we see the unlit side

B can be calculated using the 3rd equation on p345 of ESI, after calculating J and N. It can also be calculated more simply by using step 6 at the top of page 303 of MeeusII.

B’ can be calculated using the 3rd equation on p364 of ESI or by step 8 on p303 of MeeusII.

The mv correction is on p269 of MeeusII.


This routine inputs data from a file Sat2000II.dtf. For each of 14 satellites ( 2 for Mars, 4 for Jupiter, and 8 for Saturn), it inputs 13 values which are the same as elements 0 to 12 in the planet file Plan2000II.dtf discussed above. The data is stored in the array m_dblSatellites (13,12).

Also created is an array g_dblSat (7,16) for storing data about the satellite system being investigated (maximum 8 for Saturn). The data for each satellite is built up throughout the execution of the routine and is :

Element Value
0 Ecliptic Longitude
1 Ecliptic Latitude
2 Projected distance from planet (in AUs)
3 mv
4 Apparent longitude as seen from Earth (0 in front, 180 behind)
5 Ratio of radius of orbit to radius of planet
6 Difference in longitude as seen from Earth (in planetary radii)
7 Difference in latitude as seen from Earth (in planetary radii)
8 Difference in longitude as seen from Sun
9 Difference in latitude as seen from Sun
10 Data on phenomena:
1 for transit, 2 for occultation
eclipse: 3 or 23 (occultation plus eclipse)
shadow: 4 or 14 (transit plus shadow)
mutual eclipse: 80, or 84 (mutual eclipse plus shadow on planet), or 94
11 X of shadow
12 Y of shadow
13 Degrees changed per hour
14 Diameter of satellite / Diameter of planet
15 Siderial period (in days)
16 Heliocentric distance

In Chapter 43, MeeusII discusses a method of high accuracy applied to the Jovian satellites (pp 288-298). Satgraph implements the major parts of this theory and, to shorten the program, reads from two files, Jovian_Periodic_Formulae.txt and Jovian_Periodic_Coefficients.txt. To effect the calculations, an array m_dblJovSat(3,8) is created. The elements of this array will be discussed in the Jupiter section below.

Elements 3, 13, 14, and 15 are put into array g_dblSat.

For Mars, or Jupiter, or Saturn elements 5, 0, 1, and 2 are calculated . Each planet and its satellite system is treated differently as discussed below.


Satellites are treated in Chapter 12 of ESI form page 342 to page 392. Mars is from p.350 to p.354.

For Phobos and Deimos, formulae on p.351 are used to calculate the node N and the inclination J of the orbital plane on the equator of the Earth.

A periodic term in the longitude is calculated from the expression for “ l ” in the middle of p.350 and this is passed to Subroutine Longitude2, along with the other satellite details, to calculate the longitude and radius in the satellite’s orbit. Subroutine ConverttoPlane is then used to find longitude, latitude, and radius in the ecliptic plane and these are put in the array g_dblSat.


Chapter 43 of MeeusII gives periodic terms in the longitudes, the latitudes, and the radius vectors of the 4 major satellites. In total, there are 157 longitude terms, 35 latitude terms and 45 radius terms. Fortunately the coding for Satgraph can be shortened because the longitude terms are all sines of linear combinations of mean longitudes (l1,l2,l3,l4), perijoves (π1,π2,π3,π4), and nodes ( ω1, ω2, ω3, ω4) of the 4 satellites , the latitude terms are all sines of similar linear combinations and the radius terms are all cosines of the same.

So the file Jovian_Periodic_Formulae.txt contains group of 4 characters which obtain each term in the linear combination from the data on l,π, and ω stored in the array m_dblJovSat(,) and Jovian_Periodic_Coefficients.txt contains the multiplier of the sine or cosine term.

Description of m_dblJovSat(3,8)

Element 1 is 0,1,2,3 indicating the satellite Io, Europa, Ganymede, or Callisto.

For each satellite, element 2 (These are angles, stored as radians) is:

0 l mean longitudes of the satellites
1 π longitudes of the perijoves
2 ω longitudes of the ascending nodes on the equatorial plane of Jupiter
3 Σ This is the calculated sum of terms
7 B calculated latitude
8 r calculated orbital radius

The other elements are (refer to Chapter 43, p290 for definitions);

(0,4) Φ λ
(0,5) П
(0,6) 1.9927 Σ1
(1,4) Ψ
(1,5) 52.225 This is in degrees, but is converted to radians
(1,6) 1.0146Σ2
(2,4) G
(2,5) 188.37 This is in degrees, but is converted to radians
(2,6) 4.03Σ3
(3,4) G’
(3,5) 149.15 This is in degrees, but is converted to radians

Current values for elements 0, 1, and 2 are loaded from the satellite array m_dblSatellites for each of the 4 satellites.

Φ λ, Ψ, G, G’, and П are loaded according to p.290 of MeeusII.

The longitude terms for each of the 4 satellites are evaluated and summed by obtaining the argument of each sine term from file Jovian_Periodic_Formulae.txt and the multiplier from file Jovian_Periodic_Coefficients.txt. The 4 sums are Σ1, Σ2, Σ3, and Σ4.

The true longitudes L are then the sum of the crude longitudes and the Σ.

In a similar way, the latitude terms are summed for each satellite and the inverse tangent of each sum gives the corresponding latitude.

Similarly, the radius terms are summed and the radius vector Ri of each satellite calculated as on p.295.

Then, for each satellite, the longitude ( including precession), the latitude, and the radius are converted to the plane of the ecliptic and stored in elements 0, 1, and 2 of g_dblSat .


These calculations are the most complicated of all the satellite systems. They are treated in ESI on pp.364-386. They are treated in ESII on pp.356-368.

The 1999 edition of MeeusII has a new chapter on them but I have not had access to it and hence it is not included.

The first 4 satellites, Mimas, Enceladus, Tethys, and Dione are treated together and then Rhea, Titan, Hyperion and Iapetus are treated individually.

First 4 satellites ESI, pp 366-368 ESII, pp356-360

The variations in longitude, δl, are calculated as on p.367 of ESI. These, along with the satellite data in array m_dblSatellites, are passed to subroutine Longitude2 to calculate longitude and radius in the plane of the orbit and then subroutine ConverttoPlane calculates longitude, latitude, and radius in the plane of the ecliptic. All of the other satellites discussed below also call Longitude2 and ConverttoPlane after the necessary calculations.


A fuller discussion is given in the section “Technical Aspects ” below.

(Ω - Ω1)sini1 is used from ESI and the equivalent is in (6.422-1) of ESII. Degrees and radians in the latter have been converted to minutes of arc to agree with ESI and I have left both calculations in for comparison, but ESII currently supersedes ESI. The same applies to i - i1.

The periodic variation in longitude (λ ) is taken from ESII.

δП and δe are given by straightforward expressions in ESI but need to be derived from ESII (see the technical section below). This derivation gives different, but feasibly similar ball-park, values for the constants but different values for the epoch of zero argument of sine and cosine. Hence, until clarification, I will use the same argument as ESI but use the constants derived from the ESII theory.


Again, a fuller discussion is in the section Technical Aspects.

The main programming complications are :

(i) The solving of the spherical trigonometry formulae on p. 375 for Γ, Ψ and Θ ;

(ii) Solving iteratively for g on p. 373 ; and

(iii) Calculating the perturbations Δe, ΔE ( and hence ΔL), Δi , ΔΩ ,and ΔП on p.375 of ESI.

Then subroutines Longitude2 and ConverttoPlane are called.


Δe, ΔL, and ΔП are calculated as on p.378 ESI or pp. 363-4 of ESII and then subroutines Longitude2 and ConverttoPlane are called.


I am following ESII (pp. 364-6) only . The perturbations are then passed to Longitude2 and then ConverttoPlane converts to the plane of the ecliptic.

Then, for each satellite, Longitude, Latitude, Radius of orbit, and mv are stored in the array g_dblSat( ,) (see the section Satellites_Load above) in elements 0 to 3 respectively.

The elements in array g_dblSat( ,) are now used to convert data about the satellites of the current planet from ecliptic coordinates to Earth and Sun coordinates.

As seen from Earth, the longitudes and latitudes of the satellites are calculated and hence the apparent separations from the centre of the planet. Allowance is then made for differential light time and perspective effect (MeeusII, p297). The resulting values, ΔX and ΔY, are stored in array g_dblsat( , ) in elements 6 and 7 for later use by the graphics. They are also converted to polar coordinates, r and θ for checking for transit or occultation.

The same calculations are done for the satellite as seen from the Sun (for checking for eclipse or shadow) and stored in elements 8 and 9 and the radial distance is stored in element 16.

Since the planets are ellipsoids, the radius of the planet at the angle of the satellite as seen from the Earth, and from the Sun, is calculated. Then the relative diameter of the umbra of the planet’s shadow, and the thickness of the penumbra, at the distance of the satellite are calculated. It is then possible to check for phenomena - transit, occultation, eclipse, and shadow.

Data about any satellite phenomenon are stored in other arrays g_bytPhen( , ) and g_sngPhen( , ).

Transit or occultation occurs if r is less than the radius of the planet at angle θ and such is recorded in g_dblsat( ,10) and these 2 arrays with 1 for a transit and 2 for an occultation.

Also calculated (in Subroutine PhenLimits discussed below) are:

the number of minutes in the past or future the phenomenon occurred on the right;

the same on the left;

the X position on the right;

the Y position on the right;

the X position on the left;

the Y position on the left

the number of minutes to enter or exit the phenomenon.

These are stored in array g_sngPhen( , ), with the first subscript being a sequential order of the phenomenon (holding the number of the satellite) and the second subscript 0 to 6 for these values.

An eclipse is checked for and, if occurring, the phenomenon limits are again calculated and stored, along with the indicator 3 for eclipse.

If the satellite is in front of the planet and its separation from the centre less than the radius at that angle, a shadow is possible. However, it may be near the edge not seen from Earth. Also, the apparent X position (horizontal) of the shadow is different from the X position of where it would fall on the plane through the centre. (The Y position (vertical) is not significantly affected, at least in the case of Jupiter, which is more relevant than for Saturn). The apparent X and Y coordinates of the shadow are stored in g_dblSat( , ) for later use by the graphics section of SatGraph and indicator 4 stored. Limits of the shadow are not calculated.

Mention is made here the book “The Planet Jupiter - The Observer’s Handbook” by Bertrand M Peek, 2nd edition, 1981. It discusses the shadow in App. II, p. 230 with reference to the discrepancy in time of the apparent position of the centre. This is different from my calculations but shows what has to be considered. This book is a classic that was written well before the Voyager space probes revealed so much about Jupiter and its satellites, so much is outdated but it is an interesting read for any amateur observer of Jupiter. Moreover, it discusses mutual satellite phenomena , which are mentioned below.

For each pair of satellites of Jupiter and Saturn , SatGraph checks for occultations and eclipses. An occultation occurs if the apparent separation, as seen from Earth, is less than the sum of the radii.

Peek, p. 219 points out that mutual eclipses are more fascinating than occultations because the satellites may appear to be well separated. An observer may think he or she is imagining things if one fades or disappears for a short time and then emerges again. An eclipse occurs if the separation, as seen from the Sun, of the centre of the eclipsed satellite and the centre of the shadow is less than the sum of the radii of the eclipsed satellite and the shadow. Then, the data is stored in g_dblSat( , ) and in some shared variables for use by the graphics part of the program.

Graphics routines



Sub Draw

Plan View looking down on the Solar System

Firstly, this draws diagrammatically a plan view of the satellite system around the parent planet. The main VB routines used here are graphics.FillPolygon to draw a black background and a dark grey shadow thrown by the planet and graphics.DrawEllipse to draw the circular orbit of each satellite and graphics.FillElipse to draw a circular representation of the planet and each satellite. Telescopic view

It then draws the telescopic view of the planet and its satellite system. This is much more complex, as it is mostly done from mathematical calculations followed by the use of graphics.Drawline for drawing short line segments. It includes drawing bands on a planet which has an inclined equator, has rotated in longitude, and has non-zero ecliptic latitude. It also accounts for non-zero phase , rings for Saturn, and polar caps for Mars.

Firstly, it calculates the scale of the image. Unless the user chooses to zoom (done by requesting the maximum number of planetary radii to include), the scale is chosen such that all the satellites brighter than a given magnitude can fit on the screen. This is done in a subroutine called FindScale( ) which also calculates the position on the screen to put the centre of the planet.

The planet is drawn with 6 bands of colour above the equator and the same bands below the equator. A subroutine RotateVector( ) is used to rotate any 3- dimensional set of Cartesian coordinates (X, Y, and Z) by angles angZ, angY, and angX in that order to create the rotated set. The colours chosen for the bands are rather fanciful. In the case of Mars, the last band is shown as white to signify the polar cap. The latitude at which the cap starts is presumed to vary sinusoidally with amplitude 20o around mean latitude of 70o. It is assumed that the temperatures causing the caps lag 30o behind the position of the Sun.

A somewhat brute-force method is used with almost all the points on the surface of the ellipsoid being included and only those points with positive Z as seen from Earth, i.e. facing the Earth and also positive Z as seen from the Sun, being used in the drawing. This avoids drawing points on the far side of the planet or between the terminator and the edge.

The rings of Saturn are drawn in a similar way but with all points being in the plane of the equator and the radius varying from the inner edge to the outer edge of the inner ring and the same for the outer ring (only the two main rings are shown). Firstly, the shadow of the rings on the planet are drawn in dark grey and then the rings themselves are drawn and may override the shadow. Both rings are drawn in rather a brown colour with the outer ring being brighter than the inner ring. The exception is if the shadow of the planet falls on the ring, in which case it is drawn dark grey.

Apart from the bands of colour, the Great Red Spot is shown for Jupiter if it is visible. It is shown if its centre is less than 70o from the central meridian. It is assumed to cover 27o of longitude (Peek, p. 134), making it 0.15 of the major axis of Jupiter. It is also assumed to be an ellipse with vertical height 0.4 of its horizontal length and to be centred 0.4 of the major radius south of the equator. graphics.Fillellipse is used to draw it in colour brown. If it is not centred on the central meridian, its length is foreshortened accordingly.

Then, the satellites are entered onto the screen. Data about the satellites are obtained from array g_dblSat( , ) discussed above - phenomenon, mv, X and Y coordinates, X and Y coordinates of shadow. In the case of Saturn, they are only included if they exceed the limiting magnitude chosen for the particular telescope diameter. For example, Mimas, Enceladus, Hyperion, and perhaps Iapetus if its dark hemisphere is showing, may not be included. Firstly, a crude orbit is drawn in light blue. This is to give an idea of which satellite it is and how its position would change - useful for Saturn’s system with its possibly large inclination to the line of sight. The orbits are crude because they are assumed to be circles in the plane of the planet’s equator. Then, unless the satellite is occulted, it is shown in its (X,Y) position on the screen. Because the imagined orbit is crudely drawn, it is possible that the image of the satellite falls slightly off the line. If the satellite is eclipsed, it is shown in a dark grey just to show where it would be. Jupiter’s satellites are shown as white circles with their true apparent diameters. Phobos and Deimos are shown as single white points. For Saturn, the number of white points shown are proportional to the intensity, except for Titan which is shown to size, possibly with reduced brightness for each point to roughly approximate the true visual intensity. Shadows on the planet are shown as black circles of the correct apparent size at the correct positions.

If a mutual eclipse has occurred, this is shown graphically in the bottom right corner of the screen. The umbra and penumbra of the eclipsing satellite are shown, along with a red circle indicating the extent of the eclipsed satellite. The intensity of the penumbra varies from black to white in the correct manner from the inside to the outside. It is then possible to roughly estimate visually the intensity of the eclipse and thus try slightly different time to see how the red circle moves with respect to the shadow.

Technical Aspects

Planetary Perturbations

Chaper 24 of MeeusI give principal perturbations for the planets. I have included these for Jupiter, Saturn, Uranus, and Neptune. The file PlanetperturbII.txt contains the coefficients( numbers) for Jupiter and Saturn. The first 6 numbers are those for P, Q, and S on p110. There are then 36 numbers for longitude, 73 for eccentricity, 37 for perihelion, and 37 for semimajor axis of Jupiter, and the same for Saturn. Some of these numbers are zero because they correspond to non-zero terms in the other planet and coding is simplified (i.e. both planets can be processed with a single set of instructions) by covering all terms which appear in either planet. The coefficients for Uranus and Neptune are not included in the file but are coded directly into the program. VSOP87 MeusII, Chapter 31 (p.205) and Appendix II

An Excel spreadsheet, Meeus II APP II.xls has been written to implement the theory for the 6 planets Mercury to Saturn using the numbers presented by meeus in Appendix II. The calculations can be checked with MeeusII in at least 3 places. SatGraph can also be checked against the calculations done in the spreadsheet. covers the subject extensively and presents more terms than Meeus and claims greater accuracy.

(i ) Example 31.a on p. 207 of MeesII considers Venus at JD =2448976.5. This JD is entered into cell A411. Results for Venus are in column I, with L0 in row 84, L1 in row 135, L2 in row 172, L3 in row 189, L4 in row 202 and L5 in row 205. These are in agreement with MeeusII. Then L is in row 206, B in row 309 and R in row 411. These also agree with MeeusII.

(ii) Example 42.a on p. 280 considers Jupiter at JD = 2448972.50068. Results for Earth are in column N and for Jupiter are in column U of the spreadsheet.

For Earth, the results compared are

MeeusII SS
L 84o.285703 84o.28576
B 0o.000197 0o.000177
R 0.98412316 0.9841246

and for Jupiter are

MeeusII SS
L 188o.882168 188o.88221
B 1o.290464 1o.2955
R 5.44642320 5.4464430

(iii) Example 41.a on p. 275 considers Mars at JD = 2448935.500683. Results for Mars are in column Q.

For Earth, the results compared are

MeeusII SS
L 46o.843861 46o.84396
B - 0o.0001967 - 0o.000156
R 0.99041301 0.9904132

and for Mars are

MeeusII SS
L 78o.473759 78o.47373
B 0o.896321 0o.89640
R 1.5416954 1.5416535

I have made use of my own VSOP87 calculations for checking the articles by Drake and Kowal and by Nobili for the position of Neptune with respect to Jupiter in December 1612 and January 1613 recorded in Galileo’s journal. ( Refer to the Examples in the basic webpage for these references). Because the terms are truncated in AppII, the Ecliptic Longitude and Latitude calculated can be out by as much as about 5” for Jupiter and 4” for Neptune (see p 208 of MeeusII) but the accuracy is sufficient for the current task. In particular, the position of Neptune on 28 December 1612 and the surmised position by Nobili on 6 January 1613 are verified. This has been done by overriding the less-accurate calculations of Ecliptic Longitude and Ecliptic Latitude in SatGraph and proceeding to run the remainder in SatGraph. Of course, to do this, one need access to the source code.

Appearance of Mars and Jupiter

These are done in btnMarsView and btnJupiterView in SatGraph.

Refer to ESI, pp.334-337 and MeeusII, Chapter 41 for Mars and to ESI, pp338-340 and MeeusII, Chapter 42 for Jupiter. The methods used are different in the two books, so are worth comparing.

We are interested in finding De, Ds, and Longitude of Central Meridian (LCM) for each planet. In the case of Jupiter there are two LCMs.

ESI calculates these from mathematical expressions on p. 334 which use a0 and δ0 (see p. 334) and also use Δ and I which are given in Table 11.2 for Mars and Table 11.3 for Jupiter. However, SatGraph needs to calculate these and expressions on p. 333 can be used for this. Accordingly, I refer to expressions at the top of p.333 as 1(i) to 1(v) and expressions below these as 2(i) to 2(v). Similar expressions at the top of p. 334 are referred to as 3(i) to (v). To calculate:

N : 1(i)

I : 2(iii)

Δ : 2(iv)

Ω : ω from 1(iv) , Ω + ω from 2(i) and hence Ω

This concludes the calculations to create Tables 11.2 and 11.3.

Then Ds comes from the expression in the middle of p.334 and De from 3(iii).

In MarsView, SatGraph duplicates the calculations for De and Ds using MeeusII algorithm. This is done as a check.

The numerical results are the same and obviously one is redundant but they are currently left in for my benefit. They follow a series of steps in Chapter 41 and do not require the complications of finding Δ, Ω, and I. λ0 and β0 come from (41.1) on p.272 or from expressions (12.1) and (12.2) in MeeusII.

Yet another method is used for Jupiter in Chapter 42 of MeeusII and this is also done in SatGraph as a check. The results are the same but this coding is also left in.

Ae comes from dividing 2(iv) by 2(v) and hence LCM from the bottom of p.336 in ESI. This is also recalculated as per MeeusII.

JupiterView uses MeeusII, Chapter 42, calculations only and also finds P, the position angle.

Saturn’s satellite system

I do not understand the theory used to derive these expressions and I think it would be too time-consuming for me to learn it. Instead, I am merely trying to use it and to reconcile the expressions used by ESII with those of ESI.

It may be better to study Titan before Rhea because the longitude of perisaturnium of Rhea oscillates about that of Titan with a period of 38 years (ESI, p368). The amplitude of the oscillation is given as 17o.64 by ESI but it may be smaller, as discussed below.


See ESI, p368 and ESII, p360.

In both cases, the epoch used is April 0.0 1889 (1889.25, JD = 2411093), such that

d = JD - 2411093 and t = d/365.25 . However, ESI also uses year 1879.59 and ESII uses year 1890.0 (JD 2411368).

ESII uses the theory from Sinclair (1977) . See the reference below (under Titan).

Notation ESI ESII
Inclination of Saturn’s ring plane (equator) on ecliptic i1 ie
Ascending node of Saturn’s ring plane (equator) on ecliptic Ω1 Ωe
Inclination of Rhea’s Laplacian plane relative to ecliptic = i1 - 0o.0465 ie - 0o.0455
Node of Rhea’s Laplacian plane relative to ecliptic = Ω1 - 0o.0135 Ωe - 0o.0078
Longitude of perisaturnium of Rhea П1 ω
Longitude of Rhea L λ

Because of its high mass, Titan causes librations of Rhea. Hence, we require:

Longitude of perisaturnium of Titan П (p373) ωT

ESII refers to the “proper” longitude of perisaturnium of Rhea (and also the “proper” eccentricity, ep) and uses both ω p and π. See p 332 as well as p 360.

iLongitude of perisaturnium of Rhea oscillates about that of Titan in a period of 38 years (9.o5 per year according to ESI). In this calculation, both ESI and ESII appear to use the mean value of Titan, ignoring the perturbations that are taken into account when doing the calculations for Titan. ESII seems to use different values for ω0 for Rhea (276.o49 from Garcia ) than for Titan (275.o837 from Taylor and Shen). These differences probably make little difference to the accuracy of the calculations and each has probably been recently modified slightly, but they are annoying when trying to understand the reasoning. ESI also appears to use different values for each satellite and compounds the problem by using different epochs for t.

Differences in constants used ESI ESII
Starting longitude (JD 2411093) , E0 358.o3950 359.o4727
Degrees longitude/day, n 79o.6900881 79.o6900400700
N for Titan, NT 48.o5 - 0.o50t 44.o5 - 0.o5219 (t - 0.75)
= 44.o89 - 0.o5219t
θ 344.o09 - 10.o20t 356.o87 - 10.o2077t
γ0 20’.4960 = 0.o3415 0.o3305

The main differences between ESI and ESII are for eccentricity, e, and longitude of perisaturnium . ESI gives straightforward expressions, viz.

e = 0.00098 + 0,00030cos9o.5(τ - 1879.59)

Π1 = 276o.25 +0o.53t + 17o.64sin9o.5(τ - 1879.59)

( the first 2 terms should agree with the corresponding 2 terms of Π for Titan on p 373)

Alternatively, ESII requires that similar expressions be derived from the following relationships

(see Sinclair(77), p450 and ESII, p332 and p 360) :

e sinω = ef sinωT + ep sinωp

e cosω = ef cosωT + ep cosωp ,

where ωT is the longitude of pericentre of Titan, ωp is the “proper” ω for Rhea in the absence of Titan, ω is the resultant value for Rhea, ef is the forced eccentricity due to Titan, ep is the proper eccentricity, and e is the resultant eccentricity. ESII gives ef = 0.00100 and ep = 0.00021 .

The above pair of equations gives

e = (ef^2 + ep^2) ½ [1 + (2efep)/( ef^2 + ep^2)cos(ωp - ωT)]½

= 0.001022[1 + 0.4023 cos(ωp - ωT)] ½                      

≈ 0.001022[1 + 0.20 cos(ωp - ωT)]

 ≈  0.001022 + 0.000204 cos(ωp - ωT)

ESI gives e = 0.00098 + 0.00030cos(ωp - ωT)

Also, sin(ω - ωT ) = (ep/e)sin(ωp - ωT)

                           ≈  0.205 sin(ωp - ωT)

So ω - ωT ≈ 0.205 κ sin(ωp - ωT) , where κ =180/π

                          =  11o.77 sin(ωp - ωT) 

This is considerably smaller than 17o.64 used in ESI. The reason is that ESI uses ep/e = 0.306 as against 0.205.

The rate of change per year of ωp - ωT is 10o.2077 - 0o.5219 = 9o.6858 (from Sinclair and ESII) and almost the same from Garcia. This agrees closely with 9o.5 in ESI and is in agreement with the period of 38 years. However, the constant in ωp - ωT , which determines the epoch at which sin(ωp - ωT) = 0, varies greatly in the literature. Garcia uses Π0 for Rhea =105o , while Taylor and Shen, and thereby ESII, uses Π0 =305o . Neither seems to agree with τ =1879.59 in ESI.

ΦωΣΨπΩΘΔγ≈√αγδПδλ ΓκΠτ½τ


See ESI, pp373-378 and ESII, pp361-363.

In both cases, the epoch used is 1889.0 ( JD = 2411368), such that

d = JD - 2411368 and t = d/365.25 .

ESII uses the theory from Sinclair (1977), viz.

“The orbits of Tethys, Dione, Rhea, Titan and Iapetus ” by A.T. Sinclair ; Monthly Notes of the Royal Astronomical Society (1977) 180, 447-459 .

ESII is basically a copy of Sinclair (77), pp451-452, with the orbital elements taken from Taylor and Shen (1988) and the mean motions and secular rates from Garcia (1972) .

Thus, ESII uses the following values for the constants used by Sinclair (the values in brackets are used by ESI):

γ0 0o.2990 (0o.332914)
n 22o.57697385/day ( 22o.57701508/day )
N0 41o.28 ( 40o.69 )
N’ -0o.5219/year ( -0o.506/year )
ω0 275o.837 ( 276o.1283 )
ω’ 0o.5219/year ( 0o.5235/year )
a0 0.00816765AU
e0 0.028815 (0.02910)
λ0 261o.3121 (260o.4043)

Coefficients of periodic terms also vary between ESI and ESII ( taken from Sinclair(77) ). The algebraic expressions, as well as the numeric values, of many of these expressions are given by Sinclair. ESII gives only the numerical values. ESI gives algebraic expressions (p375) but evaluates them only in Example 12.11 (p 376) .

Sinclair uses the following values in evaluating the coefficients: Γ=26.o1 , ns=0.o03346/day ,

n=22o.577/day , e=0.0292 , es=0.05568 , ω’=0o.5213/year , ie=28o.072, γ0 =0o.02990

                          Sinclair                                                                                    ESI

cos(N0+N’t) in ia: κsinγ0 = 0o.2990 18.35’ = 0o.30583

sin(N0+N’t) in Ωa: κsinγ0/sinie = 0o.63551 39’ = 0o.6500

cos(2g) in e : - (1516)nsesin2 Γ/ ω’ = -0.000340 -0.000186

                          -.000184 used   (mistake in formula ?)

cos2(Ls-g) in e: (1532)(ns/n)e(1+cos Γ)^2 =0.000073 (158)( n0/n)e = 0.0000804

Note that ESI expression agrees with Sinclair if Γ taken to be 0 ( cos Γ=1)

sin(2g) in ω: (1516)nssin2 Γ/ ω’ = 0.011646 radians 22’ = 0.00640 radians (p373)

                 0.00630 radians used by Sinclair (mistake in formula?)

sin2(Ls-g) in ω: (1532)(ns/n) (1+cos Γ)^2 = 0.00250 radians (158)( n0/n) = 0.00278 radians Note that again ESI expression agrees with Sinclair if Γ taken to be 0 ( cos Γ=1)

sin(N0+N’t) in λ: κsinγ0tanie/2 = 0o.08321 4.39’ = 0o.07317

sin(ls) in λ: - ( 32)(ns/n)es(3cos 2Γ-1) = -0.000176 rads -3(ns/n) es = -0.000248 rads

 Note that again ESI expression agrees with Sinclair if  Γ taken to be 0 ( cos Γ=1)

sin(2Ls) in λ: -(34)(ns/n)sin2 Γ = -0.000215 rads -(916) (ns/n)sin2 Γ=-0.000161 rads

              Why the extra factor of 3/4 in ESI?

sin(2Ls+Ψ) in Ω: (316)(ns/n)sin Γ(1+cos Γ)/sinie = 0.0000503 rads (38)(ns/n)sin Γ/sini =0.000532 rads

Again ESI expression agrees with Sinclair if Γ taken to be 0 ( cos Γ=1) and i = ie

sin(2Ls+Ψ) in λ: +(316)(ns/n)sin Γ(1+cos Γ)tan(ie/2) = 0.000057 rads - (38)(ns/n) sin Γtan(ie/2) rads

This is obtained from the coefficient in Ω by multiplying by 2sin2(ie/2) (refer to ΔL on p375 of ESI).

Note that Sinclair has a + sign and ESI a - sign. I do not know which is correct.

cos(2Ls+ Ψ) in i: (316)(ns/n)sin Γ(1+cos Γ) =.0.000232 rads (38)(ns/n)sin Γ = 0.000247

ESI has 2 further perturbation terms in λ (in ΔE on p375), which Sinclair (and therefore ESII does not include. The values are quite small. These coefficients are :

sin2(l0- Π0) : -(94)(n0/n)e0^2 = -0.00001044 radians

sin2(l0- Π) : -(4516) (n0/n)e^2 = -0.000003534 radians


ESI, p378 and ESII, p363 agree closely in longitude and П except that ESII has a few extra minor terms.


I am using ESII, pp. 364-6 without understanding why. The mathematical expressions are quite lengthy but generally the perturbations are small.

However, Δω could possibly be large if all the sine terms are + or - 1.