def convert_tscube_old(infile, outfile):
    """Convert between old and new TSCube formats."""
    inhdulist = fits.open(infile)

    # If already in the new-style format just write and exit
    if 'DLOGLIKE_SCAN' in inhdulist['SCANDATA'].columns.names:
        if infile != outfile:
            inhdulist.writeto(outfile, clobber=True)
        return

    # Get stuff out of the input file
    nrows = inhdulist['SCANDATA']._nrows
    nebins = inhdulist['EBOUNDS']._nrows
    npts = inhdulist['SCANDATA'].data.field('NORMSCAN').shape[1] / nebins

    emin = inhdulist['EBOUNDS'].data.field('e_min') / 1E3
    emax = inhdulist['EBOUNDS'].data.field('e_max') / 1E3
    eref = np.sqrt(emin * emax)
    dnde_emin = inhdulist['EBOUNDS'].data.field('E_MIN_FL')
    dnde_emax = inhdulist['EBOUNDS'].data.field('E_MAX_FL')
    index = np.log(dnde_emin / dnde_emax) / np.log(emin / emax)

    flux = PowerLaw.eval_flux(emin, emax, [dnde_emin, index], emin)
    eflux = PowerLaw.eval_eflux(emin, emax, [dnde_emin, index], emin)
    dnde = PowerLaw.eval_dnde(np.sqrt(emin * emax), [dnde_emin, index], emin)

    ts_map = inhdulist['PRIMARY'].data.reshape((nrows))
    ok_map = inhdulist['TSMAP_OK'].data.reshape((nrows))
    n_map = inhdulist['N_MAP'].data.reshape((nrows))
    errp_map = inhdulist['ERRP_MAP'].data.reshape((nrows))
    errn_map = inhdulist['ERRN_MAP'].data.reshape((nrows))
    err_map = np.ndarray((nrows))
    m = errn_map > 0
    err_map[m] = 0.5 * (errp_map[m] + errn_map[m])
    err_map[~m] = errp_map[~m]
    ul_map = n_map + 2.0 * errp_map

    ncube = np.rollaxis(inhdulist['N_CUBE'].data,
                        0, 3).reshape((nrows, nebins))
    errpcube = np.rollaxis(
        inhdulist['ERRPCUBE'].data, 0, 3).reshape((nrows, nebins))
    errncube = np.rollaxis(
        inhdulist['ERRNCUBE'].data, 0, 3).reshape((nrows, nebins))
    tscube = np.rollaxis(inhdulist['TSCUBE'].data,
                         0, 3).reshape((nrows, nebins))
    nll_cube = np.rollaxis(
        inhdulist['NLL_CUBE'].data, 0, 3).reshape((nrows, nebins))
    ok_cube = np.rollaxis(
        inhdulist['TSCUBE_OK'].data, 0, 3).reshape((nrows, nebins))

    ul_cube = ncube + 2.0 * errpcube
    m = errncube > 0
    errcube = np.ndarray((nrows, nebins))
    errcube[m] = 0.5 * (errpcube[m] + errncube[m])
    errcube[~m] = errpcube[~m]

    norm_scan = inhdulist['SCANDATA'].data.field(
        'NORMSCAN').reshape((nrows, npts, nebins)).swapaxes(1, 2)
    nll_scan = inhdulist['SCANDATA'].data.field(
        'NLL_SCAN').reshape((nrows, npts, nebins)).swapaxes(1, 2)

    # Adjust the "EBOUNDS" hdu
    columns = inhdulist['EBOUNDS'].columns
    columns.add_col(fits.Column(name=str('e_ref'),
                                format='E', array=eref * 1E3,
                                unit='keV'))
    columns.add_col(fits.Column(name=str('ref_flux'),
                                format='D', array=flux,
                                unit='ph / (cm2 s)'))
    columns.add_col(fits.Column(name=str('ref_eflux'),
                                format='D', array=eflux,
                                unit='MeV / (cm2 s)'))
    columns.add_col(fits.Column(name=str('ref_dnde'),
                                format='D', array=dnde,
                                unit='ph / (MeV cm2 s)'))

    columns.change_name('E_MIN_FL', str('ref_dnde_e_min'))
    columns.change_unit('ref_dnde_e_min', 'ph / (MeV cm2 s)')
    columns.change_name('E_MAX_FL', str('ref_dnde_e_max'))
    columns.change_unit('ref_dnde_e_max', 'ph / (MeV cm2 s)')
    columns.change_name('NPRED', str('ref_npred'))

    hdu_e = fits.BinTableHDU.from_columns(columns, name='EBOUNDS')

    # Make the "FITDATA" hdu
    columns = fits.ColDefs([])

    columns.add_col(fits.Column(
        name=str('fit_ts'), format='E', array=ts_map))
    columns.add_col(fits.Column(
        name=str('fit_status'), format='E', array=ok_map))
    columns.add_col(fits.Column(
        name=str('fit_norm'), format='E', array=n_map))
    columns.add_col(fits.Column(
        name=str('fit_norm_err'), format='E', array=err_map))
    columns.add_col(fits.Column(
        name=str('fit_norm_errp'), format='E', array=errp_map))
    columns.add_col(fits.Column(
        name=str('fit_norm_errn'), format='E', array=errn_map))
    hdu_f = fits.BinTableHDU.from_columns(columns, name='FITDATA')

    # Make the "SCANDATA" hdu
    columns = fits.ColDefs([])

    columns.add_col(fits.Column(name=str('ts'),
                                format='%iE' % nebins, array=tscube,
                                dim=str('(%i)' % nebins)))

    columns.add_col(fits.Column(name=str('bin_status'),
                                format='%iE' % nebins, array=ok_cube,
                                dim=str('(%i)' % nebins)))

    columns.add_col(fits.Column(name=str('norm'),
                                format='%iE' % nebins, array=ncube,
                                dim=str('(%i)' % nebins)))

    columns.add_col(fits.Column(name=str('norm_ul'),
                                format='%iE' % nebins, array=ul_cube,
                                dim=str('(%i)' % nebins)))

    columns.add_col(fits.Column(name=str('norm_err'),
                                format='%iE' % nebins, array=errcube,
                                dim=str('(%i)' % nebins)))

    columns.add_col(fits.Column(name=str('norm_errp'),
                                format='%iE' % nebins, array=errpcube,
                                dim=str('(%i)' % nebins)))

    columns.add_col(fits.Column(name=str('norm_errn'),
                                format='%iE' % nebins, array=errncube,
                                dim=str('(%i)' % nebins)))

    columns.add_col(fits.Column(name=str('loglike'),
                                format='%iE' % nebins, array=nll_cube,
                                dim=str('(%i)' % nebins)))

    columns.add_col(fits.Column(name=str('norm_scan'),
                                format='%iE' % (nebins * npts),
                                array=norm_scan,
                                dim=str('(%i,%i)' % (npts, nebins))))

    columns.add_col(fits.Column(name=str('dloglike_scan'),
                                format='%iE' % (nebins * npts),
                                array=nll_scan,
                                dim=str('(%i,%i)' % (npts, nebins))))

    hdu_s = fits.BinTableHDU.from_columns(columns, name='SCANDATA')

    hdulist = fits.HDUList([inhdulist[0],
                            hdu_s,
                            hdu_f,
                            inhdulist["BASELINE"],
                            hdu_e])

    hdulist['SCANDATA'].header['UL_CONF'] = 0.95

    hdulist.writeto(outfile, clobber=True)

    return hdulist