def plot(self, **kwargs):
        """
        This will plot the data (with error) for inspection.
    
        Setting self.figures to a figure instance or list of figure instances
        will override the creation of new figures. If you specify
        a list, its length had better be at least as large as the
        number of data sets.

        kwargs will update the settings
        """
        
        # Make sure there is data to plot.
        if len(self._set_xdata)==0 or len(self._set_ydata)==0: return self
        
        # Make sure the figures is a list
        if not self.figures == None and not type(self.figures) == list:
            self.figures = [self.figures]
        
        
        # Get the trimmed and full processed data
        xts, yts, eyts = self.get_processed_data()
        xas, yas, eyas = self.get_processed_data(do_trim=False)
        
        # update settings
        for k in kwargs: self[k] = kwargs[k]

        # Calculate all studentized residuals
        if len(self.f) > 0: rt = self.studentized_residuals()
        
        # make a new figure for each data set
        for n in range(len(self._set_xdata)):
            xt = xts[n]
            xa = xas[n]
            yt = yts[n]
            ya = yas[n]
            eyt = eyts[n]
            eya = eyas[n]
            #ext = exts[n]
            eyt = eyts[n]
            
            # get the next figure
            if self.figures == None: fig = _p.figure(self['first_figure']+n)
            else: fig = self.figures[n]
                    
            # turn off interactive mode and clear the figure
            _p.ioff()
            fig.clear()
            
            # set up two axes. One for data and one for residuals.
            a1 = fig.add_subplot(211)            # Residuals
            a2 = fig.add_subplot(212, sharex=a1) # Data
            a1.set_position([0.15, 0.72, 0.75, 0.15])
            a2.set_position([0.15, 0.10, 0.75, 0.60])

            # set the scales
            a1.set_xscale(self['xscale'][n])
            a2.set_xscale(self['xscale'][n])
            a2.set_yscale(self['yscale'][n])

            # Get the function xdata
            fxa = self._get_xdata_for_function(n,xa)
            fxt = self._get_xdata_for_function(n,xt)
            
            # get the values to subtract from ydata if subtracting the background
            if self['subtract_bg'][n] and not self.bg[n] is None:

                # if we have a fit, use that for the background
                if self.results: p = self.results[0]
                else:            p = self._pguess
                
                # Get the background data
                d_ya  = self._evaluate_bg(n, xa,  p)
                d_fya = self._evaluate_bg(n, fxa, p)
                d_yt  = self._evaluate_bg(n, xt,  p)
                d_fyt = self._evaluate_bg(n, fxt, p)

            # Otherwise just make some zero arrays
            else:
                d_ya  = 0*xa
                d_fya = 0*fxa
                d_yt  = 0*xt
                d_fyt = 0*fxt

            
            
            # PLOT DATA FIRST

            # If we're supposed to, add the "all" data and function
            if self['plot_all_data'][n]:
                
                # Make it faint.
                style_data  = dict(self['style_data' ][n]); style_data ['alpha'] = 0.3
                
                if self['plot_errors'][n]: a2.errorbar(xa, ya-d_ya, eya, zorder=5, **style_data)
                else:                      a2.plot    (xa, ya-d_ya,      zorder=5, **style_data)
            
            # add the trimmed data
            if self['plot_errors'][n]:     a2.errorbar(xt, yt-d_yt, eyt, zorder=7, **self['style_data'][n])
            else:                          a2.plot(    xt, yt-d_yt,      zorder=7, **self['style_data'][n])
                           
            # Zoom on just the data for now
            _s.tweaks.auto_zoom(axes=a2, draw=False)
            
            
            
            # PLOT FUNCTIONS
            
            if n < len(self.f): # If there are any functions to plot
                
                
                # Plot the GUESS under the fit
                if self['plot_guess'][n]:
                                    
                    # FULL GUESS
                    if self['plot_all_data'][n]:
                        
                        # Make it faint.
                        style_guess = dict(self['style_guess'][n]); style_guess['alpha'] = 0.3
                
                        # FULL background GUESS
                        if self['plot_bg'][n] and self.bg[n] is not None:
                            bg_gya = self._evaluate_bg(n, fxa, self._pguess)
                            a2.plot(fxa, bg_gya-d_fya, zorder=9, **style_guess)
                        
                        # FULL guess
                        gya = self._evaluate_f (n, fxa, self._pguess)
                        a2.plot(fxa, gya-d_fya, zorder=9, **style_guess)
                        
                        # Update the trimmed datas so that the points line up
                        [fxt] = _s.fun.trim_data_uber([fxt], [fxt>=min(xt), fxt<=max(xt)])

                    # TRIMMED GUESS BACKGROUND
                    
                    # TRIMMED guess background curve
                    if self['plot_bg'][n] and self.bg[n] is not None:
                        bg_gyt = self._evaluate_bg(n, fxt, self._pguess)
                        a2.plot(fxt, bg_gyt-d_fyt, zorder=9, **self['style_guess'][n])

                    # TRIMMED main guess curve
                    gyt = self._evaluate_f (n, fxt, self._pguess)
                    a2.plot(fxt, gyt-d_fyt, zorder=9, **self['style_guess'][n])
    
                
                
                # Plot the FIT if there is one
                if not self.results == None:
    
                    # FULL FIT
                    if self['plot_all_data'][n]:
                        
                        # Make it faint.
                        style_fit = dict(self['style_fit'][n]); style_fit['alpha'] = 0.3
                
                        # FULL background fit
                        if self['plot_bg'][n] and self.bg[n] is not None:
                            bg_fya = self._evaluate_bg(n, fxa, self.results[0])
                            a2.plot(fxa, bg_fya-d_fya, zorder=10, **style_fit)
                        
                        # FULL fit
                        fya = self._evaluate_f (n, fxa, self.results[0])
                        a2.plot(fxa, fya-d_fya, zorder=10, **style_fit)
    
                        # Update the trimmed datas so that the points line up
                        [fxt] = _s.fun.trim_data_uber([fxt], [fxt>=min(xt), fxt<=max(xt)])
                            
                    # TRIMMED FIT BACKGROUND
                    if self['plot_bg'][n] and self.bg[n] is not None: 
                        bg_fyt = self._evaluate_bg(n, fxt, self.results[0])                    
                        a2.plot(fxt, bg_fyt-d_fyt, zorder=10, **self['style_fit'][n])
    
                    # TRIMMED main curve
                    fyt = self._evaluate_f(n, fxt, self.results[0])
                    a2.plot(fxt, fyt-d_fyt, zorder=10, **self['style_fit'][n])

                if self['plot_guess_zoom'][n]: _s.tweaks.auto_zoom(axes=a2, draw=False)
            
            
            # plot the residuals only if there are functions defined
            if len(self.f):
                
                # If we're supposed to also plot all the data, we have to 
                # Manually calculate the residuals. Clunky, I know.
                if self['plot_all_data'][n]:
                    
                    # Figure out what guy to use for the residuals
                    if self.results is None: 
                        p = self._pguess
                        style = style_guess
                    else:                    
                        p = self.results[0]
                        style = style_fit
                    
                    # Calculate them
                    ra = (ya-self._evaluate_f(n, xa, p))/eya
                        
                    # style_data, style_guess, and style_fit should already be faint
                    a1.errorbar(xa, ra, _n.ones(len(ra)), **style_data)
                    
                    # Put the line on top
                    a1.plot([min(xa), max(xa)], [0,0], **style)
                
                # Figure out what style to use for the line
                if self.results is None: style = self['style_guess'][n]
                else:                    style = self['style_fit'  ][n]

                # Main residuals plot
                a1.errorbar (xt, rt[n], _n.ones(len(xt)), **self['style_data'][n])
                a1.plot([min(xt), max(xt)], [0,0], **style)
                
                
            # Tidy up
            yticklabels = a1.get_yticklabels()
            for m in range(2,len(yticklabels)-2): yticklabels[m].set_visible(False)
            for m in a1.get_xticklabels(): m.set_visible(False)
            
            # Add labels to the axes
            if self['xlabel'][n] is None: a2.set_xlabel('xdata['+str(n)+']')
            else:                         a2.set_xlabel(self['xlabel'][n])
            if self['ylabel'][n] is None:
                ylabel='ydata['+str(n)+']'
                if self['subtract_bg'][n] and self.bg[n] is not None:
                    ylabel=ylabel+' - bg['+str(n)+']'
                a2.set_ylabel(ylabel)
            else:                         a2.set_ylabel(self['ylabel'][n])
            a1.set_ylabel('Studentized\nResiduals')

            
            # Assemble the title
            wrap = 80
            indent = '      '
            
            # Include the function names if available
            if n < len(self.f):
                t = _textwrap.fill('Function ('+str(n)+'/'+str(len(self._fnames)-1)+'): y = '+self._fnames[n], wrap, subsequent_indent=indent)
            else:
                t = "No functions defined. Use set_functions()."

            if len(self._cnames):
                t1 = "Constants: "
                for i in range(len(self._cnames)):
                    t1 = t1 + self._cnames[i] + "={:G}, ".format(self._constants[i])
                t = t + '\n' + _textwrap.fill(t1, wrap, subsequent_indent=indent)

            if self.results and not self.results[1] is None:
                t1 = "Fit: "
                for i in range(len(self._pnames)):
                    t1 = t1 + self._pnames[i] + "={:s}, ".format(self._format_value_error(self.results[0][i], _n.sqrt(self.results[1][i][i]), '$\pm$'))
                t1 = t1 + '$\chi^2_r$={} ({} DOF)'.format(
                        self._format_value_error(self.reduced_chi_squared(), _n.sqrt(_n.divide(2.0,self.degrees_of_freedom())), '$\pm$'), 
                        int(self.degrees_of_freedom()))
                
                t = t + '\n' + _textwrap.fill(t1, wrap, subsequent_indent=indent)

            elif self.results:
                t1 = "Fit did not converge: "
                for i in range(len(self._pnames)):
                    t1 = t1 + self._pnames[i] + "={:8G}$, "
                t = t + '\n' + _textwrap.fill(t1, wrap, subsequent_indent=indent)

            a1.set_title(t, fontsize=10, ha='left', position=(0,1))

            
            # turn back to interactive and show the plots.
            _p.ion()
            if self.figures == None:
                _p.draw()
                _p.show()
        
        # End of new figure for each data set loop
        return self
This will plot the data (with error) for inspection.
    
        Setting self.figures to a figure instance or list of figure instances
        will override the creation of new figures. If you specify
        a list, its length had better be at least as large as the
        number of data sets.

        kwargs will update the settings