docs: def York_Regression(x_segment, y_segment, x_mean, y_mean, n, lab_dc_field, steps_Arai):
      """
      input: x_segment, y_segment, x_mean, y_mean, n, lab_dc_field, steps_Arai
      output: x_err,  y_err, x_tag, y_tag, b, b_sigma, specimen_b_beta, y_intercept,
              x_intercept, x_prime, y_prime, delta_x_prime, delta_y_prime, f_Coe,
              g_Coe, g_lim, specimen_q, specimen_w, count_IZ, count_ZI, B_lab, B_anc,
              B_anc_sigma, specimen_int
      """
      x_err = x_segment - x_mean
      y_err = y_segment - y_mean
      york_b = -1* numpy.sqrt(old_div(sum(y_err**2), sum(x_err**2)) )  # averaged slope
  
      b = numpy.sign(sum(x_err * y_err)) * numpy.std(y_segment, ddof=1)/numpy.std(x_segment, ddof=1) # ddof is degrees of freedom
      if b == 0:
          york_b = 1e-10
      else:
          york_b = b
  
      york_sigma= numpy.sqrt( old_div((2 * sum(y_err**2) - 2*york_b* sum(x_err*y_err)), ( (n-2) * sum(x_err**2) )) )
      if york_sigma == 0: # prevent divide by zero
          york_sigma = 1e-10
      beta_Coe = abs(old_div(york_sigma,york_b))
      # y_T is the intercept of the extrepolated line
      # through the center of mass (see figure 7 in Coe (1978))
      y_T = y_mean - (york_b* x_mean)
      x_T = old_div((-1 * y_T), york_b)  # x intercept
      # # calculate the extrarpolated data points for f and fvds
      x_tag = old_div((y_segment - y_T ), york_b) # returns array of y points minus the y intercept, divided by slope
      y_tag = york_b*x_segment + y_T
  
      # intersect of the dashed square and the horizontal dahed line  next to delta-y-5 in figure 7, Coe (1978)
      x_prime = old_div((x_segment+x_tag), 2)
      y_prime = old_div((y_segment+y_tag), 2)
  
      delta_x_prime = abs(max(x_prime) - min(x_prime)) # TRM length of best fit line
      delta_y_prime = abs(max(y_prime) - min(y_prime)) # NRM length of best fit line
  
      f_Coe = old_div(delta_y_prime, abs(y_T))
      if delta_y_prime:
          g_Coe =  1 - (old_div(sum((y_prime[:-1]-y_prime[1:])**2), delta_y_prime ** 2))  # gap factor
      else:
          g_Coe = float('nan')
      g_lim = old_div((float(n) - 2), (float(n) - 1))
      q_Coe = abs(york_b)*f_Coe*g_Coe/york_sigma
      w_Coe = old_div(q_Coe, numpy.sqrt(n - 2))
      count_IZ = steps_Arai.count('IZ')
      count_ZI = steps_Arai.count('ZI')
      B_lab = lab_dc_field * 1e6
      B_anc = abs(york_b) * B_lab # in microtesla
      B_anc_sigma = york_sigma * B_lab
      specimen_int = -1* lab_dc_field * york_b # in tesla
      specimen_int_sigma = york_sigma * lab_dc_field
      return {'x_err': x_err, 'y_err': y_err, 'x_tag': x_tag, 'y_tag': y_tag,
              'specimen_b': york_b, 'specimen_b_sigma': york_sigma, 'specimen_b_beta': beta_Coe,
              'y_int': y_T, 'x_int': x_T, 'x_prime': x_prime, 'y_prime': y_prime,
              'delta_x_prime': delta_x_prime, 'delta_y_prime': delta_y_prime, 'specimen_f': f_Coe,
              'specimen_g': g_Coe, 'specimen_g_lim': g_lim, 'specimen_q': q_Coe, 'specimen_w': w_Coe,
              'count_IZ': count_IZ, 'count_ZI': count_ZI, 'B_lab': B_lab, 'B_anc': B_anc,
              'B_anc_sigma': B_anc_sigma, 'specimen_int': specimen_int, 'specimen_int_sigma': specimen_int_sigma}