cosinekitty/astronomy

Astronomy Engine: multi-language calculation of Sun, Moon, and planet positions. Predicts lunar phases, eclipses, transits, oppositions, conjunctions, equinoxes, solstices, rise/set times, and other events. Provides vector and angular coordinate transforms among equatorial, ecliptic, horizontal, and galactic orientations.

Astronomy Engine

Supported Programming Languages

C	Examples	Code & Docs
C#	Examples	Code & Docs
JavaScript	Examples Browser Node.js	Code & Docs
Python	Examples	Code & Docs
Kotlin / JVM	Examples Kotlin Java	Code & Docs

Overview

Astronomy Engine is a suite of open source libraries for calculating positions of the Sun, Moon, and planets, and for predicting interesting events like oppositions, conjunctions, rise and set times, lunar phases, eclipses, transits, and more.

It supports several popular programming langauges with a consistent API. Function and type names are mostly consistent across all the supported languages.

Astronomy Engine is designed to be small, fast, and accurate to within ±1 arcminute. It is based on the authoritative and well-tested models VSOP87 and NOVAS C 3.1.

These libraries are rigorously unit-tested against NOVAS, JPL Horizons, and other reliable sources of ephemeris data. Calculations are also verified to be identical among all the supported programming languages.

Features

• Provides calculations for the Sun, Moon, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, and Pluto.

• Calculates all supported objects for any calendar date and time for millennia before or after the present.

• Provides heliocentric and geocentric Cartesian vectors of all the above bodies.

• Determines apparent horizon-based positions for an observer anywhere on the Earth, given that observer's latitude, longitude, and elevation in meters. Optionally corrects for atmospheric refraction.

• Calculates rise, set, and culmination times of Sun, Moon, and planets.

• Finds civil, nautical, and astronomical twilight times (dusk and dawn).

• Finds date and time of Moon phases: new, first quarter, full, third quarter (or anywhere in between as expressed in degrees of ecliptic longitude).

• Predicts lunar and solar eclipses.

• Predicts transits of Mercury and Venus.

• Predicts lunar apogee and perigee dates, times, and distances.

• Predicts date and time of equinoxes and solstices for a given calendar year.

• Determines apparent visual magnitudes of all the supported celestial bodies.

• Predicts dates of planetary conjunctions, oppositions, and apsides.

• Predicts dates of Venus' peak visual magnitude.

• Predicts dates of maximum elongation for Mercury and Venus.

• Calculates the positions of Jupiter's four largest moons: Io, Europa, Ganymede, and Callisto.

• Allows custom simulation of the movements of user-defined small bodies, such as asteroids and comets, through the Solar System.

• Converts angular and vector coordinates among the following orientations:

  • Equatorial J2000
  • Equatorial equator-of-date
  • Ecliptic J2000
  • Topocentric Horizontal
  • Galactic (IAU 1958)

• Determines which constellation contains a given point in the sky.

• Calculates libration of the Moon.

• Calculates axis orientation and rotation angles for the Sun, Moon, and planets.

Why I Created This Thing

I have been an amateur astronomer since childhood. I still remember the amazement I felt when I saw Saturn through a backyard telescope for the first time. As a software developer, I naturally became interested in combining my love of astronomy with my computer programming skills.

In 2008, I started to learn about formulas for calculating positions of the Moon and planets. I discovered many wonderful resources, including

• Paul Schlyter's lucid and educational page How to compute planetary positions.
• Practical Astronomy with your Calculator, third edition, by Peter Duffett-Smith, Cambridge University Press. ISBN 0 521 35629 6.
• Astronomy on the Personal Computer by Oliver Montenbruck and Thomas Pfleger. ISBN-13: 978-3540672210.

I implemented algorithms based on these resources. Over time, however, I noticed that they were not quite as accurate as I would like. Their calculated positions differed from those reported by online tools like JPL Horizons and Heavens Above by large fractions of a degree in many cases.

In 2019 I renewed my interest in astronomy calculations, with the goal of creating something more accurate that could be written in JavaScript to run inside a browser. I studied how professional astronomers and space agencies did their calculations. First I looked the United States Naval Observatory's NOVAS C 3.1 library. I quickly realized it could not be ported to the browser environment, because it required very large (hundreds of megabytes) precomputed ephemeris files.

This led in turn to studying the French Bureau des Longitudes model known as VSOP87. It requires more computation but the data is much smaller, consisting of trigonometric power series coefficients. However, it was still too large to fit in a practical web page.

Furthermore, these models were extremely complicated, and far more accurate than what I needed. NOVAS, for example, performs relativistic calculations to correct for the bending of light through the gravitational fields of planets, and time dilation due to different non-intertial frames of reference! My humble needs did not require this herculean level of complexity. So I decided to create Astronomy Engine with the following engineering goals:

• Support JavaScript, C, C#, and Python with the same algorithms, and verify them to produce identical results. (Kotlin support was added in 2022.)
• No external dependencies! The code must not require anything outside the standard library for each language.
• Minified JavaScript code less than 120K. (The current size is 116485 bytes.)
• Accuracy always within 1 arcminute of results from NOVAS.
• It would be well documented, relatively easy to use, and support a wide variety of common use cases.

The solution I settled on was to truncate the VSOP87 series to make it as small as possible without exceeding the 1 arcminute error threshold. I created a code generator that converts the truncated tables into C, C#, JavaScript, and Python source code. Then I built unit tests that compare the calculations against the NOVAS C 3.1 code operating on the DE405 ephemeris and other authoritative sources, including the JPL Horizons tool. Basing on VSOP87 and verifying against independent trusted sources provides extra confidence that everything is correct.

Pluto was a special case, because VSOP87 does not include a model for it. I ended up writing a custom gravitation simulator for the major outer planets to model Pluto's orbit. The results are verified against NOVAS and the model TOP2013.

As far as I know, Astronomy Engine is the only open source solution in existence that combines very compact code for four major programming languages with such rigorous validation and testing at a reasonable accuracy threshold. The 1-arcminute accuracy is not good enough for spacecraft navigation, but it is good enough for most amateur uses, and allows the code to be much simpler, faster, and smaller.

I am committed to maintaining this project for the long term, and I am happy to answer questions about how to solve various astronomy calculation problems using Astronomy Engine. Feel free to reach out on the discussions page or submit a new issue.

Acknowledgements

astronomy's Issues

Language support: Go

Add support for the Go programming language.

Language support: C# / .NET

Especially focus on .NET Core.

Automatically validate pull requests and commits

Use a tool like Travis CI to run unit tests when commits and/or pull requests are submitted. This will help flag problems quickly so they can't slip through the cracks.

Add minified JavaScript

Provided a minified version of the generated JavaScript library. Use Google Closure compiler?

Support 1700-01-01 <= date < 2200-01-01

Mainly this will affect the Chebyshev model for Pluto. All the VSOP stuff should work fine as-is.

Host live demos of JavaScript browser examples

I was trying to use GitHub Pages just so I could host live demos of the JavaScript browser examples. It worked, but at the cost of having two different versions of this project with different Markdown renderers. I have turned off GitHub Pages. Now I need to create my own static website just for hosting the demos. Also need to automate the process of publishing there every time I hit a stable release.

Improve SearchRelativeLongitude convergence for Mercury and Mars

Mercury and Mars have high orbital eccentricities. This causes SearchRelativeLongitude to take a long time to converge, especially for Mercury. Improve the algorithm so that it converges faster. The other planets converge within 4 or 5 iterations, but Mercury takes somewhere between 20 and 100.

Determine what constellation a given (RA,DEC) is within

Given (RA,DEC) in J2000 coordinates, calculate 1875 coordinates and look up the correct constellation name. Because many users won't care about constellations, this will be an optional add-on module in the JavaScript.

Perihelion, aphelion for all planets

Calculate when any of the planets is closest to the Sun or farthest from the Sun, along with distance in AU. Will need to obtain authoritative test data: date & time of a series of apsides for every planet.

I can use the following to test apsides of Earth. Not sure how to do this for other planets; perhaps could use JPL Horizons to plot heliocentric distances of each planet with respect to time.

http://astropixels.com/ephemeris/perap2001.html

Tables not rendering correctly on github.com

Tables in documentation looks correct at:
https://cosinekitty.github.io/astronomy/source/c/

But is all messed up in:
https://github.com/cosinekitty/astronomy/tree/master/source/c/

Why don't the tables work there?

Find maximum elongation of Mercury, Venus

Calculate when the inferior planets have maximum angular separation from the Sun as seen from the center of the Earth. Report the angular separation in degrees along with the date/time of the event. Distinguish between morning sky visibility and evening sky visibility.

Add a mission statement to README.

Document why this project was created. Tell the story of how I started out trying to translate NOVAS directly to JavaScript and discovered it requires large ephemeris files.

Document project goals:

• Reliable and well-tested.
• Easier for programmers who are not professional astronomers to understand and use.
• Small and efficient code within 1 arcminute precision.
• Focused on use in amateur astronomy, calendars, etc.

Terminology: relative longitude means different things

Astronomy_SearchRelativeLongitude finds relative longitude as measured from the Sun's point of view. But astro_elongation_t::relative_longitude indicates relative longitude from the Earth's point of view. This is likely to confuse people. Consider renaming one or both of these.

Visual magnitude of Sun, Moon, and all planets. Phase of Moon and planets.

Calculate the visual magnitude of the Sun, Moon, and all planets as seen from the Earth. Geocentric calculation is fine; does not need to be micro-adjusted for a topocentric observer. Because this entails calculating distance and phase angle, include those data in the calculated results. (No phase angle for the Sun because it emits light instead of reflecting it!) Also include ring angle for Saturn in returned results.

More efficient zero-crossing search algorithm

I don't like how the Search() function does a dumb binary search for the zero crossing. Something like Newton's method or quadratic interpolation could be a lot more efficient. I don't want to sacrifice accuracy though. And there should be some empirical evidence that things are really faster with such a change: like actually timing unit tests that exercise Search.

Language support: Python 3

Port all the target code to Python 3.

Astronomy glossary

Create a document that explains astronomy terminology. Alphabetically sorted glossary. Allow deep linking from other documents.

• Observers
  • Heliocentric
  • Geocentric
  • Topocentric
• Coordinate systems
  • Equatorial
  • Ecliptic
  • Horizontal
• Different time scales (UT, TT).
  • UTC/UT1
  • Terrestrial Time
• Barycenters: SSB, EMB.
• Apsides: pericenter and apocenter.
• Equinox has two meanings:
  • A location on the Earth's orbit.
  • A time at which the Earth reaches the equinox point.
• (Search through documentation for all astronomy jargon.)

Investigate JPL Horizons calculation of refraction for negative altitude angles

JPL Horizons does not seem to calculate refraction for negative altitude angles the same way NOVAS C 3.1 does. This forced me to skip horizontal coordinate checking for tests cases of objects below the horizon. I would like to reverse engineer the JPL Horizons refraction angle formula by comparing data sets with refraction disabled (airless) and enabled. The NOVAS formula seems unpleasant to work with because it looks like there is a step discontinuity once the zenith angle goes above 91 degrees. It would be nice to see if Horizons has smooth behavior, and if so, I would like to emulate it, both for more complete unit testing and for being able to create Chebyshev approximations of objects like the Moon to have efficient rise/set/culm calculations.

Language support: C, C++

Port all the code to C with a C++ friendly wrapper.

Maximum visual magnitude of Venus

Find date/time of when Venus appears brightest in the sky. Also report the maximum magnitude value at this time. Depends on #18.

Incorrect Sun altitude

I found a case where there is a sudden jump from a negative Sun altitude to a positive, then back to negative, all over a span of 3 seconds:

1750-11-21 00:46:08.384, peak_altitude = -4.178
1750-11-21 00:46:09.384, peak_altitude = +6.752
1750-11-21 00:46:10.384, peak_altitude = -4.86

This is for the Sun with an observer at latitude +15, longitude +75, height = 0 meters.

Python Test_RiseSet seems really slow

Do profiling to figure out what's going on here. Seems much slower than it should be, compared to JavaScript and C versions.

Chebyshev coefficients slightly different between Windows and Linux

The Chebyshev coefficients generated in the Javascript code are slightly different for Pluto, when generated by Windows versus Linux. I suspect the real issue is that I'm running different versions of node:
Windows: v10.15.3
Linux: v8.15.1

Opposition, conjunction of planets

Use relative longitude between Sun and other objects to determine the times of planet conjunctions and oppositions. For Mercury and Venus, distinguish between inferior and superior conjunctions.

Eliminate use of Julian Dates

Julian Dates reduce precision by consuming so many mantissa bits. Rework to use days since 2000 (or 1900, or MJD=1857-11-17?), which is closer to the center of my eventual range of supported years (1600..2200).

SearchMaxElongation function name inconsistent with Elongation function name

The "elongation" in the name SearchMaxElongation is inconsistent with the function Elongation, which calculates an ecliptic longitude difference. It is really related to AngleFromSun.

Likewise, the type ElongationEvent is confusing. Rename whatever is needed to get consistent terminology throughout the library.

Lunar apogee, perigee

Find when moon is closest, furthest from the Earth's center. Report distance (center to center?) in kilometers.

Travis CI occasional build error trying to download ephemeris file

This is happening every now and then in the Travis CI build. Need to set a longer timeout and/or retry after errors.

$ cd generate && ./run && ./verify_clean
Ephemeris file not found: lnxp1600p2200.405
Trying to download for you from:
https://github.com/cosinekitty/ephemeris/raw/master/lnxp1600p2200.405
--2019-05-10 00:05:39--  https://github.com/cosinekitty/ephemeris/raw/master/lnxp1600p2200.405
Resolving github.com (github.com)... 192.30.253.113
Connecting to github.com (github.com)|192.30.253.113|:443... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://raw.githubusercontent.com/cosinekitty/ephemeris/master/lnxp1600p2200.405 [following]
--2019-05-10 00:05:39--  https://raw.githubusercontent.com/cosinekitty/ephemeris/master/lnxp1600p2200.405
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 151.101.0.133, 151.101.64.133, 151.101.128.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|151.101.0.133|:443... connected.
HTTP request sent, awaiting response... 503 first byte timeout
2019-05-10 00:06:40 ERROR 503: first byte timeout.
ERROR(./run): Could not download using wget: https://github.com/cosinekitty/ephemeris/raw/master/lnxp1600p2200.405
The command "cd generate && ./run && ./verify_clean" exited with 1.

Discrepancy in JavaScript and C output

Line 145936 in c_check.txt and js_check.txt have a small discrepancy:

ctest(Diff): Maximum numeric difference = 1.81899e-12, worst line number = 145936

js_check.txt

s Venus -39864.1907954055495793 -39864.1907264927140204 17.3098151195010921 -28.0371006234166060 0.3456613889446833 34.5081734386010766 -88.6913704437932324

c_check.txt

s Venus -39864.1907954055495793 -39864.1907264927140204 17.3098151195010921 -28.0371006234166096 0.3456613889446833 34.5081734385992576 -88.6913704437932608

See if this can be improved or eliminated.

rise, set, culm

Provide functions in JavaScript code to find the next rise, set, culm times for any object between a given pair of times. There may not be any (e.g. Moon never rises on some days, and some objects are circumpolar for a while at latitudes far from the equator).

Consider creating generic search functions that can find other events such as full moon, conjunctions, etc.

Travis CI: verify that generated code matches what is checked in

On Travis CI only, confirm that generated code exactly matches what was checked in by the committer. Otherwise we are not publishing the same code that the generator creates!

moon phases

Provide ability to find the next major moon phase, given a start search date: new, first quarter, full moon, third quarter.

calculate aberration and light travel correction

Include aberration in the calculation of (RA,DEC). Lack of this calculation contributes up to 1/3 of an arcminute of error.

Excessive angular error

Max angular error is reaching 1.17 arcminutes. Get this down below 1.0 arcminutes. (Start by fixing light travel time and aberration.)

horizontal coordinates

Given a location (longitude, latitude, elevation), a date and time, and equatorial coordinates (RA, DEC), calculate horizontal coordinates (azimuth, altitude) for an object. Must include refraction calculations so that rise/set times (to be implemented later) are accurate.

Investigate using finer-grained delta-t data

It may be better to use direct empirical data for delta-t values instead of calculating leap seconds.

See: https://www.usno.navy.mil/USNO/earth-orientation/eo-products/long-term

Data files:
http://maia.usno.navy.mil/ser7/historic_deltat.data
http://maia.usno.navy.mil/ser7/deltat.data
http://maia.usno.navy.mil/ser7/deltat.input
http://maia.usno.navy.mil/ser7/deltat.preds

Bright star database

Offer a database of the brighter stars (mag <= 4.0)? This will be a separate add-on module in the JavaScript version.

No need to correct Sun position for light travel time

When correcting the position of an object for light travel time, because the position is expressed in heliocentric coordinates, there is no need to correct the Sun's position because it is always at (0, 0, 0).

Create live demo right on the GitHub project page

Embed JavaScript code to demo the Astronomy code right on the GitHub project page. I looked online and apparently Markdown files are a superset of HTML, thus allowing script tags, etc. Could show the location of all the planets, the Moon's phase, stuff like that.

npm package

Create an npm package for the JavaScript version of the Astronomy library.

Improve Saturn visual magnitude calculation

It would be nice to find authoritative values for Saturn's visual magnitude to use as test data. The following paper has formulas for Saturn (and the other planets), but no test data.

https://arxiv.org/ftp/arxiv/papers/1808/1808.01973.pdf

The source code here goes with the paper, and it might be possible to adapt it to create test data:

https://sourceforge.net/projects/planetary-magnitudes/files/

add unit test based on JPL Horizons data

Create hand-crafted data file with test cases from JPL Horizons: (RA,DEC), (az,alt). Test every body from a variety of times and geographic locations. This is mainly as a sanity check that I'm not calling NOVAS functions wrong.

solstice, equinox

Find next solstice or equinox.

Enable translation of doxygen to markdown from Windows.

Currently the makedoc Linux bash script knows how to translate doxygen xml files to markdown format, but it uses xsltproc to do so. Need to add support for generating the same documentation in Windows: running doxygen and translating the resulting xml.

Run C unit tests on Windows

Add support to build C unit tests on Windows, and run them. Can create a Visual Studio solution to build the code written by generate.c. This way developers can work on Windows just as well as Linux.

Automatically trim trailing whitespace on all source files.

I encountered a problem where doxygen running on my Windows laptop causes the generated C documentation README.md to be messed up. The cause turned out to be trailing whitespace in astronomy.c. Once I cleaned that up, the markdown generated by Windows is the same as on Linux.

Plus, trailing whitespace is just annoying. I want to automatically strip trailing whitespace from all source code before generating documentation.

Figure out why c_check has lower error than js_check.

This was a real surprise:

CheckTestOutput: Verified 382306 lines of file temp/js_check.txt : max error = 0.727661 arcmin
CheckTestOutput: Verified 382306 lines of file temp/c_check.txt  : max error = 0.445261 arcmin

This makes me think there is some unnecessary error in the JS output (maybe just need more decimal places in the print statements?)

After fixing this, tighten up the error tolerances.

Tiny numeric discrepancy in Python output

I'm not sure this is just a quirk of Python floating point, or if there is a subtle bug somewhere. But there is a discrepancy between the numbers calculated by Python and those calculated by C and JavaScript:

don@spearmint:~/github/astronomy/generate $ ./ctest diff temp/{py,js}_check.txt
ctest(Diff): Maximum numeric difference = 8.90736e-11, worst line number = 33246
ERROR: Excessive error comparing files temp/py_check.txt and temp/js_check.txt
ctest exiting with 1
don@spearmint:~/github/astronomy/generate $ ./ctest diff temp/{c,js}_check.txt
ctest(Diff): Maximum numeric difference = 1.81899e-12, worst line number = 145936
ctest exiting with 0

Windows build error in novas.c

This just popped up for the first time.

novas.c(2306): warning C4701: potentially uninitialized local variable 'obl' used

Use consistent "Astronomy Engine" throughout code and documentation

Change the lingering vague names of the project to "Astronomy Engine" throughout all code and documentation.