Inputs:
    time: Numpy 1-D array of datetime.datetime OR string for dateutil.parser.parse
    glon: Numpy 2-D array of geodetic longitudes (degrees)

    Outputs:
    iyd: day of year
    utsec: seconds from midnight utc
    stl: local solar time

def datetime2gtd(time: Union[str, datetime.datetime, np.datetime64],
                 glon: Union[float, List[float], np.ndarray] = np.nan) -> Tuple[int, float, float]:
    """
    Inputs:
    time: Numpy 1-D array of datetime.datetime OR string for dateutil.parser.parse
    glon: Numpy 2-D array of geodetic longitudes (degrees)

    Outputs:
    iyd: day of year
    utsec: seconds from midnight utc
    stl: local solar time
    """
# %%
    T = np.atleast_1d(time)
    glon = np.asarray(glon)
    doy = np.empty_like(T, int)
    utsec = np.empty_like(T, float)
    stl = np.empty((T.size, *glon.shape))

    for i, t in enumerate(T):
        if isinstance(t, str):
            t = parse(t)
        elif isinstance(t, np.datetime64):
            t = t.astype(datetime.datetime)
        elif isinstance(t, (datetime.datetime, datetime.date)):
            pass
        else:
            raise TypeError('unknown time datatype {}'.format(type(t)))
# %% Day of year
        doy[i] = int(t.strftime('%j'))
# %% seconds since utc midnight
        utsec[i] = datetime2utsec(t)

        stl[i, ...] = utsec[i] / 3600. + glon / 15.

    return doy, utsec, stl