def specimens_results_magic(infile='pmag_specimens.txt', measfile='magic_measurements.txt', sampfile='er_samples.txt', sitefile='er_sites.txt', agefile='er_ages.txt', specout='er_specimens.txt', sampout='pmag_samples.txt', siteout='pmag_sites.txt', resout='pmag_results.txt', critout='pmag_criteria.txt', instout='magic_instruments.txt', plotsites=False, fmt='svg', dir_path='.', cors=[], priorities=['DA-AC-ARM', 'DA-AC-TRM'], coord='g', user='', vgps_level='site', do_site_intensity=True, DefaultAge=["none"], avg_directions_by_sample=False, avg_intensities_by_sample=False, avg_all_components=False, avg_by_polarity=False, skip_directions=False, skip_intensities=False, use_sample_latitude=False, use_paleolatitude=False, use_criteria='default'):
    """
    Writes magic_instruments, er_specimens, pmag_samples, pmag_sites, pmag_criteria, and pmag_results. The data used to write this is obtained by reading a pmag_speciemns, a magic_measurements, a er_samples, a er_sites, a er_ages.
    @param -> infile: path from the WD to the pmag speciemns table
    @param -> measfile: path from the WD to the magic measurement file
    @param -> sampfile: path from the WD to the er sample file
    @param -> sitefile: path from the WD to the er sites data file
    @param -> agefile: path from the WD to the er ages data file
    @param -> specout: path from the WD to the place to write the er specimens data file
    @param -> sampout: path from the WD to the place to write the pmag samples data file
    @param -> siteout: path from the WD to the place to write the pmag sites data file
    @param -> resout: path from the WD to the place to write the pmag results data file
    @param -> critout: path from the WD to the place to write the pmag criteria file
    @param -> instout: path from th WD to the place to write the magic instruments file
    @param -> documentation incomplete if you know more about the purpose of the parameters in this function and it's side effects please extend and complete this string
    """
    # initialize some variables
    plotsites = False  # cannot use draw_figs from within ipmag
    Comps = []  # list of components
    version_num = pmag.get_version()
    args = sys.argv
    model_lat_file = ""
    Dcrit, Icrit, nocrit = 0, 0, 0
    corrections = []
    nocorrection = ['DA-NL', 'DA-AC', 'DA-CR']

    # do some data adjustments
    for cor in cors:
        nocorrection.remove('DA-' + cor)
        corrections.append('DA-' + cor)

    for p in priorities:
        if not p.startswith('DA-AC-'):
            p = 'DA-AC-' + p

    # translate coord into coords
    if coord == 's':
        coords = ['-1']
    if coord == 'g':
        coords = ['0']
    if coord == 't':
        coords = ['100']
    if coord == 'b':
        coords = ['0', '100']

    if vgps_level == 'sample':
        vgps = 1  # save sample level VGPS/VADMs
    else:
        vgps = 0  # site level

    if do_site_intensity:
        nositeints = 0
    else:
        nositeints = 1

    # chagne these all to True/False instead of 1/0

    if not skip_intensities:
        # set model lat and
        if use_sample_latitude and use_paleolatitude:
            print("you should set a paleolatitude file OR use present day lat - not both")
            return False
        elif use_sample_latitude:
            get_model_lat = 1
        elif use_paleolatitude:
            get_model_lat = 2
            try:
                model_lat_file = dir_path + '/' + args[ind + 1]
                get_model_lat = 2
                mlat = open(model_lat_file, 'r')
                ModelLats = []
                for line in mlat.readlines():
                    ModelLat = {}
                    tmp = line.split()
                    ModelLat["er_site_name"] = tmp[0]
                    ModelLat["site_model_lat"] = tmp[1]
                    ModelLat["er_sample_name"] = tmp[0]
                    ModelLat["sample_lat"] = tmp[1]
                    ModelLats.append(ModelLat)
                mlat.clos()
            except:
                print("use_paleolatitude option requires a valid paleolatitude file")
        else:
            get_model_lat = 0  # skips VADM calculation entirely

    if plotsites and not skip_directions:  # plot by site - set up plot window
        EQ = {}
        EQ['eqarea'] = 1
        # define figure 1 as equal area projection
        pmagplotlib.plot_init(EQ['eqarea'], 5, 5)
        # I don't know why this has to be here, but otherwise the first plot
        # never plots...
        pmagplotlib.plot_net(EQ['eqarea'])
        pmagplotlib.draw_figs(EQ)

    infile = os.path.join(dir_path, infile)
    measfile = os.path.join(dir_path, measfile)
    instout = os.path.join(dir_path, instout)
    sampfile = os.path.join(dir_path, sampfile)
    sitefile = os.path.join(dir_path, sitefile)
    agefile = os.path.join(dir_path, agefile)
    specout = os.path.join(dir_path, specout)
    sampout = os.path.join(dir_path, sampout)
    siteout = os.path.join(dir_path, siteout)
    resout = os.path.join(dir_path, resout)
    critout = os.path.join(dir_path, critout)

    if use_criteria == 'none':
        Dcrit, Icrit, nocrit = 1, 1, 1  # no selection criteria
        crit_data = pmag.default_criteria(nocrit)
    elif use_criteria == 'default':
        crit_data = pmag.default_criteria(nocrit)  # use default criteria
    elif use_criteria == 'existing':
        crit_data, file_type = pmag.magic_read(
            critout)  # use pmag_criteria file
        print("Acceptance criteria read in from ", critout)
    accept = {}
    for critrec in crit_data:
        for key in list(critrec.keys()):
            # need to migrate specimen_dang to specimen_int_dang for intensity
            # data using old format
            if 'IE-SPEC' in list(critrec.keys()) and 'specimen_dang' in list(critrec.keys()) and 'specimen_int_dang' not in list(critrec.keys()):
                critrec['specimen_int_dang'] = critrec['specimen_dang']
                del critrec['specimen_dang']
# need to get rid of ron shaars sample_int_sigma_uT
            if 'sample_int_sigma_uT' in list(critrec.keys()):
                critrec['sample_int_sigma'] = '%10.3e' % (
                    eval(critrec['sample_int_sigma_uT']) * 1e-6)
            if key not in list(accept.keys()) and critrec[key] != '':
                accept[key] = critrec[key]

    if use_criteria == 'default':
        pmag.magic_write(critout, [accept], 'pmag_criteria')
        print("\n Pmag Criteria stored in ", critout, '\n')

# now we're done slow dancing

    # read in site data - has the lats and lons
    SiteNFO, file_type = pmag.magic_read(sitefile)
    # read in site data - has the lats and lons
    SampNFO, file_type = pmag.magic_read(sampfile)
    # find all the sites with height info.
    height_nfo = pmag.get_dictitem(SiteNFO, 'site_height', '', 'F')
    if agefile:
        AgeNFO, file_type = pmag.magic_read(
            agefile)  # read in the age information
    # read in specimen interpretations
    Data, file_type = pmag.magic_read(infile)
    # retrieve specimens with intensity data
    IntData = pmag.get_dictitem(Data, 'specimen_int', '', 'F')
    comment, orient = "", []
    samples, sites = [], []
    for rec in Data:  # run through the data filling in missing keys and finding all components, coordinates available
        # fill in missing fields, collect unique sample and site names
        if 'er_sample_name' not in list(rec.keys()):
            rec['er_sample_name'] = ""
        elif rec['er_sample_name'] not in samples:
            samples.append(rec['er_sample_name'])
        if 'er_site_name' not in list(rec.keys()):
            rec['er_site_name'] = ""
        elif rec['er_site_name'] not in sites:
            sites.append(rec['er_site_name'])
        if 'specimen_int' not in list(rec.keys()):
            rec['specimen_int'] = ''
        if 'specimen_comp_name' not in list(rec.keys()) or rec['specimen_comp_name'] == "":
            rec['specimen_comp_name'] = 'A'
        if rec['specimen_comp_name'] not in Comps:
            Comps.append(rec['specimen_comp_name'])
        rec['specimen_tilt_correction'] = rec['specimen_tilt_correction'].strip(
            '\n')
        if "specimen_tilt_correction" not in list(rec.keys()):
            rec["specimen_tilt_correction"] = "-1"  # assume sample coordinates
        if rec["specimen_tilt_correction"] not in orient:
            # collect available coordinate systems
            orient.append(rec["specimen_tilt_correction"])
        if "specimen_direction_type" not in list(rec.keys()):
            # assume direction is line - not plane
            rec["specimen_direction_type"] = 'l'
        if "specimen_dec" not in list(rec.keys()):
            # if no declination, set direction type to blank
            rec["specimen_direction_type"] = ''
        if "specimen_n" not in list(rec.keys()):
            rec["specimen_n"] = ''  # put in n
        if "specimen_alpha95" not in list(rec.keys()):
            rec["specimen_alpha95"] = ''  # put in alpha95
        if "magic_method_codes" not in list(rec.keys()):
            rec["magic_method_codes"] = ''
     # start parsing data into SpecDirs, SpecPlanes, SpecInts
    SpecInts, SpecDirs, SpecPlanes = [], [], []
    samples.sort()  # get sorted list of samples and sites
    sites.sort()
    if not skip_intensities:  # don't skip intensities
        # retrieve specimens with intensity data
        IntData = pmag.get_dictitem(Data, 'specimen_int', '', 'F')
        if nocrit == 0:  # use selection criteria
            for rec in IntData:  # do selection criteria
                kill = pmag.grade(rec, accept, 'specimen_int')
                if len(kill) == 0:
                    # intensity record to be included in sample, site
                    # calculations
                    SpecInts.append(rec)
        else:
            SpecInts = IntData[:]  # take everything - no selection criteria
# check for required data adjustments
        if len(corrections) > 0 and len(SpecInts) > 0:
            for cor in corrections:
                # only take specimens with the required corrections
                SpecInts = pmag.get_dictitem(
                    SpecInts, 'magic_method_codes', cor, 'has')
        if len(nocorrection) > 0 and len(SpecInts) > 0:
            for cor in nocorrection:
                # exclude the corrections not specified for inclusion
                SpecInts = pmag.get_dictitem(
                    SpecInts, 'magic_method_codes', cor, 'not')
# take top priority specimen of its name in remaining specimens (only one
# per customer)
        PrioritySpecInts = []
        specimens = pmag.get_specs(SpecInts)  # get list of uniq specimen names
        for spec in specimens:
            # all the records for this specimen
            ThisSpecRecs = pmag.get_dictitem(
                SpecInts, 'er_specimen_name', spec, 'T')
            if len(ThisSpecRecs) == 1:
                PrioritySpecInts.append(ThisSpecRecs[0])
            elif len(ThisSpecRecs) > 1:  # more than one
                prec = []
                for p in priorities:
                    # all the records for this specimen
                    ThisSpecRecs = pmag.get_dictitem(
                        SpecInts, 'magic_method_codes', p, 'has')
                    if len(ThisSpecRecs) > 0:
                        prec.append(ThisSpecRecs[0])
                PrioritySpecInts.append(prec[0])  # take the best one
        SpecInts = PrioritySpecInts  # this has the first specimen record
    if not skip_directions:  # don't skip directions
        # retrieve specimens with directed lines and planes
        AllDirs = pmag.get_dictitem(Data, 'specimen_direction_type', '', 'F')
        # get all specimens with specimen_n information
        Ns = pmag.get_dictitem(AllDirs, 'specimen_n', '', 'F')
        if nocrit != 1:  # use selection criteria
            for rec in Ns:  # look through everything with specimen_n for "good" data
                kill = pmag.grade(rec, accept, 'specimen_dir')
                if len(kill) == 0:  # nothing killed it
                    SpecDirs.append(rec)
        else:  # no criteria
            SpecDirs = AllDirs[:]  # take them all
# SpecDirs is now the list of all specimen directions (lines and planes)
# that pass muster
    # list of all sample data and list of those that pass the DE-SAMP criteria
    PmagSamps, SampDirs = [], []
    PmagSites, PmagResults = [], []  # list of all site data and selected results
    SampInts = []
    for samp in samples:  # run through the sample names
        if avg_directions_by_sample:  # average by sample if desired
            # get all the directional data for this sample
            SampDir = pmag.get_dictitem(SpecDirs, 'er_sample_name', samp, 'T')
            if len(SampDir) > 0:  # there are some directions
                for coord in coords:  # step through desired coordinate systems
                    # get all the directions for this sample
                    CoordDir = pmag.get_dictitem(
                        SampDir, 'specimen_tilt_correction', coord, 'T')
                    if len(CoordDir) > 0:  # there are some with this coordinate system
                        if not avg_all_components:  # look component by component
                            for comp in Comps:
                                # get all directions from this component
                                CompDir = pmag.get_dictitem(
                                    CoordDir, 'specimen_comp_name', comp, 'T')
                                if len(CompDir) > 0:  # there are some
                                    # get a sample average from all specimens
                                    PmagSampRec = pmag.lnpbykey(
                                        CompDir, 'sample', 'specimen')
                                    # decorate the sample record
                                    PmagSampRec["er_location_name"] = CompDir[0]['er_location_name']
                                    PmagSampRec["er_site_name"] = CompDir[0]['er_site_name']
                                    PmagSampRec["er_sample_name"] = samp
                                    PmagSampRec["er_citation_names"] = "This study"
                                    PmagSampRec["er_analyst_mail_names"] = user
                                    PmagSampRec['magic_software_packages'] = version_num
                                    if CompDir[0]['specimen_flag'] == 'g':
                                        PmagSampRec['sample_flag'] = 'g'
                                    else:
                                        PmagSampRec['sample_flag'] = 'b'
                                    if nocrit != 1:
                                        PmagSampRec['pmag_criteria_codes'] = "ACCEPT"
                                    if agefile != "":
                                        PmagSampRec = pmag.get_age(
                                            PmagSampRec, "er_site_name", "sample_inferred_", AgeNFO, DefaultAge)
                                    site_height = pmag.get_dictitem(
                                        height_nfo, 'er_site_name', PmagSampRec['er_site_name'], 'T')
                                    if len(site_height) > 0:
                                        # add in height if available
                                        PmagSampRec["sample_height"] = site_height[0]['site_height']
                                    PmagSampRec['sample_comp_name'] = comp
                                    PmagSampRec['sample_tilt_correction'] = coord
                                    PmagSampRec['er_specimen_names'] = pmag.get_list(
                                        CompDir, 'er_specimen_name')  # get a list of the specimen names used
                                    PmagSampRec['magic_method_codes'] = pmag.get_list(
                                        CompDir, 'magic_method_codes')  # get a list of the methods used
                                    if nocrit != 1:  # apply selection criteria
                                        kill = pmag.grade(
                                            PmagSampRec, accept, 'sample_dir')
                                    else:
                                        kill = []
                                    if len(kill) == 0:
                                        SampDirs.append(PmagSampRec)
                                        if vgps == 1:  # if sample level VGP info desired, do that now
                                            PmagResRec = pmag.getsampVGP(
                                                PmagSampRec, SiteNFO)
                                            if PmagResRec != "":
                                                PmagResults.append(PmagResRec)
#                                   print(PmagSampRec)
                                    PmagSamps.append(PmagSampRec)
                        if avg_all_components:  # average all components together  basically same as above
                            PmagSampRec = pmag.lnpbykey(
                                CoordDir, 'sample', 'specimen')
                            PmagSampRec["er_location_name"] = CoordDir[0]['er_location_name']
                            PmagSampRec["er_site_name"] = CoordDir[0]['er_site_name']
                            PmagSampRec["er_sample_name"] = samp
                            PmagSampRec["er_citation_names"] = "This study"
                            PmagSampRec["er_analyst_mail_names"] = user
                            PmagSampRec['magic_software_packages'] = version_num
                            if all(i['specimen_flag'] == 'g' for i in CoordDir):
                                PmagSampRec['sample_flag'] = 'g'
                            else:
                                PmagSampRec['sample_flag'] = 'b'
                            if nocrit != 1:
                                PmagSampRec['pmag_criteria_codes'] = ""
                            if agefile != "":
                                PmagSampRec = pmag.get_age(
                                    PmagSampRec, "er_site_name", "sample_inferred_", AgeNFO, DefaultAge)
                            site_height = pmag.get_dictitem(
                                height_nfo, 'er_site_name', site, 'T')
                            if len(site_height) > 0:
                                # add in height if available
                                PmagSampRec["sample_height"] = site_height[0]['site_height']
                            PmagSampRec['sample_tilt_correction'] = coord
                            PmagSampRec['sample_comp_name'] = pmag.get_list(
                                CoordDir, 'specimen_comp_name')  # get components used
                            PmagSampRec['er_specimen_names'] = pmag.get_list(
                                CoordDir, 'er_specimen_name')  # get specimne names averaged
                            PmagSampRec['magic_method_codes'] = pmag.get_list(
                                CoordDir, 'magic_method_codes')  # assemble method codes
                            if nocrit != 1:  # apply selection criteria
                                kill = pmag.grade(
                                    PmagSampRec, accept, 'sample_dir')
                                if len(kill) == 0:  # passes the mustard
                                    SampDirs.append(PmagSampRec)
                                    if vgps == 1:
                                        PmagResRec = pmag.getsampVGP(
                                            PmagSampRec, SiteNFO)
                                        if PmagResRec != "":
                                            PmagResults.append(PmagResRec)
                            else:  # take everything
                                SampDirs.append(PmagSampRec)
                                if vgps == 1:
                                    PmagResRec = pmag.getsampVGP(
                                        PmagSampRec, SiteNFO)
                                    if PmagResRec != "":
                                        PmagResults.append(PmagResRec)
                            PmagSamps.append(PmagSampRec)
        if avg_intensities_by_sample:  # average by sample if desired
            # get all the intensity data for this sample
            SampI = pmag.get_dictitem(SpecInts, 'er_sample_name', samp, 'T')
            if len(SampI) > 0:  # there are some
                # get average intensity stuff
                PmagSampRec = pmag.average_int(SampI, 'specimen', 'sample')
                # decorate sample record
                PmagSampRec["sample_description"] = "sample intensity"
                PmagSampRec["sample_direction_type"] = ""
                PmagSampRec['er_site_name'] = SampI[0]["er_site_name"]
                PmagSampRec['er_sample_name'] = samp
                PmagSampRec['er_location_name'] = SampI[0]["er_location_name"]
                PmagSampRec["er_citation_names"] = "This study"
                PmagSampRec["er_analyst_mail_names"] = user
                if agefile != "":
                    PmagSampRec = pmag.get_age(
                        PmagSampRec, "er_site_name", "sample_inferred_", AgeNFO, DefaultAge)
                site_height = pmag.get_dictitem(
                    height_nfo, 'er_site_name', PmagSampRec['er_site_name'], 'T')
                if len(site_height) > 0:
                    # add in height if available
                    PmagSampRec["sample_height"] = site_height[0]['site_height']
                PmagSampRec['er_specimen_names'] = pmag.get_list(
                    SampI, 'er_specimen_name')
                PmagSampRec['magic_method_codes'] = pmag.get_list(
                    SampI, 'magic_method_codes')
                if nocrit != 1:  # apply criteria!
                    kill = pmag.grade(PmagSampRec, accept, 'sample_int')
                    if len(kill) == 0:
                        PmagSampRec['pmag_criteria_codes'] = "ACCEPT"
                        SampInts.append(PmagSampRec)
                        PmagSamps.append(PmagSampRec)
                    else:
                        PmagSampRec = {}  # sample rejected
                else:  # no criteria
                    SampInts.append(PmagSampRec)
                    PmagSamps.append(PmagSampRec)
                    PmagSampRec['pmag_criteria_codes'] = ""
                if vgps == 1 and get_model_lat != 0 and PmagSampRec != {}:
                    if get_model_lat == 1:  # use sample latitude
                        PmagResRec = pmag.getsampVDM(PmagSampRec, SampNFO)
                        # get rid of the model lat key
                        del(PmagResRec['model_lat'])
                    elif get_model_lat == 2:  # use model latitude
                        PmagResRec = pmag.getsampVDM(PmagSampRec, ModelLats)
                        if PmagResRec != {}:
                            PmagResRec['magic_method_codes'] = PmagResRec['magic_method_codes'] + ":IE-MLAT"
                    if PmagResRec != {}:
                        PmagResRec['er_specimen_names'] = PmagSampRec['er_specimen_names']
                        PmagResRec['er_sample_names'] = PmagSampRec['er_sample_name']
                        PmagResRec['pmag_criteria_codes'] = 'ACCEPT'
                        PmagResRec['average_int_sigma_perc'] = PmagSampRec['sample_int_sigma_perc']
                        PmagResRec['average_int_sigma'] = PmagSampRec['sample_int_sigma']
                        PmagResRec['average_int_n'] = PmagSampRec['sample_int_n']
                        PmagResRec['vadm_n'] = PmagSampRec['sample_int_n']
                        PmagResRec['data_type'] = 'i'
                        PmagResults.append(PmagResRec)
    if len(PmagSamps) > 0:
        # fill in missing keys from different types of records
        TmpSamps, keylist = pmag.fillkeys(PmagSamps)
        # save in sample output file
        pmag.magic_write(sampout, TmpSamps, 'pmag_samples')
        print(' sample averages written to ', sampout)

#
# create site averages from specimens or samples as specified
#
    for site in sites:
        for coord in coords:
            if not avg_directions_by_sample:
                key, dirlist = 'specimen', SpecDirs  # if specimen averages at site level desired
            if avg_directions_by_sample:
                key, dirlist = 'sample', SampDirs  # if sample averages at site level desired
            # get all the sites with  directions
            tmp = pmag.get_dictitem(dirlist, 'er_site_name', site, 'T')
            # use only the last coordinate if avg_all_components==False
            tmp1 = pmag.get_dictitem(tmp, key + '_tilt_correction', coord, 'T')
            # fish out site information (lat/lon, etc.)
            sd = pmag.get_dictitem(SiteNFO, 'er_site_name', site, 'T')
            if len(sd) > 0:
                sitedat = sd[0]
                if not avg_all_components:  # do component wise averaging
                    for comp in Comps:
                        # get all components comp
                        siteD = pmag.get_dictitem(
                            tmp1, key + '_comp_name', comp, 'T')
                        # remove bad data from means
                        quality_siteD = []
                        # remove any records for which specimen_flag or sample_flag are 'b'
                        # assume 'g' if flag is not provided
                        for rec in siteD:
                            spec_quality = rec.get('specimen_flag', 'g')
                            samp_quality = rec.get('sample_flag', 'g')
                            if (spec_quality == 'g') and (samp_quality == 'g'):
                                quality_siteD.append(rec)
                        siteD = quality_siteD
                        if len(siteD) > 0:  # there are some for this site and component name
                            # get an average for this site
                            PmagSiteRec = pmag.lnpbykey(siteD, 'site', key)
                            # decorate the site record
                            PmagSiteRec['site_comp_name'] = comp
                            PmagSiteRec["er_location_name"] = siteD[0]['er_location_name']
                            PmagSiteRec["er_site_name"] = siteD[0]['er_site_name']
                            PmagSiteRec['site_tilt_correction'] = coord
                            PmagSiteRec['site_comp_name'] = pmag.get_list(
                                siteD, key + '_comp_name')
                            if avg_directions_by_sample:
                                PmagSiteRec['er_sample_names'] = pmag.get_list(
                                    siteD, 'er_sample_name')
                            else:
                                PmagSiteRec['er_specimen_names'] = pmag.get_list(
                                    siteD, 'er_specimen_name')
    # determine the demagnetization code (DC3,4 or 5) for this site
                            AFnum = len(pmag.get_dictitem(
                                siteD, 'magic_method_codes', 'LP-DIR-AF', 'has'))
                            Tnum = len(pmag.get_dictitem(
                                siteD, 'magic_method_codes', 'LP-DIR-T', 'has'))
                            DC = 3
                            if AFnum > 0:
                                DC += 1
                            if Tnum > 0:
                                DC += 1
                            PmagSiteRec['magic_method_codes'] = pmag.get_list(
                                siteD, 'magic_method_codes') + ':' + 'LP-DC' + str(DC)
                            PmagSiteRec['magic_method_codes'].strip(":")
                            if plotsites:
                                print(PmagSiteRec['er_site_name'])
                                # plot and list the data
                                pmagplotlib.plot_site(
                                    EQ['eqarea'], PmagSiteRec, siteD, key)
                                pmagplotlib.draw_figs(EQ)
                            PmagSites.append(PmagSiteRec)
                else:  # last component only
                    # get the last orientation system specified
                    siteD = tmp1[:]
                    if len(siteD) > 0:  # there are some
                        # get the average for this site
                        PmagSiteRec = pmag.lnpbykey(siteD, 'site', key)
                        # decorate the record
                        PmagSiteRec["er_location_name"] = siteD[0]['er_location_name']
                        PmagSiteRec["er_site_name"] = siteD[0]['er_site_name']
                        PmagSiteRec['site_comp_name'] = comp
                        PmagSiteRec['site_tilt_correction'] = coord
                        PmagSiteRec['site_comp_name'] = pmag.get_list(
                            siteD, key + '_comp_name')
                        PmagSiteRec['er_specimen_names'] = pmag.get_list(
                            siteD, 'er_specimen_name')
                        PmagSiteRec['er_sample_names'] = pmag.get_list(
                            siteD, 'er_sample_name')
                        AFnum = len(pmag.get_dictitem(
                            siteD, 'magic_method_codes', 'LP-DIR-AF', 'has'))
                        Tnum = len(pmag.get_dictitem(
                            siteD, 'magic_method_codes', 'LP-DIR-T', 'has'))
                        DC = 3
                        if AFnum > 0:
                            DC += 1
                        if Tnum > 0:
                            DC += 1
                        PmagSiteRec['magic_method_codes'] = pmag.get_list(
                            siteD, 'magic_method_codes') + ':' + 'LP-DC' + str(DC)
                        PmagSiteRec['magic_method_codes'].strip(":")
                        if not avg_directions_by_sample:
                            PmagSiteRec['site_comp_name'] = pmag.get_list(
                                siteD, key + '_comp_name')
                        if plotsites:
                            pmagplotlib.plot_site(
                                EQ['eqarea'], PmagSiteRec, siteD, key)
                            pmagplotlib.draw_figs(EQ)
                        PmagSites.append(PmagSiteRec)
            else:
                print('site information not found in er_sites for site, ',
                      site, ' site will be skipped')
    for PmagSiteRec in PmagSites:  # now decorate each dictionary some more, and calculate VGPs etc. for results table
        PmagSiteRec["er_citation_names"] = "This study"
        PmagSiteRec["er_analyst_mail_names"] = user
        PmagSiteRec['magic_software_packages'] = version_num
        if agefile != "":
            PmagSiteRec = pmag.get_age(
                PmagSiteRec, "er_site_name", "site_inferred_", AgeNFO, DefaultAge)
        PmagSiteRec['pmag_criteria_codes'] = 'ACCEPT'
        if 'site_n_lines' in list(PmagSiteRec.keys()) and 'site_n_planes' in list(PmagSiteRec.keys()) and PmagSiteRec['site_n_lines'] != "" and PmagSiteRec['site_n_planes'] != "":
            if int(PmagSiteRec["site_n_planes"]) > 0:
                PmagSiteRec["magic_method_codes"] = PmagSiteRec['magic_method_codes'] + ":DE-FM-LP"
            elif int(PmagSiteRec["site_n_lines"]) > 2:
                PmagSiteRec["magic_method_codes"] = PmagSiteRec['magic_method_codes'] + ":DE-FM"
            kill = pmag.grade(PmagSiteRec, accept, 'site_dir')
            if len(kill) == 0:
                PmagResRec = {}  # set up dictionary for the pmag_results table entry
                PmagResRec['data_type'] = 'i'  # decorate it a bit
                PmagResRec['magic_software_packages'] = version_num
                PmagSiteRec['site_description'] = 'Site direction included in results table'
                PmagResRec['pmag_criteria_codes'] = 'ACCEPT'
                dec = float(PmagSiteRec["site_dec"])
                inc = float(PmagSiteRec["site_inc"])
                if 'site_alpha95' in list(PmagSiteRec.keys()) and PmagSiteRec['site_alpha95'] != "":
                    a95 = float(PmagSiteRec["site_alpha95"])
                else:
                    a95 = 180.
                sitedat = pmag.get_dictitem(SiteNFO, 'er_site_name', PmagSiteRec['er_site_name'], 'T')[
                    0]  # fish out site information (lat/lon, etc.)
                lat = float(sitedat['site_lat'])
                lon = float(sitedat['site_lon'])
                plon, plat, dp, dm = pmag.dia_vgp(
                    dec, inc, a95, lat, lon)  # get the VGP for this site
                if PmagSiteRec['site_tilt_correction'] == '-1':
                    C = ' (spec coord) '
                if PmagSiteRec['site_tilt_correction'] == '0':
                    C = ' (geog. coord) '
                if PmagSiteRec['site_tilt_correction'] == '100':
                    C = ' (strat. coord) '
                PmagResRec["pmag_result_name"] = "VGP Site: " + \
                    PmagSiteRec["er_site_name"]  # decorate some more
                PmagResRec["result_description"] = "Site VGP, coord system = " + \
                    str(coord) + ' component: ' + comp
                PmagResRec['er_site_names'] = PmagSiteRec['er_site_name']
                PmagResRec['pmag_criteria_codes'] = 'ACCEPT'
                PmagResRec['er_citation_names'] = 'This study'
                PmagResRec['er_analyst_mail_names'] = user
                PmagResRec["er_location_names"] = PmagSiteRec["er_location_name"]
                if avg_directions_by_sample:
                    PmagResRec["er_sample_names"] = PmagSiteRec["er_sample_names"]
                else:
                    PmagResRec["er_specimen_names"] = PmagSiteRec["er_specimen_names"]
                PmagResRec["tilt_correction"] = PmagSiteRec['site_tilt_correction']
                PmagResRec["pole_comp_name"] = PmagSiteRec['site_comp_name']
                PmagResRec["average_dec"] = PmagSiteRec["site_dec"]
                PmagResRec["average_inc"] = PmagSiteRec["site_inc"]
                PmagResRec["average_alpha95"] = PmagSiteRec["site_alpha95"]
                PmagResRec["average_n"] = PmagSiteRec["site_n"]
                PmagResRec["average_n_lines"] = PmagSiteRec["site_n_lines"]
                PmagResRec["average_n_planes"] = PmagSiteRec["site_n_planes"]
                PmagResRec["vgp_n"] = PmagSiteRec["site_n"]
                PmagResRec["average_k"] = PmagSiteRec["site_k"]
                PmagResRec["average_r"] = PmagSiteRec["site_r"]
                PmagResRec["average_lat"] = '%10.4f ' % (lat)
                PmagResRec["average_lon"] = '%10.4f ' % (lon)
                if agefile != "":
                    PmagResRec = pmag.get_age(
                        PmagResRec, "er_site_names", "average_", AgeNFO, DefaultAge)
                site_height = pmag.get_dictitem(
                    height_nfo, 'er_site_name', site, 'T')
                if len(site_height) > 0:
                    PmagResRec["average_height"] = site_height[0]['site_height']
                PmagResRec["vgp_lat"] = '%7.1f ' % (plat)
                PmagResRec["vgp_lon"] = '%7.1f ' % (plon)
                PmagResRec["vgp_dp"] = '%7.1f ' % (dp)
                PmagResRec["vgp_dm"] = '%7.1f ' % (dm)
                PmagResRec["magic_method_codes"] = PmagSiteRec["magic_method_codes"]
                if '0' in PmagSiteRec['site_tilt_correction'] and "DA-DIR-GEO" not in PmagSiteRec['magic_method_codes']:
                    PmagSiteRec['magic_method_codes'] = PmagSiteRec['magic_method_codes'] + ":DA-DIR-GEO"
                if '100' in PmagSiteRec['site_tilt_correction'] and "DA-DIR-TILT" not in PmagSiteRec['magic_method_codes']:
                    PmagSiteRec['magic_method_codes'] = PmagSiteRec['magic_method_codes'] + ":DA-DIR-TILT"
                PmagSiteRec['site_polarity'] = ""
                if avg_by_polarity:  # assign polarity based on angle of pole lat to spin axis - may want to re-think this sometime
                    angle = pmag.angle([0, 0], [0, (90 - plat)])
                    if angle <= 55.:
                        PmagSiteRec["site_polarity"] = 'n'
                    if angle > 55. and angle < 125.:
                        PmagSiteRec["site_polarity"] = 't'
                    if angle >= 125.:
                        PmagSiteRec["site_polarity"] = 'r'
                PmagResults.append(PmagResRec)
    if avg_by_polarity:
        # find the tilt corrected data
        crecs = pmag.get_dictitem(
            PmagSites, 'site_tilt_correction', '100', 'T')
        if len(crecs) < 2:
            # if there aren't any, find the geographic corrected data
            crecs = pmag.get_dictitem(
                PmagSites, 'site_tilt_correction', '0', 'T')
        if len(crecs) > 2:  # if there are some,
            comp = pmag.get_list(crecs, 'site_comp_name').split(':')[
                0]  # find the first component
            # fish out all of the first component
            crecs = pmag.get_dictitem(crecs, 'site_comp_name', comp, 'T')
            precs = []
            for rec in crecs:
                precs.append({'dec': rec['site_dec'], 'inc': rec['site_inc'],
                              'name': rec['er_site_name'], 'loc': rec['er_location_name']})
            # calculate average by polarity
            polpars = pmag.fisher_by_pol(precs)
            # hunt through all the modes (normal=A, reverse=B, all=ALL)
            for mode in list(polpars.keys()):
                PolRes = {}
                PolRes['er_citation_names'] = 'This study'
                PolRes["pmag_result_name"] = "Polarity Average: Polarity " + mode
                PolRes["data_type"] = "a"
                PolRes["average_dec"] = '%7.1f' % (polpars[mode]['dec'])
                PolRes["average_inc"] = '%7.1f' % (polpars[mode]['inc'])
                PolRes["average_n"] = '%i' % (polpars[mode]['n'])
                PolRes["average_r"] = '%5.4f' % (polpars[mode]['r'])
                PolRes["average_k"] = '%6.0f' % (polpars[mode]['k'])
                PolRes["average_alpha95"] = '%7.1f' % (
                    polpars[mode]['alpha95'])
                PolRes['er_site_names'] = polpars[mode]['sites']
                PolRes['er_location_names'] = polpars[mode]['locs']
                PolRes['magic_software_packages'] = version_num
                PmagResults.append(PolRes)

    if not skip_intensities and nositeints != 1:
        for site in sites:  # now do intensities for each site
            if plotsites:
                print(site)
            if not avg_intensities_by_sample:
                key, intlist = 'specimen', SpecInts  # if using specimen level data
            if avg_intensities_by_sample:
                key, intlist = 'sample', PmagSamps  # if using sample level data
            # get all the intensities  for this site
            Ints = pmag.get_dictitem(intlist, 'er_site_name', site, 'T')
            if len(Ints) > 0:  # there are some
                # get average intensity stuff for site table
                PmagSiteRec = pmag.average_int(Ints, key, 'site')
                # get average intensity stuff for results table
                PmagResRec = pmag.average_int(Ints, key, 'average')
                if plotsites:  # if site by site examination requested - print this site out to the screen
                    for rec in Ints:
                        print(rec['er_' + key + '_name'], ' %7.1f' %
                              (1e6 * float(rec[key + '_int'])))
                    if len(Ints) > 1:
                        print('Average: ', '%7.1f' % (
                            1e6 * float(PmagResRec['average_int'])), 'N: ', len(Ints))
                        print('Sigma: ', '%7.1f' % (
                            1e6 * float(PmagResRec['average_int_sigma'])), 'Sigma %: ', PmagResRec['average_int_sigma_perc'])
                    input('Press any key to continue\n')
                er_location_name = Ints[0]["er_location_name"]
                # decorate the records
                PmagSiteRec["er_location_name"] = er_location_name
                PmagSiteRec["er_citation_names"] = "This study"
                PmagResRec["er_location_names"] = er_location_name
                PmagResRec["er_citation_names"] = "This study"
                PmagSiteRec["er_analyst_mail_names"] = user
                PmagResRec["er_analyst_mail_names"] = user
                PmagResRec["data_type"] = 'i'
                if not avg_intensities_by_sample:
                    PmagSiteRec['er_specimen_names'] = pmag.get_list(
                        Ints, 'er_specimen_name')  # list of all specimens used
                    PmagResRec['er_specimen_names'] = pmag.get_list(
                        Ints, 'er_specimen_name')
                PmagSiteRec['er_sample_names'] = pmag.get_list(
                    Ints, 'er_sample_name')  # list of all samples used
                PmagResRec['er_sample_names'] = pmag.get_list(
                    Ints, 'er_sample_name')
                PmagSiteRec['er_site_name'] = site
                PmagResRec['er_site_names'] = site
                PmagSiteRec['magic_method_codes'] = pmag.get_list(
                    Ints, 'magic_method_codes')
                PmagResRec['magic_method_codes'] = pmag.get_list(
                    Ints, 'magic_method_codes')
                kill = pmag.grade(PmagSiteRec, accept, 'site_int')
                if nocrit == 1 or len(kill) == 0:
                    b, sig = float(PmagResRec['average_int']), ""
                    if(PmagResRec['average_int_sigma']) != "":
                        sig = float(PmagResRec['average_int_sigma'])
                    # fish out site direction
                    sdir = pmag.get_dictitem(
                        PmagResults, 'er_site_names', site, 'T')
                    # get the VDM for this record using last average
                    # inclination (hope it is the right one!)
                    if len(sdir) > 0 and sdir[-1]['average_inc'] != "":
                        inc = float(sdir[0]['average_inc'])
                        # get magnetic latitude using dipole formula
                        mlat = pmag.magnetic_lat(inc)
                        # get VDM with magnetic latitude
                        PmagResRec["vdm"] = '%8.3e ' % (pmag.b_vdm(b, mlat))
                        PmagResRec["vdm_n"] = PmagResRec['average_int_n']
                        if 'average_int_sigma' in list(PmagResRec.keys()) and PmagResRec['average_int_sigma'] != "":
                            vdm_sig = pmag.b_vdm(
                                float(PmagResRec['average_int_sigma']), mlat)
                            PmagResRec["vdm_sigma"] = '%8.3e ' % (vdm_sig)
                        else:
                            PmagResRec["vdm_sigma"] = ""
                    mlat = ""  # define a model latitude
                    if get_model_lat == 1:  # use present site latitude
                        mlats = pmag.get_dictitem(
                            SiteNFO, 'er_site_name', site, 'T')
                        if len(mlats) > 0:
                            mlat = mlats[0]['site_lat']
                    # use a model latitude from some plate reconstruction model
                    # (or something)
                    elif get_model_lat == 2:
                        mlats = pmag.get_dictitem(
                            ModelLats, 'er_site_name', site, 'T')
                        if len(mlats) > 0:
                            PmagResRec['model_lat'] = mlats[0]['site_model_lat']
                        mlat = PmagResRec['model_lat']
                    if mlat != "":
                        # get the VADM using the desired latitude
                        PmagResRec["vadm"] = '%8.3e ' % (
                            pmag.b_vdm(b, float(mlat)))
                        if sig != "":
                            vdm_sig = pmag.b_vdm(
                                float(PmagResRec['average_int_sigma']), float(mlat))
                            PmagResRec["vadm_sigma"] = '%8.3e ' % (vdm_sig)
                            PmagResRec["vadm_n"] = PmagResRec['average_int_n']
                        else:
                            PmagResRec["vadm_sigma"] = ""
                    # fish out site information (lat/lon, etc.)
                    sitedat = pmag.get_dictitem(
                        SiteNFO, 'er_site_name', PmagSiteRec['er_site_name'], 'T')
                    if len(sitedat) > 0:
                        sitedat = sitedat[0]
                        PmagResRec['average_lat'] = sitedat['site_lat']
                        PmagResRec['average_lon'] = sitedat['site_lon']
                    else:
                        PmagResRec['average_lon'] = 'UNKNOWN'
                        PmagResRec['average_lon'] = 'UNKNOWN'
                    PmagResRec['magic_software_packages'] = version_num
                    PmagResRec["pmag_result_name"] = "V[A]DM: Site " + site
                    PmagResRec["result_description"] = "V[A]DM of site"
                    PmagResRec["pmag_criteria_codes"] = "ACCEPT"
                    if agefile != "":
                        PmagResRec = pmag.get_age(
                            PmagResRec, "er_site_names", "average_", AgeNFO, DefaultAge)
                    site_height = pmag.get_dictitem(
                        height_nfo, 'er_site_name', site, 'T')
                    if len(site_height) > 0:
                        PmagResRec["average_height"] = site_height[0]['site_height']
                    PmagSites.append(PmagSiteRec)
                    PmagResults.append(PmagResRec)
    if len(PmagSites) > 0:
        Tmp, keylist = pmag.fillkeys(PmagSites)
        pmag.magic_write(siteout, Tmp, 'pmag_sites')
        print(' sites written to ', siteout)
    else:
        print("No Site level table")
    if len(PmagResults) > 0:
        TmpRes, keylist = pmag.fillkeys(PmagResults)
        pmag.magic_write(resout, TmpRes, 'pmag_results')
        print(' results written to ', resout)
    else:
        print("No Results level table")