standish_module – fortran-astrodynamics-toolkit

Approximate positions of the major planets.

Reference

E.M. Standish, Solar System Dynamics Group, JPL/Caltech, "Keplerian Elements for Approximate Positions of the Major Planets"

History

Original version copyright 2018 https://github.com/CumuloEpsilon BSD License.
Jacob Williams, extensive refactoring with some modifications, and integration into the FAT ephemeris module.

Original license

   Copyright (c) 2018, CumuloEpsilon
   All rights reserved.

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions are met:

   * Redistributions of source code must retain the above copyright notice, this
     list of conditions and the following disclaimer.

   * Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
   FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
   DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
   SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
   CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Uses

Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Used by

Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Variables

Type	Visibility	Attributes		Name		Initial
real(kind=wp),	private,	parameter	::	obliquity	=	23.43928_wp	obliquity at J2000 [deg]
real(kind=wp),	private,	parameter	::	s_sobl	=	sin(obliquity*deg2rad)	sin of j2000 obliquity
real(kind=wp),	private,	parameter	::	s_cobl	=	cos(obliquity*deg2rad)	cos of j2000 obliquity
real(kind=wp),	private,	parameter	::	epoch	=	2451545.0_wp	Julian date of J2000 epoch
type(ephem),	private,	parameter	::	eph1	=	ephem(1, 'Keplerian Elements Valid 1800AD-2050AD', [2378497.0_wp, 2470172.0_wp], reshape([0.38709927_wp, 0.20563594_wp, 0.12225995_wp, 4.4025989_wp, 1.3518935_wp, 0.84353095_wp, 3.70000009e-07_wp, 1.90600003e-05_wp, -1.03803286e-04_wp, 2608.7903_wp, 2.80085020e-03_wp, -2.18760967e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.72333568_wp, 6.77671982e-03_wp, 5.92482723e-02_wp, 3.1761343_wp, 2.2968962_wp, 1.3383157_wp, 3.90000014e-06_wp, -4.10700013e-05_wp, -1.37689030e-05_wp, 1021.3286_wp, 4.68322469e-05_wp, -4.84667765e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0000026_wp, 1.67112295e-02_wp, -2.67209913e-07_wp, 1.7534375_wp, 1.7966015_wp, 0.0_wp, 5.62000014e-06_wp, -4.39200012e-05_wp, -2.25962198e-04_wp, 628.30756_wp, 5.64218918e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.5237104_wp, 9.33941007e-02_wp, 3.22832055e-02_wp, -7.94723779e-02_wp, -0.41789517_wp, 0.86497712_wp, 1.84700002e-05_wp, 7.88199977e-05_wp, -1.41918135e-04_wp, 334.06131_wp, 7.75643345e-03_wp, -5.10636950e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 5.2028871_wp, 4.83862385e-02_wp, 2.27660220e-02_wp, 0.60033119_wp, 0.25706047_wp, 1.7536005_wp, -1.16069998e-04_wp, -1.32529996e-04_wp, -3.20641411e-05_wp, 52.966312_wp, 3.70929041e-03_wp, 3.57253314e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 9.5366764_wp, 5.38617894e-02_wp, 4.33887430e-02_wp, 0.87186599_wp, 1.6161553_wp, 1.9837835_wp, -1.25059998e-03_wp, -5.09909994e-04_wp, 3.37911442e-05_wp, 21.336540_wp, -7.31244357e-03_wp, -5.03838016e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 19.189165_wp, 4.72574383e-02_wp, 1.34850740e-02_wp, 5.4670362_wp, 2.9837148_wp, 1.2918390_wp, -1.96175999e-03_wp, -4.39700016e-05_wp, -4.24008576e-05_wp, 7.4784222_wp, 7.12186471e-03_wp, 7.40122399e-04_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 30.069923_wp, 8.59048031e-03_wp, 3.08930874e-02_wp, -0.96202600_wp, 0.78478318_wp, 2.3000686_wp, 2.62910005e-04_wp, 5.10499995e-05_wp, 6.17357864e-06_wp, 3.8128369_wp, -5.62719675e-03_wp, -8.87786155e-05_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 39.482117_wp, 0.24882729_wp, 0.29914966_wp, 4.1700983_wp, 3.9107401_wp, 1.9251670_wp, -3.15960002e-04_wp, 5.17000008e-05_wp, 8.40899645e-07_wp, 2.5343544_wp, -7.09117157e-04_wp, -2.06556579e-04_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp], [16, 9]))	The first ephemeris data table (this is standish's table 1). keplerian elements valid 1800 ad - 2050 ad
type(ephem),	private,	parameter	::	eph2	=	ephem(2, 'Keplerian Elements Valid 3000BC-3000AD', [625674.0_wp, 2816788.0_wp], reshape([0.38709843_wp, 0.20563661_wp, 0.12227069_wp, 4.4026222_wp, 1.3518922_wp, 0.84368551_wp, 0.0_wp, 2.12300001e-05_wp, -1.03002007e-04_wp, 2608.7903_wp, 2.78205727e-03_wp, -2.13177688e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.72332102_wp, 6.76399004e-03_wp, 5.93023673e-02_wp, 3.1761451_wp, 2.2997777_wp, 1.3381896_wp, -2.60000007e-07_wp, -5.10700011e-05_wp, 7.59113527e-06_wp, 1021.3286_wp, 9.91285546e-04_wp, -4.76024114e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0000002_wp, 1.67316291e-02_wp, -9.48516663e-06_wp, 1.7534785_wp, 1.7964685_wp, -8.92317668e-02_wp, -2.99999989e-08_wp, -3.66099994e-05_wp, -2.33381579e-04_wp, 628.30762_wp, 5.54932002e-03_wp, -4.21040738e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.5237124_wp, 9.33651105e-02_wp, 3.23203318e-02_wp, -7.97289312e-02_wp, -0.41743821_wp, 0.86765921_wp, 9.69999974e-07_wp, 9.14900011e-05_wp, -1.26493964e-04_wp, 334.06125_wp, 7.89301097e-03_wp, -4.68663359e-03_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 5.2024803_wp, 4.85358983e-02_wp, 2.26650927e-02_wp, 0.59925520_wp, 0.24914493_wp, 1.7504400_wp, -2.86400009e-05_wp, 1.80260002e-04_wp, -5.63216017e-05_wp, 52.969063_wp, 3.17635899e-03_wp, 2.27322499e-03_wp, -2.17328397e-06_wp, 1.05837814e-03_wp, -6.21955749e-03_wp, 0.66935557_wp, 9.5414991_wp, 5.55082485e-02_wp, 4.35327180e-02_wp, 0.87398607_wp, 1.6207365_wp, 1.9833919_wp, -3.06500006e-05_wp, -3.20440013e-04_wp, 7.88834659e-05_wp, 21.329931_wp, 9.45610274e-03_wp, -4.36594151e-03_wp, 4.52022823e-06_wp, -2.34475732e-03_wp, 1.52402408e-02_wp, 0.66935557_wp, 19.187979_wp, 4.68574017e-02_wp, 1.34910680e-02_wp, 5.4838729_wp, 3.0095420_wp, 1.2908891_wp, -2.04550000e-04_wp, -1.54999998e-05_wp, -3.14429781e-05_wp, 7.4786506_wp, 1.61739404e-03_wp, 1.00176642e-03_wp, 1.01806800e-05_wp, -1.70574244e-02_wp, 3.08735552e-03_wp, 0.13387112_wp, 30.069527_wp, 8.95438995e-03_wp, 3.08932904e-02_wp, 5.3096914_wp, 0.81474739_wp, 2.3001058_wp, 6.44699976e-05_wp, 8.17999990e-06_wp, 3.90953755e-06_wp, 3.8129361_wp, 1.76267436e-04_wp, -1.05819658e-04_wp, -7.21658762e-06_wp, 1.19286822e-02_wp, -1.77369907e-03_wp, 0.13387112_wp, 39.486862_wp, 0.24885239_wp, 0.29916763_wp, 4.1707320_wp, 3.9112310_wp, 1.9251275_wp, 4.49750992e-03_wp, 6.01600004e-05_wp, 8.74410020e-08_wp, 2.5338767_wp, -1.69092222e-04_wp, -1.41368364e-04_wp, -2.20386923e-04_wp, 0.0_wp, 0.0_wp, 0.0_wp], [16, 9]))	The first ephemeris data table (this is standish's table 2). keplerian elements valid 3000 bc - 3000 ad
type(ephem),	private,	parameter, dimension(2)	::	eph_set	=	[eph1, eph2]	the set of ephem options that are available.

Derived Types

type, public, extends(ephemeris_class) :: standish_ephemeris

Standish ephemeris class for computing the approximate positions of the major planets.

Type-Bound Procedures

procedure, public :: get_rv => standish_rv_func

type, private, extends(base_class) :: ephem

an ephemeris defined using a date range and a set of elements from the reference.

Components

Type	Visibility	Attributes		Name		Initial
integer,	public		::	id	=	0	a unique ID code that distinguishes a variable from other variables of the same type.
character(len=name_len),	public		::	name	=	''	the variable name
real(kind=wp),	public,	dimension(2)	::	jd_range	=	zero	valid julian date range
real(kind=wp),	public,	dimension (16, 9)	::	o	=	zero	keplerian elements terms

Type-Bound Procedures

generic, public :: operator(==) => base_class_equal
generic, public :: operator(/=) => base_class_not_equal

Functions

private pure function spice_id_to_standish_id(spice_id) result(old_id)

Author: Jacob Williams
Date: 3/4/2018

Convert the NAIF SPICE ID code to the one used by the standish ephemeris. Returns 0 if the body was not found.

Arguments

Type	Intent	Optional	Attributes		Name
integer,	intent(in)			::	spice_id	the ID code used by SPICE

Return Value integer

the ID code used by this module

private pure function kepler(ma, ec)

solve kepler's equation ma = ea + ec*sin(ea)

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in)			::	ma	eccentricity
real(kind=wp),	intent(in)			::	ec	mean anomaly in rad

Return Value real(kind=wp)

private pure function tbl(jd) result(itbl)

Determine which data set to use for highest accuracy for the given julian date.

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in)			::	jd	julian date (eg 2451545.0)

Return Value integer

Which data set to use:

itbl=1 jd in range of table 1 (1800ad-2050ad) - highest accuracy
itbl=2 jd outside range of table 1 but in range of table 2 (3000bc-3000ad)
itbl=0 3000bc3000ad julian date out of range for ephemeris.

Subroutines

private subroutine standish_rv_func(me, et, targ, obs, rv, status_ok)

Author: Jacob Williams
Date: 3/4/2018

Return the state of the targ body relative to the obs body, in the inertial frame [ICRF].

Arguments

Type	Intent	Attributes		Name
class(standish_ephemeris),	intent(inout)		::	me
real(kind=wp),	intent(in)		::	et	ephemeris time [sec]
type(celestial_body),	intent(in)		::	targ	target body
type(celestial_body),	intent(in)		::	obs	observer body
real(kind=wp),	intent(out),	dimension(6)	::	rv	state of targ w.r.t. obs [km, km/s]
logical,	intent(out)		::	status_ok	true if there were no problems

private pure subroutine helio(np, jd, p, itbl)

For planet np and julian date jd and using using table itbl, return j2000 ecliptic position (au) and velocity (au/yr). in cartesian coordinates (p(1)-p(6)).

Arguments

Type	Intent	Attributes		Name
integer,	intent(in)		::	np	planet 1-9
real(kind=wp),	intent(in)		::	jd	julian date
real(kind=wp),	intent(out),	dimension(6)	::	p	position (au)/velocity (au/yr)
integer,	intent(out)		::	itbl	table used or error if zero

private pure subroutine calcelements(np, jd, itbl, z)

Calculate current elements z(jd) for planet j from jpl data

Arguments

Type	Intent	Attributes		Name
integer,	intent(in)		::	np	planet number (1-9)
real(kind=wp),	intent(in)		::	jd	julian date
integer,	intent(in)		::	itbl	which table to use (1-2)
real(kind=wp),	intent(out),	dimension(8)	::	z	elements for `jd` Read more…

private pure subroutine el2op(z, po)

heliocentric coordinates for orbital plane from elements

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in),		dimension(8)	::	z	elements [a,e,i,l,w,o,ma,ea]
real(kind=wp),	intent(out),		dimension(6)	::	po	coordinates and velocities

private pure subroutine op2ec(z, po, pe)

heliocentric coordinates j2000 ecliptic plane from orbital plane

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(8)	::	z	elements a,e,i,l,w,o,ma,ea
real(kind=wp),	intent(in),	dimension(6)	::	po	orbital plane coordinates
real(kind=wp),	intent(out),	dimension(6)	::	pe	j2000 ecliptic plane coordinates

private pure subroutine ec2eq(pe, pq)

converts cartesian heliocentric j2000 ecliptic to equatorial

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in),		dimension(6)	::	pe	ecliptic
real(kind=wp),	intent(out),		dimension(6)	::	pq	equatorial

private pure subroutine eq2ec(pq, pe)

converts cartesian heliocentric equatorial to ecliptic

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in),		dimension(6)	::	pq	equatorial
real(kind=wp),	intent(out),		dimension(6)	::	pe	ecliptic

private pure subroutine sphere(x, y, z, rho, theta, phi)

cartesian to spherical coordinates (angles in radians)

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	x	x = r cos(phi) cos (theta)
real(kind=wp),	intent(in)	::	y	y = r cos(phi) sin(theta)
real(kind=wp),	intent(in)	::	z	z = r sin(phi)
real(kind=wp),	intent(out)	::	rho	distance
real(kind=wp),	intent(out)	::	theta	longitude
real(kind=wp),	intent(out)	::	phi	latitude

public subroutine standish_module_test()

Author: Jacob Williams
Date: 3/4/2018

Test routine for the standish_module routines.

Arguments

None

standish_module Module

Contents

Variables

Derived Types

Functions

Subroutines

See also

Reference

History

Original license

Uses

Used by

Variables

Derived Types

type, public, extends(ephemeris_class) :: standish_ephemeris

Type-Bound Procedures

type, private, extends(base_class) :: ephem

Components

Type-Bound Procedures

Functions

private pure function spice_id_to_standish_id(spice_id) result(old_id)

Arguments

Return Value integer

private pure function kepler(ma, ec)

Arguments

Return Value real(kind=wp)

private pure function tbl(jd) result(itbl)

Arguments

Return Value integer

Subroutines

private subroutine standish_rv_func(me, et, targ, obs, rv, status_ok)

Arguments

private pure subroutine helio(np, jd, p, itbl)

Arguments

private pure subroutine calcelements(np, jd, itbl, z)

Arguments

private pure subroutine el2op(z, po)

Arguments

private pure subroutine op2ec(z, po, pe)

Arguments

private pure subroutine ec2eq(pe, pq)

Arguments

private pure subroutine eq2ec(pq, pe)

Arguments

private pure subroutine sphere(x, y, z, rho, theta, phi)

Arguments

public subroutine standish_module_test()

Arguments