def fill_livetime_hist(skydir, tab_sc, tab_gti, zmax, costh_edges):
    """Generate a sequence of livetime distributions at the sky
    positions given by ``skydir``.  The output of the method are two
    NxM arrays containing a sequence of histograms for N sky positions
    and M incidence angle bins where the bin edges are defined by
    ``costh_edges``.  This method uses the same algorithm as
    `gtltcube` with the exception that SC time intervals are assumed
    to be aligned with GTIs.

    Parameters
    ----------
    skydir : `~astropy.coordinates.SkyCoord`    
        Vector of sky directions for which livetime histograms will be
        accumulated.

    tab_sc : `~astropy.table.Table`    
        Spacecraft table.  Must contain the following columns: START,
        STOP, LIVETIME, RA_SCZ, DEC_SZ, RA_ZENITH, DEC_ZENITH.

    tab_gti : `~astropy.table.Table`
        Table of good time intervals (GTIs).

    zmax : float
        Zenith cut.

    costh_edges : `~numpy.ndarray`
        Incidence angle bin edges in cos(angle).

    Returns
    -------
    lt : `~numpy.ndarray`
        Array of livetime histograms.

    lt_wt : `~numpy.ndarray`    
        Array of histograms of weighted livetime (livetime x livetime
        fraction).
    """

    if len(tab_gti) == 0:
        shape = (len(costh_edges) - 1, len(skydir))
        return (np.zeros(shape), np.zeros(shape))

    m = (tab_sc['START'] < tab_gti['STOP'][-1])
    m &= (tab_sc['STOP'] > tab_gti['START'][0])
    tab_sc = tab_sc[m]

    cos_zmax = np.cos(np.radians(zmax))
    sc_t0 = np.array(tab_sc['START'].data)
    sc_t1 = np.array(tab_sc['STOP'].data)
    sc_live = np.array(tab_sc['LIVETIME'].data)
    sc_lfrac = sc_live / (sc_t1 - sc_t0)

    sc_xyz = angle_to_cartesian(np.radians(tab_sc['RA_SCZ'].data),
                                np.radians(tab_sc['DEC_SCZ'].data))
    zn_xyz = angle_to_cartesian(np.radians(tab_sc['RA_ZENITH'].data),
                                np.radians(tab_sc['DEC_ZENITH'].data))

    tab_gti_t0 = np.array(tab_gti['START'].data)
    tab_gti_t1 = np.array(tab_gti['STOP'].data)

    # Index of the closest GTI interval
    idx = np.digitize(sc_t0, tab_gti_t0) - 1

    # start/stop time of closest GTI interval
    gti_t0 = np.zeros_like(sc_t0)
    gti_t1 = np.zeros_like(sc_t1)
    gti_t0[idx >= 0] = tab_gti_t0[idx[idx >= 0]]
    gti_t1[idx >= 0] = tab_gti_t1[idx[idx >= 0]]

    nbin = len(costh_edges) - 1

    lt = np.zeros((nbin,) + skydir.shape)
    lt_wt = np.zeros((nbin,) + skydir.shape)

    m0 = (idx >= 0) & (sc_t0 >= gti_t0) & (sc_t1 <= gti_t1)

    xyz = angle_to_cartesian(skydir.ra.rad, skydir.dec.rad)

    for i, t in enumerate(xyz):
        cos_sep = utils.dot_prod(t, sc_xyz)
        cos_zn = utils.dot_prod(t, zn_xyz)
        m = m0 & (cos_zn > cos_zmax) & (cos_sep > 0.0)
        bins = np.digitize(cos_sep[m], bins=costh_edges) - 1
        bins = np.clip(bins, 0, nbin - 1)
        lt[:, i] = np.bincount(bins, weights=sc_live[m], minlength=nbin)
        lt_wt[:, i] = np.bincount(bins, weights=sc_live[m] * sc_lfrac[m],
                                  minlength=nbin)

    return lt, lt_wt
Generate a sequence of livetime distributions at the sky
    positions given by ``skydir``.  The output of the method are two
    NxM arrays containing a sequence of histograms for N sky positions
    and M incidence angle bins where the bin edges are defined by
    ``costh_edges``.  This method uses the same algorithm as
    `gtltcube` with the exception that SC time intervals are assumed
    to be aligned with GTIs.

    Parameters
    ----------
    skydir : `~astropy.coordinates.SkyCoord`    
        Vector of sky directions for which livetime histograms will be
        accumulated.

    tab_sc : `~astropy.table.Table`    
        Spacecraft table.  Must contain the following columns: START,
        STOP, LIVETIME, RA_SCZ, DEC_SZ, RA_ZENITH, DEC_ZENITH.

    tab_gti : `~astropy.table.Table`
        Table of good time intervals (GTIs).

    zmax : float
        Zenith cut.

    costh_edges : `~numpy.ndarray`
        Incidence angle bin edges in cos(angle).

    Returns
    -------
    lt : `~numpy.ndarray`
        Array of livetime histograms.

    lt_wt : `~numpy.ndarray`    
        Array of histograms of weighted livetime (livetime x livetime
        fraction).