Different/better ways to calculate the shadow region passage of the satellite and the effect of orbital elements, other than shadow cone technique?

The NASA Technical Publication Method for the calculation of spacecraft umbra and penumbra shadow terminator points; TP-3547 describes the following:

A method for calculating orbital shadow terminator points is presented. The current method employs the use of an iterative process which is used for an accurate determination of shadow points. This calculation methodology is required since orbital perturbation effects can introduce large errors when a spacecraft orbits a planet in a high altitude and/or highly elliptical orbit. To compensate for the required iteration methodology, all reference frame change definitions and calculations are performed with quaternions. Quaternion algebra significantly reduces the computational time required for the accurate determination of shadow terminator points.

Are there any other, different, better and/or alternative ways to calculate the the shadow region passage of the satellite and its effect over orbital elements?

• @uhoh Thank you for your help. I'm expecting an equation-based answer to use for my orbital propagator. Commented Feb 28, 2022 at 3:28

TL;DR: No, some kind of iteration is required the find the exact entrance or exit of the penumbra.

To calculate whether or not a point on an orbit is in the shadow, one does not need any iteration: the shadow cones from your referenced paper explains it well (for an HTML version of a similar algorithm, one may refer to the MathSpec of Nyx and Basilisk).

However, this algorithm requires the knowledge of the spacecraft state, and further, you wish you find the exact location of where the shadow region starts. The problem is that the point on an orbit isn't only "in the shadow" or "in the light": there's a third state, the "penumbra". Entry/exit of the penumbra is the definition of the shadow region (as the solar panels will generate less energy, etc.). As explained in the MathSpecs linked above, this calculation corresponds to a percentage of illumination from the shadow bodies. Depending on the percentage of illumination, a spacecraft will be subjected to more or less solar radiation pressure (Nyx MathSpec), which in turn will affect its orbital elements as described in your reference.

The method in your reference is quite clever as it analyzes the change in orbital elements directly. A more typical method (in my experience), is to simply propagate an orbit and analyze the influence of different orbital dynamics on the shape of that orbit.

In fact, that's how one would handle that in Nyx, cf. this custom example. Propagate a spacecraft forward in time, generate an interpolated trajectory, and search for a specific event along that trajectory using a Brent solver:

Min event: [Earth J2000] 2022-02-28T00:42:13 UTC        sma = 6828.135682 km    ecc = 0.001000  inc = 63.400003 deg     raan = 135.000001 deg   aop = 90.000132 deg     ta = 222.218854 deg => Penumbra 0.00%
Min event - 1 s: [Earth J2000] 2022-02-28T00:42:12 UTC  sma = 6828.135689 km    ecc = 0.001000  inc = 63.400003 deg     raan = 135.000001 deg   aop = 90.000106 deg     ta = 222.154863 deg => Penumbra 6.16%
Min event + 1 s: [Earth J2000] 2022-02-28T00:42:14 UTC  sma = 6828.135674 km    ecc = 0.001000  inc = 63.400003 deg     raan = 135.000001 deg   aop = 90.000158 deg     ta = 222.282845 deg => Umbra

Max event: [Earth J2000] 2022-02-28T00:23:24 UTC        sma = 6828.137397 km    ecc = 0.001000  inc = 63.400003 deg     raan = 135.000001 deg   aop = 90.010516 deg     ta = 149.960681 deg => Visibilis
Max event - 1 s: [Earth J2000] 2022-02-28T00:23:23 UTC  sma = 6828.137265 km    ecc = 0.001000  inc = 63.400003 deg     raan = 135.000001 deg   aop = 90.010071 deg     ta = 149.897126 deg => Visibilis
Max event + 1 s: [Earth J2000] 2022-02-28T00:23:25 UTC  sma = 6828.137526 km    ecc = 0.001000  inc = 63.400003 deg     raan = 135.000001 deg   aop = 90.010950 deg     ta = 150.024248 deg => Visibilis