docs: def domean(data, start, end, calculation_type):
      """
      Gets average direction using Fisher or principal component analysis (line
      or plane) methods
  
      Parameters
      ----------
      data : nest list of data: [[treatment,dec,inc,int,quality],...]
      start : step being used as start of fit (often temperature minimum)
      end : step being used as end of fit (often temperature maximum)
      calculation_type : string describing type of calculation to be made
      'DE-BFL' (line), 'DE-BFL-A' (line-anchored), 'DE-BFL-O' (line-with-origin),
      'DE-BFP' (plane), 'DE-FM' (Fisher mean)
  
      Returns
      -------
      mpars : dictionary with the keys "specimen_n","measurement_step_min",
      "measurement_step_max","specimen_mad","specimen_dec","specimen_inc"
      """
      mpars = {}
      datablock = []
      start0, end0 = start, end
      # indata = [rec.append('g') if len(rec)<6 else rec for rec in indata] #
      # this statement doesn't work!
      indata = []
      for rec in data:
          if len(rec) < 6:
              rec.append('g')
          indata.append(rec)
      if indata[start0][5] == 'b':
          print("Can't select 'bad' point as start for PCA")
      flags = [x[5] for x in indata]
      bad_before_start = flags[:start0].count('b')
      bad_in_mean = flags[start0:end0 + 1].count('b')
      start = start0 - bad_before_start
      end = end0 - bad_before_start - bad_in_mean
      datablock = [x for x in indata if x[5] == 'g']
      if indata[start0] != datablock[start]:
          print('problem removing bad data in pmag.domean start of datablock shifted:\norigional: %d\nafter removal: %d' % (
              start0, indata.index(datablock[start])))
      if indata[end0] != datablock[end]:
          print('problem removing bad data in pmag.domean end of datablock shifted:\norigional: %d\nafter removal: %d' % (
              end0, indata.index(datablock[end])))
      mpars["calculation_type"] = calculation_type
      rad = old_div(np.pi, 180.)
      if end > len(datablock) - 1 or end < start:
          end = len(datablock) - 1
      control, data, X, Nrec = [], [], [], float(end - start + 1)
      cm = [0., 0., 0.]
  #
  #  get cartesian coordinates
  #
      fdata = []
      for k in range(start, end + 1):
          if calculation_type == 'DE-BFL' or calculation_type == 'DE-BFL-A' or calculation_type == 'DE-BFL-O':  # best-fit line
              data = [datablock[k][1], datablock[k][2], datablock[k][3]]
          else:
              data = [datablock[k][1], datablock[k][2], 1.0]  # unit weight
          fdata.append(data)
          cart = dir2cart(data)
          X.append(cart)
      if calculation_type == 'DE-BFL-O':  # include origin as point
          X.append([0., 0., 0.])
          # pass
      if calculation_type == 'DE-FM':  # for fisher means
          fpars = fisher_mean(fdata)
          mpars["specimen_direction_type"] = 'l'
          mpars["specimen_dec"] = fpars["dec"]
          mpars["specimen_inc"] = fpars["inc"]
          mpars["specimen_alpha95"] = fpars["alpha95"]
          mpars["specimen_n"] = fpars["n"]
          mpars["specimen_r"] = fpars["r"]
          mpars["measurement_step_min"] = indata[start0][0]
          mpars["measurement_step_max"] = indata[end0][0]
          mpars["center_of_mass"] = cm
          mpars["specimen_dang"] = -1
          return mpars
  #
  #   get center of mass for principal components (DE-BFL or DE-BFP)
  #
      for cart in X:
          for l in range(3):
              cm[l] += old_div(cart[l], Nrec)
      mpars["center_of_mass"] = cm
  
  #
  #   transform to center of mass (if best-fit line)
  #
      if calculation_type != 'DE-BFP':
          mpars["specimen_direction_type"] = 'l'
      if calculation_type == 'DE-BFL' or calculation_type == 'DE-BFL-O':  # not for planes or anchored lines
          for k in range(len(X)):
              for l in range(3):
                  X[k][l] = X[k][l] - cm[l]
      else:
          mpars["specimen_direction_type"] = 'p'
  
  #
  #   put in T matrix
  #
      T = np.array(Tmatrix(X))
  #
  #   get sorted evals/evects
  #
      t, V = tauV(T)
      if t == []:
          mpars["specimen_direction_type"] = "Error"
          print("Error in calculation")
          return mpars
      v1, v3 = V[0], V[2]
      if t[2] < 0:
          t[2] = 0  # make positive
      if calculation_type == 'DE-BFL-A':
          Dir, R = vector_mean(fdata)
          mpars["specimen_direction_type"] = 'l'
          mpars["specimen_dec"] = Dir[0]
          mpars["specimen_inc"] = Dir[1]
          mpars["specimen_n"] = len(fdata)
          mpars["measurement_step_min"] = indata[start0][0]
          mpars["measurement_step_max"] = indata[end0][0]
          mpars["center_of_mass"] = cm
          s1 = np.sqrt(t[0])
          MAD = old_div(np.arctan(old_div(np.sqrt(t[1] + t[2]), s1)), rad)
          if np.iscomplexobj(MAD):
              MAD = MAD.real
          # I think this is how it is done - i never anchor the "PCA" - check
          mpars["specimen_mad"] = MAD
          return mpars
      if calculation_type != 'DE-BFP':
          #
          #   get control vector for principal component direction
          #
          rec = [datablock[start][1], datablock[start][2], datablock[start][3]]
          P1 = dir2cart(rec)
          rec = [datablock[end][1], datablock[end][2], datablock[end][3]]
          P2 = dir2cart(rec)
  #
  #   get right direction along principal component
  ##
          for k in range(3):
              control.append(P1[k] - P2[k])
          # changed by rshaar
          # control is taken as the center of mass
          # control=cm
  
          dot = 0
          for k in range(3):
              dot += v1[k] * control[k]
          if dot < -1:
              dot = -1
          if dot > 1:
              dot = 1
          if np.arccos(dot) > old_div(np.pi, 2.):
              for k in range(3):
                  v1[k] = -v1[k]
  #   get right direction along principal component
  #
          s1 = np.sqrt(t[0])
          Dir = cart2dir(v1)
          MAD = old_div(np.arctan(old_div(np.sqrt(t[1] + t[2]), s1)), rad)
          if np.iscomplexobj(MAD):
              MAD = MAD.real
      if calculation_type == "DE-BFP":
          Dir = cart2dir(v3)
          MAD = old_div(
              np.arctan(np.sqrt(old_div(t[2], t[1]) + old_div(t[2], t[0]))), rad)
          if np.iscomplexobj(MAD):
              MAD = MAD.real
  #
  #   get angle with  center of mass
  #
      CMdir = cart2dir(cm)
      Dirp = [Dir[0], Dir[1], 1.]
      dang = angle(CMdir, Dirp)
      mpars["specimen_dec"] = Dir[0]
      mpars["specimen_inc"] = Dir[1]
      mpars["specimen_mad"] = MAD
      # mpars["specimen_n"]=int(Nrec)
      mpars["specimen_n"] = len(X)
      mpars["specimen_dang"] = dang[0]
      mpars["measurement_step_min"] = indata[start0][0]
      mpars["measurement_step_max"] = indata[end0][0]
      return mpars