Index: pr_inversion/c_extensions/Cinvertor.c
===================================================================
--- pr_inversion/c_extensions/Cinvertor.c	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/c_extensions/Cinvertor.c	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
@@ -222,4 +222,19 @@
 }
 
+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);	
+}
+
+/**
+ * Gets the maximum distance
+ */
+static PyObject * get_alpha(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.alpha);	
+}
+
 /**
  * Gets the number of x points
@@ -246,6 +261,6 @@
 /**
  * Function to call to evaluate the residuals
- * @param args: input q or [q,phi]
- * @return: function value
+ * @param args: input parameters
+ * @return: list of residuals
  */
 static PyObject * residuals(Cinvertor *self, PyObject *args) {
@@ -257,5 +272,6 @@
 	double residual, diff;
 	// Regularization factor
-	double lmult = 0.0;
+	double regterm = 0.0;
+	double tmp = 0.0;
 	
 	PyObject *data_obj;
@@ -266,14 +282,18 @@
 	OUTVECTOR(data_obj,pars,npars);
 		
-	//printf("npars: %i", 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);
     
     for(i=0; i<self->params.npoints; i++) {
-    	//printf("\n%i: ", 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]);
-    	//printf("%i %g\n", i, residual);
+    	tmp = residual;
+    	
+    	// regularization term
+    	residual += self->params.alpha * regterm;
+    	
     	if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
     	    PyErr_SetString(CinvertorError, 
@@ -286,4 +306,52 @@
 	return 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;
+	PyObject* temp;
+	double *res;
+	int i;
+	double residual, diff;
+	// Regularization factor
+	double regterm = 0.0;
+	double tmp = 0.0;
+	
+	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);
+
+    
+    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;
+}
 
 /**
@@ -323,7 +391,31 @@
 }
 
+/**
+ * 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);
+	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);	
+}
+
 
 static PyMethodDef Cinvertor_methods[] = {
 		   {"residuals", (PyCFunction)residuals, METH_VARARGS, "Get the list of residuals"},
+		   {"pr_residuals", (PyCFunction)pr_residuals, METH_VARARGS, "Get the list of residuals"},
 		   {"set_x", (PyCFunction)set_x, METH_VARARGS, ""},
 		   {"get_x", (PyCFunction)get_x, METH_VARARGS, ""},
@@ -334,4 +426,6 @@
 		   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""},
 		   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""},
+		   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, ""},
+		   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, ""},
 		   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, ""},
 		   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, ""},
@@ -339,4 +433,5 @@
 		   {"iq", (PyCFunction)get_iq, METH_VARARGS, ""},
 		   {"pr", (PyCFunction)get_pr, METH_VARARGS, ""},
+		   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, ""},
 		   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""},
    
Index: pr_inversion/c_extensions/invertor.c
===================================================================
--- pr_inversion/c_extensions/invertor.c	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/c_extensions/invertor.c	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
@@ -40,5 +40,4 @@
 double ortho(double d_max, int n, double r) {
 	return 2.0*r*sin(pi*n*r/d_max);
-	
 }
 
@@ -76,5 +75,5 @@
  */
 double pr(double *pars, double d_max, int n_c, double r) {
-    double sum = 0.0; 
+    double sum = 0.0;  
 	int i;
     for (i=0; i<n_c; i++) {
@@ -84,3 +83,53 @@
 }
 
+/**
+ * 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;
+    }
+    *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) {
+    double sum = 0.0; 
+    double r;
+    double deriv;
+	int i;
+    for (i=0; i<25; i++) {
+    	r = d_max/25.0*i;
+    	deriv = dprdr(pars, d_max, n_c, r);
+        sum += deriv*deriv;
+    }
+    return sum/25.0*d_max;
+}
+
+
Index: pr_inversion/c_extensions/invertor.h
===================================================================
--- pr_inversion/c_extensions/invertor.h	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/c_extensions/invertor.h	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
@@ -10,6 +10,7 @@
     int npoints;    
     int ny;    
-    int nerr;    
-} Invertor_params;
+    int nerr;  
+    double alpha;
+} Invertor_params; 
 
 void invertor_dealloc(Invertor_params *pars);
@@ -24,4 +25,8 @@
 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);
+void pr_err(double *pars, double *err, double d_max, int n_c, 
+		double r, double *pr_value, double *pr_value_err);
 
 #endif
Index: pr_inversion/invertor.py
===================================================================
--- pr_inversion/invertor.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/invertor.py	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
@@ -3,4 +3,7 @@
 
 class Invertor(Cinvertor):
+    
+    ## Chisqr of the last computation
+    chisqr = 0
     
     def __init__(self):
@@ -14,4 +17,6 @@
         """
         if   name=='x':
+            if 0.0 in value:
+                raise ValueError, "Invertor: one of your q-values is zero. Delete that entry before proceeding"
             return self.set_x(value)
         elif name=='y':
@@ -21,4 +26,6 @@
         elif name=='d_max':
             return self.set_dmax(value)
+        elif name=='alpha':
+            return self.set_alpha(value)
             
         return Cinvertor.__setattr__(self, name, value)
@@ -45,4 +52,6 @@
         elif name=='d_max':
             return self.get_dmax()
+        elif name=='alpha':
+            return self.get_alpha()
         elif name in self.__dict__:
             return self.__dict__[name]
@@ -62,5 +71,51 @@
         out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True)
         
+        # Compute chi^2
+        res = self.residuals(out)
+        chisqr = 0
+        for i in range(len(res)):
+            chisqr += res[i]
+        
+        self.chi2 = chisqr
+        
         return out, cov_x
+    
+    def pr_fit(self, nfunc=5):
+        """
+            Perform inversion to P(r)
+        """
+        from scipy import optimize
+        
+        # First, check that the current data is valid
+        if self.is_valid()<=0:
+            raise RuntimeError, "Invertor.invert: Data arrays are of different length"
+        
+        p = numpy.ones(nfunc)
+        out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, full_output=1, warning=True)
+        
+        # Compute chi^2
+        res = self.pr_residuals(out)
+        chisqr = 0
+        for i in range(len(res)):
+            chisqr += res[i]
+        
+        self.chisqr = chisqr
+        
+        return out, cov_x
+    
+    def pr_err(self, c, c_cov, r):
+        import math
+        c_err = numpy.zeros(len(c))
+        for i in range(len(c)):
+            try:
+                c_err[i] = math.sqrt(math.fabs(c_cov[i][i]))
+            except:
+                import sys
+                print sys.exc_value
+                print "oups", c_cov[i][i]
+                c_err[i] = c[i]
+
+        return self.get_pr_err(c, c_err, r)
+        
     
 if __name__ == "__main__":
Index: pr_inversion/test/sphere_80.txt
===================================================================
--- pr_inversion/test/sphere_80.txt	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
+++ pr_inversion/test/sphere_80.txt	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
@@ -0,0 +1,100 @@
+<X>   <Y>
+0.0010101  2.14186e+008
+0.0020202  2.13348e+008
+0.0030303  2.11958e+008
+0.0040404  2.10025e+008
+0.00505051  2.07561e+008
+0.00606061  2.04584e+008
+0.00707071  2.01112e+008
+0.00808081  1.97169e+008
+0.00909091  1.92781e+008
+0.010101  1.87975e+008
+0.0111111  1.82783e+008
+0.0121212  1.77237e+008
+0.0131313  1.71373e+008
+0.0141414  1.65226e+008
+0.0151515  1.58833e+008
+0.0161616  1.52233e+008
+0.0171717  1.45465e+008
+0.0181818  1.38567e+008
+0.0191919  1.31578e+008
+0.020202  1.24536e+008
+0.0212121  1.17479e+008
+0.0222222  1.10442e+008
+0.0232323  1.03461e+008
+0.0242424  9.65692e+007
+0.0252525  8.97973e+007
+0.0262626  8.31748e+007
+0.0272727  7.67288e+007
+0.0282828  7.04839e+007
+0.0292929  6.4462e+007
+0.030303  5.86828e+007
+0.0313131  5.31627e+007
+0.0323232  4.79158e+007
+0.0333333  4.29533e+007
+0.0343434  3.82836e+007
+0.0353535  3.39126e+007
+0.0363636  2.98434e+007
+0.0373737  2.60768e+007
+0.0383838  2.2611e+007
+0.0393939  1.9442e+007
+0.040404  1.65639e+007
+0.0414141  1.39687e+007
+0.0424242  1.16467e+007
+0.0434343  9.58662e+006
+0.0444444  7.77602e+006
+0.0454545  6.20118e+006
+0.0464646  4.84751e+006
+0.0474747  3.69971e+006
+0.0484848  2.74197e+006
+0.0494949  1.95815e+006
+0.0505051  1.33203e+006
+0.0515152  847402
+0.0525253  488319
+0.0535354  239201
+0.0545455  84994.7
+0.0555556  11292.7
+0.0565657  4445.63
+0.0575758  51653.5
+0.0585859  141041
+0.059596  261717
+0.0606061  403816
+0.0616162  558522
+0.0626263  718084
+0.0636364  875805
+0.0646465  1.02603e+006
+0.0656566  1.1641e+006
+0.0666667  1.28634e+006
+0.0676768  1.38998e+006
+0.0686869  1.47309e+006
+0.069697  1.53454e+006
+0.0707071  1.57391e+006
+0.0717172  1.59141e+006
+0.0727273  1.5878e+006
+0.0737374  1.5643e+006
+0.0747475  1.52256e+006
+0.0757576  1.46449e+006
+0.0767677  1.39227e+006
+0.0777778  1.30823e+006
+0.0787879  1.21477e+006
+0.079798  1.11433e+006
+0.0808081  1.00933e+006
+0.0818182  902079
+0.0828283  794770
+0.0838384  689421
+0.0848485  587856
+0.0858586  491670
+0.0868687  402219
+0.0878788  320609
+0.0888889  247687
+0.089899  184050
+0.0909091  130047
+0.0919192  85790.3
+0.0929293  51176.4
+0.0939394  25900.5
+0.0949495  9481.12
+0.0959596  1283.93
+0.0969697  547.497
+0.0979798  6409.98
+0.0989899  17935.9
+0.1  34142.4
Index: pr_inversion/test/utest_invertor.py
===================================================================
--- pr_inversion/test/utest_invertor.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/test/utest_invertor.py	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
@@ -8,5 +8,5 @@
 
 
-import unittest, math, numpy
+import unittest, math, numpy, pylab
 from sans.pr.invertor import Invertor
         
@@ -21,5 +21,5 @@
         self.x_in = numpy.ones(self.ntest)
         for i in range(self.ntest):
-            self.x_in[i] = 1.0*i
+            self.x_in[i] = 1.0*(i+1)
 
 
@@ -31,4 +31,12 @@
         self.invertor.d_max = value
         self.assertEqual(self.invertor.d_max, value)
+        
+    def testset_alpha(self):
+        """
+            Set and read alpha
+        """
+        value = 15.0
+        self.invertor.alpha = value
+        self.assertEqual(self.invertor.alpha, value)
         
     def testset_x_1(self):
@@ -113,4 +121,160 @@
         self.assertEqual(self.invertor.test_no_data, None)
         
+    def test_inversion(self):
+        """
+            Test an inversion for which we know the answer
+        """
+        x, y, err = load("sphere_80.txt")
+
+        # Choose the right d_max...
+        self.invertor.d_max = 160.0
+        # Set a small alpha
+        self.invertor.alpha = 1e-7
+        # Set data
+        self.invertor.x   = x
+        self.invertor.y   = y
+        self.invertor.err = err
+        # Perform inversion
+        out, cov = self.invertor.invert(10)
+        # This is a very specific case
+        # We should make sure it always passes
+        self.assertTrue(self.invertor.chi2/len(x)<200.00)
+        
+        # Check the computed P(r) with the theory
+        # for shpere of radius 80
+        x = pylab.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0)
+        y = numpy.zeros(len(x))
+        dy = numpy.zeros(len(x))
+        y_true = numpy.zeros(len(x))
+
+        sum = 0.0
+        sum_true = 0.0
+        for i in range(len(x)):
+            #y[i] = self.invertor.pr(out, x[i])
+            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i])
+            sum += y[i]
+            if x[i]<80.0:
+                y_true[i] = pr_theory(x[i], 80.0)
+            else:
+                y_true[i] = 0
+            sum_true += y_true[i]
+            
+        y = y/sum*self.invertor.d_max/len(x)
+        dy = dy/sum*self.invertor.d_max/len(x)
+        y_true = y_true/sum_true*self.invertor.d_max/len(x)
+        
+        chi2 = 0.0
+        for i in range(len(x)):
+            res = (y[i]-y_true[i])/dy[i]
+            chi2 += res*res
+            
+        try:
+            self.assertTrue(chi2/51.0<10.0)
+        except:
+            print "chi2 =", chi2/51.0
+            raise
+            
+    def test_q_zero(self):
+        """
+            Test error condition where a point has q=0
+        """
+        x, y, err = load("sphere_80.txt")
+        x[0] = 0.0
+        
+        # Choose the right d_max...
+        self.invertor.d_max = 160.0
+        # Set a small alpha
+        self.invertor.alpha = 1e-7
+        # Set data
+        def doit():
+            self.invertor.x   = x
+        self.assertRaises(ValueError, doit)
+        
+                            
+    def test_q_neg(self):
+        """
+            Test error condition where a point has q<0
+        """
+        x, y, err = load("sphere_80.txt")
+        x[0] = -0.2
+        
+        # Choose the right d_max...
+        self.invertor.d_max = 160.0
+        # Set a small alpha
+        self.invertor.alpha = 1e-7
+        # Set data
+        self.invertor.x   = x
+        self.invertor.y   = y
+        self.invertor.err = err
+        # Perform inversion
+        out, cov = self.invertor.invert(4)
+        
+        try:
+            self.assertTrue(self.invertor.chi2>0)
+        except:
+            print "Chi2 =", self.invertor.chi2
+            raise
+                            
+    def test_Iq_zero(self):
+        """
+            Test error condition where a point has q<0
+        """
+        x, y, err = load("sphere_80.txt")
+        y[0] = 0.0
+        
+        # Choose the right d_max...
+        self.invertor.d_max = 160.0
+        # Set a small alpha
+        self.invertor.alpha = 1e-7
+        # Set data
+        self.invertor.x   = x
+        self.invertor.y   = y
+        self.invertor.err = err
+        # Perform inversion
+        out, cov = self.invertor.invert(4)
+        
+        try:
+            self.assertTrue(self.invertor.chi2>0)
+        except:
+            print "Chi2 =", self.invertor.chi2
+            raise
+                            
+        
+
+def pr_theory(r, R):
+    """
+       P(r) for a sphere
+    """
+    if r<=2*R:
+        return 12.0* ((0.5*r/R)**2) * ((1.0-0.5*r/R)**2) * ( 2.0 + 0.5*r/R )
+    else:
+        return 0.0
+    
+def load(path = "sphere_60_q0_2.txt"):
+    import numpy, math, sys
+    # Read the data from the data file
+    data_x   = numpy.zeros(0)
+    data_y   = numpy.zeros(0)
+    data_err = numpy.zeros(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()
+                x = float(toks[0])
+                y = float(toks[1])
+                data_x = numpy.append(data_x, x)
+                data_y = numpy.append(data_y, y)
+                # Set the error of the first point to 5%
+                # to make theory look like data
+                scale = 0.1/math.sqrt(data_x[0])
+                data_err = numpy.append(data_err, scale*math.sqrt(y))
+            except:
+                pass
+               
+    return data_x, data_y, data_err 
+       
 if __name__ == '__main__':
     unittest.main()
