Comparison between the is_sunlit() function and JPL Horizons

[ejmontiel]

I'm a recent user of Skyfield and it's been an excellent tool! I've been looking at satellite illuminations and was using JPL's Horizon app as a check. I was working through your ISS example from "Bluffton":

and when I compare it to JPL Horizons, it suggests that the ISS should be illuminated for the first set events

When I explored this further with the Hubble Space Telescope and stayed in the Geocentric frame, I also found a similar discrepancy

and there's a strong offset between where Skyfield and Horizons say the maximum and minimum illumination should occur. So, I went looking at the source code, specifically the "intersect_line_and_sphere" in geometry.py. I noticed your return of the quadratic roots

and was wondering what guarantees the minus_b variable to actually be -b in the quadratic formula. So, I did the obvious thing of trying an explicit call of "-minus_b" in the return

to see what kind of result I would get. When I re-run my code looking at the illumination for HST and compare it to Horizons

and was happy to find an agreement between the two methods (I'll note that around 30% illumination in Horizions is where Skyfield makes its binary call of sunlit or not sunlit). I spent some time trying to see if anyone else has ever raised and mainly wondering if I'm doing something wrong (with the Skyfield code as it comes downloaded, not my temporary modification). Apologies if this is more "discussion" then "issue".

[brandon-rhodes]

Could you re-examine your HORIZONS output and look for the section, probably down near the bottom, where it defines what the Illu% column means, and paste that definition here for us? I would suspect that the illumination percent designed for a planet or the Moon, and has nothing to do with whether an object is in the Earth's shadow.

[ejmontiel]

Of course:

--Illuminated fraction--
Percent of target objects' assumed circular disk illuminated by Sun
(phase), as seen by observer.
  Units: PERCENT
  Labels:  Illu%

I agree that Horizons has no idea about an object's size based on a TLE. I'm hoping it would make a satellite-sized circle to try and illuminate.

[ejmontiel]

I would suspect that the illumination percent designed for a planet or the Moon, and has nothing to do with whether an object is in the Earth's shadow.

I believe I'm starting to understand what you mean by this. I went to a couple of ISS tracking visualizers tonight. I happened to be able to watch the ISS as it came out of Earth's shadow traveling from over Africa towards central Asia, and then back into the shadow over the Pacific heading towards South America. The way Skyfield would naturally calculate intersect_line_and_sphere (e.g. how you coded it) for is_sunlit() agrees with what I watched. While both Horizons and my modification to intersect_line_and_sphere disagrees with that behavior.

This got me thinking about what Horizons returns for lunar eclipses (spoiler: Illu% stays near 100%, which we know it shouldn't). However, this did lead me to the correct observer table settings I need to be returning:

12. Angular separation/visibility
The angular separation between the center of the target object and the center of the (remote) primary body it revolves around, as seen by the observer, with target visibility code. The observer cannot be on the primary body.

Visibility codes (refers to limb-to-limb):

    /t = Transiting primary body disk    /O = Occulted by primary body disk
    /p = Partial umbral eclipse          /P = Occulted partial umbral eclipse
    /u = Total umbral eclipse            /U = Occulted total umbral eclipse
    /- = Target is the primary body      /* = None of above ("free and clear")

The radius of both primary and target body is taken to be the equatorial value (maximum, given a triaxial shape). Atmospheric effects and oblateness aspect are _NOT_ currently considered.  Light-time is considered.

  Labels:   ang-sep/v

Now going back to the Skyfield example with observing the ISS from Bluffton:

*******************************************************************************
 Date__(UT)__HR:MN:SC.fff     Azi____(a-app)___Elev      Illu%   ang-sep/v
**************************************************************************
2014-Jan-23 06:25:37.000  m  202.975934  30.240492   90.03142  433496.4/u
2014-Jan-23 06:26:58.000  e  138.572074   54.332623   92.03801  520110.2/u
2014-Jan-23 06:28:19.000  m   75.874248  29.900247   62.09970  431472.9/u

2014-Jan-23 12:54:56.000 *m  308.761366  30.384440   92.23595  432954.9/*
2014-Jan-23 12:56:27.000 *m   50.773339  85.566749   48.64725  631598.4/*
2014-Jan-23 12:57:58.000 *m  123.896314  29.617829    6.97462  431005.5/*

The first three noted as "/u = Total umbral eclipse", while the second three are "/* = None of above ("free and clear")" matches your example.

I'm sorry for not assuming the problem could be with Horizons or, in this case, my own user error in what I was asking for from it. Many thanks for taking the time to respond to my question, and for writing and maintaining this awesome software!