Sunrise / sunset

Orekit usage
1

Hello,
I am trying to calculate the visibility of the sun from a point on Earth as follows (using the python interface v 12.1.2 and the data downloaded from the orekit repo):

point = GeodeticPoint(radians(48.843391), radians(2.390224), 0.0)
sun = CelestialBodyFactory.getSun()
earth = CelestialBodyFactory.getEarth()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth_ellipsoid = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                                   Constants.WGS84_EARTH_FLATTENING,
                                   itrf)
topo_rf = TopocentricFrame(earth_ellipsoid, point, "Topocentric")
now = AbsoluteDate(2024, 10, 18, 16, 49, 0.0, pop.utc)
sun_el = topo_rf.getElevation(
    ExtendedPVCoordinatesProvider.cast_(sun).getPosition(now, earth.bodyOrientedFrame),
    pop.earth.bodyOrientedFrame,
    now
)
print(now)
np.rad2deg(sun_el)

With the numbers defined above I get that the elevation is 0 at 18:49 GMT+2 (=16:49 UTC) whereas any internet calculator at the same coordinates gives a sunset at 18:56 GMT+2, i.e. about 6 minutes after.

Where is the difference in calculation in your opinion? Do I need a better model of the Earth than the WGS84? Is the data I am using imprecise? Am I doing anything wrong in the computation?

Thanks!

2

There are three main things I can think about :

  • your computation and internet calculators may use different Earth Orientation Parameters (i.e. the content of orekit-data folder for example)
  • here you compute direction of the center of Sun and ignore completely refraction, which is quite big at the horizon
  • Sun setting is not at 0°, the civil dawn/dusk is 6° below horizon, the nautical dawn/dusk is 12° below horizon and the astronomical dawn/dusk is 12° below horizon

I would say the first reason is probably not the culprit, it would most probably give a difference less than 0.9 second (because UT1-UTC is guaranteed to be in a ±0.9s range) ; the second reason, refraction, may be more important and the third reason is the biggest difference : it takes a few minutes to go from 0° to -6°/-12°/-18°.

3

Thanks luc.
I honestly doubt the internet calculators consider refraction and calculating dusk at -6/-12/-18 degrees below the horizon adds time in the order of hours rather than minutes.

Your suggestion about the centre of the Sun looks the best. I just tried considering the apparent radius of the Sun to calculate dusk and this clearly improves the error. To be precise I updated the calculation as

sun_distance = ExtendedPVCoordinatesProvider.cast_(sun).getPosition(date, earth.bodyOrientedFrame).norm
sun_angular_radius = 2*np.rad2deg(np.arctan(Constants.SUN_RADIUS/sun_distance))
new_sun_el = sun_el+sun_angular_radius

I see new_sun_el crosses the 0 mark about 5 minutes later (with test case) this way, reducing the difference wrt internet tools to about 1 minute.
I am guessing the remaining error is about the altitude which I cannot set nor see on online tools and which I am ignoring.

1 Like