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``.

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