docs: def RV_RMSgrid(self,ts,res=20,mres=None,Pres=None,
                     conf=0.95,measured_rms=None,drv=0,
                     plot=True,fig=None,contour=True,sigma=1):
          """Writes a grid of RV RMS values, assuming observations at given times.
  
          Caveat Emptor: Written a long time ago, and
          hasn't really been tested.
  
          :param ts:
              Times of observations
  
          :param res, mres, Pres: (optional)
              Resolution of grids.  ``res`` relates to both mass and period;
              otherwise ``mres`` and ``Pres`` can be set separately.
  
          :param conf: (optional)
              Confidence at which to exclude regions.  Used if ``measured_rms``
              is ``None``.
  
          :param measured_rms: (optional)
              Measured RV RMS, if applicable [not sure exactly how this is used]
  
          :param drv: (optional)
              Uncertainties in RV to simulate.
  
          :param plot: (optional)
              Whether to plot result.
  
          :param fig: (optional)
              Passed to :func:`plotutils.setfig`.
  
          :param contour: (optional)
              Whether to plot contours.
  
          :param sigma: (optional)
              Level at which to exclude, based on ``measured_rms``.
  
          """
          RVs = self.RV_timeseries(ts)
          RVs += rand.normal(size=np.size(RVs)).reshape(RVs.shape)*drv
          rms = RVs.std(axis=0)
  
          if mres is None:
              mres = res
          if Pres is None:
              Pres = res
  
          mbins = np.linspace(self.M2.min(),self.M2.max(),mres+1)
          Pbins = np.logspace(np.log10(self.P.min()),np.log10(self.P.max()),Pres+1)
          logPbins = np.log10(Pbins)
  
          mbin_centers = (mbins[:-1] + mbins[1:])/2.
          logPbin_centers = (logPbins[:-1] + logPbins[1:])/2.
  
          minds = np.digitize(self.M2,mbins)
          Pinds = np.digitize(self.P,Pbins)
  
          pctiles = np.zeros((mres,Pres))
          ns = np.zeros((mres,Pres))
  
          for i in np.arange(mres):
              for j in np.arange(Pres):
                  w = np.where((minds==i+1) & (Pinds==j+1))
                  these = rms[w]
                  n = size(these)
                  ns[i,j] = n
                  if measured_rms is not None:
                      pctiles[i,j] = (these > sigma*measured_rms).sum()/float(n)
                  else:
                      inds = np.argsort(these)
                      pctiles[i,j] = these[inds][int((1-conf)*n)]
  
          Ms,logPs = np.meshgrid(mbin_centers,logPbin_centers)
  
          if plot:
              setfig(fig)
  
              if contour:
                  mbin_centers = (mbins[:-1] + mbins[1:])/2.
                  logPbins = np.log10(Pbins)
                  logPbin_centers = (logPbins[:-1] + logPbins[1:])/2.
                  if measured_rms is not None:
                      levels = [0.68,0.95,0.99]
                  else:
                      levels = np.arange(0,20,2)
                  c = plt.contour(logPbin_centers,mbin_centers,pctiles,levels=levels,colors='k')
                  plt.clabel(c, fontsize=10, inline=1)
  
              else:
                  extent = [np.log10(self.P.min()),np.log10(self.P.max()),self.M2.min(),self.M2.max()]
                  im = plt.imshow(pctiles,cmap='Greys',extent=extent,aspect='auto')
  
                  fig = plt.gcf()
                  ax = plt.gca()
  
  
                  if measured_rms is None:
                      cbarticks = np.arange(0,21,2)
                  else:
                      cbarticks = np.arange(0,1.01,0.1)
                  cbar = fig.colorbar(im, ticks=cbarticks)
  
              plt.xlabel('Log P')
              plt.ylabel('M2')
  
          return mbins,Pbins,pctiles,ns