docs: def fields(fieldfile, only_header=False, only_istep=False):
      """Extract fields data.
  
      Args:
          fieldfile (:class:`pathlib.Path`): path of the binary field file.
          only_header (bool): when True (and :data:`only_istep` is False), only
              :data:`header` is returned.
          only_istep (bool): when True, only :data:`istep` is returned.
  
      Returns:
          depends on flags.: :obj:`int`: istep
              If :data:`only_istep` is True, this function returns the time step
              at which the binary file was written.
          :obj:`dict`: header
              Else, if :data:`only_header` is True, this function returns a dict
              containing the header informations of the binary file.
          :class:`numpy.array`: fields
              Else, this function returns the tuple :data:`(header, fields)`.
              :data:`fields` is an array of scalar fields indexed by variable,
              x-direction, y-direction, z-direction, block.
      """
      # something to skip header?
      if not fieldfile.is_file():
          return None
      header = {}
      with fieldfile.open('rb') as fid:
          readbin = partial(_readbin, fid)
          magic = readbin()
          if magic > 8000:  # 64 bits
              magic -= 8000
              readbin()  # need to read 4 more bytes
              readbin = partial(readbin, file64=True)
  
          # check nb components
          nval = 1
          if magic > 400:
              nval = 4
          elif magic > 300:
              nval = 3
  
          magic %= 100
  
          # extra ghost point in horizontal direction
          header['xyp'] = int(magic >= 9 and nval == 4)
  
          # total number of values in relevant space basis
          # (e1, e2, e3) = (theta, phi, radius) in spherical geometry
          #              = (x, y, z)            in cartesian geometry
          header['nts'] = readbin(nwords=3)
  
          # number of blocks, 2 for yinyang or cubed sphere
          header['ntb'] = readbin() if magic >= 7 else 1
  
          # aspect ratio
          header['aspect'] = readbin('f', 2)
  
          # number of parallel subdomains
          header['ncs'] = readbin(nwords=3)  # (e1, e2, e3) space
          header['ncb'] = readbin() if magic >= 8 else 1  # blocks
  
          # r - coordinates
          # rgeom[0:self.nrtot+1, 0] are edge radial position
          # rgeom[0:self.nrtot, 1] are cell-center radial position
          if magic >= 2:
              header['rgeom'] = readbin('f', header['nts'][2] * 2 + 1)
          else:
              header['rgeom'] = np.array(range(0, header['nts'][2] * 2 + 1))\
                  * 0.5 / header['nts'][2]
          header['rgeom'] = np.resize(header['rgeom'], (header['nts'][2] + 1, 2))
  
          header['rcmb'] = readbin('f') if magic >= 7 else None
  
          header['ti_step'] = readbin() if magic >= 3 else 0
          if only_istep:
              return header['ti_step']
          header['ti_ad'] = readbin('f') if magic >= 3 else 0
          header['erupta_total'] = readbin('f') if magic >= 5 else 0
          header['bot_temp'] = readbin('f') if magic >= 6 else 1
  
          if magic >= 4:
              header['e1_coord'] = readbin('f', header['nts'][0])
              header['e2_coord'] = readbin('f', header['nts'][1])
              header['e3_coord'] = readbin('f', header['nts'][2])
          else:
              # could construct them from other info
              raise ParsingError(fieldfile,
                                 'magic >= 4 expected to get grid geometry')
  
          if only_header:
              return header
  
          # READ FIELDS
          # number of points in (e1, e2, e3) directions PER CPU
          npc = header['nts'] // header['ncs']
          # number of blocks per cpu
          nbk = header['ntb'] // header['ncb']
          # number of values per 'read' block
          npi = (npc[0] + header['xyp']) * (npc[1] + header['xyp']) * npc[2] * \
              nbk * nval
  
          header['scalefac'] = readbin('f') if nval > 1 else 1
  
          flds = np.zeros((nval,
                           header['nts'][0] + header['xyp'],
                           header['nts'][1] + header['xyp'],
                           header['nts'][2],
                           header['ntb']))
  
          # loop over parallel subdomains
          for icpu in product(range(header['ncb']),
                              range(header['ncs'][2]),
                              range(header['ncs'][1]),
                              range(header['ncs'][0])):
              # read the data for one CPU
              data_cpu = readbin('f', npi) * header['scalefac']
  
              # icpu is (icpu block, icpu z, icpu y, icpu x)
              # data from file is transposed to obtained a field
              # array indexed with (x, y, z, block), as in StagYY
              flds[:,
                   icpu[3] * npc[0]:(icpu[3] + 1) * npc[0] + header['xyp'],  # x
                   icpu[2] * npc[1]:(icpu[2] + 1) * npc[1] + header['xyp'],  # y
                   icpu[1] * npc[2]:(icpu[1] + 1) * npc[2],  # z
                   icpu[0] * nbk:(icpu[0] + 1) * nbk  # block
                   ] = np.transpose(data_cpu.reshape(
                       (nbk, npc[2], npc[1] + header['xyp'],
                        npc[0] + header['xyp'], nval)))
      return header, flds