docs: def int_pars(x, y, vds, **kwargs):
      """
       calculates York regression and Coe parameters (with Tauxe Fvds)
      """
      # first do linear regression a la York
      # do Data Model 3 way:
      if 'version' in list(kwargs.keys()) and kwargs['version'] == 3:
          n_key = 'int_n_measurements'
          b_key = 'int_b'
          sigma_key = 'int_b_sigma'
          f_key = 'int_f'
          fvds_key = 'int_fvds'
          g_key = 'int_g'
          q_key = 'int_q'
          b_beta_key = 'int_b_beta'
  
      else:  # version 2
          n_key = 'specimen_int_n'
          b_key = 'specimen_b'
          sigma_key = 'specimen_b_sigma'
          f_key = 'specimen_f'
          fvds_key = 'specimen_fvds'
          g_key = 'specimen_g'
          q_key = 'specimen_q'
          b_beta_key = 'specimen_b_beta'
  
      xx, yer, xer, xyer, yy, xsum, ysum, xy = 0., 0., 0., 0., 0., 0., 0., 0.
      xprime, yprime = [], []
      pars = {}
      pars[n_key] = len(x)
      n = float(len(x))
      if n <= 2:
          print("shouldn't be here at all!")
          return pars, 1
      for i in range(len(x)):
          xx += x[i]**2.
          yy += y[i]**2.
          xy += x[i] * y[i]
          xsum += x[i]
          ysum += y[i]
      xsig = np.sqrt(old_div((xx - (old_div(xsum**2., n))), (n - 1.)))
      ysig = np.sqrt(old_div((yy - (old_div(ysum**2., n))), (n - 1.)))
      sum = 0
      for i in range(int(n)):
          yer += (y[i] - old_div(ysum, n))**2.
          xer += (x[i] - old_div(xsum, n))**2.
          xyer += (y[i] - old_div(ysum, n)) * (x[i] - old_div(xsum, n))
      slop = -np.sqrt(old_div(yer, xer))
      pars[b_key] = slop
      s1 = 2. * yer - 2. * slop * xyer
      s2 = (n - 2.) * xer
      sigma = np.sqrt(old_div(s1, s2))
      pars[sigma_key] = sigma
      s = old_div((xy - (xsum * ysum / n)), (xx - old_div((xsum**2.), n)))
      r = old_div((s * xsig), ysig)
      pars["specimen_rsc"] = r**2.
      ytot = abs(old_div(ysum, n) - slop * xsum / n)
      for i in range(int(n)):
          xprime.append(old_div((slop * x[i] + y[i] - ytot), (2. * slop)))
          yprime.append((old_div((slop * x[i] + y[i] - ytot), 2.)) + ytot)
      sumdy, dy = 0, []
      dyt = abs(yprime[0] - yprime[int(n) - 1])
      for i in range((int(n) - 1)):
          dy.append(abs(yprime[i + 1] - yprime[i]))
          sumdy += dy[i]**2.
      f = old_div(dyt, ytot)
      pars[f_key] = f
      pars["specimen_ytot"] = ytot
      ff = old_div(dyt, vds)
      pars[fvds_key] = ff
      ddy = (old_div(1., dyt)) * sumdy
      g = 1. - old_div(ddy, dyt)
      pars[g_key] = g
      q = abs(slop) * f * g / sigma
      pars[q_key] = q
      pars[b_beta_key] = old_div(-sigma, slop)
      return pars, 0