2
$\begingroup$

I'm needing to accurately find the sun and moon position in ECEF and I'm wanting to make use of JPL's Horisons VECTORS to do this.

What data should I export from the web interface and how should I convert it?


Background

I need the accurate positions to find 3d angles of incoming light onto slopped surfaces and also if the view of a point on that surface of the sun or moon is obstructed by another location on that surface that i have mapped in ECEF.

There are around 500 surfaces about 20sqkm each derived from LiDAR sources located at various parts of the earth.


What worked for me

I ended up following @Alfonso Gonzalez's advice and getting the position directly from SPICE, and I'm able to get the same results as the Horizons OBSERVER data but doing the math in the ITRF93/ECEF reference frame.

This is the python code I'm using:

import spiceypy as spice

spice.furnsh('earth_000101_210617_210326.bpc')
spice.furnsh('naif0012.tls')
spice.furnsh('de432s.bsp')

et = spice.str2et('2021-03-23T00:00:00')

sun_naif_id = spice.bods2c('SUN')
earth_naif_id = spice.bods2c('EARTH')

sun, lt = spice.spkezp(sun_naif_id, et, 'ITRF93', 'LT+S', earth_naif_id)
sun *= 1000

print("Sun Position at 2021-03-23 00:00 UTC in ECEF\n  ", sun)

Outputs:

Sun Position at 2021-03-23 00:00 UTC in ECEF
   [-1.49005330e+11 -4.33523011e+09  2.67037253e+09]

As test I tried calculating the sun elevation angle for Adelaide, Australia at 2020-03-23 00:00:00 UTC, and got 36.060167 and Horizons OBSERVER gives 36.0602 for the same time and location.

import numpy as np

adelaide = np.array([-3927127.15169169, 3462288.85408112, -3630705.11577511])
adelaide_z = np.array([-4542168.38388766, 4004530.12124128, -4203154.84360054])

ba = adelaide_z-adelaide
bc = sun-adelaide
cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
z_angle = np.arccos(cosine_angle)

print("Sun Elevation for Adelaide 2021-03-23 00:00 UTC\n  ", 90 - np.degrees(z_angle))

Outputs:

Sun Elevation for Adelaide 2021-03-23 00:00 UTC
   36.06016701724842
$\endgroup$
6
  • $\begingroup$ My guess is it's a problem with the reference frame, if I was to create rotation matrix by comparing locations of the same point on the earth then applying that to the sun/moon vectors, would that work? $\endgroup$ – complistic Mar 26 at 20:57
  • 1
    $\begingroup$ Your conversion to ECEF is flat-out wrong. $\endgroup$ – David Hammen Mar 27 at 0:32
  • 1
    $\begingroup$ Do you really need the 3D position X, Y, Z in space in ECEF, or just the apparent position on the celestial sphere RA/Dec or locally in Alt/Az given a geographic lat/lon on the body in question? I ask because both of those are readily available in Horizons. See for example what I did in this answer. I ask because this might be an instance of an X, Y problem rather than an X, Y, Z problem (pun intended) $\endgroup$ – uhoh Mar 27 at 0:43
  • 1
    $\begingroup$ @uhoh, I need to find the angles from the normals of a 3d surface for both sun and moon. And also if a point on that surface view of the sun or moon is obstructed by another point on that surface. I have ~500 20sqkm surfaces (8 billion faces) located on various parts on the earth. So while I could get the RA/Dec Alt/Az from the OBSERVER data and convert that to an angle in ECEF, I would need to do that for ~200000 locations (every 1km) and interpolate that to get half-decent results. $\endgroup$ – complistic Mar 27 at 12:19
  • 1
    $\begingroup$ Okay I got ya. Yes at only 380,000 km the direction of the Moon can change by up to +/- 1° depending on where you are on Earth, so you do need the distances, not just the normals, unless you don't want to worry about that effect. If you want to assume that the Earth's axis points in a fixed direction then all you need is the rate of Earth's rotation and some simple spherical trigonometry. I know there are standard formulae for these, and some even for the precession of the Earth's axis. I'll look around, including on this site. Thanks for the feedback! $\endgroup$ – uhoh Mar 27 at 12:31
2
$\begingroup$

enter image description here You could do all this analysis using SPICE, including getting ephemeris data from SPICE kernels instead of Horizons. I'm not sure which programming language you're using, but you can get the toolkit directly in C, Fortran, MATLAB, and Java (https://naif.jpl.nasa.gov/naif/toolkit.html) or there's a 3rd party Python wrapper that I use all the time called SpiceyPy (https://spiceypy.readthedocs.io/en/master/documentation.html)

First, if you want to get the ephemeris from SPICE, you'll need an ephemeris kernel, which can be found here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/ de432s.bsp will work for what you're doing. And then you'll need a leapseconds kernel, which can be found here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/ naif0012.tls will work.

In order to convert from inertial to ECEF vectors, there are several options for Earth orientation kernels, but I like to use earth_200101_990628_predict.bpc most of the time, which can be found here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/

Once you have all these, you can use the SPICE toolkit to get the ECEF vectors. Here is how I did it in Python for the groundtrack plot:

enter image description here

If you want to use the data you got from Horizons and convert it to ECEF, you can use SPICE's pxform function to obtain the rotation matrices at each time step. That would look something like this: enter image description here

If you want to do this in python, you can install spiceypy using pip:

$ python3 -m pip install spiceypy

The conversion to latitude / longitude coordinates is out of the scope of the question, but if you're interested, that bf2latlong function is using SPICE's reclat function, which converts from rectangular to latitudinal coordinates. https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/reclat_c.html

Also, SPICE has functions for calculating illumination angles and local solar time, which you could also use for your analysis: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ilumin_c.html

$\endgroup$
5
  • 1
    $\begingroup$ Doing everything in one place using established code is definitely the better way to do it. +1 $\endgroup$ – uhoh Mar 27 at 16:05
  • $\begingroup$ I'll have a look into getting the data this way, I did look into the SPICE and SOFA C libraries but it was not so obvious. If I get a python version to work I'll have a look at the source to see how to implement it in C. $\endgroup$ – complistic Mar 27 at 17:31
  • $\begingroup$ @complistic Spiceypy has an example that goes step by step on how to plot the Cassini trajectory using SPICE kernels, I think you'll find that useful: spiceypy.readthedocs.io/en/master/exampleone.html Also if you go to the "orientation models" section of pck00010.tpc, it explains how they are calculating the rotations to ECEF. You can find that file in the same place the earth orientation kernel is: naif.jpl.nasa.gov/pub/naif/generic_kernels/pck $\endgroup$ – Alfonso Gonzalez Mar 27 at 18:11
  • 1
    $\begingroup$ @AlfonsoGonzalez, This seems to be close, I'm getting 36.0546 for the elevation (Horizons gives 36.0602) but this dependency might be due to how I'm calculating Adelaide's location and zenith, I'll see if using SPICE rather then Proj4 to do the WGS84 to ITRF93 conversion helps. $\endgroup$ – complistic Mar 27 at 19:22
  • 1
    $\begingroup$ I found the correct setting to get the same result as Horizons, I needed to use earth_000101_210617_210326.bpc with the LT+S option to account for stellar aberration. $\endgroup$ – complistic Mar 27 at 20:57
1
$\begingroup$

This is a simple answer that is not super-accurate but it will get you started.

Hopefully someone will chime in with corrections/refinements!


If we work in ecliptic coordinates then for a spherical Earth of radius $r_0$ the position of a location is just

\begin{align} x & = r_0 \cos\phi \cos(\lambda + \theta)\\ y & = r_0 \cos\phi \sin(\lambda + \theta)\\ z & = r_0 \sin\phi \\ \end{align}

where $\phi$ and $\lambda$ are latitude and longitude, $\theta$ is the rotational angle of the Earth.

Presumably $x, y, z$ will be arrays in your program, of the same size and dimensionality as $\phi$ and $\lambda$.

If you want to use an oblate Earth and include height above the WGS84 geoid (i.e. altitude) the see this answer.

One approximation here is that the precession of Earth's axis is not taken into account. Below you'll see that the epoch for the Earth's equator is J2000 and it's moved a tiny bit since then.

From Earth rotation angle:

$$\theta(t_U) = 2 \pi (0.779\,057\,273\,2640+1.002\,737\,811\,911\,354\,48t_{U})$$

where $t_U$ is the Julian UT1 date − 2451545.0.

For this approximation you can use any Julian date converter your programming environment provides. For example the Julian date of 27-Mar-2021 12:00 UTC is about 2459301.00. Note that JD crosses zero integer value at noon, not midnight, so as not to confuse poor astronomers in the middle of the night.

Please be sure to ask Horizons to give you the positions of the Sun, Moon, and Earth in

reference plane: Earth mean equator and equinox of reference epoch
reference system: ICRF/J2000.0

per this answer. That looks like this:

reference plane: Earth mean equator and equinox of reference epoch reference system: ICRF/J2000.0

Then subtract the Sun's position from your $x, y, z$ array of points on the Earth to get vectors at each point that point towards the Sun, and similar for the Moon. You should see as much as a +/- 1° shift in those vectors for the Moon at opposite sides of the Earth because the Moon is so close.

Then once you have the vectors from the faces to the Sun or to the Moon, you can normalize them.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.