docs: def dftphotom(cfg):
      """Run the discrete-Fourier-transform photometry algorithm.
  
      See the module-level documentation and the output of ``casatask dftphotom
      --help`` for help. All of the algorithm configuration is specified in the
      *cfg* argument, which is an instance of :class:`Config`.
  
      """
      tb = util.tools.table()
      ms = util.tools.ms()
      me = util.tools.measures()
  
      # Read stuff in. Even if the weight values don't have their
      # absolute scale set correctly, we can still use them to set the
      # relative weighting of the data points.
      #
      # datacol is (ncorr, nchan, nchunk)
      # flag is (ncorr, nchan, nchunk)
      # weight is (ncorr, nchunk)
      # uvw is (3, nchunk)
      # time is (nchunk)
      # axis_info.corr_axis is (ncorr)
      # axis_info.freq_axis.chan_freq is (nchan, 1) [for now?]
      #
      # Note that we apply msselect() again when reading the data because
      # selectinit() is broken, but the invocation here is good because it
      # affects the results from ms.range() and friends.
  
      if ':' in (cfg.spw or ''):
          warn('it looks like you are attempting to select channels within one or more spws')
          warn('this is NOT IMPLEMENTED; I will average over the whole spw instead')
  
      ms.open(b(cfg.vis))
      totrows = ms.nrow()
      ms_sels = dict((n, cfg.get(n)) for n in util.msselect_keys
                     if cfg.get(n) is not None)
      ms.msselect(b(ms_sels))
  
      rangeinfo = ms.range(b'data_desc_id field_id'.split())
      ddids = rangeinfo['data_desc_id']
      fields = rangeinfo['field_id']
      colnames = [cfg.datacol] + 'flag weight time axis_info'.split()
      rephase = (cfg.rephase is not None)
  
      if fields.size != 1:
          # I feel comfortable making this a fatal error, even if we're
          # not rephasing.
          die('selected data should contain precisely one field; got %d', fields.size)
  
      if rephase:
          fieldid = fields[0]
          tb.open(b(os.path.join(cfg.vis, 'FIELD')))
          phdirinfo = tb.getcell(b'PHASE_DIR', fieldid)
          tb.close()
  
          if phdirinfo.shape[1] != 1:
              die('trying to rephase but target field (#%d) has a '
                  'time-variable phase center, which I can\'t handle', fieldid)
          ra0, dec0 = phdirinfo[:,0] # in radians.
  
          # based on intflib/pwflux.py, which was copied from
          # hex/hex-lib-calcgainerr:
  
          dra = cfg.rephase[0] - ra0
          dec = cfg.rephase[1]
          l = np.sin(dra) * np.cos(dec)
          m = np.sin(dec) * np.cos(dec0) - np.cos(dra) * np.cos(dec) * np.sin(dec0)
          n = np.sin(dec) * np.sin(dec0) + np.cos(dra) * np.cos(dec) * np.cos(dec0)
          n -= 1 # makes the work below easier
          lmn = np.asarray([l, m, n])
          colnames.append('uvw')
  
          # Also need this although 99% of the time `ddid` and `spwid` are the same
          tb.open(b(os.path.join(cfg.vis, 'DATA_DESCRIPTION')))
          ddspws = np.asarray(tb.getcol(b'SPECTRAL_WINDOW_ID'))
          tb.close()
  
      tbins = {}
      colnames = b(colnames)
  
      for ddindex, ddid in enumerate(ddids):
          # Starting in CASA 4.6, selectinit(ddid) stopped actually filtering
          # your data to match the specified DDID! What garbage. Work around
          # with our own filtering.
          ms_sels['taql'] = 'DATA_DESC_ID == %d' % ddid
          ms.msselect(b(ms_sels))
  
          ms.selectinit(ddid)
          if cfg.polarization is not None:
              ms.selectpolarization(b(cfg.polarization.split(',')))
          ms.iterinit(maxrows=4096)
          ms.iterorigin()
  
          while True:
              cols = ms.getdata(items=colnames)
  
              if rephase:
                  # With appropriate spw/DDID selection, `freqs` has shape
                  # (nchan, 1). Convert to m^-1 so we can multiply against UVW
                  # directly.
                  freqs = cols['axis_info']['freq_axis']['chan_freq']
                  assert freqs.shape[1] == 1, 'internal inconsistency, chan_freq??'
                  freqs = freqs[:,0] * util.INVERSE_C_MS
  
              for i in range(cols['time'].size): # all records
                  time = cols['time'][i]
                  # get out of UTC as fast as we can! For some reason
                  # giving 'unit=s' below doesn't do what one might hope it would.
                  # CASA can convert to a variety of timescales; TAI is probably
                  # the safest conversion in terms of being helpful while remaining
                  # close to the fundamental data, but TT is possible and should
                  # be perfectly precise for standard applications.
                  mq = me.epoch(b'utc', b({'value': time / 86400., 'unit': 'd'}))
                  mjdtt = me.measure(b(mq), b'tt')['m0']['value']
  
                  tdata = tbins.get(mjdtt, None)
                  if tdata is None:
                      tdata = tbins[mjdtt] = [0., 0., 0., 0., 0]
  
                  if rephase:
                      uvw = cols['uvw'][:,i]
                      ph = np.exp((0-2j) * np.pi * np.dot(lmn, uvw) * freqs)
  
                  for j in range(cols['flag'].shape[0]): # all polns
                      # We just average together all polarizations right now!
                      # (Not actively, but passively by just iterating over them.)
                      data = cols[cfg.datacol][j,:,i]
                      flags = cols['flag'][j,:,i]
  
                      # XXXXX casacore is currently (ca. 2012) broken and
                      # returns the raw weights from the dataset rather than
                      # applying the polarization selection. Fortunately all of
                      # our weights are the same, and you can never fetch more
                      # pol types than the dataset has, so this bit works
                      # despite the bug.
  
                      w = np.where(~flags)[0]
                      if not w.size:
                          continue # all flagged
  
                      if rephase:
                          data *= ph
  
                      d = data[w].mean()
                      # account for flagged parts. 90% sure this is the
                      # right thing to do:
                      wt = cols['weight'][j,i] * w.size / data.size
                      wd = wt * d
                      # note a little bit of a hack here to encode real^2 and
                      # imag^2 separately:
                      wd2 = wt * (d.real**2 + (1j) * d.imag**2)
  
                      tdata[0] += wd
                      tdata[1] += wd2
                      tdata[2] += wt
                      tdata[3] += wt**2
                      tdata[4] += 1
  
              if not ms.iternext():
                  break
  
          ms.reset() # reset selection filter so we can get next DDID
  
      ms.close()
  
      # Could gain some efficiency by using a better data structure than a dict().
      smjd = sorted(six.iterkeys(tbins))
      cfg.format.header(cfg)
  
      for mjd in smjd:
          wd, wd2, wt, wt2, n = tbins[mjd]
          if n < 3: # not enough data for meaningful statistics
              continue
  
          dtmin = 1440 * (mjd - smjd[0])
          r_sc = wd.real / wt * cfg.datascale
          i_sc = wd.imag / wt * cfg.datascale
          r2_sc = wd2.real / wt * cfg.datascale**2
          i2_sc = wd2.imag / wt * cfg.datascale**2
  
          if cfg.believeweights:
              ru_sc = wt**-0.5 * cfg.datascale
              iu_sc = wt**-0.5 * cfg.datascale
          else:
              rv_sc = r2_sc - r_sc**2 # variance among real/imag msmts
              iv_sc = i2_sc - i_sc**2
              ru_sc = np.sqrt(rv_sc * wt2) / wt # uncert in mean real/img values
              iu_sc = np.sqrt(iv_sc * wt2) / wt
  
          mag = np.sqrt(r_sc**2 + i_sc**2)
          umag = np.sqrt(r_sc**2 * ru_sc**2 + i_sc**2 * iu_sc**2) / mag
          cfg.format.row(cfg, mjd, dtmin, r_sc, ru_sc, i_sc, iu_sc, mag, umag, n)