docs: 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")