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.
Introduction
SatGraph program
Main.vb
The source code discussed here can be obtained from the file main.txt
frmPlanet Form
Shared
Form1_Load
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.txt are input 17 values for each of the 8 planets. The contents of this file and the corresponding array, m_dblPlanets(7 ,16) are described below in the section “Routine btnCalculateSolar”.
. It inputs data from any of the 8 files VSOP87D.xxx.txt, one for each planet.
Subroutine Longitude2
This is called in btnCalculateSolar to calculate the longitude and distance from the Sun of the Earth and the planets. Now that the VSOP87 theory has been used in Satgraph, Longitude2 is used only for the final calculation for the planet (and it may not be necessary) and even this has been overridden by the subroutine VSOP87.
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 semi-major 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.
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.
The change in Dynamic time can be entered (in seconds) or roughly estimated by the program. It is added to the Julian Date.
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(7,16) are used by this routine. Hence, the data are described here. The array covers 8 planets.
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 of equator per century |
VSOP87. The subroutine VSOP87 calculates the ecliptic longitude and latitude of the Earth.
Subroutine GeoCoords is used to calculate geocentric longitude, latitude, projected radius, and “ radius” i.e. distance , as seen from Earth. The distance is then used to calculate the delay due to the finite speed of light, or “light time”.
The next two paragraphs are obsolete. Instead VSOP87 is used as discussed below.
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 semi-major 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 again called. This has been left in (at least for now, as it outputs values which are still used elsewhere).
SubroutineVSOP87 is called, followed by subroutine GeoCoords.
Subroutine VSOP87
The VSOP87 theory is discussed in Meeus II, Chapter 31.
I will find more-original references later.
VSOP87 seems to be the theory of choice now. It is used in the Astronomical software program STELLARIUM and Jay Tranter supplies information on facilitating its use.
The ‘raw’ files, which contain hundreds of thousands of numbers can be obtained at :
https://github.com/cdtk/VSOP87/blob/master/VSOP87x.ppp where x is A, B, C, or D and ppp is mer, ven, ear, mar, jup, sat, ura, or nep.
I am using VSOP87D files.
Programming can be tricky as it involves essentially descrambling a single string containing all the data for the particular planet.
(To be completed!!!!!!!!!!!!)
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.
frmLocalCoords_Load
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.
btnMarsView
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 which (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.
btnJupiterView
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 set up (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.
btnSaturnsRings
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.
frmSatellites
Satellites_Load
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 of longitude 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, Jupiter, or Saturn, elements 5, 0, 1, and 2 are calculated . Each planet and its satellite system is treated differently as discussed below.
Mars
Satellites are treated in Chapter 12 of ESI from 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.
Jupiter
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);
Current values for elements 0, 1, and 2 are loaded from the satellite array m_dblSatellites for each of the 4 satellites.
| (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 |
Φ, λ, Ψ, 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 .
Saturn
These calculations are the most complicated of all the satellite systems. They are treated in ESI on pp.366-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.
Rhea
ESI, p368
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.
Titan
ESI, p373
The main programming complications are :
(i) The solving of the spherical trigonometry formulae on p. 375 of ESI 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.
Hyperion
ESI, p378
Δe, ΔL, and ΔП are calculated as on p.378 ESI or pp. 363-4 of ESII and then subroutines Longitude2 and ConverttoPlane are called.
Iapetus
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.
Satellite Phenomena
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 in array element g_dblsat( ,10)
If the satellite is in front of the planet and its separation from the centre is 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 in array element g_dblsat( ,10) . Limits of the shadow are not calculated.
Mention is made here of 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.
Mutual satellite phenomena
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 of the 2 satellites.
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 occurrence is recorded in g_dblSat( , ) and the data is stored in some shared variables for use by the graphics part of the program.
Graphics routines
Graphics.vb
GraphicEllipse.vb
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 is 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. Also, if the rings are deemed to be invisible (because the Sun is shining on the other side from that seen from Earth), they are 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 70 degrees from the central meridian. It is assumed to cover 27degrees 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.
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 a slightly different time to see how the red circle moves with respect to the shadow.
Technical Aspects
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. The satellite “number” ( 1 to 4 for Jupiter,1 to 8 for Saturn) is shown above the satellite. 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.
VSOP87
Chapter 31 of Meeus II concerns the planetary theory VSOP87 and Appendix II contains the most important terms for this theory.
Satgraph includes the full list of terms given in the theory.
The relevant files can be found at
https://github.com/cdtk/VSOP87/blob/master/VSOP87D.ppp, where ppp is
mer,ven,ear,mar,jup,sat,ura, or nep.
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.
http:www.neoprogrammics.com/VSOP87/source_code_generator_tool/index.php 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.
Meeus uses the complete VSOP87 theory.
For Earth, the results compared are
| Item | MeeusII / Spreadsheet | |
|---|---|---|
| L | 84o.285703 degrees | 84o.28576 |
| B | 0o.000197 degrees | 0o.000177 |
| R | 0o.98412316 | 0o.9841246 |
and for Jupiter are
| Item | MeeusII / Spreadsheet | |
|---|---|---|
| L | 181o.882168 | 181o.88221 |
| B | 1o.290464 | 1o.2955 |
| R | 5o.44642320 | 5o.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
| Item | MeeusII / Spreadsheet | |
|---|---|---|
| L | 46o.843861 degrees | 46o.84396 |
| B | – 0o.000167 degrees | 0o.000156 |
| R | 0o.99041301 | 0o.9904132 |
and for Mars are
| Item | MeeusII / Spreadsheet | |
|---|---|---|
| L | 78o.473759 | 78o.47373 |
| B | 0o.896321 | 0o.89640 |
| R | 1o.5416954 | 1o.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.
Rhea
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 wp and π. See p332 as well as p360.
Longitude 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 | 358o.3950 | 359o.4727 |
| Degrees longitude/day, n | 79o.6900881 | 79o.6900400700 |
| N for Titan, NT | 48o.5 – 0o.50t | 44o.5 – 0o.5219 (t – 0.75) = 44o.89 – 0.o5219t |
| θ | 344o.09 – 10o.20t | 356o.87 – 10o.2077t |
| γ0 | 20′.49/60 = 0o.3415 | 0o.3305 |
The main differences between ESI and ESII are for eccentricity, e, and longitude of perisaturnium .
ESI gives straightforward expressions, viz.
e = 0.00098 + 0.00030cos9.5(τ – 1879.59)
Π1 = 276o.25 +0o.53t + 17o.64sin9.5(τ – 1879.59)
( the first 2 terms should agree with the corresponding 2 terms of Π for Titan on p373)
Alternatively, ESII requires that similar expressions be derived from the following relationships
(see Sinclair(77), p450 and ESII, p332 and p360) :
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 = (ef2 + ep2)½ [1 + (2ef.ep)/( ef2 + ep2)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/π
= 11.77 sin(ωp - ωT) degrees
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.
ΦωΣΨπΩΘΔγ≈√αγδПδλ ΓκΠτ½τ
Titan
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 | 0.2990o | (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: Γ=26o.1 , ns=0o.03346/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 : – (15⁄16)nsesin2Γω’ = -0.000340 -0.000186
-.000184 used (mistake in formula ?)
cos2(Ls-g) in e: (15⁄32)(ne/n)e(1+cos Γ)2 =0.000073 (15⁄8)( n0/n)e = 0.0000804
Note that ESI expression agrees with Sinclair if Γ taken to be 0 (cos Γ=1)
sin(2g) in ω: (15⁄16)nssin2Γ/ω’ = 0.011646 radians 22′ = 0.00640 rads (p373)
0.00630 radians used by Sinclair (mistake in formula?)
sin2(Ls-g) in ω: (15⁄32)(ns/n) (1+cos Γ)2 = 0.00250 rads (15⁄8)(n0/n) = 0.00278 rads 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 λ: – ( 3⁄2)(ns/n)es(3cos2Γ-1) = -0.000176rads -3(ns/n) es= -0.000248rads
Note that again ESI expression agrees with Sinclair if Γ taken to be 0 ( cos Γ=1)
sin(2Ls) in λ: -(3⁄4)(ns/n)sin2Γ = -0.000215rads -(9⁄16) (ns/n)sin2Γ=-0.000161rads
Why the extra factor of 3/4 in ESI?
sin(2Ls+Ψ) in Ω: (3⁄16)(ns/n)sinΓ(1+cos Γ)/sinie = 0.0000503 rads (3⁄8)(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 λ: +(3⁄16)(ns/n)sin Γ(1+cos Γ)tan(ie/2) = 0.000057 rads – (3⁄8)(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: (3⁄16)(ns/n)sin Γ(1+cos Γ) =.0.000232 rads (3⁄8)(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) : -(9⁄4)(n0/n)e02 = -0.00001044 radians
sin2(l0– Π) : -(45⁄16) (n0/n)e2 = -0.000003534 radians
Hyperion
ESI, p378 and ESII, p363 agree closely in longitude and П except that ESII has a few extra minor terms.
Iapetus
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.
¼καβγδ½¾νχ