Index: sansmodels/src/c_models/libfunc.c
===================================================================
--- sansmodels/src/c_models/libfunc.c	(revision f8644dd1e4253bc00e1fcb3452ed2071cf21f61e)
+++ sansmodels/src/c_models/libfunc.c	(revision 5da3cc57a0d2568f3d27414a6c82006204b81eca)
@@ -11,102 +11,68 @@
 //used in Si func
 
-int factorial(int i)
+int factorial(int i) {
+
+	int k, j;
+	if (i<2){
+		return 1;
+	}
+
+	k=1;
+
+	for(j=1;j<i;j++) {
+		k=k*(j+1);
+	}
+
+	return k;
+
+}
+
+
+
+// Used in pearl nec model
+
+// Sine integral function: approximated within 1%!!!
+
+// integral of sin(x)/x up to namx term nmax=6 looks the best.
+
+double Si(double x)
 
 {
-
-	int k, j;
-
-	if (i<2){
-
-		return 1;
-
-	}
-
-	k=1;
-
-	for(j=1;j<i;j++)
-
-	{
-
-		k=k*(j+1);
-
-	}
-
-	return k;
-
-}
-
-
-
-// Used in pearl nec model
-
-// Sine integral function: approximated within 1%!!!
-
-// integral of sin(x)/x up to namx term nmax=6 looks the best.
-
-double Si(double x)
-
-{
-
 	int i;
-
 	int nmax=6;
-
 	double out;
-
 	long double power;
-
 	double pi = 4.0*atan(1.0);
 
 	if (x >= pi*6.2/4.0){
-
-
-
 		double out_sin = 0.0;
-
 		double out_cos = 0.0;
-
 		out = pi/2.0;
 
-		for (i=0; i<nmax-2; i+=1){
-
+		for (i=0; i<nmax-2; i+=1) {
 			out_cos += pow(-1.0, i) * (double)factorial(2*i) / pow(x, 2*i+1);
-
 			out_sin += pow(-1.0, i) * (double)factorial(2*i+1) / pow(x, 2*i+2);
-
 		}
 
 		out -= cos(x) * out_cos;
-
 		out -= sin(x) * out_sin;
-
 		return out;
-
 	}
 
 	out = 0.0;
 
-	for (i=0; i<nmax; i+=1)
-
-	{
-
-		if (i==0){
-
+	for (i=0; i<nmax; i+=1)	{
+		if (i==0) {
 			out += x;
-
 			continue;
-
 		}
 
 		power = pow(x,(2 * i + 1));
-
 		out += (double)pow(-1, i) * power / ((2.0 * (double)i + 1.0) * (double)factorial(2 * i + 1));
 
 		//printf ("Si=%g %g %d\n", x, out, i);
-
 	}
 
 	return out;
-
 }
 
@@ -138,2 +104,102 @@
 }
 
+/** Modifications below by kieranrcampbell@gmail.com
+    Institut Laue-Langevin, July 2012
+**/
+
+/**
+   Implements eq 6.2.5 (small gamma) of Numerical Recipes in C, essentially
+   the incomplete gamma function multiplied by the gamma function.
+   Required for implementation of fast error function (erf)
+**/
+
+
+#define ITMAX 100
+#define EPS 3.0e-7
+#define FPMIN 1.0e-30
+
+void gser(float *gamser, float a, float x, float *gln) {
+  int n;
+  float sum,del,ap;
+
+  *gln = gamln(a);
+  if(x <= 0.0) {
+    if (x < 0.0) printf("Error: x less than 0 in routine gser");
+    *gamser = 0.0;
+    return;
+  } else {
+    ap = a;
+    del = sum = 1.0/a;
+    
+    for(n=1;n<=ITMAX;n++) {
+      ++ap;
+      del *= x/ap;
+      sum += del;
+      if(fabs(del) < fabs(sum)*EPS) {
+	*gamser = sum * exp(-x + a * log(x) - (*gln));
+	return;
+      }
+    }
+    printf("a too large, ITMAX too small in routine gser");
+    return;
+
+  }
+
+
+}
+
+/**
+   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction 
+   representation 
+**/
+
+void gcf(float *gammcf, float a, float x, float *gln) {
+  int i;
+  float an,b,c,d,del,h;
+
+  *gln = gamln(a);
+  b = x+1.0-a;
+  c = 1.0/FPMIN;
+  d = 1.0/b;
+  h=d;
+  for (i=1;i <= ITMAX; i++) {
+    an = -i*(i-a);
+    b += 2.0;
+    d = an*d + b;
+    if (fabs(d) < FPMIN) d = FPMIN;
+    c = b+an/c;
+    if (fabs(c) < FPMIN) c = FPMIN;
+    d = 1.0/d;
+    del = d*c;
+    h += del;
+    if (fabs(del-1.0) < EPS) break;
+  }
+  if (i > ITMAX) printf("a too large, ITMAX too small in gcf");
+  *gammcf = exp(-x+a*log(x)-(*gln))*h;
+  return;
+}
+
+
+/**
+   Represents incomplete error function, P(a,x)
+**/
+float gammp(float a, float x) {
+  float gamser,gammcf,gln;
+  if(x < 0.0 || a <= 0.0) printf("Invalid arguments in routine gammp");
+  if (x < (a+1.0)) {
+    gser(&gamser,a,x,&gln);
+    return gamser;
+  } else {
+    gcf(&gammcf,a,x,&gln);
+    return 1.0 - gammcf;
+  }
+}
+
+/** 
+    Implementation of the error function, erf(x)
+**/
+
+float erff(float x) {
+  return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
+}
+
Index: sansmodels/src/python_wrapper/WrapperGenerator.py
===================================================================
--- sansmodels/src/python_wrapper/WrapperGenerator.py	(revision fa6db8bddd1335d532cec83bac89eb720c235a65)
+++ sansmodels/src/python_wrapper/WrapperGenerator.py	(revision 5da3cc57a0d2568f3d27414a6c82006204b81eca)
@@ -92,4 +92,6 @@
         # Model category
         self.category = None
+        # Whether model belongs to multifunc
+        self.is_multifunc = False
         ## output directory for wrappers
         self.output_dir = output_dir
@@ -228,5 +230,14 @@
                     raise ValueError, "Could not parse file %s" % self.file
 
-
+            # is_multifunc
+            key = "[MULTIPLICITY_INFO]"
+            if line.count(key) > 0:
+                self.is_multifunc = True
+                try:
+                    index = line.index(key)
+                    toks = line[index:].split("=")
+                    self.multiplicity_info = toks[1].lstrip().rstrip()
+                except:
+                    raise ValueError, "Could not parse file %s" % self.file
 
             # Catch struct name
@@ -331,4 +342,14 @@
             newline = self.replaceToken(newline, 
                                         "[MODELSTRUCT]", self.structName)
+
+            # Sort model initialization based on multifunc
+            if(self.is_multifunc):
+                line = "int level = 1;\nPyArg_ParseTuple(args,\"i\",&level);\n"
+                line += "self->model = new " + self.pythonClass + "(level);"
+            else:
+                line = "self->model = new " + self.pythonClass + "();"
+    
+            newline = self.replaceToken(newline,"[INITIALIZE_MODEL]",
+                                            line)
             
             # Dictionary initialization
@@ -465,4 +486,15 @@
             newline = self.replaceToken(newline, 
                                         "[PAR_DETAILS]", self.details)
+            
+            # Call base constructor
+            if self.is_multifunc:
+                newline = self.replaceToken(newline,"[CALL_CPYTHON_INIT]",
+                                            'C' + self.pythonClass + ".__init__(self,multfactor)\n\tself.is_multifunc = True")
+                newline = self.replaceToken(newline,"[MULTIPLICITY_INFO]",self.multiplicity_info)
+            else:
+                newline = self.replaceToken(newline,"[CALL_CPYTHON_INIT]",
+                                            'C' + self.pythonClass + ".__init__(self)\n\tself.is_multifunc = False")
+                newline = self.replaceToken(newline,"[MULTIPLICITY_INFO]","None")
+
            
             # fixed list  details
Index: sansmodels/src/python_wrapper/classTemplate.txt
===================================================================
--- sansmodels/src/python_wrapper/classTemplate.txt	(revision e08bd5b099f35eabce0c3e5a23f6752b42dad397)
+++ sansmodels/src/python_wrapper/classTemplate.txt	(revision 5da3cc57a0d2568f3d27414a6c82006204b81eca)
@@ -86,6 +86,7 @@
         self->params = PyDict_New();
         self->dispersion = PyDict_New();
-        self->model = new [CMODEL]();
-        
+
+	[INITIALIZE_MODEL]
+
         [INITDICTIONARY]
          
Index: sansmodels/src/python_wrapper/modelTemplate.txt
===================================================================
--- sansmodels/src/python_wrapper/modelTemplate.txt	(revision fa6db8bddd1335d532cec83bac89eb720c235a65)
+++ sansmodels/src/python_wrapper/modelTemplate.txt	(revision 5da3cc57a0d2568f3d27414a6c82006204b81eca)
@@ -43,5 +43,5 @@
     """
         
-    def __init__(self):
+    def __init__(self,multfactor=1):
         """ Initialization """
         self.__dict__ = {}
@@ -50,6 +50,7 @@
         BaseComponent.__init__(self)
         #apply([CPYTHONCLASS].__init__, (self,)) 
-        [CPYTHONCLASS].__init__(self)
-        
+
+	[CALL_CPYTHON_INIT]
+		        
         ## Name of the model
         self.name = "[PYTHONCLASS]"
@@ -70,4 +71,6 @@
 
         self.category = [CATEGORY]
+        self.multiplicity_info = [MULTIPLICITY_INFO]
+        
 
     def __setstate__(self, state):
