Index: src/sans/pr/Doxyfile
===================================================================
--- src/sans/pr/Doxyfile	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/Doxyfile	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,263 @@
+# Doxyfile 1.5.1-p1
+
+#---------------------------------------------------------------------------
+# Project related configuration options
+#---------------------------------------------------------------------------
+PROJECT_NAME           = sanspr
+PROJECT_NUMBER         = 
+OUTPUT_DIRECTORY       = doc
+CREATE_SUBDIRS         = NO
+OUTPUT_LANGUAGE        = English
+USE_WINDOWS_ENCODING   = NO
+BRIEF_MEMBER_DESC      = YES
+REPEAT_BRIEF           = YES
+ABBREVIATE_BRIEF       = "The $name class" \
+                         "The $name widget" \
+                         "The $name file" \
+                         is \
+                         provides \
+                         specifies \
+                         contains \
+                         represents \
+                         a \
+                         an \
+                         the
+ALWAYS_DETAILED_SEC    = NO
+INLINE_INHERITED_MEMB  = NO
+FULL_PATH_NAMES        = YES
+STRIP_FROM_PATH        = 
+STRIP_FROM_INC_PATH    = 
+SHORT_NAMES            = NO
+JAVADOC_AUTOBRIEF      = NO
+MULTILINE_CPP_IS_BRIEF = NO
+DETAILS_AT_TOP         = NO
+INHERIT_DOCS           = YES
+SEPARATE_MEMBER_PAGES  = NO
+TAB_SIZE               = 4
+ALIASES                = 
+OPTIMIZE_OUTPUT_FOR_C  = NO
+OPTIMIZE_OUTPUT_JAVA   = YES
+BUILTIN_STL_SUPPORT    = NO
+DISTRIBUTE_GROUP_DOC   = NO
+SUBGROUPING            = YES
+#---------------------------------------------------------------------------
+# Build related configuration options
+#---------------------------------------------------------------------------
+EXTRACT_ALL            = NO
+EXTRACT_PRIVATE        = NO
+EXTRACT_STATIC         = NO
+EXTRACT_LOCAL_CLASSES  = YES
+EXTRACT_LOCAL_METHODS  = NO
+HIDE_UNDOC_MEMBERS     = NO
+HIDE_UNDOC_CLASSES     = NO
+HIDE_FRIEND_COMPOUNDS  = NO
+HIDE_IN_BODY_DOCS      = NO
+INTERNAL_DOCS          = NO
+CASE_SENSE_NAMES       = YES
+HIDE_SCOPE_NAMES       = NO
+SHOW_INCLUDE_FILES     = YES
+INLINE_INFO            = YES
+SORT_MEMBER_DOCS       = YES
+SORT_BRIEF_DOCS        = NO
+SORT_BY_SCOPE_NAME     = NO
+GENERATE_TODOLIST      = YES
+GENERATE_TESTLIST      = YES
+GENERATE_BUGLIST       = YES
+GENERATE_DEPRECATEDLIST= YES
+ENABLED_SECTIONS       = 
+MAX_INITIALIZER_LINES  = 30
+SHOW_USED_FILES        = YES
+SHOW_DIRECTORIES       = YES
+FILE_VERSION_FILTER    = 
+#---------------------------------------------------------------------------
+# configuration options related to warning and progress messages
+#---------------------------------------------------------------------------
+QUIET                  = NO
+WARNINGS               = YES
+WARN_IF_UNDOCUMENTED   = YES
+WARN_IF_DOC_ERROR      = YES
+WARN_NO_PARAMDOC       = NO
+WARN_FORMAT            = "$file:$line: $text"
+WARN_LOGFILE           = 
+#---------------------------------------------------------------------------
+# configuration options related to the input files
+#---------------------------------------------------------------------------
+INPUT                  = .
+FILE_PATTERNS          = *.c \
+                         *.cc \
+                         *.cxx \
+                         *.cpp \
+                         *.c++ \
+                         *.d \
+                         *.java \
+                         *.ii \
+                         *.ixx \
+                         *.ipp \
+                         *.i++ \
+                         *.inl \
+                         *.h \
+                         *.hh \
+                         *.hxx \
+                         *.hpp \
+                         *.h++ \
+                         *.idl \
+                         *.odl \
+                         *.cs \
+                         *.php \
+                         *.php3 \
+                         *.inc \
+                         *.m \
+                         *.mm \
+                         *.dox \
+                         *.py
+RECURSIVE              = YES
+EXCLUDE                = build \
+                         dist \
+                         test
+EXCLUDE_SYMLINKS       = NO
+EXCLUDE_PATTERNS       = 
+EXAMPLE_PATH           = 
+EXAMPLE_PATTERNS       = *
+EXAMPLE_RECURSIVE      = NO
+IMAGE_PATH             = doc
+INPUT_FILTER           = 
+FILTER_PATTERNS        = 
+FILTER_SOURCE_FILES    = NO
+#---------------------------------------------------------------------------
+# configuration options related to source browsing
+#---------------------------------------------------------------------------
+SOURCE_BROWSER         = NO
+INLINE_SOURCES         = NO
+STRIP_CODE_COMMENTS    = YES
+REFERENCED_BY_RELATION = NO
+REFERENCES_RELATION    = NO
+REFERENCES_LINK_SOURCE = YES
+USE_HTAGS              = NO
+VERBATIM_HEADERS       = YES
+#---------------------------------------------------------------------------
+# configuration options related to the alphabetical class index
+#---------------------------------------------------------------------------
+ALPHABETICAL_INDEX     = NO
+COLS_IN_ALPHA_INDEX    = 5
+IGNORE_PREFIX          = 
+#---------------------------------------------------------------------------
+# configuration options related to the HTML output
+#---------------------------------------------------------------------------
+GENERATE_HTML          = YES
+HTML_OUTPUT            = html
+HTML_FILE_EXTENSION    = .html
+HTML_HEADER            = 
+HTML_FOOTER            = 
+HTML_STYLESHEET        = 
+HTML_ALIGN_MEMBERS     = YES
+GENERATE_HTMLHELP      = NO
+CHM_FILE               = 
+HHC_LOCATION           = 
+GENERATE_CHI           = NO
+BINARY_TOC             = NO
+TOC_EXPAND             = NO
+DISABLE_INDEX          = NO
+ENUM_VALUES_PER_LINE   = 4
+GENERATE_TREEVIEW      = NO
+TREEVIEW_WIDTH         = 250
+#---------------------------------------------------------------------------
+# configuration options related to the LaTeX output
+#---------------------------------------------------------------------------
+GENERATE_LATEX         = NO
+LATEX_OUTPUT           = latex
+LATEX_CMD_NAME         = latex
+MAKEINDEX_CMD_NAME     = makeindex
+COMPACT_LATEX          = NO
+PAPER_TYPE             = a4wide
+EXTRA_PACKAGES         = 
+LATEX_HEADER           = 
+PDF_HYPERLINKS         = NO
+USE_PDFLATEX           = YES
+LATEX_BATCHMODE        = NO
+LATEX_HIDE_INDICES     = NO
+#---------------------------------------------------------------------------
+# configuration options related to the RTF output
+#---------------------------------------------------------------------------
+GENERATE_RTF           = NO
+RTF_OUTPUT             = rtf
+COMPACT_RTF            = NO
+RTF_HYPERLINKS         = NO
+RTF_STYLESHEET_FILE    = 
+RTF_EXTENSIONS_FILE    = 
+#---------------------------------------------------------------------------
+# configuration options related to the man page output
+#---------------------------------------------------------------------------
+GENERATE_MAN           = NO
+MAN_OUTPUT             = man
+MAN_EXTENSION          = .3
+MAN_LINKS              = NO
+#---------------------------------------------------------------------------
+# configuration options related to the XML output
+#---------------------------------------------------------------------------
+GENERATE_XML           = NO
+XML_OUTPUT             = xml
+XML_SCHEMA             = 
+XML_DTD                = 
+XML_PROGRAMLISTING     = YES
+#---------------------------------------------------------------------------
+# configuration options for the AutoGen Definitions output
+#---------------------------------------------------------------------------
+GENERATE_AUTOGEN_DEF   = NO
+#---------------------------------------------------------------------------
+# configuration options related to the Perl module output
+#---------------------------------------------------------------------------
+GENERATE_PERLMOD       = NO
+PERLMOD_LATEX          = NO
+PERLMOD_PRETTY         = YES
+PERLMOD_MAKEVAR_PREFIX = 
+#---------------------------------------------------------------------------
+# Configuration options related to the preprocessor   
+#---------------------------------------------------------------------------
+ENABLE_PREPROCESSING   = YES
+MACRO_EXPANSION        = NO
+EXPAND_ONLY_PREDEF     = NO
+SEARCH_INCLUDES        = YES
+INCLUDE_PATH           = 
+INCLUDE_FILE_PATTERNS  = 
+PREDEFINED             = 
+EXPAND_AS_DEFINED      = 
+SKIP_FUNCTION_MACROS   = YES
+#---------------------------------------------------------------------------
+# Configuration::additions related to external references   
+#---------------------------------------------------------------------------
+TAGFILES               = 
+GENERATE_TAGFILE       = 
+ALLEXTERNALS           = NO
+EXTERNAL_GROUPS        = YES
+PERL_PATH              = /usr/bin/perl
+#---------------------------------------------------------------------------
+# Configuration options related to the dot tool   
+#---------------------------------------------------------------------------
+CLASS_DIAGRAMS         = YES
+HIDE_UNDOC_RELATIONS   = YES
+HAVE_DOT               = YES
+CLASS_GRAPH            = YES
+COLLABORATION_GRAPH    = NO
+GROUP_GRAPHS           = YES
+UML_LOOK               = YES
+TEMPLATE_RELATIONS     = YES
+INCLUDE_GRAPH          = YES
+INCLUDED_BY_GRAPH      = YES
+CALL_GRAPH             = YES
+CALLER_GRAPH           = YES
+GRAPHICAL_HIERARCHY    = YES
+DIRECTORY_GRAPH        = YES
+DOT_IMAGE_FORMAT       = png
+DOT_PATH               = 
+DOTFILE_DIRS           = 
+MAX_DOT_GRAPH_WIDTH    = 1024
+MAX_DOT_GRAPH_HEIGHT   = 1024
+MAX_DOT_GRAPH_DEPTH    = 1000
+DOT_TRANSPARENT        = NO
+DOT_MULTI_TARGETS      = NO
+GENERATE_LEGEND        = YES
+DOT_CLEANUP            = YES
+#---------------------------------------------------------------------------
+# Configuration::additions related to the search engine   
+#---------------------------------------------------------------------------
+SEARCHENGINE           = NO
Index: src/sans/pr/__init__.py
===================================================================
--- src/sans/pr/__init__.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/__init__.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,106 @@
+"""
+    P(r) inversion for SANS
+"""
+## \mainpage P(r) inversion for SANS
+#
+# \section intro_sec Introduction
+# This module provides calculations to transform scattering intensity data 
+# I(q) into distance distribution function P(r). A description of the 
+# technique can be found elsewhere [1-5]. The module is useable as a 
+# standalone application but its functionality is meant to be presented 
+# to end-users through the user interface developed as part of the SANS 
+# flagship application. 
+#
+# Procedure: We will follow the procedure of Moore [1].
+# 
+# [1] P.B. Moore, J.Appl. Cryst (1980) 13, 168-175.
+#
+# [2] O. Glatter, J.Appl. Cryst (1977) 10, 415-421.
+#
+# [3] D.I. Svergun, J.Appl. Cryst (1991) 24, 485-492.
+#
+# [4] D.I. Svergun, J.Appl. Cryst (1992) 25, 495-503.
+#
+# [5] S. Hansen and J. Skov Pedersen, J.Appl. Cryst (1991) 24, 541-548.
+#
+## \subsection class Class Diagram:
+# The following shows a partial class diagram with the main attributes
+# and methods of the invertor.
+#
+# \image html architecture.png
+#
+# \section install_sec Installation
+#
+# \subsection obtain Obtaining the Code
+#
+# The code is available here:
+# \verbatim
+#$ svn co svn://danse.us/sans/pr_inversion
+# \endverbatim
+#
+# \subsection depends External Dependencies
+# scipy, numpy
+#
+# \subsection build Building the code
+# The standard python package can be built with distutils.
+# \verbatim
+#$ python setup.py build
+#$ python setup.py install
+# \endverbatim
+#
+#
+# \subsection Tutorial
+# To create an inversion object:
+# \verbatim
+#from sans.pr.invertor import Invertor
+#    invertor = Invertor()
+# \endverbatim
+#
+# To set the maximum distance between any two points:
+# \verbatim
+#    invertor.d_max = 160.0
+# \endverbatim
+#
+# To set the regularization constant:
+# \verbatim
+#    invertor.alpha = 0.0007
+# \endverbatim
+#
+# To set the q, I(q) and error on I(q):
+# \verbatim
+#    invertor.x = q_vector
+#    invertor.y = Iq_vector
+#    invertor.err = dIq_vector
+# \endverbatim
+#
+# To perform the inversion. In this example, we choose
+# a P(r) expension wit 10 base functions.
+# \verbatim
+#    c_out, c_cov = invertor.invert(10)
+# \endverbatim
+# The c_out and c_cov are the set of coefficients and the covariance 
+# matrix for those coefficients, respectively.
+#
+# To get P(r):
+# \verbatim
+#    r = 10.0
+#    pr = invertor.pr(c_out, r)
+# \endverbatim
+# Alternatively, one can get P(r) with the error on P(r):
+# \verbatim
+#    r = 10.0
+#    pr, dpr = invertor.pr_err(c_out, c_cov, r)
+# \endverbatim
+#
+# To get the output I(q) from the set of coefficients found:
+# \verbatim
+#    q = 0.001
+#    iq = invertor.iq(c_out, q)
+# \endverbatim
+#
+# Examples are available as unit tests under sans.pr_inversion.test.
+#
+# \section help_sec Contact Info
+# Code and Documentation produced as part of the DANSE project.
+
+__author__ = 'University of Tennessee'
Index: src/sans/pr/c_extensions/Cinvertor.c
===================================================================
--- src/sans/pr/c_extensions/Cinvertor.c	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/c_extensions/Cinvertor.c	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,1141 @@
+/**
+ * C implementation of the P(r) inversion
+ * Cinvertor is the base class for the Invertor class
+ * and provides the underlying computations.
+ *
+ */
+#include <Python.h>
+#include "structmember.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+#include "invertor.h"
+
+
+/// Error object for raised exceptions
+PyObject * CinvertorError;
+
+#define INVECTOR(obj,buf,len)										\
+    do { \
+        int err = PyObject_AsReadBuffer(obj, (const void **)(&buf), &len); \
+        if (err < 0) return NULL; \
+        len /= sizeof(*buf); \
+    } while (0)
+
+#define OUTVECTOR(obj,buf,len) \
+    do { \
+        int err = PyObject_AsWriteBuffer(obj, (void **)(&buf), &len); \
+        if (err < 0) return NULL; \
+        len /= sizeof(*buf); \
+    } while (0)
+
+
+// Class definition
+/**
+ * C implementation of the P(r) inversion
+ * Cinvertor is the base class for the Invertor class
+ * and provides the underlying computations.
+ *
+ */
+typedef struct {
+    PyObject_HEAD
+    /// Internal data structure
+    Invertor_params params;
+} Cinvertor;
+
+
+static void
+Cinvertor_dealloc(Cinvertor* self)
+{
+    invertor_dealloc(&(self->params));
+
+    self->ob_type->tp_free((PyObject*)self);
+
+}
+
+static PyObject *
+Cinvertor_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    Cinvertor *self;
+
+    self = (Cinvertor *)type->tp_alloc(type, 0);
+
+    return (PyObject *)self;
+}
+
+static int
+Cinvertor_init(Cinvertor *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {
+    	// Create parameters
+    	invertor_init(&(self->params));
+    }
+    return 0;
+}
+
+static PyMemberDef Cinvertor_members[] = {
+    //{"params", T_OBJECT, offsetof(Cinvertor, params), 0,
+    // "Parameters"},
+    {NULL}  /* Sentinel */
+};
+
+const char set_x_doc[] =
+	"Function to set the x data\n"
+	"Takes an array of doubles as input.\n"
+	" @return: number of entries found";
+
+/**
+ * Function to set the x data
+ * Takes an array of doubles as input
+ * Returns the number of entries found
+ */
+static PyObject * set_x(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+	int i;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,data,ndata);
+
+	free(self->params.x);
+	self->params.x = (double*) malloc(ndata*sizeof(double));
+
+	if(self->params.x==NULL) {
+	    PyErr_SetString(CinvertorError,
+	    	"Cinvertor.set_x: problem allocating memory.");
+		return NULL;
+	}
+
+	for (i=0; i<ndata; i++) {
+		self->params.x[i] = data[i];
+	}
+
+	//self->params.x = data;
+	self->params.npoints = ndata;
+	return Py_BuildValue("i", self->params.npoints);
+}
+
+const char get_x_doc[] =
+	"Function to get the x data\n"
+	"Takes an array of doubles as input.\n"
+	" @return: number of entries found";
+
+static PyObject * get_x(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+    int i;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj, data, ndata);
+
+	// Check that the input array is large enough
+	if (ndata < self->params.npoints) {
+	    PyErr_SetString(CinvertorError,
+	    	"Cinvertor.get_x: input array too short for data.");
+		return NULL;
+	}
+
+	for(i=0; i<self->params.npoints; i++){
+		data[i] = self->params.x[i];
+	}
+
+	return Py_BuildValue("i", self->params.npoints);
+}
+
+const char set_y_doc[] =
+	"Function to set the y data\n"
+	"Takes an array of doubles as input.\n"
+	" @return: number of entries found";
+
+/**
+ * Function to set the y data
+ * Takes an array of doubles as input
+ * Returns the number of entries found
+ */
+static PyObject * set_y(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+	int i;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,data,ndata);
+
+	free(self->params.y);
+	self->params.y = (double*) malloc(ndata*sizeof(double));
+
+	if(self->params.y==NULL) {
+	    PyErr_SetString(CinvertorError,
+	    	"Cinvertor.set_y: problem allocating memory.");
+		return NULL;
+	}
+
+	for (i=0; i<ndata; i++) {
+		self->params.y[i] = data[i];
+	}
+
+	//self->params.y = data;
+	self->params.ny = ndata;
+	return Py_BuildValue("i", self->params.ny);
+}
+
+const char get_y_doc[] =
+	"Function to get the y data\n"
+	"Takes an array of doubles as input.\n"
+	" @return: number of entries found";
+
+static PyObject * get_y(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+    int i;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj, data, ndata);
+
+	// Check that the input array is large enough
+	if (ndata < self->params.ny) {
+	    PyErr_SetString(CinvertorError,
+	    	"Cinvertor.get_y: input array too short for data.");
+		return NULL;
+	}
+
+	for(i=0; i<self->params.ny; i++){
+		data[i] = self->params.y[i];
+	}
+
+	return Py_BuildValue("i", self->params.npoints);
+}
+
+const char set_err_doc[] =
+	"Function to set the err data\n"
+	"Takes an array of doubles as input.\n"
+	" @return: number of entries found";
+
+/**
+ * Function to set the x data
+ * Takes an array of doubles as input
+ * Returns the number of entries found
+ */
+static PyObject * set_err(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+	int i;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,data,ndata);
+
+	free(self->params.err);
+	self->params.err = (double*) malloc(ndata*sizeof(double));
+
+	if(self->params.err==NULL) {
+	    PyErr_SetString(CinvertorError,
+	    	"Cinvertor.set_err: problem allocating memory.");
+		return NULL;
+	}
+
+	for (i=0; i<ndata; i++) {
+		self->params.err[i] = data[i];
+	}
+
+	//self->params.err = data;
+	self->params.nerr = ndata;
+	return Py_BuildValue("i", self->params.nerr);
+}
+
+const char get_err_doc[] =
+	"Function to get the err data\n"
+	"Takes an array of doubles as input.\n"
+	" @return: number of entries found";
+
+static PyObject * get_err(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+    int i;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj, data, ndata);
+
+	// Check that the input array is large enough
+	if (ndata < self->params.nerr) {
+	    PyErr_SetString(CinvertorError,
+	    	"Cinvertor.get_err: input array too short for data.");
+		return NULL;
+	}
+
+	for(i=0; i<self->params.nerr; i++){
+		data[i] = self->params.err[i];
+	}
+
+	return Py_BuildValue("i", self->params.npoints);
+}
+
+const char is_valid_doc[] =
+	"Check the validity of the stored data\n"
+	" @return: Returns the number of points if it's all good, -1 otherwise";
+
+/**
+ * Check the validity of the stored data
+ * Returns the number of points if it's all good, -1 otherwise
+ */
+static PyObject * is_valid(Cinvertor *self, PyObject *args) {
+	if(self->params.npoints==self->params.ny &&
+			self->params.npoints==self->params.nerr) {
+		return Py_BuildValue("i", self->params.npoints);
+	} else {
+		return Py_BuildValue("i", -1);
+	}
+}
+
+const char set_has_bck_doc[] =
+	"Sets background flag\n";
+
+/**
+ * Sets the maximum distance
+ */
+static PyObject * set_has_bck(Cinvertor *self, PyObject *args) {
+	int has_bck;
+
+	if (!PyArg_ParseTuple(args, "i", &has_bck)) return NULL;
+	self->params.has_bck = has_bck;
+	return Py_BuildValue("i", self->params.has_bck);
+}
+
+const char get_has_bck_doc[] =
+	"Gets background flag\n";
+
+/**
+ * Gets the maximum distance
+ */
+static PyObject * get_has_bck(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("i", self->params.has_bck);
+}
+
+const char set_dmax_doc[] =
+	"Sets the maximum distance\n";
+
+/**
+ * Sets the maximum distance
+ */
+static PyObject * set_dmax(Cinvertor *self, PyObject *args) {
+	double d_max;
+
+	if (!PyArg_ParseTuple(args, "d", &d_max)) return NULL;
+	self->params.d_max = d_max;
+	return Py_BuildValue("d", self->params.d_max);
+}
+
+const char get_dmax_doc[] =
+	"Gets the maximum distance\n";
+
+/**
+ * Gets the maximum distance
+ */
+static PyObject * get_dmax(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.d_max);
+}
+
+const char set_slit_height_doc[] =
+	"Sets the slit height in units of q [A-1]\n";
+
+/**
+ * Sets the slit height
+ */
+static PyObject * set_slit_height(Cinvertor *self, PyObject *args) {
+	double slit_height;
+
+	if (!PyArg_ParseTuple(args, "d", &slit_height)) return NULL;
+	self->params.slit_height = slit_height;
+	return Py_BuildValue("d", self->params.slit_height);
+}
+
+const char get_slit_height_doc[] =
+	"Gets the slit height\n";
+
+/**
+ * Gets the slit height
+ */
+static PyObject * get_slit_height(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.slit_height);
+}
+
+const char set_slit_width_doc[] =
+	"Sets the slit width in units of q [A-1]\n";
+
+/**
+ * Sets the slit width
+ */
+static PyObject * set_slit_width(Cinvertor *self, PyObject *args) {
+	double slit_width;
+
+	if (!PyArg_ParseTuple(args, "d", &slit_width)) return NULL;
+	self->params.slit_width = slit_width;
+	return Py_BuildValue("d", self->params.slit_width);
+}
+
+const char get_slit_width_doc[] =
+	"Gets the slit width\n";
+
+/**
+ * Gets the slit width
+ */
+static PyObject * get_slit_width(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.slit_width);
+}
+
+
+const char set_qmin_doc[] =
+	"Sets the minimum q\n";
+
+/**
+ * Sets the minimum q
+ */
+static PyObject * set_qmin(Cinvertor *self, PyObject *args) {
+	double q_min;
+
+	if (!PyArg_ParseTuple(args, "d", &q_min)) return NULL;
+	self->params.q_min = q_min;
+	return Py_BuildValue("d", self->params.q_min);
+}
+
+const char get_qmin_doc[] =
+	"Gets the minimum q\n";
+
+/**
+ * Gets the minimum q
+ */
+static PyObject * get_qmin(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.q_min);
+}
+
+const char set_qmax_doc[] =
+	"Sets the maximum q\n";
+
+/**
+ * Sets the maximum q
+ */
+static PyObject * set_qmax(Cinvertor *self, PyObject *args) {
+	double q_max;
+
+	if (!PyArg_ParseTuple(args, "d", &q_max)) return NULL;
+	self->params.q_max = q_max;
+	return Py_BuildValue("d", self->params.q_max);
+}
+
+const char get_qmax_doc[] =
+	"Gets the maximum q\n";
+
+/**
+ * Gets the maximum q
+ */
+static PyObject * get_qmax(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.q_max);
+}
+
+const char set_alpha_doc[] =
+	"Sets the alpha parameter\n";
+
+static PyObject * set_alpha(Cinvertor *self, PyObject *args) {
+	double alpha;
+
+	if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL;
+	self->params.alpha = alpha;
+	return Py_BuildValue("d", self->params.alpha);
+}
+
+const char get_alpha_doc[] =
+	"Gets the alpha parameter\n";
+
+/**
+ * Gets the maximum distance
+ */
+static PyObject * get_alpha(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.alpha);
+}
+
+const char get_nx_doc[] =
+	"Gets the number of x points\n";
+
+/**
+ * Gets the number of x points
+ */
+static PyObject * get_nx(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("i", self->params.npoints);
+}
+
+const char get_ny_doc[] =
+	"Gets the number of y points\n";
+
+/**
+ * Gets the number of y points
+ */
+static PyObject * get_ny(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("i", self->params.ny);
+}
+
+const char get_nerr_doc[] =
+	"Gets the number of err points\n";
+
+/**
+ * Gets the number of error points
+ */
+static PyObject * get_nerr(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("i", self->params.nerr);
+}
+
+
+const char residuals_doc[] =
+	"Function to call to evaluate the residuals\n"
+	"for P(r) inversion\n"
+	" @param args: input parameters\n"
+	" @return: list of residuals";
+
+/**
+ * Function to call to evaluate the residuals
+ * @param args: input parameters
+ * @return: list of residuals
+ */
+static PyObject * residuals(Cinvertor *self, PyObject *args) {
+	double *pars;
+	PyObject* residuals;
+	int i;
+	double residual, diff;
+	// Regularization factor
+	double regterm = 0.0;
+	double tmp = 0.0;
+	// Number of slices in regularization term estimate
+	int nslice = 25;
+
+	PyObject *data_obj;
+	Py_ssize_t npars;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+
+	OUTVECTOR(data_obj,pars,npars);
+
+    // PyList of residuals
+	// Should create this list only once and refill it
+    residuals = PyList_New(self->params.npoints);
+
+    regterm = reg_term(pars, self->params.d_max, npars, nslice);
+
+    for(i=0; i<self->params.npoints; i++) {
+    	diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]);
+    	residual = diff*diff / (self->params.err[i]*self->params.err[i]);
+    	tmp = residual;
+
+    	// regularization term
+    	residual += self->params.alpha * regterm;
+
+    	if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
+    	    PyErr_SetString(CinvertorError,
+    	    	"Cinvertor.residuals: error setting residual.");
+    		return NULL;
+    	};
+
+    }
+
+	return residuals;
+}
+
+const char pr_residuals_doc[] =
+	"Function to call to evaluate the residuals\n"
+	"for P(r) minimization (for testing purposes)\n"
+	" @param args: input parameters\n"
+	" @return: list of residuals";
+
+/**
+ * Function to call to evaluate the residuals
+ * for P(r) minimization (for testing purposes)
+ * @param args: input parameters
+ * @return: list of residuals
+ */
+static PyObject * pr_residuals(Cinvertor *self, PyObject *args) {
+	double *pars;
+	PyObject* residuals;
+	int i;
+	double residual, diff;
+	// Regularization factor
+	double regterm = 0.0;
+	double tmp = 0.0;
+	// Number of slices in regularization term estimate
+	int nslice = 25;
+
+	PyObject *data_obj;
+	Py_ssize_t npars;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+
+	OUTVECTOR(data_obj,pars,npars);
+
+	// Should create this list only once and refill it
+    residuals = PyList_New(self->params.npoints);
+
+    regterm = reg_term(pars, self->params.d_max, npars, nslice);
+
+
+    for(i=0; i<self->params.npoints; i++) {
+    	diff = self->params.y[i] - pr(pars, self->params.d_max, npars, self->params.x[i]);
+    	residual = diff*diff / (self->params.err[i]*self->params.err[i]);
+    	tmp = residual;
+
+    	// regularization term
+    	residual += self->params.alpha * regterm;
+
+    	if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
+    	    PyErr_SetString(CinvertorError,
+    	    	"Cinvertor.residuals: error setting residual.");
+    		return NULL;
+    	};
+
+    }
+
+	return residuals;
+}
+
+const char get_iq_doc[] =
+	"Function to call to evaluate the scattering intensity\n"
+	" @param args: c-parameters, and q\n"
+	" @return: I(q)";
+
+/**
+ * Function to call to evaluate the scattering intensity
+ * @param args: c-parameters, and q
+ * @return: I(q)
+ */
+static PyObject * get_iq(Cinvertor *self, PyObject *args) {
+	double *pars;
+	double q, iq_value;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+
+	if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	iq_value = iq(pars, self->params.d_max, npars, q);
+	return Py_BuildValue("f", iq_value);
+}
+
+const char get_iq_smeared_doc[] =
+	"Function to call to evaluate the scattering intensity.\n"
+	"The scattering intensity is slit-smeared."
+	" @param args: c-parameters, and q\n"
+	" @return: I(q)";
+
+/**
+ * Function to call to evaluate the scattering intensity
+ * The scattering intensity is slit-smeared.
+ * @param args: c-parameters, and q
+ * @return: I(q)
+ */
+static PyObject * get_iq_smeared(Cinvertor *self, PyObject *args) {
+	double *pars;
+	double q, iq_value;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+
+	if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	iq_value = iq_smeared(pars, self->params.d_max, npars,
+							self->params.slit_height, self->params.slit_width,
+							q, 21);
+	return Py_BuildValue("f", iq_value);
+}
+
+const char get_pr_doc[] =
+	"Function to call to evaluate P(r)\n"
+	" @param args: c-parameters and r\n"
+	" @return: P(r)";
+
+/**
+ * Function to call to evaluate P(r)
+ * @param args: c-parameters and r
+ * @return: P(r)
+ */
+static PyObject * get_pr(Cinvertor *self, PyObject *args) {
+	double *pars;
+	double r, pr_value;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+
+	if (!PyArg_ParseTuple(args, "Od", &data_obj, &r)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	pr_value = pr(pars, self->params.d_max, npars, r);
+	return Py_BuildValue("f", pr_value);
+}
+
+const char get_pr_err_doc[] =
+	"Function to call to evaluate P(r) with errors\n"
+	" @param args: c-parameters and r\n"
+	" @return: (P(r),dP(r))";
+
+/**
+ * Function to call to evaluate P(r) with errors
+ * @param args: c-parameters and r
+ * @return: P(r)
+ */
+static PyObject * get_pr_err(Cinvertor *self, PyObject *args) {
+	double *pars;
+	double *pars_err;
+	double pr_err_value;
+	double r, pr_value;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	PyObject *err_obj;
+	Py_ssize_t npars2;
+
+	if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	if (err_obj == Py_None) {
+		pr_value = pr(pars, self->params.d_max, npars, r);
+		pr_err_value = 0.0;
+	} else {
+		OUTVECTOR(err_obj,pars_err,npars2);
+		pr_err(pars, pars_err, self->params.d_max, npars, r, &pr_value, &pr_err_value);
+	}
+	return Py_BuildValue("ff", pr_value, pr_err_value);
+}
+
+const char basefunc_ft_doc[] =
+	"Returns the value of the nth Fourier transofrmed base function\n"
+	" @param args: c-parameters, n and q\n"
+	" @return: nth Fourier transformed base function, evaluated at q";
+
+static PyObject * basefunc_ft(Cinvertor *self, PyObject *args) {
+	double d_max, q;
+	int n;
+
+	if (!PyArg_ParseTuple(args, "did", &d_max, &n, &q)) return NULL;
+	return Py_BuildValue("f", ortho_transformed(d_max, n, q));
+
+}
+
+const char oscillations_doc[] =
+	"Returns the value of the oscillation figure of merit for\n"
+	"the given set of coefficients. For a sphere, the oscillation\n"
+	"figure of merit is 1.1.\n"
+	" @param args: c-parameters\n"
+	" @return: oscillation figure of merit";
+
+static PyObject * oscillations(Cinvertor *self, PyObject *args) {
+	double *pars;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	double oscill, norm;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	oscill = reg_term(pars, self->params.d_max, npars, 100);
+	norm   = int_p2(pars, self->params.d_max, npars, 100);
+	return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max );
+
+}
+
+const char get_peaks_doc[] =
+	"Returns the number of peaks in the output P(r) distrubution\n"
+	"for the given set of coefficients.\n"
+	" @param args: c-parameters\n"
+	" @return: number of P(r) peaks";
+
+static PyObject * get_peaks(Cinvertor *self, PyObject *args) {
+	double *pars;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	int count;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	count = npeaks(pars, self->params.d_max, npars, 100);
+
+	return Py_BuildValue("i", count );
+
+}
+
+const char get_positive_doc[] =
+	"Returns the fraction of P(r) that is positive over\n"
+	"the full range of r for the given set of coefficients.\n"
+	" @param args: c-parameters\n"
+	" @return: fraction of P(r) that is positive";
+
+static PyObject * get_positive(Cinvertor *self, PyObject *args) {
+	double *pars;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	double fraction;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	fraction = positive_integral(pars, self->params.d_max, npars, 100);
+
+	return Py_BuildValue("f", fraction );
+
+}
+
+const char get_pos_err_doc[] =
+	"Returns the fraction of P(r) that is 1 standard deviation\n"
+	"above zero over the full range of r for the given set of coefficients.\n"
+	" @param args: c-parameters\n"
+	" @return: fraction of P(r) that is positive";
+
+static PyObject * get_pos_err(Cinvertor *self, PyObject *args) {
+	double *pars;
+	double *pars_err;
+	PyObject *data_obj;
+	PyObject *err_obj;
+	Py_ssize_t npars;
+	Py_ssize_t npars2;
+	double fraction;
+
+	if (!PyArg_ParseTuple(args, "OO", &data_obj, &err_obj)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+	OUTVECTOR(err_obj,pars_err,npars2);
+
+	fraction = positive_errors(pars, pars_err, self->params.d_max, npars, 51);
+
+	return Py_BuildValue("f", fraction );
+
+}
+
+const char get_rg_doc[] =
+	"Returns the value of the radius of gyration Rg.\n"
+	" @param args: c-parameters\n"
+	" @return: Rg";
+
+static PyObject * get_rg(Cinvertor *self, PyObject *args) {
+	double *pars;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	double value;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	value = rg(pars, self->params.d_max, npars, 101);
+
+	return Py_BuildValue("f", value );
+
+}
+
+const char get_iq0_doc[] =
+	"Returns the value of I(q=0).\n"
+	" @param args: c-parameters\n"
+	" @return: I(q=0)";
+
+static PyObject * get_iq0(Cinvertor *self, PyObject *args) {
+	double *pars;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	double value;
+
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+
+	value = 4.0*acos(-1.0)*int_pr(pars, self->params.d_max, npars, 101);
+
+	return Py_BuildValue("f", value );
+
+}
+
+/**
+ * Check whether a q-value is within acceptabel limits
+ * Return 1 if accepted, 0 if rejected.
+ */
+int accept_q(Cinvertor *self, double q) {
+    if (self->params.q_min>0 && q<self->params.q_min) return 0;
+    if (self->params.q_max>0 && q>self->params.q_max) return 0;
+    return 1;
+}
+
+const char get_matrix_doc[] =
+	"Returns A matrix and b vector for least square problem.\n"
+	" @param nfunc: number of base functions\n"
+	" @param nr: number of r-points used when evaluating reg term.\n"
+	" @param a: A array to fill\n"
+	" @param b: b vector to fill\n"
+	" @return: 0";
+
+static PyObject * get_matrix(Cinvertor *self, PyObject *args) {
+	double *a;
+	double *b;
+	PyObject *a_obj;
+	PyObject *b_obj;
+	Py_ssize_t n_a;
+	Py_ssize_t n_b;
+	// Number of bins for regularization term evaluation
+	int nr, nfunc;
+	int i, j, i_r;
+	double r, sqrt_alpha, pi;
+	double tmp;
+	int offset;
+
+	if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &b_obj)) return NULL;
+	OUTVECTOR(a_obj,a,n_a);
+	OUTVECTOR(b_obj,b,n_b);
+
+	assert(n_b>=nfunc);
+	assert(n_a>=nfunc*(nr+self->params.npoints));
+
+	sqrt_alpha = sqrt(self->params.alpha);
+	pi = acos(-1.0);
+	offset = (self->params.has_bck==1) ? 0 : 1;
+
+    for (j=0; j<nfunc; j++) {
+        for (i=0; i<self->params.npoints; i++) {
+            if (self->params.err[i]==0.0) {
+              PyErr_SetString(CinvertorError,
+                "Cinvertor.get_matrix: Some I(Q) points have no error.");
+              return NULL;
+            }
+            if (accept_q(self, self->params.x[i])){
+                if (self->params.has_bck==1 && j==0) {
+                    a[i*nfunc+j] = 1.0/self->params.err[i];
+                } else {
+                	if (self->params.slit_width>0 || self->params.slit_height>0) {
+                		a[i*nfunc+j] = ortho_transformed_smeared(self->params.d_max,
+                				j+offset, self->params.slit_height, self->params.slit_width,
+                				self->params.x[i], 21)/self->params.err[i];
+                	} else {
+                		a[i*nfunc+j] = ortho_transformed(self->params.d_max, j+offset, self->params.x[i])/self->params.err[i];
+                	}
+            	}
+            }
+        }
+        for (i_r=0; i_r<nr; i_r++){
+            if (self->params.has_bck==1 && j==0) {
+                a[(i_r+self->params.npoints)*nfunc+j] = 0.0;
+            } else {
+	            r = self->params.d_max/nr*i_r;
+	            tmp = pi*(j+offset)/self->params.d_max;
+	            a[(i_r+self->params.npoints)*nfunc+j] = sqrt_alpha * 1.0/nr*self->params.d_max*2.0*
+	            (2.0*pi*(j+offset)/self->params.d_max*cos(pi*(j+offset)*r/self->params.d_max) +
+	            tmp*tmp*r * sin(pi*(j+offset)*r/self->params.d_max));
+	        }
+        }
+    }
+
+    for (i=0; i<self->params.npoints; i++) {
+        if (accept_q(self, self->params.x[i])){
+            b[i] = self->params.y[i]/self->params.err[i];
+         }
+	}
+
+	return Py_BuildValue("i", 0);
+
+}
+
+const char get_invcov_matrix_doc[] =
+	" Compute the inverse covariance matrix, defined as inv_cov = a_transposed x a.\n"
+	" @param nfunc: number of base functions\n"
+	" @param nr: number of r-points used when evaluating reg term.\n"
+	" @param a: A array to fill\n"
+	" @param inv_cov: inverse covariance array to be filled\n"
+	" @return: 0";
+
+static PyObject * get_invcov_matrix(Cinvertor *self, PyObject *args) {
+	double *a;
+	PyObject *a_obj;
+	Py_ssize_t n_a;
+	double *inv_cov;
+	PyObject *cov_obj;
+	Py_ssize_t n_cov;
+	int nr, nfunc;
+	int i, j, k;
+
+	if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &cov_obj)) return NULL;
+	OUTVECTOR(a_obj,a,n_a);
+	OUTVECTOR(cov_obj,inv_cov,n_cov);
+
+	assert(n_cov>=nfunc*nfunc);
+	assert(n_a>=nfunc*(nr+self->params.npoints));
+
+	for (i=0; i<nfunc; i++) {
+		for (j=0; j<nfunc; j++) {
+			inv_cov[i*nfunc+j] = 0.0;
+			for (k=0; k<nr+self->params.npoints; k++) {
+				inv_cov[i*nfunc+j] += a[k*nfunc+i]*a[k*nfunc+j];
+			}
+		}
+	}
+	return Py_BuildValue("i", 0);
+}
+
+const char get_reg_size_doc[] =
+	" Compute the covariance matrix, defined as inv_cov = a_transposed x a.\n"
+	" @param nfunc: number of base functions\n"
+	" @param nr: number of r-points used when evaluating reg term.\n"
+	" @param a: A array to fill\n"
+	" @param inv_cov: inverse covariance array to be filled\n"
+	" @return: 0";
+
+static PyObject * get_reg_size(Cinvertor *self, PyObject *args) {
+	double *a;
+	PyObject *a_obj;
+	Py_ssize_t n_a;
+	int nr, nfunc;
+	int i, j;
+	double sum_sig, sum_reg;
+
+	if (!PyArg_ParseTuple(args, "iiO", &nfunc, &nr, &a_obj)) return NULL;
+	OUTVECTOR(a_obj,a,n_a);
+
+	assert(n_a>=nfunc*(nr+self->params.npoints));
+
+    sum_sig = 0.0;
+    sum_reg = 0.0;
+    for (j=0; j<nfunc; j++){
+        for (i=0; i<self->params.npoints; i++){
+            if (accept_q(self, self->params.x[i])==1)
+                sum_sig += (a[i*nfunc+j])*(a[i*nfunc+j]);
+        }
+        for (i=0; i<nr; i++){
+            sum_reg += (a[(i+self->params.npoints)*nfunc+j])*(a[(i+self->params.npoints)*nfunc+j]);
+        }
+    }
+    return Py_BuildValue("ff", sum_sig, sum_reg);
+}
+
+const char eeeget_qmin_doc[] = "\
+This is a multiline doc string.\n\
+\n\
+This is the second line.";
+const char eeeset_qmin_doc[] =
+	"This is a multiline doc string.\n"
+	"\n"
+	"This is the second line.";
+
+static PyMethodDef Cinvertor_methods[] = {
+		   {"residuals", (PyCFunction)residuals, METH_VARARGS, residuals_doc},
+		   {"pr_residuals", (PyCFunction)pr_residuals, METH_VARARGS, pr_residuals_doc},
+		   {"set_x", (PyCFunction)set_x, METH_VARARGS, set_x_doc},
+		   {"get_x", (PyCFunction)get_x, METH_VARARGS, get_x_doc},
+		   {"set_y", (PyCFunction)set_y, METH_VARARGS, set_y_doc},
+		   {"get_y", (PyCFunction)get_y, METH_VARARGS, get_y_doc},
+		   {"set_err", (PyCFunction)set_err, METH_VARARGS, set_err_doc},
+		   {"get_err", (PyCFunction)get_err, METH_VARARGS, get_err_doc},
+		   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, set_dmax_doc},
+		   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, get_dmax_doc},
+		   {"set_qmin", (PyCFunction)set_qmin, METH_VARARGS, set_qmin_doc},
+		   {"get_qmin", (PyCFunction)get_qmin, METH_VARARGS, get_qmin_doc},
+		   {"set_qmax", (PyCFunction)set_qmax, METH_VARARGS, set_qmax_doc},
+		   {"get_qmax", (PyCFunction)get_qmax, METH_VARARGS, get_qmax_doc},
+		   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, set_alpha_doc},
+		   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, get_alpha_doc},
+		   {"set_slit_width", (PyCFunction)set_slit_width, METH_VARARGS, set_slit_width_doc},
+		   {"get_slit_width", (PyCFunction)get_slit_width, METH_VARARGS, get_slit_width_doc},
+		   {"set_slit_height", (PyCFunction)set_slit_height, METH_VARARGS, set_slit_height_doc},
+		   {"get_slit_height", (PyCFunction)get_slit_height, METH_VARARGS, get_slit_height_doc},
+		   {"set_has_bck", (PyCFunction)set_has_bck, METH_VARARGS, set_has_bck_doc},
+		   {"get_has_bck", (PyCFunction)get_has_bck, METH_VARARGS, get_has_bck_doc},
+		   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, get_nx_doc},
+		   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, get_ny_doc},
+		   {"get_nerr", (PyCFunction)get_nerr, METH_VARARGS, get_nerr_doc},
+		   {"iq", (PyCFunction)get_iq, METH_VARARGS, get_iq_doc},
+		   {"iq_smeared", (PyCFunction)get_iq_smeared, METH_VARARGS, get_iq_smeared_doc},
+		   {"pr", (PyCFunction)get_pr, METH_VARARGS, get_pr_doc},
+		   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, get_pr_err_doc},
+		   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, is_valid_doc},
+		   {"basefunc_ft", (PyCFunction)basefunc_ft, METH_VARARGS, basefunc_ft_doc},
+		   {"oscillations", (PyCFunction)oscillations, METH_VARARGS, oscillations_doc},
+		   {"get_peaks", (PyCFunction)get_peaks, METH_VARARGS, get_peaks_doc},
+		   {"get_positive", (PyCFunction)get_positive, METH_VARARGS, get_positive_doc},
+		   {"get_pos_err", (PyCFunction)get_pos_err, METH_VARARGS, get_pos_err_doc},
+		   {"rg", (PyCFunction)get_rg, METH_VARARGS, get_rg_doc},
+		   {"iq0", (PyCFunction)get_iq0, METH_VARARGS, get_iq0_doc},
+		   {"_get_matrix", (PyCFunction)get_matrix, METH_VARARGS, get_matrix_doc},
+		   {"_get_invcov_matrix", (PyCFunction)get_invcov_matrix, METH_VARARGS, get_invcov_matrix_doc},
+		   {"_get_reg_size", (PyCFunction)get_reg_size, METH_VARARGS, get_reg_size_doc},
+
+   {NULL}
+};
+
+static PyTypeObject CinvertorType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "Cinvertor",             /*tp_name*/
+    sizeof(Cinvertor),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)Cinvertor_dealloc, /*tp_dealloc*/
+    0,                         /*tp_print*/
+    0,                         /*tp_getattr*/
+    0,                         /*tp_setattr*/
+    0,                         /*tp_compare*/
+    0,                         /*tp_repr*/
+    0,                         /*tp_as_number*/
+    0,                         /*tp_as_sequence*/
+    0,                         /*tp_as_mapping*/
+    0,                         /*tp_hash */
+    0,                         /*tp_call*/
+    0,                         /*tp_str*/
+    0,                         /*tp_getattro*/
+    0,                         /*tp_setattro*/
+    0,                         /*tp_as_buffer*/
+    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
+    "Cinvertor objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    Cinvertor_methods,             /* tp_methods */
+    Cinvertor_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)Cinvertor_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    Cinvertor_new,                 /* tp_new */
+};
+
+
+static PyMethodDef module_methods[] = {
+    {NULL}
+};
+
+/**
+ * Function used to add the model class to a module
+ * @param module: module to add the class to
+ */
+void addCinvertor(PyObject *module) {
+	PyObject *d;
+
+    if (PyType_Ready(&CinvertorType) < 0)
+        return;
+
+    Py_INCREF(&CinvertorType);
+    PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType);
+
+    d = PyModule_GetDict(module);
+    CinvertorError = PyErr_NewException("sans.pr.invertor.Cinvertor.InvertorError", PyExc_RuntimeError, NULL);
+    PyDict_SetItemString(d, "CinvertorError", CinvertorError);
+}
+
+
+#ifndef PyMODINIT_FUNC	/* declarations for DLL import/export */
+#define PyMODINIT_FUNC void
+#endif
+PyMODINIT_FUNC
+initpr_inversion(void)
+{
+    PyObject* m;
+
+    m = Py_InitModule3("pr_inversion", module_methods,
+                       "C extension module for inversion to P(r).");
+
+    addCinvertor(m);
+}
Index: src/sans/pr/c_extensions/__init__.py
===================================================================
--- src/sans/pr/c_extensions/__init__.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/c_extensions/__init__.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,3 @@
+"""
+    C extensions to provide the P(r) inversion computations.
+"""
Index: src/sans/pr/c_extensions/invertor.c
===================================================================
--- src/sans/pr/c_extensions/invertor.c	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/c_extensions/invertor.c	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,315 @@
+#include <math.h>
+#include "invertor.h"
+#include <memory.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+double pi = 3.1416;
+
+/**
+ * Deallocate memory
+ */
+void invertor_dealloc(Invertor_params *pars) {
+	free(pars->x);
+	free(pars->y);
+	free(pars->err);
+}
+
+void invertor_init(Invertor_params *pars) {
+	pars->d_max = 180;
+	pars->q_min = -1.0;
+	pars->q_max = -1.0;
+	pars->has_bck = 0;
+}
+
+
+/**
+ * P(r) of a sphere, for test purposes
+ *
+ * @param R: radius of the sphere
+ * @param r: distance, in the same units as the radius
+ * @return: P(r)
+ */
+double pr_sphere(double R, double r) {
+    if (r <= 2.0*R) {
+        return 12.0* pow(0.5*r/R, 2.0) * pow(1.0-0.5*r/R, 2.0) * ( 2.0 + 0.5*r/R );
+    } else {
+        return 0.0;
+    }
+}
+
+/**
+ * Orthogonal functions:
+ * B(r) = 2r sin(pi*nr/d)
+ *
+ */
+double ortho(double d_max, int n, double r) {
+	return 2.0*r*sin(pi*n*r/d_max);
+}
+
+/**
+ * Fourier transform of the nth orthogonal function
+ *
+ */
+double ortho_transformed(double d_max, int n, double q) {
+	return 8.0*pow(pi, 2.0)/q * d_max * n * pow(-1.0, n+1)
+		*sin(q*d_max) / ( pow(pi*n, 2.0) - pow(q*d_max, 2.0) );
+}
+
+/**
+ * Slit-smeared Fourier transform of the nth orthogonal function.
+ * Smearing follows Lake, Acta Cryst. (1967) 23, 191.
+ */
+double ortho_transformed_smeared(double d_max, int n, double height, double width, double q, int npts) {
+	double sum, y, z;
+	int i, j, n_height, n_width;
+	double count_w;
+	double fnpts;
+	sum = 0.0;
+	fnpts = (float)npts-1.0;
+
+	// Check for zero slit size
+	n_height = (height>0) ? npts : 1;
+	n_width  = (width>0)  ? npts : 1;
+
+	count_w = 0.0;
+
+	for(j=0; j<n_height; j++) {
+		if(height>0){
+			z = height/fnpts*(float)j;
+		} else {
+			z = 0.0;
+		}
+
+		for(i=0; i<n_width; i++) {
+			if(width>0){
+				y = -width/2.0+width/fnpts*(float)i;
+			} else {
+				y = 0.0;
+			}
+			if (((q-y)*(q-y)+z*z)<=0.0) continue;
+			count_w += 1.0;
+			sum += ortho_transformed(d_max, n, sqrt((q-y)*(q-y)+z*z));
+		}
+	}
+	return sum/count_w;
+}
+
+/**
+ * First derivative in of the orthogonal function dB(r)/dr
+ *
+ */
+double ortho_derived(double d_max, int n, double r) {
+    return 2.0*sin(pi*n*r/d_max) + 2.0*r*cos(pi*n*r/d_max);
+}
+
+/**
+ * Scattering intensity calculated from the expansion.
+ */
+double iq(double *pars, double d_max, int n_c, double q) {
+    double sum = 0.0;
+	int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * ortho_transformed(d_max, i+1, q);
+    }
+    return sum;
+}
+
+/**
+ * Scattering intensity calculated from the expansion,
+ * slit-smeared.
+ */
+double iq_smeared(double *pars, double d_max, int n_c, double height, double width, double q, int npts) {
+    double sum = 0.0;
+	int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * ortho_transformed_smeared(d_max, i+1, height, width, q, npts);
+    }
+    return sum;
+}
+
+/**
+ * P(r) calculated from the expansion.
+ */
+double pr(double *pars, double d_max, int n_c, double r) {
+    double sum = 0.0;
+	int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * ortho(d_max, i+1, r);
+    }
+    return sum;
+}
+
+/**
+ * P(r) calculated from the expansion, with errors
+ */
+void pr_err(double *pars, double *err, double d_max, int n_c,
+		double r, double *pr_value, double *pr_value_err) {
+    double sum = 0.0;
+    double sum_err = 0.0;
+    double func_value;
+	int i;
+    for (i=0; i<n_c; i++) {
+    	func_value = ortho(d_max, i+1, r);
+        sum += pars[i] * func_value;
+        //sum_err += err[i]*err[i]*func_value*func_value;
+        sum_err += err[i*n_c+i]*func_value*func_value;
+    }
+    *pr_value = sum;
+    if (sum_err>0) {
+    	*pr_value_err = sqrt(sum_err);
+    } else {
+    	*pr_value_err = sum;
+    }
+}
+
+/**
+ * dP(r)/dr calculated from the expansion.
+ */
+double dprdr(double *pars, double d_max, int n_c, double r) {
+    double sum = 0.0;
+	int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * 2.0*(sin(pi*(i+1)*r/d_max) + pi*(i+1)*r/d_max * cos(pi*(i+1)*r/d_max));
+    }
+    return sum;
+}
+
+/**
+ * regularization term calculated from the expansion.
+ */
+double reg_term(double *pars, double d_max, int n_c, int nslice) {
+    double sum = 0.0;
+    double r;
+    double deriv;
+	int i;
+    for (i=0; i<nslice; i++) {
+    	r = d_max/(1.0*nslice)*i;
+    	deriv = dprdr(pars, d_max, n_c, r);
+        sum += deriv*deriv;
+    }
+    return sum/(1.0*nslice)*d_max;
+}
+
+/**
+ * regularization term calculated from the expansion.
+ */
+double int_p2(double *pars, double d_max, int n_c, int nslice) {
+    double sum = 0.0;
+    double r;
+    double value;
+	int i;
+    for (i=0; i<nslice; i++) {
+    	r = d_max/(1.0*nslice)*i;
+    	value = pr(pars, d_max, n_c, r);
+        sum += value*value;
+    }
+    return sum/(1.0*nslice)*d_max;
+}
+
+/**
+ * Integral of P(r)
+ */
+double int_pr(double *pars, double d_max, int n_c, int nslice) {
+    double sum = 0.0;
+    double r;
+    double value;
+	int i;
+    for (i=0; i<nslice; i++) {
+    	r = d_max/(1.0*nslice)*i;
+    	value = pr(pars, d_max, n_c, r);
+        sum += value;
+    }
+    return sum/(1.0*nslice)*d_max;
+}
+
+/**
+ * Get the number of P(r) peaks.
+ */
+int npeaks(double *pars, double d_max, int n_c, int nslice) {
+    double r;
+    double value;
+	int i;
+	double previous = 0.0;
+	double slope    = 0.0;
+	int count = 0;
+    for (i=0; i<nslice; i++) {
+    	r = d_max/(1.0*nslice)*i;
+    	value = pr(pars, d_max, n_c, r);
+    	if (previous<=value){
+    		//if (slope<0) count += 1;
+    		slope = 1;
+    	} else {
+    		//printf("slope -1");
+    		if (slope>0) count += 1;
+    		slope = -1;
+    	}
+    	previous = value;
+    }
+    return count;
+}
+
+/**
+ * Get the fraction of the integral of P(r) over the whole range
+ * of r that is above zero.
+ * A valid P(r) is define as being positive for all r.
+ */
+double positive_integral(double *pars, double d_max, int n_c, int nslice) {
+    double r;
+    double value;
+	int i;
+	double sum_pos = 0.0;
+	double sum = 0.0;
+
+    for (i=0; i<nslice; i++) {
+    	r = d_max/(1.0*nslice)*i;
+    	value = pr(pars, d_max, n_c, r);
+    	if (value>0.0) sum_pos += value;
+    	sum += fabs(value);
+    }
+    return sum_pos/sum;
+}
+
+/**
+ * Get the fraction of the integral of P(r) over the whole range
+ * of r that is at least one sigma above zero.
+ */
+double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice) {
+    double r;
+	int i;
+	double sum_pos = 0.0;
+	double sum = 0.0;
+	double pr_val;
+	double pr_val_err;
+
+    for (i=0; i<nslice; i++) {
+    	r = d_max/(1.0*nslice)*i;
+    	pr_err(pars, err, d_max, n_c, r, &pr_val, &pr_val_err);
+    	if (pr_val>pr_val_err) sum_pos += pr_val;
+    	sum += fabs(pr_val);
+
+
+    }
+    return sum_pos/sum;
+}
+
+/**
+ * R_g radius of gyration calculation
+ *
+ * R_g**2 = integral[r**2 * p(r) dr] /  (2.0 * integral[p(r) dr])
+ */
+double rg(double *pars, double d_max, int n_c, int nslice) {
+    double sum_r2 = 0.0;
+    double sum    = 0.0;
+    double r;
+    double value;
+	int i;
+    for (i=0; i<nslice; i++) {
+    	r = d_max/(1.0*nslice)*i;
+    	value = pr(pars, d_max, n_c, r);
+        sum += value;
+        sum_r2 += r*r*value;
+    }
+    return sqrt(sum_r2/(2.0*sum));
+}
+
Index: src/sans/pr/c_extensions/invertor.h
===================================================================
--- src/sans/pr/c_extensions/invertor.h	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/c_extensions/invertor.h	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,59 @@
+#if !defined(invertor_h)
+#define invertor_h
+
+/**
+ * Internal data structure for P(r) inversion
+ */
+typedef struct {
+	/// Maximum distance between any two points in the system
+    double d_max;
+    /// q data
+    double *x;
+    /// I(q) data
+    double *y;
+    /// dI(q) data
+    double *err;
+    /// Number of q points
+    int npoints;
+    /// Number of I(q) points
+    int ny;
+    /// Number of dI(q) points
+    int nerr;
+    /// Alpha value
+    double alpha;
+    /// Minimum q to include in inversion
+    double q_min;
+    /// Maximum q to include in inversion
+    double q_max;
+    /// Flag for whether or not to evalute a constant background while inverting
+    int has_bck;
+    /// Slit height in units of q [A-1]
+    double slit_height;
+    /// Slit width in units of q [A-1]
+    double slit_width;
+} Invertor_params;
+
+void invertor_dealloc(Invertor_params *pars);
+
+void invertor_init(Invertor_params *pars);
+
+double pr_sphere(double R, double r);
+double ortho(double d_max, int n, double r);
+double ortho_transformed(double d_max, int n, double q);
+double ortho_derived(double d_max, int n, double r);
+double iq(double *c, double d_max, int n_c, double q);
+double pr(double *c, double d_max, int n_c, double r);
+double dprdr(double *pars, double d_max, int n_c, double r);
+double reg_term(double *pars, double d_max, int n_c, int nslice);
+double int_p2(double *pars, double d_max, int n_c, int nslice);
+void pr_err(double *pars, double *err, double d_max, int n_c,
+		double r, double *pr_value, double *pr_value_err);
+int npeaks(double *pars, double d_max, int n_c, int nslice);
+double positive_integral(double *pars, double d_max, int n_c, int nslice);
+double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice);
+double rg(double *pars, double d_max, int n_c, int nslice);
+double int_pr(double *pars, double d_max, int n_c, int nslice);
+double ortho_transformed_smeared(double d_max, int n, double heigth, double width, double q, int npts);
+double iq_smeared(double *pars, double d_max, int n_c, double height, double width, double q, int npts);
+
+#endif
Index: src/sans/pr/distance_explorer.py
===================================================================
--- src/sans/pr/distance_explorer.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/distance_explorer.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,107 @@
+
+################################################################################
+#This software was developed by the University of Tennessee as part of the
+#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+#project funded by the US National Science Foundation. 
+#
+#See the license text in license.txt
+#
+#copyright 2009, University of Tennessee
+################################################################################
+
+"""
+Module to explore the P(r) inversion results for a range
+of D_max value. User picks a number of points and a range of
+distances, then get a series of outputs as a function of D_max
+over that range.
+"""
+import sys
+
+
+class Results:
+    """
+    Class to hold the inversion output parameters
+    as a function of D_max
+    """
+    def __init__(self):
+        """
+        Initialization. Create empty arrays
+        and dictionary of labels.
+        """
+        # Array of output for each inversion
+        self.chi2 = []
+        self.osc = []
+        self.pos = []
+        self.pos_err = []
+        self.rg = []
+        self.iq0 = []
+        self.bck = []
+        self.d_max = []
+        ## List of errors found during the last exploration
+        self.errors = []
+        
+        
+class DistExplorer(object):
+    """
+    The explorer class
+    """
+    
+    def __init__(self, pr_state):
+        """
+        Initialization.
+        
+        :param pr_state: sans.pr.invertor.Invertor object
+        
+        """
+        self.pr_state = pr_state
+        self._default_min = 0.8 * self.pr_state.d_max
+        self._default_max = 1.2 * self.pr_state.d_max
+        
+    def __call__(self, dmin=None, dmax=None, npts=10):
+        """
+        Compute the outputs as a function of D_max.
+        
+        :param dmin: minimum value for D_max
+        :param dmax: maximum value for D_max
+        :param npts: number of points for D_max
+        
+        """
+        # Take care of the defaults if needed
+        if dmin is None:
+            dmin = self._default_min
+            
+        if dmax is None:
+            dmax = self._default_max
+            
+        # Results object to store the computation outputs.
+        results = Results()
+        
+        # Loop over d_max values
+        for i in range(npts):
+            d = dmin + i * (dmax - dmin) / (npts - 1.0)
+            self.pr_state.d_max = d
+            try:
+                out, cov = self.pr_state.invert(self.pr_state.nfunc)
+            
+                # Store results
+                iq0 = self.pr_state.iq0(out)
+                rg = self.pr_state.rg(out)
+                pos = self.pr_state.get_positive(out)
+                pos_err = self.pr_state.get_pos_err(out, cov)
+                osc = self.pr_state.oscillations(out)
+            
+                results.d_max.append(self.pr_state.d_max)
+                results.bck.append(self.pr_state.background)
+                results.chi2.append(self.pr_state.chi2)
+                results.iq0.append(iq0)
+                results.rg.append(rg)
+                results.pos.append(pos)
+                results.pos_err.append(pos_err)
+                results.osc.append(osc)
+            except:
+                # This inversion failed, skip this D_max value
+                msg = "ExploreDialog: inversion failed for "
+                msg += "D_max=%s\n %s" % (str(d), sys.exc_value)
+                results.errors.append(msg)
+            
+        return results
Index: src/sans/pr/invertor.py
===================================================================
--- src/sans/pr/invertor.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/invertor.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,739 @@
+"""
+Module to perform P(r) inversion.
+The module contains the Invertor class.
+"""
+
+import numpy
+import sys
+import math
+import time
+import copy
+import os
+import re
+from numpy.linalg import lstsq
+from scipy import optimize
+from sans.pr.core.pr_inversion import Cinvertor
+
+def help():
+    """
+    Provide general online help text
+    Future work: extend this function to allow topic selection
+    """
+    info_txt  = "The inversion approach is based on Moore, J. Appl. Cryst. "
+    info_txt += "(1980) 13, 168-175.\n\n"
+    info_txt += "P(r) is set to be equal to an expansion of base functions "
+    info_txt += "of the type "
+    info_txt += "phi_n(r) = 2*r*sin(pi*n*r/D_max). The coefficient of each "
+    info_txt += "base functions "
+    info_txt += "in the expansion is found by performing a least square fit "
+    info_txt += "with the "
+    info_txt += "following fit function:\n\n"
+    info_txt += "chi**2 = sum_i[ I_meas(q_i) - I_th(q_i) ]**2/error**2 +"
+    info_txt += "Reg_term\n\n"
+    info_txt += "where I_meas(q) is the measured scattering intensity and "
+    info_txt += "I_th(q) is "
+    info_txt += "the prediction from the Fourier transform of the P(r) "
+    info_txt += "expansion. "
+    info_txt += "The Reg_term term is a regularization term set to the second"
+    info_txt += " derivative "
+    info_txt += "d**2P(r)/dr**2 integrated over r. It is used to produce "
+    info_txt += "a smooth P(r) output.\n\n"
+    info_txt += "The following are user inputs:\n\n"
+    info_txt += "   - Number of terms: the number of base functions in the P(r)"
+    info_txt += " expansion.\n\n"
+    info_txt += "   - Regularization constant: a multiplicative constant "
+    info_txt += "to set the size of "
+    info_txt += "the regularization term.\n\n"
+    info_txt += "   - Maximum distance: the maximum distance between any "
+    info_txt += "two points in the system.\n"
+     
+    return info_txt
+    
+
+class Invertor(Cinvertor):
+    """
+    Invertor class to perform P(r) inversion
+    
+    The problem is solved by posing the problem as  Ax = b,
+    where x is the set of coefficients we are looking for.
+    
+    Npts is the number of points.
+    
+    In the following i refers to the ith base function coefficient.
+    The matrix has its entries j in its first Npts rows set to
+        A[j][i] = (Fourier transformed base function for point j)
+        
+    We them choose a number of r-points, n_r, to evaluate the second
+    derivative of P(r) at. This is used as our regularization term.
+    For a vector r of length n_r, the following n_r rows are set to
+        A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2,
+        evaluated at r[j])
+        
+    The vector b has its first Npts entries set to
+        b[j] = (I(q) observed for point j)
+        
+    The following n_r entries are set to zero.
+    
+    The result is found by using scipy.linalg.basic.lstsq to invert
+    the matrix and find the coefficients x.
+    
+    Methods inherited from Cinvertor:
+    - get_peaks(pars): returns the number of P(r) peaks
+    - oscillations(pars): returns the oscillation parameters for the output P(r)
+    - get_positive(pars): returns the fraction of P(r) that is above zero
+    - get_pos_err(pars): returns the fraction of P(r) that is 1-sigma above zero
+    """
+    ## Chisqr of the last computation
+    chi2  = 0
+    ## Time elapsed for last computation
+    elapsed = 0
+    ## Alpha to get the reg term the same size as the signal
+    suggested_alpha = 0
+    ## Last number of base functions used
+    nfunc = 10
+    ## Last output values
+    out = None
+    ## Last errors on output values
+    cov = None
+    ## Background value
+    background = 0
+    ## Information dictionary for application use
+    info = {}
+    
+    def __init__(self):
+        Cinvertor.__init__(self)
+        
+    def __setstate__(self, state):
+        """
+        restore the state of invertor for pickle
+        """
+        (self.__dict__, self.alpha, self.d_max,
+         self.q_min, self.q_max,
+         self.x, self.y,
+         self.err, self.has_bck,
+         self.slit_height, self.slit_width) = state
+   
+    def __reduce_ex__(self, proto):
+        """
+        Overwrite the __reduce_ex__
+        """
+
+        state = (self.__dict__,
+                 self.alpha, self.d_max,
+                 self.q_min, self.q_max,
+                 self.x, self.y,
+                 self.err, self.has_bck,
+                 self.slit_height, self.slit_width,
+                 )
+        return (Invertor, tuple(), state, None, None)
+    
+    def __setattr__(self, name, value):
+        """
+        Set the value of an attribute.
+        Access the parent class methods for
+        x, y, err, d_max, q_min, q_max and alpha
+        """
+        if   name == 'x':
+            if 0.0 in value:
+                msg = "Invertor: one of your q-values is zero. "
+                msg += "Delete that entry before proceeding"
+                raise ValueError, msg
+            return self.set_x(value)
+        elif name == 'y':
+            return self.set_y(value)
+        elif name == 'err':
+            value2 = abs(value)
+            return self.set_err(value2)
+        elif name == 'd_max':
+            return self.set_dmax(value)
+        elif name == 'q_min':
+            if value == None:
+                return self.set_qmin(-1.0)
+            return self.set_qmin(value)
+        elif name == 'q_max':
+            if value == None:
+                return self.set_qmax(-1.0)
+            return self.set_qmax(value)
+        elif name == 'alpha':
+            return self.set_alpha(value)
+        elif name == 'slit_height':
+            return self.set_slit_height(value)
+        elif name == 'slit_width':
+            return self.set_slit_width(value)
+        elif name == 'has_bck':
+            if value == True:
+                return self.set_has_bck(1)
+            elif value == False:
+                return self.set_has_bck(0)
+            else:
+                raise ValueError, "Invertor: has_bck can only be True or False"
+            
+        return Cinvertor.__setattr__(self, name, value)
+    
+    def __getattr__(self, name):
+        """
+        Return the value of an attribute
+        """
+        #import numpy
+        if name == 'x':
+            out = numpy.ones(self.get_nx())
+            self.get_x(out)
+            return out
+        elif name == 'y':
+            out = numpy.ones(self.get_ny())
+            self.get_y(out)
+            return out
+        elif name == 'err':
+            out = numpy.ones(self.get_nerr())
+            self.get_err(out)
+            return out
+        elif name == 'd_max':
+            return self.get_dmax()
+        elif name == 'q_min':
+            qmin = self.get_qmin()
+            if qmin < 0:
+                return None
+            return qmin
+        elif name == 'q_max':
+            qmax = self.get_qmax()
+            if qmax < 0:
+                return None
+            return qmax
+        elif name == 'alpha':
+            return self.get_alpha()
+        elif name == 'slit_height':
+            return self.get_slit_height()
+        elif name == 'slit_width':
+            return self.get_slit_width()
+        elif name == 'has_bck':
+            value = self.get_has_bck()
+            if value == 1:
+                return True
+            else:
+                return False
+        elif name in self.__dict__:
+            return self.__dict__[name]
+        return None
+    
+    def clone(self):
+        """
+        Return a clone of this instance
+        """
+        #import copy
+        
+        invertor = Invertor()
+        invertor.chi2    = self.chi2
+        invertor.elapsed = self.elapsed
+        invertor.nfunc   = self.nfunc
+        invertor.alpha   = self.alpha
+        invertor.d_max   = self.d_max
+        invertor.q_min   = self.q_min
+        invertor.q_max   = self.q_max
+        
+        invertor.x = self.x
+        invertor.y = self.y
+        invertor.err = self.err
+        invertor.has_bck = self.has_bck
+        invertor.slit_height = self.slit_height
+        invertor.slit_width = self.slit_width
+        
+        invertor.info = copy.deepcopy(self.info)
+        
+        return invertor
+    
+    def invert(self, nfunc=10, nr=20):
+        """
+        Perform inversion to P(r)
+        
+        The problem is solved by posing the problem as  Ax = b,
+        where x is the set of coefficients we are looking for.
+        
+        Npts is the number of points.
+        
+        In the following i refers to the ith base function coefficient.
+        The matrix has its entries j in its first Npts rows set to
+            A[i][j] = (Fourier transformed base function for point j)
+            
+        We them choose a number of r-points, n_r, to evaluate the second
+        derivative of P(r) at. This is used as our regularization term.
+        For a vector r of length n_r, the following n_r rows are set to
+            A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j])
+            
+        The vector b has its first Npts entries set to
+            b[j] = (I(q) observed for point j)
+            
+        The following n_r entries are set to zero.
+        
+        The result is found by using scipy.linalg.basic.lstsq to invert
+        the matrix and find the coefficients x.
+        
+        :param nfunc: number of base functions to use.
+        :param nr: number of r points to evaluate the 2nd derivative at for the reg. term.
+        :return: c_out, c_cov - the coefficients with covariance matrix
+        
+        """
+        # Reset the background value before proceeding
+        self.background = 0.0
+        return self.lstsq(nfunc, nr=nr)
+    
+    def iq(self, out, q):
+        """
+        Function to call to evaluate the scattering intensity
+        
+        :param args: c-parameters, and q
+        :return: I(q)
+        
+        """
+        return Cinvertor.iq(self, out, q) + self.background
+    
+    def invert_optimize(self, nfunc=10, nr=20):
+        """
+        Slower version of the P(r) inversion that uses scipy.optimize.leastsq.
+        
+        This probably produce more reliable results, but is much slower.
+        The minimization function is set to
+        sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term,
+        where the reg_term is given by Svergun: it is the integral of
+        the square of the first derivative
+        of P(r), d(P(r))/dr, integrated over the full range of r.
+        
+        :param nfunc: number of base functions to use.
+        :param nr: number of r points to evaluate the 2nd derivative at
+            for the reg. term.
+        
+        :return: c_out, c_cov - the coefficients with covariance matrix
+        
+        """
+        self.nfunc = nfunc
+        # First, check that the current data is valid
+        if self.is_valid() <= 0:
+            msg = "Invertor.invert: Data array are of different length"
+            raise RuntimeError, msg
+        
+        p = numpy.ones(nfunc)
+        t_0 = time.time()
+        out, cov_x, _, _, _ = optimize.leastsq(self.residuals,
+                                                            p, full_output=1)
+        
+        # Compute chi^2
+        res = self.residuals(out)
+        chisqr = 0
+        for i in range(len(res)):
+            chisqr += res[i]
+        
+        self.chi2 = chisqr
+
+        # Store computation time
+        self.elapsed = time.time() - t_0
+        
+        if cov_x is None:
+            cov_x = numpy.ones([nfunc, nfunc])
+            cov_x *= math.fabs(chisqr)
+        return out, cov_x
+    
+    def pr_fit(self, nfunc=5):
+        """
+        This is a direct fit to a given P(r). It assumes that the y data
+        is set to some P(r) distribution that we are trying to reproduce
+        with a set of base functions.
+        
+        This method is provided as a test.
+        """
+        # First, check that the current data is valid
+        if self.is_valid() <= 0:
+            msg = "Invertor.invert: Data arrays are of different length"
+            raise RuntimeError, msg
+        
+        p = numpy.ones(nfunc)
+        t_0 = time.time()
+        out, cov_x, _, _, _ = optimize.leastsq(self.pr_residuals, p,
+                                                            full_output=1)
+        
+        # Compute chi^2
+        res = self.pr_residuals(out)
+        chisqr = 0
+        for i in range(len(res)):
+            chisqr += res[i]
+        
+        self.chisqr = chisqr
+        
+        # Store computation time
+        self.elapsed = time.time() - t_0
+
+        return out, cov_x
+    
+    def pr_err(self, c, c_cov, r):
+        """
+        Returns the value of P(r) for a given r, and base function
+        coefficients, with error.
+        
+        :param c: base function coefficients
+        :param c_cov: covariance matrice of the base function coefficients
+        :param r: r-value to evaluate P(r) at
+        
+        :return: P(r)
+        
+        """
+        return self.get_pr_err(c, c_cov, r)
+       
+    def _accept_q(self, q):
+        """
+        Check q-value against user-defined range
+        """
+        if not self.q_min == None and q < self.q_min:
+            return False
+        if not self.q_max == None and q > self.q_max:
+            return False
+        return True
+       
+    def lstsq(self, nfunc=5, nr=20):
+        """
+        The problem is solved by posing the problem as  Ax = b,
+        where x is the set of coefficients we are looking for.
+        
+        Npts is the number of points.
+        
+        In the following i refers to the ith base function coefficient.
+        The matrix has its entries j in its first Npts rows set to
+            A[i][j] = (Fourier transformed base function for point j)
+            
+        We them choose a number of r-points, n_r, to evaluate the second
+        derivative of P(r) at. This is used as our regularization term.
+        For a vector r of length n_r, the following n_r rows are set to
+            A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2,
+            evaluated at r[j])
+            
+        The vector b has its first Npts entries set to
+            b[j] = (I(q) observed for point j)
+            
+        The following n_r entries are set to zero.
+        
+        The result is found by using scipy.linalg.basic.lstsq to invert
+        the matrix and find the coefficients x.
+        
+        :param nfunc: number of base functions to use.
+        :param nr: number of r points to evaluate the 2nd derivative at
+            for the reg. term.
+
+        If the result does not allow us to compute the covariance matrix,
+        a matrix filled with zeros will be returned.
+
+        """
+        # Note: To make sure an array is contiguous:
+        # blah = numpy.ascontiguousarray(blah_original)
+        # ... before passing it to C
+        
+        if self.is_valid() < 0:
+            msg = "Invertor: invalid data; incompatible data lengths."
+            raise RuntimeError, msg
+        
+        self.nfunc = nfunc
+        # a -- An M x N matrix.
+        # b -- An M x nrhs matrix or M vector.
+        npts = len(self.x)
+        nq   = nr
+        sqrt_alpha = math.sqrt(math.fabs(self.alpha))
+        if sqrt_alpha < 0.0:
+            nq = 0
+
+        # If we need to fit the background, add a term
+        if self.has_bck == True:
+            nfunc_0 = nfunc
+            nfunc += 1
+
+        a = numpy.zeros([npts + nq, nfunc])
+        b = numpy.zeros(npts + nq)
+        err = numpy.zeros([nfunc, nfunc])
+        
+        # Construct the a matrix and b vector that represent the problem
+        t_0 = time.time()
+        try:
+            self._get_matrix(nfunc, nq, a, b)
+        except:
+            raise RuntimeError, "Invertor: could not invert I(Q)\n  %s" % sys.exc_value
+             
+        # Perform the inversion (least square fit)
+        c, chi2, _, _ = lstsq(a, b)
+        # Sanity check
+        try:
+            float(chi2)
+        except:
+            chi2 = -1.0
+        self.chi2 = chi2
+                
+        inv_cov = numpy.zeros([nfunc, nfunc])
+        # Get the covariance matrix, defined as inv_cov = a_transposed * a
+        self._get_invcov_matrix(nfunc, nr, a, inv_cov)
+                    
+        # Compute the reg term size for the output
+        sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a)
+                    
+        if math.fabs(self.alpha) > 0:
+            new_alpha = sum_sig / (sum_reg / self.alpha)
+        else:
+            new_alpha = 0.0
+        self.suggested_alpha = new_alpha
+        
+        try:
+            cov = numpy.linalg.pinv(inv_cov)
+            err = math.fabs(chi2 / float(npts - nfunc)) * cov
+        except:
+            # We were not able to estimate the errors
+            # Return an empty error matrix
+            pass
+            
+        # Keep a copy of the last output
+        if self.has_bck == False:
+            self.background = 0
+            self.out = c
+            self.cov = err
+        else:
+            self.background = c[0]
+            
+            err_0 = numpy.zeros([nfunc, nfunc])
+            c_0 = numpy.zeros(nfunc)
+            
+            for i in range(nfunc_0):
+                c_0[i] = c[i+1]
+                for j in range(nfunc_0):
+                    err_0[i][j] = err[i+1][j+1]
+                    
+            self.out = c_0
+            self.cov = err_0
+            
+        # Store computation time
+        self.elapsed = time.time() - t_0
+        
+        return self.out, self.cov
+        
+    def estimate_numterms(self, isquit_func=None):
+        """
+        Returns a reasonable guess for the
+        number of terms
+        
+        :param isquit_func: reference to thread function to call to
+                            check whether the computation needs to
+                            be stopped.
+        
+        :return: number of terms, alpha, message
+        
+        """
+        from num_term import Num_terms
+        estimator = Num_terms(self.clone())
+        try:
+            return estimator.num_terms(isquit_func)
+        except:
+            # If we fail, estimate alpha and return the default
+            # number of terms
+            best_alpha, _, _ = self.estimate_alpha(self.nfunc)
+            return self.nfunc, best_alpha, "Could not estimate number of terms"
+                    
+    def estimate_alpha(self, nfunc):
+        """
+        Returns a reasonable guess for the
+        regularization constant alpha
+        
+        :param nfunc: number of terms to use in the expansion.
+        
+        :return: alpha, message, elapsed
+        
+        where alpha is the estimate for alpha,
+        message is a message for the user,
+        elapsed is the computation time
+        """
+        #import time
+        try:
+            pr = self.clone()
+            
+            # T_0 for computation time
+            starttime = time.time()
+            elapsed = 0
+            
+            # If the current alpha is zero, try
+            # another value
+            if pr.alpha <= 0:
+                pr.alpha = 0.0001
+                 
+            # Perform inversion to find the largest alpha
+            out, _ = pr.invert(nfunc)
+            elapsed = time.time() - starttime
+            initial_alpha = pr.alpha
+            initial_peaks = pr.get_peaks(out)
+    
+            # Try the inversion with the estimated alpha
+            pr.alpha = pr.suggested_alpha
+            out, _ = pr.invert(nfunc)
+    
+            npeaks = pr.get_peaks(out)
+            # if more than one peak to start with
+            # just return the estimate
+            if npeaks > 1:
+                #message = "Your P(r) is not smooth,
+                #please check your inversion parameters"
+                message = None
+                return pr.suggested_alpha, message, elapsed
+            else:
+                
+                # Look at smaller values
+                # We assume that for the suggested alpha, we have 1 peak
+                # if not, send a message to change parameters
+                alpha = pr.suggested_alpha
+                best_alpha = pr.suggested_alpha
+                found = False
+                for i in range(10):
+                    pr.alpha = (0.33)**(i+1) * alpha
+                    out, _ = pr.invert(nfunc)
+                    
+                    peaks = pr.get_peaks(out)
+                    if peaks > 1:
+                        found = True
+                        break
+                    best_alpha = pr.alpha
+                    
+                # If we didn't find a turning point for alpha and
+                # the initial alpha already had only one peak,
+                # just return that
+                if not found and initial_peaks == 1 and \
+                    initial_alpha < best_alpha:
+                    best_alpha = initial_alpha
+                    
+                # Check whether the size makes sense
+                message = ''
+                
+                if not found:
+                    message = None
+                elif best_alpha >= 0.5 * pr.suggested_alpha:
+                    # best alpha is too big, return a
+                    # reasonable value
+                    message  = "The estimated alpha for your system is too "
+                    message += "large. "
+                    message += "Try increasing your maximum distance."
+                
+                return best_alpha, message, elapsed
+    
+        except:
+            message = "Invertor.estimate_alpha: %s" % sys.exc_value
+            return 0, message, elapsed
+    
+    def to_file(self, path, npts=100):
+        """
+        Save the state to a file that will be readable
+        by SliceView.
+        
+        :param path: path of the file to write
+        :param npts: number of P(r) points to be written
+        
+        """
+        file = open(path, 'w')
+        file.write("#d_max=%g\n" % self.d_max)
+        file.write("#nfunc=%g\n" % self.nfunc)
+        file.write("#alpha=%g\n" % self.alpha)
+        file.write("#chi2=%g\n" % self.chi2)
+        file.write("#elapsed=%g\n" % self.elapsed)
+        file.write("#qmin=%s\n" % str(self.q_min))
+        file.write("#qmax=%s\n" % str(self.q_max))
+        file.write("#slit_height=%g\n" % self.slit_height)
+        file.write("#slit_width=%g\n" % self.slit_width)
+        file.write("#background=%g\n" % self.background)
+        if self.has_bck == True:
+            file.write("#has_bck=1\n")
+        else:
+            file.write("#has_bck=0\n")
+        file.write("#alpha_estimate=%g\n" % self.suggested_alpha)
+        if not self.out == None:
+            if len(self.out) == len(self.cov):
+                for i in range(len(self.out)):
+                    file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]),
+                                                    str(self.cov[i][i])))
+        file.write("<r>  <Pr>  <dPr>\n")
+        r = numpy.arange(0.0, self.d_max, self.d_max/npts)
+        
+        for r_i in r:
+            (value, err) = self.pr_err(self.out, self.cov, r_i)
+            file.write("%g  %g  %g\n" % (r_i, value, err))
+    
+        file.close()
+     
+    def from_file(self, path):
+        """
+        Load the state of the Invertor from a file,
+        to be able to generate P(r) from a set of
+        parameters.
+        
+        :param path: path of the file to load
+        
+        """
+        #import os
+        #import re
+        if os.path.isfile(path):
+            try:
+                fd = open(path, 'r')
+                
+                buff = fd.read()
+                lines = buff.split('\n')
+                for line in lines:
+                    if line.startswith('#d_max='):
+                        toks = line.split('=')
+                        self.d_max = float(toks[1])
+                    elif line.startswith('#nfunc='):
+                        toks = line.split('=')
+                        self.nfunc = int(toks[1])
+                        self.out = numpy.zeros(self.nfunc)
+                        self.cov = numpy.zeros([self.nfunc, self.nfunc])
+                    elif line.startswith('#alpha='):
+                        toks = line.split('=')
+                        self.alpha = float(toks[1])
+                    elif line.startswith('#chi2='):
+                        toks = line.split('=')
+                        self.chi2 = float(toks[1])
+                    elif line.startswith('#elapsed='):
+                        toks = line.split('=')
+                        self.elapsed = float(toks[1])
+                    elif line.startswith('#alpha_estimate='):
+                        toks = line.split('=')
+                        self.suggested_alpha = float(toks[1])
+                    elif line.startswith('#qmin='):
+                        toks = line.split('=')
+                        try:
+                            self.q_min = float(toks[1])
+                        except:
+                            self.q_min = None
+                    elif line.startswith('#qmax='):
+                        toks = line.split('=')
+                        try:
+                            self.q_max = float(toks[1])
+                        except:
+                            self.q_max = None
+                    elif line.startswith('#slit_height='):
+                        toks = line.split('=')
+                        self.slit_height = float(toks[1])
+                    elif line.startswith('#slit_width='):
+                        toks = line.split('=')
+                        self.slit_width = float(toks[1])
+                    elif line.startswith('#background='):
+                        toks = line.split('=')
+                        self.background = float(toks[1])
+                    elif line.startswith('#has_bck='):
+                        toks = line.split('=')
+                        if int(toks[1]) == 1:
+                            self.has_bck = True
+                        else:
+                            self.has_bck = False
+            
+                    # Now read in the parameters
+                    elif line.startswith('#C_'):
+                        toks = line.split('=')
+                        p = re.compile('#C_([0-9]+)')
+                        m = p.search(toks[0])
+                        toks2 = toks[1].split('+-')
+                        i = int(m.group(1))
+                        self.out[i] = float(toks2[0])
+                        
+                        self.cov[i][i] = float(toks2[1])
+            
+            except:
+                msg = "Invertor.from_file: corrupted file\n%s" % sys.exc_value
+                raise RuntimeError, msg
+        else:
+            msg = "Invertor.from_file: '%s' is not a file" % str(path)
+            raise RuntimeError, msg
Index: src/sans/pr/num_term.py
===================================================================
--- src/sans/pr/num_term.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/num_term.py	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,204 @@
+import math
+import numpy
+import copy
+from sans.pr.invertor import Invertor
+
+
+class Num_terms():
+    """
+    """
+    def __init__(self, invertor):
+        """
+        """
+        self.invertor = invertor
+        self.nterm_min = 10
+        self.nterm_max = len(self.invertor.x)
+        if self.nterm_max > 50:
+            self.nterm_max = 50
+        self.isquit_func = None
+         
+        self.osc_list = []
+        self.err_list = []
+        self.alpha_list = []
+        self.mess_list = []
+         
+        self.dataset = []
+     
+    def is_odd(self, n):
+        """
+        """
+        return bool(n % 2)
+
+    def sort_osc(self):
+        """
+        """
+        #import copy
+        osc = copy.deepcopy(self.dataset)
+        lis = []
+        for i in range(len(osc)):
+            osc.sort()
+            re = osc.pop(0)
+            lis.append(re)
+        return lis
+           
+    def median_osc(self):
+        """
+        """
+        osc = self.sort_osc()
+        dv = len(osc)
+        med = float(dv) / 2.0
+        odd = self.is_odd(dv)
+        medi = 0
+        for i in range(dv):
+            if odd == True:
+                medi = osc[int(med)]
+            else:
+                medi = osc[int(med) - 1]
+        return medi
+
+    def get0_out(self):
+        """
+        """
+        inver = self.invertor
+        self.osc_list = []
+        self.err_list = []
+        self.alpha_list = []
+        for k in range(self.nterm_min, self.nterm_max, 1):
+            if self.isquit_func != None:
+                self.isquit_func()
+            best_alpha, message, _ = inver.estimate_alpha(k)
+            inver.alpha = best_alpha
+            inver.out, inver.cov = inver.lstsq(k)
+            osc = inver.oscillations(inver.out)
+            err = inver.get_pos_err(inver.out, inver.cov)
+            if osc > 10.0:
+                break
+            self.osc_list.append(osc)
+            self.err_list.append(err)
+            self.alpha_list.append(inver.alpha)
+            self.mess_list.append(message)
+         
+        new_osc1 = []
+        new_osc2 = []
+        new_osc3 = []
+        flag9 = False
+        flag8 = False
+        flag7 = False
+        for i in range(len(self.err_list)):
+            if self.err_list[i] <= 1.0 and self.err_list[i] >= 0.9:
+                new_osc1.append(self.osc_list[i])
+                flag9 = True
+            if self.err_list[i] < 0.9 and self.err_list[i] >= 0.8:
+                new_osc2.append(self.osc_list[i])
+                flag8 = True
+            if self.err_list[i] < 0.8 and self.err_list[i] >= 0.7:
+                new_osc3.append(self.osc_list[i])
+                flag7 = True
+                 
+        if flag9 == True:
+            self.dataset = new_osc1
+        elif flag8 == True:
+            self.dataset = new_osc2
+        else:
+            self.dataset = new_osc3
+         
+        return self.dataset
+        
+    def ls_osc(self):
+        """
+        """
+        # Generate data
+        ls_osc = self.get0_out()
+        med = self.median_osc()
+        
+        #TODO: check 1
+        ls_osc = self.dataset
+        ls = []
+        for i in range(len(ls_osc)):
+            if int(med) == int(ls_osc[i]):
+                ls.append(ls_osc[i])
+        return ls
+
+    def compare_err(self):
+        """
+        """
+        ls = self.ls_osc()
+        #print "ls", ls
+        nt_ls = []
+        for i in range(len(ls)):
+            r = ls[i]
+            n = self.osc_list.index(r) + 10
+            #er = self.err_list[n]
+            #nt = self.osc_list.index(r) + 10
+            nt_ls.append(n)
+        #print "nt list", nt_ls
+        return nt_ls
+
+    def num_terms(self, isquit_func=None):
+        """
+        """
+        try:
+            self.isquit_func = isquit_func
+            nts = self.compare_err()
+            div = len(nts)
+            tem = float(div)/2.0
+            odd = self.is_odd(div)
+            if odd == True:
+                nt = nts[int(tem)]
+            else:
+                nt = nts[int(tem) - 1]
+            return nt, self.alpha_list[nt - 10], self.mess_list[nt-10]
+        except:
+            return self.nterm_min, self.alpha_list[10], self.mess_list[10]
+
+
+#For testing
+def load(path):
+    # Read the data from the data file
+    data_x   = numpy.zeros(0)
+    data_y   = numpy.zeros(0)
+    data_err = numpy.zeros(0)
+    scale    = None
+    min_err  = 0.0
+    if not path == None:
+        input_f = open(path,'r')
+        buff    = input_f.read()
+        lines   = buff.split('\n')
+        for line in lines:
+            try:
+                toks = line.split()
+                test_x = float(toks[0])
+                test_y = float(toks[1])
+                if len(toks) > 2:
+                    err = float(toks[2])
+                else:
+                    if scale == None:
+                        scale = 0.05 * math.sqrt(test_y)
+                        #scale = 0.05/math.sqrt(y)
+                        min_err = 0.01 * y
+                    err = scale * math.sqrt(test_y) + min_err
+                    #err = 0
+                    
+                data_x = numpy.append(data_x, test_x)
+                data_y = numpy.append(data_y, test_y)
+                data_err = numpy.append(data_err, err)
+            except:
+                pass
+               
+    return data_x, data_y, data_err
+
+
+if __name__ == "__main__":
+    i = Invertor()
+    x, y, erro = load("test/Cyl_A_D102.txt")
+    i.d_max = 102.0
+    i.nfunc = 10
+    #i.q_max = 0.4
+    #i.q_min = 0.07
+    i.x = x
+    i.y = y
+    i.err = erro
+    #i.out, i.cov = i.lstsq(10)
+    # Testing estimator
+    est = Num_terms(i)
+    print est.num_terms()
Index: src/sans/pr/release_notes.txt
===================================================================
--- src/sans/pr/release_notes.txt	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
+++ src/sans/pr/release_notes.txt	(revision 5777106e8a4a0a62fe2ee0fbaf7c0aa19c0bb92b)
@@ -0,0 +1,64 @@
+Release Notes
+=============
+
+Package name: pr_inversion 0.2
+
+1- 	Version 0.2.2
+	- Release date: 7/4/2009
+	- Minor fixes
+	  
+ 	Version 0.2.1
+	- Release date: 7/8/2008
+	- Minor change to clone method
+	  
+	Version 0.2.0
+	
+	- Release date: 7/4/2008
+	- Number of term estimator
+	- Faster implementation
+
+	Version 0.1.0
+
+	- Release date: 5/28/2008
+	- Provide python module sans.pr
+
+2- Downloading and Installing
+
+	2.1- System Requirements:
+		- Python version >= 2.4 should be running on the system
+
+	2.2- Installing:
+		- Get the code from svn://danse.us/sans/releases/pr_inversion-0.1.0
+		- The following modules are required:
+			* matplotlib
+			* numpy
+			* scipy
+
+3- Known Issues
+
+	3.1- All systems:
+		- None
+
+	3.2- Windows:
+		- None
+
+	3.3- Linux:
+		- None
+
+4- Troubleshooting
+
+	- None
+
+5- Frequently Asked Questions
+
+	- None
+
+6- Other Resources
+
+	- See: http://danse.chem.utk.edu/prview.html
+
+	
+
+	
+
+
