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
