A bug in fraction_illuminated() on Venus?

[sgbmzm]

A bug in fraction_illuminated() on Venus?

It worked well in the skyfield for the moon using fraction_illuminated(). But when I tried to do it on Venus, the results were completely wrong and weird. I could not in any way use fraction_illuminated() on Venus and get a sensible result. (After all, it is known that Venus also has different appearances, and different illuminations, just like the moon has). I use Stellarium software. The illumination percentage of Venus for today September 24, 2022 is 99 percent. Whereas for the date August 11, 2023 the lighting percentage is 0.9. If I use skyfield according to the example given on the skyfield examples page (the example is for the moon, and I changed it for Venus) the result is that today, September 24, 2022, the illumination percentage of Venus is 0.4, while on August 11, 2023, the illumination percentage is 0.5. What is the problem? And how should it be solved? I didn't check for Mercury, but it also has different illuminated, like the moon has. Thank you.

from skyfield.api import load
from skyfield.framelib import ecliptic_frame

ts = load.timescale()
t = ts.utc(2022, 9, 24, 15, 36)

eph = load('de440s.bsp')
sun, venus, earth = eph['sun'], eph['venus'], eph['earth']

e = earth.at(t)
s = e.observe(sun).apparent()
v = e.observe(venus).apparent()

percent = 100.0 * v.fraction_illuminated(sun)

print('Percent illuminated: {0:.1f}%'.format(percent))

Percent illuminated: 0.4%

[Bernmeister]

I put together a test script to compare PyEphem with Skyfield:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from skyfield.api import load, wgs84
import datetime, ephem, skyfield

latitude = -33.8
longitude = 151.2
elevation = 0

observer = ephem.Observer()
observer.lat = str( latitude )
observer.lon = str( longitude )
observer.elev = elevation
observer.date = ephem.Date( datetime.datetime.utcnow() )

def observePyEphem( observer, body ):
    body.compute( observer )
    print( body.name, body.az, body.alt, body.phase )

observePyEphem( observer, ephem.Moon() )
observePyEphem( observer, ephem.Venus() )

def obsverveSkyfield( locationAtNow, body, sun ):
    apparent = locationAtNow.observe( body ).apparent()
    alt, az, distance = apparent.altaz()
    illumination = apparent.fraction_illuminated( sun ) * 100
    print( ' '.join( body.target_name.split()[ 1 : ] ).title(), az, alt, illumination )

ephemeris = load( "de421.bsp" )

timeScale = load.timescale( builtin = True )
now = timeScale.now()

location = wgs84.latlon( latitude, longitude, elevation )
locationAtNow = ( ephemeris[ "earth" ] + location ).at( now )

obsverveSkyfield( locationAtNow, ephemeris[ "moon" ], ephemeris[ "sun" ] )
obsverveSkyfield( locationAtNow, ephemeris[ "venus" ], ephemeris[ "sun" ] )

Got a similar result:

Moon 27:19:25.9 44:24:26.8 1.0348552465438843
Venus 32:36:33.2 47:40:15.9 99.17821502685547
Moon 27deg 19' 24.8" 44deg 23' 28.3" 1.0747593036094294
Venus 32deg 36' 33.3" 47deg 39' 23.7" 0.4221313414979222

Moon looks good, but Venus is off, unless I've blundered. Dunno about other planets/bodies.

[brandon-rhodes]

I've re-thought and rewritten the fraction-illuminated logic and this time have tested it against JPL HORIZONS results. See if installing the development version of Skyfield fixes things!

To install the development version of this project, run:

pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip

[Bernmeister]

Moon 343:42:46.3 47:26:39.2 0.8711615204811096
Venus 344:31:52.5 51:43:30.4 99.18321228027344
Moon 343deg 42' 44.5" 47deg 25' 46.0" 0.9617487871472319
Venus 344deg 31' 51.7" 51deg 42' 45.1" 99.18072805488858

Venus now looks good; Moon is a little off though (unless the issue is within PyEphem).

[Bernmeister]

I modified the script so that both PyEphem and Skyfield use the same date/time:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from skyfield.api import load, wgs84
import datetime, ephem, skyfield

latitude = -33.8
longitude = 151.2
elevation = 0

utcNow = datetime.datetime.utcnow()
print( utcNow )

observer = ephem.Observer()
observer.lat = str( latitude )
observer.lon = str( longitude )
observer.elev = elevation
observer.date = ephem.Date( utcNow )

def observePyEphem( observer, body ):
    body.compute( observer )
    print( body.name, body.az, body.alt, body.phase )

observePyEphem( observer, ephem.Moon() )
observePyEphem( observer, ephem.Venus() )

def obsverveSkyfield( locationAtNow, body, sun ):
    apparent = locationAtNow.observe( body ).apparent()
    alt, az, distance = apparent.altaz()
    illumination = apparent.fraction_illuminated( sun ) * 100
    print( ' '.join( body.target_name.split()[ 1 : ] ).title(), az, alt, illumination )

ephemeris = load( "de421.bsp" )

timeScale = load.timescale( builtin = True )
now = timeScale.utc( utcNow.year, utcNow.month, utcNow.day, utcNow.hour, utcNow.minute, utcNow.second )

location = wgs84.latlon( latitude, longitude, elevation )
locationAtNow = ( ephemeris[ "earth" ] + location ).at( now )

obsverveSkyfield( locationAtNow, ephemeris[ "moon" ], ephemeris[ "sun" ] )
obsverveSkyfield( locationAtNow, ephemeris[ "venus" ], ephemeris[ "sun" ] )

Getting the results:

2022-09-25 04:15:49.510725
Moon 305:35:14.2 31:39:07.6 0.7093373537063599
Venus 303:29:20.1 34:50:32.0 99.18865966796875
Moon 305deg 35' 20.6" 31deg 37' 39.1" 0.8262983964746395
Venus 303deg 29' 27.1" 34deg 49' 15.3" 0.41677288119522093

Comparing with Horizons (and this is my first time using so fingers crossed I've pressed the correct buttons):

...
Center geodetic : 151.283300,-33.916600,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
...
2022-Sep-25 04:15 *m  11 38 26.39 +07 01 22.0   -4.991   5.882    0.82994  0.00258709285924   0.2092554   10.4320 /L  169.5410   25.315206   127.74826   14.800949         n.a.     n.a.

The illumination % is 0.82994 which agrees with Skyfield. @brandon-rhodes I presume you'll bounce this over to PyEphem...

[sgbmzm]

Thank you very much everyone. Note that this will not work if the body is the sun itself. I think the illuminated should be 100, but the result is 50. So maybe should add a error message that arrived in case you are trying to calculate the percentage of the sun's illumination, or alternatively set it to 100 permanently.

percent = 100.0 * eph['earth'].at(time).observe(eph[body]).apparent().fraction_illuminated(eph['sun'])

[brandon-rhodes]

I think the illuminated should be 100, but the result is 50.

But whichever angle you look at the Sun, you only see half of its illuminated surface, so 50% is kind of hilariously correct, right? And so I might keep that result because it’s so much fun 🙂

[brandon-rhodes]

brandon-rhodes I presume you'll bounce this over to PyEphem...

I'm not sure I have time at this point to dive into XEphem's old C code and figure out how its calculations differ—since it uses old Meeus formulas for many items, and tends to start with geocentric positions and then use sines and cosines to adjust for where the observer is on the Earth's surface, it often strikes me as non-obvious how its calculations should be compared with Skyfield's. (It also uses VSOP87 instead of a modern ephemeris, but I would be surprised if that is a big enough difference to affect phase and illuminated fraction?)

But if anyone is interested in diving into the XEphem / PyEphem code, I always welcome refinements to its calculations!

[hidp123]

I think the illuminated should be 100, but the result is 50.

But whichever angle you look at the Sun, you only see half of its illuminated surface, so 50% is kind of hilariously correct, right? And so I might keep that result because it’s so much fun 🙂

When we say the moon is 100% illuminated then it's still half the moon 😊

[sgbmzm]

The sun is illuminated at 100 percent of its circumference. On the other hand, the moon is always lit only on 50 percent of its circumference: only on the side facing the sun. When the side of the moon that faces the earth, faces the sun, we see the entire illuminated part which is about 50 percent of the moon's circumference, and we call it 100% illuminated.

So the sun has 200% illuminated. But we only see 100 of what is there, so 100 out of 200 is 50 percent.

So everyone is right, but the question is what will the user understand.