How does one use SPICE to calculate sunrise and sunset?

I've seen it mentioned here that SPICE can be used to calculate sunrise and sunset. How?

According to the official documentation, SPICE is used to calculate:

• Time system conversions
• Positions of spacecraft and natural bodies
• Orientation of spacecraft and instruments
• Reference frame transformations
• Illumination angles

I know sunrise and sunset can be calculated from the Sun's coordinates and the sunrise equation.

I was wondering if SPICE had a specific algorithm for this purpose?

• This is an interesting and legitimate question; using SPICE you can specify a location on many solar system bodies using a latitude and longitude and calculate an altitude and azimuth of another body, and at least for Earth you can correct for atmospheric refraction, and I believe SPICE will help you to search for a time when altitude is zero (sunrise/sunset). Here's an example of doing it with JPL's Horizons website which I think itself uses SPICE. But I wonder if in addition to your question as asked, you would be interested in also trying
– uhoh
May 15 at 3:14
• …something that mighrt be a lot easier? If you would like to dabble in Python, then check rhodesmill.org/skyfield/almanac.html Skyfield uses the same JPL Development Ephemerides that Horizons and SPICE do, but is a heck of a lot easier and more intuitive to use because it's a modern Python implementation rather than something originally designed for career JPL mission planners using mainframes in the 1980's. What a difference 40 years makes!
– uhoh
May 15 at 3:15

I'm no SPICE expert but here are some potential solutions (unless of course you want to try Skyfield's Almanac methods!)

Firsrt possibility, but all these may not work

this answer links to this comment links to SPICE Tutorials, Updated December 11, 2019 links to the 69 slide presentation SPICE Geometry Finder (GF) Subsystem; Searching for times when specified geometric conditions occur dated January 2020.

• Much SPICE software computes a geometry parameter at a given time, t, i.e. x = f(t). Example: on 2011 MAR 30 14:57:08, what is the spacecraft’s altitude above Mars?
• The Geometry Finder subsystem does the inverse: it finds times when specified geometric events occur

GF provides two primary types of event-finding APIs

• Boolean: a geometric condition (an event) is true or false – Numeric: a geometric quantity has a given value, is within a given range or has achieved a local or global maximum or minimum

In your case you would specify a position on a body, then look for times when elevation was 0 or 0+/- small.

Another possibility, this seems to address exactly what you are looking for

• https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/spkcpo_c.html

Examples

The numerical results shown for these examples may differ across platforms. The results depend on the SPICE kernels used as input, the compiler and supporting libraries, and the machine
specific arithmetic implementation.

1. Compute apparent solar azimuth and elevation as seen from a specified surface point on the earth.

In this example we'll use the location of the DSN station
DSS-14 as our surface point.

We'll perform the solar azimuth and elevation computation two ways:

• Using a station frame kernel to provide the
specification of a topocentric reference frame centered at DSS-14.

• Computing inline the transformation from the earth-fixed, earth-centered frame ITRF93 to a topocentric frame centered at DSS-14.

Note that results of the two computations will differ slightly. There are three sources of the differences:

1. The station position is time-dependent due to tectonic plate motion, and epochs of the station positions used to specify the axes of the topocentric frame are different in the two cases. This gives rise to different orientations of the frame's axes relative to the frame ITRF93.

2. The two computations use different earth radii; this results in computation of different geodetic latitudes of the station. This difference also affects the topocentric frame orientation relative to ITRF93.

3. The station movement between ET and the epoch at which the DSS-14_TOPO frame is specified contributes a very small offset---on the order of 10 cm---to the station-sun position vector, expressed in the ITRF93 frame.

and there's a lot more there.

It seems you can do the same thing in Python using SpicyPy

Find the azimuth and elevation of the apparent position of the Moon as seen from the DSN station DSS-13 by the following steps:

1. Find the apparent position vector of the Moon relative to the DSN station DSS-13 in the topocentric reference frame DSS-13_TOPO at epoch ET. Use light time and stellar aberration corrections.

For this step, you'll need to have loaded a station SPK file providing geocentric station position vectors, as well as a frame kernel specifying topocentric reference frames centered at the respective DSN stations. (Other kernels will be needed as well; you must choose these.)

1. Convert the position vector to latitudinal coordinates. Use the routine spiceypy.reclat for this computation.

2. Compute the Moon's azimuth and elevation as follows: azimuth is the negative of topocentric longitude and lies within the range 0-360 degrees; elevation is equal to the topocentric latitude. Display the results in degrees.

• Thanks for the detailed answer May 18 at 2:49