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?
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)
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))
Sun Elevation for Adelaide 2021-03-23 00:00 UTC 36.06016701724842