How to use the fitsio.FITSHDR function in fitsio

plt.clf()
                img = tim.getImage()
                mn,mx = [np.percentile(img,p) for p in [25,99]]
                dimshow(img, vmin=mn, vmax=mx)
                xx,yy = [],[]
                for src in srcs:
                    x,y = tim.getWcs().positionToPixel(src.getPosition())
                    xx.append(x)
                    yy.append(y)
                ax = plt.axis()
                plt.plot(xx, yy, 'r+')
                plt.axis(ax)
                plt.savefig('tim-%s%i.png' % (band, ifield))

            # Populate the FITS header we will put in the output file.
            hdr = fitsio.FITSHDR()
            hdr.add_record(dict(name='RUN', value=run, comment='SDSS run'))
            hdr.add_record(dict(name='CAMCOL', value=camcol, comment='SDSS camcol'))
            hdr.add_record(dict(name='FIELD', value=field, comment='SDSS field'))
            hdr.add_record(dict(name='BAND', value=band, comment='SDSS band'))

            # Copy from input "frame" header
            orighdr = tinfo['hdr']
            for key in ['NMGY']:
                hdr.add_record(dict(name=key, value=orighdr[key],
                                    comment=orighdr.get_comment(key)))

            hdr.add_record(dict(name='CALIB', value=tim.sdss_calib,
                                comment='Mean "calibvec" value for this image'))
            hdr.add_record(dict(name='SKY', value=tim.sdss_sky,
                                comment='SDSS sky estimate at image center'))
            hdr.add_record(dict(name='GAIN', value=tim.sdss_gain,

for (key, value, comment) in addToHeader:
                hdr.add_record(dict(name=key, value=value, comment=comment))
                
            # Convert tim's WCS to a local Tan.
            cd = tim.wcs.cdAtPixel(x, y)
            tanwcs = Tan(ra, dec, fx, fy,
                         cd[0,0], cd[0,1], cd[1,0], cd[1,1],
                         x1-x0, y1-y0)
            tim.wcs = ConstantFitsWcs(tanwcs)

            imghdr = fitsio.FITSHDR()
            imghdr.add_record(dict(name='DATATYPE', value='IMAGE',
                                   comment='SDSS image pixels in nanomaggies'))
            tanwcs.add_to_header(imghdr)

            ivhdr = fitsio.FITSHDR()
            ivhdr.add_record(dict(name='DATATYPE', value='INVVAR',
                                   comment='SDSS image pixel inverse variance in 1/nmgy^2'))
            tanwcs.add_to_header(ivhdr)
            
            fn = stamp_pattern % band
            print('writing', fn)

            F = fitsio.FITS(fn, mode='rw', clobber=True)
            tim.toFits(F, primheader=hdr,
                       imageheader=imghdr, invvarheader=ivhdr)
                       
            continue
        
            
            # Resample to common grid
            th,tw = tim.shape

for iband,band in enumerate(bands):
        i = allbands.index(band)
        TT.decam_apflux[:,i,:] = AP.get('apflux_img_%s' % band)
        TT.decam_apflux_ivar[:,i,:] = AP.get('apflux_img_ivar_%s' % band)
        TT.decam_apflux_resid[:,i,:] = AP.get('apflux_resid_%s' % band)

    cat.thawAllRecursive()
    hdr = None
    T2,hdr = prepare_fits_catalog(cat, invvars, TT, hdr, bands, fs,
                                  allbands=allbands)

    # mod
    T2.ra += (T2.ra <   0) * 360.
    T2.ra -= (T2.ra > 360) * 360.

    primhdr = fitsio.FITSHDR()
    for r in version_header.records():
        primhdr.add_record(r)

    primhdr.add_record(dict(name='ALLBANDS', value=allbands,
                            comment='Band order in array values'))

    for i,ap in enumerate(apertures_arcsec):
        primhdr.add_record(dict(name='APRAD%i' % i, value=ap,
                                comment='Aperture radius, in arcsec'))

    bits = CP_DQ_BITS.values()
    bits.sort()
    bitmap = dict((v,k) for k,v in CP_DQ_BITS.items())
    for i in range(16):
        bit = 1<

hdrs = []
    #s = ''
    hdr = None
    for line in f.readlines():
        #print('Read %i chars' % len(line), ': ' + line[:25]+'...')
        line = line.strip()
        #print('Stripped to %i' % len(line))
        #line = line + ' ' * (80 - len(line))
        #assert(len(line) == 80)
        #s += line
        if line.startswith('END'):
            #print('Found END')
            hdr = None
            continue
        if hdr is None:
            hdr = fitsio.FITSHDR()
            hdrs.append(hdr)
            #print('Started header number', len(hdrs))
        hdr.add_record(line)
    return hdrs

def write_fits(self, filename, hdr=None):
        import fitsio
        
        tt = type(self)
        sky_type = '%s.%s' % (tt.__module__, tt.__name__)
        if hdr is None:
            hdr = fitsio.FITSHDR()
        hdr.add_record(dict(name='SKY', value=sky_type,
                            comment='Sky class'))
        self.toFitsHeader(hdr, prefix='SKY_')

        fitsio.write(filename, None, header=hdr, clobber=True)

#version_hdr.add_record(dict(name='CPFILE', value=im.imgfn, comment='DECam comm.pipeline file'))
    version_hdr.add_record(dict(name='CPFILE', value=fname, comment='DECam comm.pipeline file'))
    version_hdr.add_record(dict(name='CPHDU', value=im.hdu, comment='DECam comm.pipeline ext'))
    version_hdr.add_record(dict(name='CAMERA', value='DECam', comment='Dark Energy Camera'))
    version_hdr.add_record(dict(name='EXPNUM', value=im.expnum, comment='DECam exposure num'))
    version_hdr.add_record(dict(name='CCDNAME', value=im.extname, comment='DECam CCD name'))
    version_hdr.add_record(dict(name='FILTER', value=tim.band, comment='Bandpass of this image'))
    version_hdr.add_record(dict(name='EXPOSURE', value='decam-%s-%s' % (im.expnum, im.extname), comment='Name of this image'))
    
    keys = ['TELESCOP','OBSERVAT','OBS-LAT','OBS-LONG','OBS-ELEV',
            'INSTRUME']
    for key in keys:
        if key in tim.primhdr:
            version_hdr.add_record(dict(name=key, value=tim.primhdr[key]))

    hdr = fitsio.FITSHDR()

    units = {'mjd':'sec', 'exptime':'sec', 'flux':'nanomaggy',
             'flux_ivar':'1/nanomaggy^2'}
    columns = F.get_columns()
    for i,col in enumerate(columns):
        if col in units:
            hdr.add_record(dict(name='TUNIT%i' % (i+1), value=units[col]))

    outdir = os.path.dirname(outfn)
    if len(outdir):
        trymakedirs(outdir)
    fitsio.write(outfn, None, header=version_hdr, clobber=True)
    F.writeto(outfn, header=hdr, append=True)
    print 'Wrote', outfn

    print 'Finished forced phot:', Time()-t0

plt.ylabel('Tractor W1 mag')
            plt.title('WISE catalog vs Tractor forced photometry')
            plt.axis([cathi,lo,cathi,lo])
            ps.savefig()

        print 'Tile', tile.coadd_id, 'band', wband, 'took', Time()-tb0

    T.cut(inbounds)

    T.delete_column('psfflux')
    T.delete_column('cmodelflux')
    T.delete_column('devflux')
    T.delete_column('expflux')
    T.treated_as_pointsource = T.treated_as_pointsource.astype(np.uint8)

    hdr = fitsio.FITSHDR()
    hdr.add_record(dict(name='SEQ_VER', value=version['Revision'],
                        comment='SVN revision'))
    hdr.add_record(dict(name='SEQ_URL', value=version['URL'], comment='SVN URL'))
    hdr.add_record(dict(name='SEQ_DATE', value=datetime.datetime.now().isoformat(),
                        comment='forced phot run time'))
    hdr.add_record(dict(name='SEQ_NNEG', value=opt.nonneg, comment='non-negative?'))
    hdr.add_record(dict(name='SEQ_SKY', value=opt.sky, comment='fit sky?'))
    for band in bands:
        minsig = getattr(opt, 'minsig%i' % band)
        hdr.add_record(dict(name='SEQ_MNS%i' % band, value=minsig,
                            comment='min surf brightness in sig, band %i' % band))
    hdr.add_record(dict(name='SEQ_BL', value=opt.blocks, comment='image blocks'))
    hdr.add_record(dict(name='SEQ_CERE', value=opt.ceres, comment='use Ceres?'))
    hdr.add_record(dict(name='SEQ_ERRF', value=opt.errfrac, comment='error flux fraction'))
    if opt.ceres:
        hdr.add_record(dict(name='SEQ_CEBL', value=opt.ceresblock,

def main(opt, cs82field):
    t0 = Time()
    
    bands = opt.bands

    if opt.plots:
        ps = PlotSequence(opt.prefix)
    else:
        ps = None
    
    version = get_git_version()
    print('git version info:', version)

    hdr = fitsio.FITSHDR()
    hdr.add_record(dict(name='PHO_VER', value=version['commit'],
                        comment='cs82.py code version (git commit)'))
    if 'describe' in version:
        hdr.add_record(dict(name='PHO_GIT', value=version['describe'],
                            comment='cs82.py code version (git describe)'))
    hdr.add_record(dict(name='PHO_DATE', value=datetime.datetime.now().isoformat(),
                        comment='cs82.py photometry run time'))
    hdr.add_record(dict(name='PHO_MACH', value=socket.getfqdn(),
                        comment='machine where phot was run'))

    T,F = getTables(cs82field, enclosed=False)

    masks = get_cs82_masks(cs82field)
    
    sdss = DR9(basedir='data/unzip')
    if opt.local:

# Running out of memory can cause failure to converge
                # and term status = 2.
                # Fail completely in this case.
                if term != 0:
                    raise RuntimeError(
                        'Ceres terminated with status %i' % term)

            if wantims:
                ims0 = R.ims0
                ims1 = R.ims1
            IV, fs = R.IV, R.fitstats

            if save_fits:
                import fitsio
                (dat, mod, ie, chi, roi) = ims1[0]
                wcshdr = fitsio.FITSHDR()
                tim.wcs.wcs.add_to_header(wcshdr)

                tag = 'fit-%s-w%i' % (tile.coadd_id, band)
                fitsio.write('%s-data.fits' %
                             tag, dat, clobber=True, header=wcshdr)
                fitsio.write('%s-mod.fits' % tag,  mod,
                             clobber=True, header=wcshdr)
                fitsio.write('%s-chi.fits' % tag,  chi,
                             clobber=True, header=wcshdr)

            if get_models:
                (dat, mod, ie, chi, roi) = ims1[0]
                models[(tile.coadd_id, band)] = (mod, tim.roi)

            if ps:
                tag = '%s W%i' % (tile.coadd_id, band)

def getStandardFitsHeader(self, header=None):
        if header is None:
            import fitsio
            hdr = fitsio.FITSHDR()
        else:
            hdr = header
        psf = self.getPsf()
        wcs = self.getWcs()
        sky = self.getSky()
        pcal = self.getPhotoCal()
        psf.toStandardFitsHeader(hdr)
        wcs.toStandardFitsHeader(hdr)
        sky.toStandardFitsHeader(hdr)
        pcal.toStandardFitsHeader(hdr)
        return hdr