Generate a Day/Night boundary polygon

[hamiltoncj]

Earlier I asked a question similar to this and got an answer that at the time I thought worked but it didn't. I have been trying to incorporate some skyfield capabilities into QGIS (https://qgis.org). See the discussion: #621. The one algorithm I am missing is to be able to plot the day/night boundary for civil twilight, nautical twilight, astronomical twilight, and night. I can find where the sun is at its azimuth with this code shared by @JoshPaterson .

from skyfield.api import load, wgs84
import numpy as np
from skyfield.positionlib import Geocentric

eph = load('de421.bsp')
earth = eph['earth'] # vector from solar system barycenter to geocenter
sun = eph['sun'] # vector from solar system barycenter to sun
geocentric_sun = sun - earth # vector from geocenter to sun

ts = load.timescale()
t = ts.utc(2021, 5, 19, 12, 26)

sun_subpoint = wgs84.subpoint(geocentric_sun.at(t)) # subpoint method requires a geocentric position

print('subpoint latitude: ', sun_subpoint.latitude.degrees)
print('subpoint longitude: ', sun_subpoint.longitude.degrees)

I am guessing that the way to solve this is to start at the sun's azimuth coordinate on earth, then travel a certain distance in every angle and that would give the area being illuminated by the sun. How do I know where that illumination boundary point is located at starting from the solar azimuth traveling in some direction or angle?

The original code suggested seems to be close for the solstice, but not for equinox dates. This is what the original suggested code generated for various dates.

For the summer solstice it should really look something like this using the lightest boundary.

For the fall equinox it should be the lightest boundary in this image.

Can anyone help out? Math is not my strong point especially celestial math.

Thanks!