Astro::Coord::ECI::VSOP87D - Implement the VSOP87D planetary position model.
package My_Object; use base qw{ Astro::Coord::ECI }; use Astro::Coord::ECI::VSOP87D qw{ :mixin }; sub __model_definition { return [ # VSOP87D coeficients as digested from the distro by # tools/decode ]; }
This Perl module implements the VSOP87D model of planetary motion in a manner that is consistent with the Astro::Coord::ECI hierarchy.
The VSOP87 models calculate the positions of the planets through Neptune as a series of cosine terms. They were created by Pierre Bretagnon and Gerard Francou of the Bureau des Longitudes in Paris France. Depending on the model you select, you can calculate either instantaneous orbital parameters, Heliocentric ecliptic position in either spherical or Cartesian coordinates and either current or J2000.0 equinox, or barycentric ecliptic position.
The models are stated to be within 1 second of arc for varying times around J2000 (about noon January 1 2000 UT), as follows:
4000 years for Mercury, Venus, Earth-Moon barycenter, Mars 2000 years for Jupiter and Saturn 6000 years for Uranus and Neptune
The actual model coefficients are available from http://cdsarc.u-strasbg.fr/viz-bin/qcat?VI/81/, as is additional information the accuracy of the model, and a reference implementation in Fortran.
VSOP87D provides Heliocentric spherical coordinates referred to the current equinox. This is consistent with the existing members of the Astro::Coord::ECI hierarchy, but more importantly Jean Meeus' Astronomical Algorithms provides worked examples for the Sun and Venus.
The VSOP87D model itself is simple enough to implement, and occupies about 50 lines of Perl. Testing shows that the Perl solution for a given body and time is within 5e-9 of the Fortran reference implementation over selected years in the range -4000 to 2300. The units of that 5e-9 vary, being radians for longitude and latitude, AU for radius, and radians/day and AU/day for the corresponding velocities.
5e-9
One of the advantages of the VSOP87 models over the previous VSOP82 models is that they are easily truncated when less than full precision is required. Meeus has his own truncation, which is the default one. In addition, this code provides truncation 'none', which uses the whole model, and a mechanism to define custom truncations.
'none'
As the examples make clear, the devil is in the details, and a significant amount of work needs to be done to get from Heliocentric geometric ecliptic coordinates to Geocentric apparent equatorial coordinates. Meeus' examples are embodied in t/meeus_25b.t and t/meeus_33a.t.
For the Sun, the Perl implementation agrees pretty well with the worked example, differing by 0.001 seconds of right ascension (not seconds of arc), 0.007 seconds of declination (seconds of arc) and 1e-7 AU in radius.
0.001
0.007
1e-7
For Venus, the Perl implementation agrees less well with the worked example, differing by 0.14 seconds of right ascension, 0.07 seconds of declination, and 5e-8 AU in radius. I have gone over the calculation step-by-step, and as of this writing most of the difference seems to come in at the computation of Geocentric ecliptic position from its Heliocentric Cartesian ecliptic position. This is simply a square root and two atan2() calls, and if it is subtly wrong some way I can not see it.
0.14
0.07
5e-8
atan2()
Because these are the only two worked examples I found in Meeus, I felt the need to sanity-check the rest of the distribution in other ways. These tests are contained in the t/planet_*.t test files.
Rise and set times came from the United States Naval Observatory (USNO) Rise/Set/Transit Times for Major Solar System Bodies and Bright Stars page at http://aa.usno.navy.mil/data/docs/mrst.php, and are given to the nearest minute. These models almost always reproduce the USNO's times in the cases tested. I have no idea how the USNO data are calculated or to what accuracy.
Conjunctions and so forth came from Guy Ottewell's Astronomical Calendar at http://www.universalworkshop.com/, and are given to the nearest hour. These models usually give the same. The 2018 edition also gives JDE dates for the phenomena to three decimal places, or an accuracy of about 90 seconds. These models typically are within about 10 minutes of the Ottewell figures. Again I have no idea how the Astronomical Calendar data are calculated or to what accuracy.
This specific module is not intended to be used directly. Instead, subclass a member of the Astro::Coord::ECI hierarchy and import methods from this module. In most cases this means whatever the :mixin tag provides, though for the Sun the :sun tag should be used instead. You must also provide a __model_definition() method that returns the model parameters you wish to implement. For the Sun, a __model() method is also required, since the design of this distribution hangs the Earth's __model_definition() information on the Sun.
:mixin
:sun
But the end user does not do this either. Instead this distribution provides three such subclasses:
This is a subclass of Astro::Coord::ECI::Sun. The __model() method (private to this distribution) simply returns zeroes. But this class carries the model parameters for the Earth, which are subtracted from whatever __model() returns to get Geocentric coordinates.
This is a subclass of Astro::Coord::ECI, which pulls in everything relevant from this module (i.e. imports :mixin), and provides almanac methods appropriate to inferior planets. Mercury and Venus are subclassed from this.
This is a subclass of Astro::Coord::ECI, which pulls in everything relevant from this module (i.e. imports :mixin), and provides almanac methods appropriate to superior planets. Mars, Jupiter, Saturn, Uranus, and Neptune are subclassed from this.
This module does not implement a class, and is better regarded as a mixin. It supports the following methods, which are public and exportable unless documented otherwise.
This method reports, creates, and deletes model cutoff definitions.
The first argument is the name of the model cutoff. If this is the only argument, a reference to a hash defining the named model cutoff is returned. This return is a deep clone of the actual definition.
If the second argument is undef, the named model cutoff is deleted. If the model cutoff does not exist, the call does nothing. It is an error to try to delete built-in cutoffs 'none' and 'Meeus'.
undef
'Meeus'
If the second argument is a reference to a hash, this defines or redefines a model cutoff. The keys to the hash are the names of VSOP87D series ('L0' through 'L5', 'B0' through 'B5', and 'R0' through 'R5'), and the value of each key is the number of terms of that series to use. If one of the keys is omitted or has a false value, that series is not used.
'L0'
'L5'
'B0'
'B5'
'R0'
'R5'
If the second argument is a scalar, it is expected to be a number, and a model cutoff is generated consisting of all terms whose coefficient ('A' in Meeus' terminology) is equal to or greater than the number.
'A'
If the second argument is a code reference, this code is expected to return a reference to a valid model cutoff hash as described two paragraphs previously. Its arguments are the individual series hashes, extracted from the model. Each hash will have the following keys:
The name of the series (e.g. 'L0').
An array reference containing the terms of the series. Each term is a reference to an array containing in order, in Meeus' terms, values 'A', 'B', and 'C'.
'B'
'C'
This method is exportable, either by name or via the :mixin or :sun tags.
This method returns the geometric longitude of the body. This is after correction for aberration and light time (for bodies other than the Sun) and conversion to FK5.
This method is exportable, either by name or via the :sun tag.
my ( $delta_psi, $delta_epsilon ) = $self->nutation( $dynamical_time, $cutoff );
This method calculates the nutation in ecliptic longitude ($delta_psi) and latitude ($delta_epsilon) at the given dynamical time in seconds since the epoch (i.e. Perl time), according to the IAU 1980 model.
$delta_psi
$delta_epsilon
The $time argument is optional, and defaults to the object's current dynamical time.
$time
The $cutoff argument is optional; if specified as a number larger than 0, terms whose amplitudes are smaller than the nutation cutoff (in milli arc seconds) are ignored. The Meeus version of the algorithm is specified by a value of 3. The default is specified by the nutation_cutoff attribute.
$cutoff
0
3
nutation_cutoff
The model itself is the IAU 1980 nutation model. Later models exist, but this was chosen because of the desire to be compatible with Meeus' examples. The implementation itself actually comes from Meeus chapter 22. The model parameters were not transcribed from that source, however, but were taken from the source IAU C reference implementation of the algorithm, src/nut80.c, with the minimum modifications necessary to make the C code into Perl code. This file is contained in http://www.iausofa.org/2018_0130_C/sofa_c-20180130.tar.gz.
$epsilon = $self->obliquity( $time );
This method calculates the obliquity of the ecliptic in radians at the given dynamical time. If the time is omitted or specified as undef, the current dynamical time of the object is used.
The algorithm is equation 22.3 from Jean Meeus' "Astronomical Algorithms", 2nd Edition, Chapter 22, pages 143ff.
say 'Order from Sun: ', $self->order();
This method returns the order of the body from the Sun, with the Sun itself being 0. The number 3 is skipped, since that would represent the Earth.
$self->period()
This method returns the sidereal period of the object, calculated from the coefficient of its first L1 term.
L1
The algorithm is the author's, and is a first approximation. That is. it is just the tropical period plus however long it takes the object to cover the amount of precession during the tropical year.
$self->synodic_period()
This method returns the synodic period of the object -- that is to say the mean interval between oppositions or conjunctions of superior planets or between corresponding conjunctions of inferior planets.
This method is exportable, either by name or via the :mixin tag.
$self->time_set()
This method is not normally called by the user. It is called by Astro::Coord::ECI to compute the position once the time has been set.
It returns the invocant.
$self->year()
This method returns the length of the tropical year of the object, calculated from the coefficient of its first L1 term.
This method is private to this distribution, and may be changed or revoked without notice at any time. Documentation is for the benefit of the author.
This method returns the hash containing VSOP87D-specific attributes, creating it if necessary.
my @state_vector = $class->__model( $time_dynamical, ... );
This static method is private to this package, and may be changed or revoked without notice at any time. Documentation is for the benefit of the author.
This static method executes the VSOP87D model returned by $class->__model_definition(), and returns the computed state vector. The components of the state vector are Heliocentric ecliptic longitude and latitude in radians, radius in AU, and the associated velocities in radians/day and AU/day.
$class->__model_definition()
Additional optional arguments can be specified as name/value pairs. The following are defined:
If this Boolean argument is true, whatever debug data the author finds useful at the moment is written to standard error. Only false values of this argument are supported.
This is the model cutoff definition hash to use in the computation. If unspecified the whole model is used.
This method is exportable, either by name or via the :mixin tag. It is not exported by the :sun tag.
my $model = $class->__model_definition( 'model' )
This static method returns model-related information. The argument describes the information to return. The following arguments are valid:
This argument returns the default value of 'model_cutoff_definition' for a new object.
'model_cutoff_definition'
This argument returns the terms of the VSOP87D model, expressed as nested array references.
The top level is indexed on the zero-based coordinate index (i.e. one less than the indices documented in vsop87.txt. The next level is indexed on the exponent of T represented, and is expected to be in the range 0-5.
T
0-5
The third level is a hash describing the individual series. It contains the following keys:
This is the name of the individual series (e.g. 'L0'.
This is a reference to an array containing the terms of the series. Each term is a reference to an array containing quantities A, B, and C to be plugged into the equation A * cos( B + C * T) * T ** n.
A * cos( B + C * T) * T ** n
T is dynamical time in Julian millennia since J2000.0.
This is the sidereal period in seconds, and is returned by period().
This is the tropical period, in seconds, and is returned by year().
This method is not exportable. It is expected that each derived class will implement its own version of this method.
Astro::Coord::ECI
https://www.caglow.com/info/compute/vsop87
The author wishes to acknowledge and thank the following individuals and organizations.
Pierre Bretagnon and Gerard Francou of the Bureau des Longitudes, Paris, France, whose VSOP87 solutions of planetary motion form the basis of these modules. Both the coefficients for the various VSOP87 models and a reference implementation in FORTRAN are available from http://cdsarc.u-strasbg.fr/viz-bin/qcat?VI/81/. Without these, this module would not exist.
Jean Meeus, whose book "Astronomical Algorithms" (second edition) contained invaluable worked examples for the Sun and Venus, which eventually became t/meeus_25b.t and t/meeus_33a.t. Without these I would never have been able to reduce the Heliocentric ecliptic coordinates generated by the VSOP87D models to Geocentric equatorial coordinates.
The International Astronomical Union (IAU), whose SOFA software collection provided both coefficients and a reference implementation for the nutation calculation. Formal credit is given them under COPYRIGHT AND LICENSE.
Guy Ottewell, whose Astronomical Calendar has long been part of my life. Paper publication of this ceased with the 2016 edition, but it continues on line at http://www.universalworkshop.com/, and formed the basis of the conjunction portion of the t/planet_*.t.
The United States Naval Observatory, whose Data Services, available at http://aa.usno.navy.mil/data/index.php, have long been my reference for daily events involving astronomical bodies.
Support is by the author. Please file bug reports at https://github.com/trwyant/perl-Astro-Coord-ECI-VSOP87D/issues, or in electronic mail to the author.
Thomas R. Wyant, III wyant at cpan dot org
Copyright (C) 2018-2021 by Thomas R. Wyant, III
This program is free software; you can redistribute it and/or modify it under the same terms as Perl 5.10.0. For more details, see the full text of the licenses in the directory LICENSES.
This program is distributed in the hope that it will be useful, but without any warranty; without even the implied warranty of merchantability or fitness for a particular purpose.
Software Routines from the IAU SOFA Collection were used. Copyright (C) International Astronomical Union Standards of Fundamental Astronomy (http://www.iausofa.org)
To install Astro::Coord::ECI::VSOP87D, copy and paste the appropriate command in to your terminal.
cpanm
cpanm Astro::Coord::ECI::VSOP87D
CPAN shell
perl -MCPAN -e shell install Astro::Coord::ECI::VSOP87D
For more information on module installation, please visit the detailed CPAN module installation guide.