def ztwilight(alt) :
    """Estimate twilight contribution to zenith sky brightness, in magnitudes
    per square arcsecond.  

    Evaluates a polynomial approximation to observational data (see source for 
    reference) of zenith sky brightness (blue) as a function of the sun's elevation
    from -0.9 degrees to -18 degrees.  For reference, 3 mag is roughly Nautical
    twilight and looks 'pretty dark'; something like 10 mag is about the maximum
    for broadband sky flats in many cases.

    Parameters
    ----------
    alt : Angle
        Sun's elevation.  Meaningful range is -0.9 to -18 degrees.    

    Returns 
        float.  If the sun is up, returns 20; if the sun below -18, returns 0.

    """
# Given an Angle alt in the range -0.9 to -18 degrees, 
# evaluates a polynomial expansion for the approximate brightening
# in magnitudes of the zenith in twilight compared to its 
# value at full night, as function of altitude of the sun (in degrees).
# To get this expression I looked in Meinel, A.,
# & Meinel, M., "Sunsets, Twilight, & Evening Skies", Cambridge U.
# Press, 1983; there's a graph on p. 38 showing the decline of 
# zenith twilight.  I read points off this graph and fit them with a
# polynomial.
# Comparison with Ashburn, E. V. 1952, JGR, v.57, p.85 shows that this
# is a good fit to his B-band measurements.  

    if alt.deg > -0.9 : return 20.  # flag for sun up, also not grossly wrong.
    if alt.deg < -18. : return 0.   # fully dark, no contrib to zenith skyglow.

    y = (-1. * alt.deg - 9.0) / 9.0  # my polynomial's argument
    return ((2.0635175 * y + 1.246602) * y - 9.4084495)*y + 6.132725