docs: def loadtxt(fname, dtype="float", delimiter="\t", usecols=None, comments="#"):
      r"""
      Load unyt_arrays with unit information from a text file. Each row in the
      text file must have the same number of values.
  
      Parameters
      ----------
      fname : str
          Filename to read.
      dtype : data-type, optional
          Data-type of the resulting array; default: float.
      delimiter : str, optional
          The string used to separate values.  By default, this is any
          whitespace.
      usecols : sequence, optional
          Which columns to read, with 0 being the first.  For example,
          ``usecols = (1,4,5)`` will extract the 2nd, 5th and 6th columns.
          The default, None, results in all columns being read.
      comments : str, optional
          The character used to indicate the start of a comment;
          default: '#'.
  
      Examples
      --------
      >>> temp, velx = loadtxt(
      ...    "sphere.dat", usecols=(1,2), delimiter="\t")  # doctest: +SKIP
      """
      f = open(fname, "r")
      next_one = False
      units = []
      num_cols = -1
      for line in f.readlines():
          words = line.strip().split()
          if len(words) == 0:
              continue
          if line[0] == comments:
              if next_one:
                  units = words[1:]
              if len(words) == 2 and words[1] == "Units":
                  next_one = True
          else:
              # Here we catch the first line of numbers
              col_words = line.strip().split(delimiter)
              for word in col_words:
                  float(word)
              num_cols = len(col_words)
              break
      f.close()
      if len(units) != num_cols:
          units = ["dimensionless"] * num_cols
      arrays = np.loadtxt(
          fname,
          dtype=dtype,
          comments=comments,
          delimiter=delimiter,
          converters=None,
          unpack=True,
          usecols=usecols,
          ndmin=0,
      )
      if len(arrays.shape) < 2:
          arrays = [arrays]
      if usecols is not None:
          units = [units[col] for col in usecols]
      ret = tuple([unyt_array(arr, unit) for arr, unit in zip(arrays, units)])
      if len(ret) == 1:
          return ret[0]
      return ret