docs: def r_division(onarray, offarray, rarray, mode='mean'):
      """Apply R division.
  
      Args:
          onarray (decode.array): Decode array of on-point observations.
          offarray (decode.array): Decode array of off-point observations.
          rarray (decode.array): Decode array of R observations.
          mode (str): Method for the selection of nominal R value.
              'mean': Mean.
              'median': Median.
  
      Returns:
          onarray_cal (decode.array): Calibrated array of on-point observations.
          offarray_cal (decode.array): Calibrated array of off-point observations.
      """
      logger = getLogger('decode.models.r_division')
      logger.info('mode')
      logger.info('{}'.format(mode))
  
      offid = np.unique(offarray.scanid)
      onid  = np.unique(onarray.scanid)
      rid   = np.unique(rarray.scanid)
  
      onarray   = onarray.copy() # Xarray
      onvalues  = onarray.values
      onscanid  = onarray.scanid.values
      offarray  = offarray.copy() # Xarray
      offvalues = offarray.values
      offscanid = offarray.scanid.values
      rarray    = rarray.copy() # Xarray
      rvalues   = rarray.values
      rscanid   = rarray.scanid.values
      for i in onid:
          rleftid  = np.searchsorted(rid, i) - 1
          rrightid = np.searchsorted(rid, i)
  
          if rleftid == -1:
              Xr   = rvalues[rscanid == rid[rrightid]]
              Xr_m = getattr(np, 'nan'+mode)(Xr, axis=0)
          elif rrightid == len(rid):
              Xr   = rvalues[rscanid == rid[rleftid]]
              Xr_m = getattr(np, 'nan'+mode)(Xr, axis=0)
          else:
              Xr_l  = rvalues[rscanid == rid[rleftid]]
              Xr_r  = rvalues[rscanid == rid[rrightid]]
              Xr_m  = getattr(np, 'nan'+mode)(np.vstack([Xr_l, Xr_r]), axis=0)
          onvalues[onscanid == i] /= Xr_m
  
      for j in offid:
          rleftid  = np.searchsorted(rid, j) - 1
          rrightid = np.searchsorted(rid, j)
  
          Xoff_m = getattr(np, 'nan'+mode)(offvalues[offscanid == j], axis=0)
          if rleftid == -1:
              Xr   = rvalues[rscanid == rid[rrightid]]
              Xr_m = getattr(np, 'nan'+mode)(Xr, axis=0)
          elif rrightid == len(rid):
              Xr   = rvalues[rscanid == rid[rleftid]]
              Xr_m = getattr(np, 'nan'+mode)(Xr, axis=0)
          else:
              Xr_l  = rvalues[rscanid == rid[rleftid]]
              Xr_r  = rvalues[rscanid == rid[rrightid]]
              Xr_m  = getattr(np, 'nan'+mode)(np.vstack([Xr_l, Xr_r]), axis=0)
          offvalues[offscanid == j] /= Xr_m
  
      Xon_rdiv    = dc.full_like(onarray, onarray)
      Xoff_rdiv   = dc.full_like(offarray, offarray)
      Xonoff_rdiv = dc.concat([Xon_rdiv, Xoff_rdiv], dim='t')
      Xonoff_rdiv_sorted = Xonoff_rdiv[np.argsort(Xonoff_rdiv.time.values)]
  
      scantype    = Xonoff_rdiv_sorted.scantype.values
      newscanid   = np.cumsum(np.hstack([False, scantype[1:] != scantype[:-1]]))
      onmask      = np.in1d(Xonoff_rdiv_sorted.scanid, onid)
      offmask     = np.in1d(Xonoff_rdiv_sorted.scanid, offid)
      Xon_rdiv    = Xonoff_rdiv_sorted[onmask]
      Xoff_rdiv   = Xonoff_rdiv_sorted[offmask]
      Xon_rdiv.coords.update({'scanid': ('t', newscanid[onmask])})
      Xoff_rdiv.coords.update({'scanid': ('t', newscanid[offmask])})
  
      return Xon_rdiv, Xoff_rdiv