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