Index: src/sas/sascalc/calculator/c_extensions/__init__.py
===================================================================
--- src/sas/sascalc/calculator/c_extensions/__init__.py	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/__init__.py	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,3 @@
+"""
+    C extensions to provide the sas_gen computations.
+"""
Index: src/sas/sascalc/calculator/c_extensions/libfunc.c
===================================================================
--- src/sas/sascalc/calculator/c_extensions/libfunc.c	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/libfunc.c	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,332 @@
+// by jcho
+
+#include <math.h>
+
+#include "libfunc.h"
+
+#include <stdio.h>
+
+
+
+//used in Si func
+
+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 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) {
+			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) {
+			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;
+}
+
+
+
+double sinc(double x)
+{
+  if (x==0.0){
+    return 1.0;
+  }
+  return sin(x)/x;
+}
+
+
+double gamln(double xx) {
+
+    double x,y,tmp,ser;
+    static double cof[6]={76.18009172947146,-86.50532032941677,
+    24.01409824083091,-1.231739572450155,
+    0.1208650973866179e-2,-0.5395239384953e-5};
+    int j;
+
+    y=x=xx;
+    tmp=x+5.5;
+    tmp -= (x+0.5)*log(tmp);
+    ser=1.000000000190015;
+    for (j=0;j<=5;j++) ser += cof[j]/++y;
+    return -tmp+log(2.5066282746310005*ser/x);
+}
+
+// calculate magnetic sld and return total sld
+// bn : contrast (not just sld of the layer)
+// m0: max mag of M; mtheta: angle from x-z plane;
+// mphi: angle (anti-clock-wise)of x-z projection(M) from x axis
+// spinfraci: the fraction of UP among UP+Down (before sample)
+// spinfracf: the fraction of UP among UP+Down (after sample and before detector)
+// spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis
+// Note: all angles are in degrees.
+polar_sld cal_msld(int isangle, double qx, double qy, double bn,
+				double m01, double mtheta1, double mphi1,
+				double spinfraci, double spinfracf, double spintheta)
+{
+	//locals
+	double q_x = qx;
+	double q_y = qy;
+	double sld = bn;
+	int is_angle = isangle;
+	double pi = 4.0*atan(1.0);
+	double s_theta = spintheta * pi/180.0;
+	double m_max = m01;
+	double m_phi = mphi1;
+	double m_theta = mtheta1;
+	double in_spin = spinfraci;
+	double out_spin = spinfracf;
+
+	double m_perp = 0.0;
+	double m_perp_z = 0.0;
+	double m_perp_y = 0.0;
+	double m_perp_x = 0.0;
+	double m_sigma_x = 0.0;
+	double m_sigma_z = 0.0;
+	double m_sigma_y = 0.0;
+	//double b_m = 0.0;
+	double q_angle = 0.0;
+	double mx = 0.0;
+	double my = 0.0;
+	double mz = 0.0;
+	polar_sld p_sld;
+	p_sld.uu = sld;
+	p_sld.dd = sld;
+	p_sld.re_ud = 0.0;
+	p_sld.im_ud = 0.0;
+	p_sld.re_du = 0.0;
+	p_sld.im_du = 0.0;
+
+	//No mag means no further calculation
+	if (isangle>0){
+		if (m_max < 1.0e-32){
+			p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+			p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+			return p_sld;
+		}
+	}
+	else{
+		if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){
+			p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+			p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+			return p_sld;
+		}
+	}
+
+	//These are needed because of the precision of inputs
+	if (in_spin < 0.0) in_spin = 0.0;
+	if (in_spin > 1.0) in_spin = 1.0;
+	if (out_spin < 0.0) out_spin = 0.0;
+	if (out_spin > 1.0) out_spin = 1.0;
+
+	if (q_x == 0.0) q_angle = pi / 2.0;
+	else q_angle = atan(q_y/q_x);
+	if (q_y < 0.0 && q_x < 0.0) q_angle -= pi;
+	else if (q_y > 0.0 && q_x < 0.0) q_angle += pi;
+
+	q_angle = pi/2.0 - q_angle;
+	if (q_angle > pi) q_angle -= 2.0 * pi;
+	else if (q_angle < -pi) q_angle += 2.0 * pi;
+
+	if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){
+		m_perp = 0.0;
+		}
+	else {
+		m_perp = m_max;
+		}
+	if (is_angle > 0){
+		m_phi *= pi/180.0;
+		m_theta *= pi/180.0;
+		mx = m_perp * cos(m_theta) * cos(m_phi);
+		my = m_perp * sin(m_theta);
+		mz = -(m_perp * cos(m_theta) * sin(m_phi));
+	}
+	else{
+		mx = m_perp;
+		my = m_phi;
+		mz = m_theta;
+	}
+	//ToDo: simplify these steps
+	// m_perp1 -m_perp2
+	m_perp_x = (mx) *  cos(q_angle);
+	m_perp_x -= (my) * sin(q_angle);
+	m_perp_y = m_perp_x;
+	m_perp_x *= cos(-q_angle);
+	m_perp_y *= sin(-q_angle);
+	m_perp_z = mz;
+
+	m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta));
+	m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta));
+	m_sigma_z = (m_perp_z);
+
+	//Find b
+	p_sld.uu -= m_sigma_x;
+	p_sld.dd += m_sigma_x;
+	p_sld.re_ud = m_sigma_y;
+	p_sld.re_du = m_sigma_y;
+	p_sld.im_ud = m_sigma_z;
+	p_sld.im_du = -m_sigma_z;
+
+	p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+	p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+
+	p_sld.re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.re_ud;
+	p_sld.im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.im_ud;
+	p_sld.re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.re_du;
+	p_sld.im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.im_du;
+
+	return p_sld;
+}
+
+
+/** 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: src/sas/sascalc/calculator/c_extensions/libfunc.h
===================================================================
--- src/sas/sascalc/calculator/c_extensions/libfunc.h	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/libfunc.h	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,31 @@
+#if !defined(o_h)
+#define libfunc_h
+typedef struct {
+	double uu;
+	double dd;
+	double re_ud;
+	double im_ud;
+	double re_du;
+	double im_du;
+} polar_sld;
+
+int factorial(int i);
+
+double Si(double x);
+
+double sinc(double x);
+
+double gamln(double x);
+
+polar_sld cal_msld(int isangle, double qx, double qy, double bn, double m01, double mtheta1, 
+			double mphi1, double spinfraci, double spinfracf, double spintheta);
+
+void gser(float *gamser, float a, float x, float *gln);
+
+void gcf(float *gammcf, float a, float x, float *gln);
+
+float gammp(float a,float x);
+
+float erff(float x);
+
+#endif
Index: src/sas/sascalc/calculator/c_extensions/librefl.c
===================================================================
--- src/sas/sascalc/calculator/c_extensions/librefl.c	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/librefl.c	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,331 @@
+// The original code, of which work was not DANSE funded,
+// was provided by J. Cho.
+// And modified to fit sansmodels/sansview: JC
+
+#include <math.h>
+#include "librefl.h"
+#include <stdio.h>
+#include <stdlib.h>
+#if defined(_MSC_VER)
+#include "winFuncs.h"
+#endif
+
+complex cassign(real, imag)
+	double real, imag;
+{
+	complex x;
+	x.re = real;
+	x.im = imag;
+	return x;
+}
+
+
+complex cplx_add(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re + y.re;
+	z.im = x.im + y.im;
+	return z;
+}
+
+complex rcmult(x,y)
+	double x;
+    complex y;
+{
+	complex z;
+	z.re = x*y.re;
+	z.im = x*y.im;
+	return z;
+}
+
+complex cplx_sub(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re - y.re;
+	z.im = x.im - y.im;
+	return z;
+}
+
+
+complex cplx_mult(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = x.re*y.re - x.im*y.im;
+	z.im = x.re*y.im + x.im*y.re;
+	return z;
+}
+
+complex cplx_div(x,y)
+	complex x,y;
+{
+	complex z;
+	z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
+	z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
+	return z;
+}
+
+complex cplx_exp(b)
+	complex b;
+{
+	complex z;
+	double br,bi;
+	br=b.re;
+	bi=b.im;
+	z.re = exp(br)*cos(bi);
+	z.im = exp(br)*sin(bi);
+	return z;
+}
+
+
+complex cplx_sqrt(z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
+	complex z;
+{
+	complex c;
+	double zr,zi,x,y,r,w;
+
+	zr=z.re;
+	zi=z.im;
+
+	if (zr==0.0 && zi==0.0)
+	{
+    c.re=0.0;
+	c.im=0.0;
+    return c;
+	}
+	else
+	{
+		x=fabs(zr);
+		y=fabs(zi);
+		if (x>y)
+		{
+			r=y/x;
+			w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
+		}
+		else
+		{
+			r=x/y;
+			w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
+		}
+		if (zr >=0.0)
+		{
+			c.re=w;
+			c.im=zi/(2.0*w);
+		}
+		else
+		{
+			c.im=(zi >= 0) ? w : -w;
+			c.re=zi/(2.0*c.im);
+		}
+		return c;
+	}
+}
+
+complex cplx_cos(b)
+	complex b;
+{
+	complex zero,two,z,i,bi,negbi;
+	zero = cassign(0.0,0.0);
+	two = cassign(2.0,0.0);
+	i = cassign(0.0,1.0);
+	bi = cplx_mult(b,i);
+	negbi = cplx_sub(zero,bi);
+	z = cplx_div(cplx_add(cplx_exp(bi),cplx_exp(negbi)),two);
+	return z;
+}
+
+// normalized and modified erf
+//   |
+// 1 +                __  - - - -
+//   |             _
+//	 |            _
+//   |        __
+// 0 + - - -
+//   |-------------+------------+--
+//   0           center       n_sub    --->
+//                                     ind
+//
+// n_sub = total no. of bins(or sublayers)
+// ind = x position: 0 to max
+// nu = max x to integration
+double err_mod_func(double n_sub, double ind, double nu)
+{
+  double center, func;
+  if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+
+	//ind = (n_sub-1.0)/2.0-1.0 +ind;
+	center = n_sub/2.0;
+	// transform it so that min(ind) = 0
+	ind -= center;
+	// normalize by max limit
+	ind /= center;
+	// divide by sqrt(2) to get Gaussian func
+	nu /= sqrt(2.0);
+	ind *= nu;
+	// re-scale and normalize it so that max(erf)=1, min(erf)=0
+	func = erf(ind)/erf(nu)/2.0;
+	// shift it by +0.5 in y-direction so that min(erf) = 0
+	func += 0.5;
+
+	return func;
+}
+double linearfunc(double n_sub, double ind, double nu)
+{
+  double bin_size, func;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	func = ind;
+
+	return func;
+}
+// use the right hand side from the center of power func
+double power_r(double n_sub, double ind, double nu)
+{
+  double bin_size,func;
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	func = pow(ind, nu);
+
+	return func;
+}
+// use the left hand side from the center of power func
+double power_l(double n_sub, double ind, double nu)
+{
+  double bin_size, func;
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	func = 1.0-pow((1.0-ind),nu);
+
+	return func;
+}
+// use 1-exp func from x=0 to x=1
+double exp_r(double n_sub, double ind, double nu)
+{
+  double bin_size, func;
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	// modify func so that func(0) =0 and func(max)=1
+	func = 1.0-exp(-nu*ind);
+	// normalize by its max
+	func /= (1.0-exp(-nu));
+
+	return func;
+}
+
+// use the left hand side mirror image of exp func
+double exp_l(double n_sub, double ind, double nu)
+{
+  double bin_size, func;
+	if (nu == 0.0)
+		nu = 1e-14;
+	if (n_sub == 0.0)
+		n_sub = 1.0;
+
+	bin_size = 1.0/n_sub;  //size of each sub-layer
+	// rescale
+	ind *= bin_size;
+	// modify func
+	func = exp(-nu*(1.0-ind))-exp(-nu);
+	// normalize by its max
+	func /= (1.0-exp(-nu));
+
+	return func;
+}
+
+// To select function called
+// At nu = 0 (singular point), call line function
+double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r)
+{
+	double sld_i, func;
+	// this condition protects an error from the singular point
+	if (nu == 0.0){
+		nu = 1e-13;
+	}
+	// select func
+	switch(fun_type){
+		case 1 :
+			func = power_r(n_sub, i, nu);
+			break;
+		case 2 :
+			func = power_l(n_sub, i, nu);
+			break;
+		case 3 :
+			func = exp_r(n_sub, i, nu);
+			break;
+		case 4 :
+			func = exp_l(n_sub, i, nu);
+			break;
+		case 5 :
+			func = linearfunc(n_sub, i, nu);
+			break;
+		default:
+			func = err_mod_func(n_sub, i, nu);
+			break;
+	}
+	// compute sld
+	if (sld_r>sld_l){
+		sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
+	}
+	else if (sld_r<sld_l){
+		func = 1.0-func;
+		sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
+	}
+	else{
+		sld_i = sld_r;
+	}
+	return sld_i;
+}
+
+
+// used by refl.c
+double interfunc(int fun_type, double n_sub, double i, double sld_l, double sld_r)
+{
+	double sld_i, func;
+	switch(fun_type){
+		case 0 :
+			func = err_mod_func(n_sub, i, 2.5);
+			break;
+		default:
+			func = linearfunc(n_sub, i, 1.0);
+			break;
+	}
+	if (sld_r>sld_l){
+		sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
+	}
+	else if (sld_r<sld_l){
+		func = 1.0-func;
+		sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
+	}
+	else{
+		sld_i = sld_r;
+	}
+	return sld_i;
+}
Index: src/sas/sascalc/calculator/c_extensions/librefl.h
===================================================================
--- src/sas/sascalc/calculator/c_extensions/librefl.h	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/librefl.h	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,37 @@
+#if !defined(librefl_h)
+#define librefl_h
+
+typedef struct {
+	double re;
+	double im;
+} complex;
+
+typedef struct {
+	complex a;
+	complex b;
+	complex c;
+	complex d;
+} matrix;
+
+complex cassign(double real, double imag);
+
+complex cplx_add(complex x,complex y);
+
+complex rcmult(double x,complex y);
+
+complex cplx_sub(complex x,complex y);
+
+complex cplx_mult(complex x,complex y);
+
+complex cplx_div(complex x,complex y);
+
+complex cplx_exp(complex b);
+
+complex cplx_sqrt(complex z);
+
+complex cplx_cos(complex b);
+
+double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r);
+double interfunc(int fun_type, double n_sub, double i, double sld_l, double sld_r);
+double linePq(int fun_type, double thick, double sld_in, double sld_out,double r, double q);
+#endif
Index: src/sas/sascalc/calculator/c_extensions/sld2i.cpp
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i.cpp	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/sld2i.cpp	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,210 @@
+/**
+Computes the (magnetic) scattering form sld (n and m) profile
+ */
+#include "sld2i.hh"
+#include <stdio.h>
+#include <math.h>
+using namespace std;
+extern "C" {
+	#include "libfunc.h"
+	#include "librefl.h"
+}
+/**
+ * Constructor for GenI
+ *
+ * binning
+ * //@param qx: array of Qx values
+ * //@param qy: array of Qy values
+ * //@param qz: array of Qz values
+ * @param x: array of x values
+ * @param y: array of y values
+ * @param z: array of z values
+ * @param sldn: array of sld n
+ * @param mx: array of sld mx
+ * @param my: array of sld my
+ * @param mz: array of sld mz
+ * @param in_spin: ratio of up spin in Iin
+ * @param out_spin: ratio of up spin in Iout
+ * @param s_theta: angle (from x-axis) of the up spin in degree
+ */
+GenI :: GenI(int npix, double* x, double* y, double* z, double* sldn,
+			double* mx, double* my, double* mz, double* voli,
+			double in_spin, double out_spin,
+			double s_theta) {
+	//this->qx_val = qx;
+	//this->qy_val = qy;
+	//this->qz_val = qz;
+	this->n_pix = npix;
+	this->x_val = x;
+	this->y_val = y;
+	this->z_val = z;
+	this->sldn_val = sldn;
+	this->mx_val = mx;
+	this->my_val = my;
+	this->mz_val = mz;
+	this->vol_pix = voli;
+	this->inspin = in_spin;
+	this->outspin = out_spin;
+	this->stheta = s_theta;
+};
+
+/**
+ * Compute 2D anisotropic
+ */
+void GenI :: genicomXY(int npoints, double *qx, double *qy, double *I_out){
+	//npoints is given negative for angular averaging 
+	// Assumes that q doesn't have qz component and sld_n is all real
+	//double q = 0.0;
+	//double Pi = 4.0*atan(1.0);
+	polar_sld b_sld;
+	double qr = 0.0;
+	complex iqr = cassign(0.0, 0.0);
+	complex ephase = cassign(0.0, 0.0);
+	complex comp_sld = cassign(0.0, 0.0);
+
+	complex sumj_uu;
+	complex sumj_ud;
+	complex sumj_du;
+	complex sumj_dd;
+	complex temp_fi;
+
+	double count = 0.0;
+	//check if this computation is for averaging
+
+	//Assume that pixel volumes are given in vol_pix in A^3 unit
+	//int x_size = 0; //in Ang
+	//int y_size = 0; //in Ang
+	//int z_size = 0; //in Ang
+	
+	// Loop over q-values and multiply apply matrix
+	
+	for(int i=0; i<npoints; i++){
+		//I_out[i] = 0.0;
+		sumj_uu = cassign(0.0, 0.0);
+		sumj_ud = cassign(0.0, 0.0);
+		sumj_du = cassign(0.0, 0.0);
+		sumj_dd = cassign(0.0, 0.0);		
+		//printf ("%d ", i);
+		//q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]);
+
+		for(int j=0; j<n_pix; j++){
+			if (sldn_val[j]!=0.0||mx_val[j]!=0.0||my_val[j]!=0.0||mz_val[j]!=0.0)
+			{	
+				//anisotropic
+				temp_fi = cassign(0.0, 0.0);
+				b_sld = cal_msld(0, qx[i], qy[i], sldn_val[j],
+							 mx_val[j], my_val[j], mz_val[j],
+			 				 inspin, outspin, stheta);
+				qr = (qx[i]*x_val[j] + qy[i]*y_val[j]);
+				iqr = cassign(0.0, qr);
+				ephase = cplx_exp(iqr);
+				
+				//Let's multiply pixel(atomic) volume here
+				ephase = rcmult(vol_pix[j], ephase);
+				//up_up
+				if (inspin > 0.0 && outspin > 0.0){
+					comp_sld = cassign(b_sld.uu, 0.0);
+					temp_fi = cplx_mult(comp_sld, ephase);
+					sumj_uu = cplx_add(sumj_uu, temp_fi);
+				}
+				//down_down
+				if (inspin < 1.0 && outspin < 1.0){
+					comp_sld = cassign(b_sld.dd, 0.0);
+					temp_fi = cplx_mult(comp_sld, ephase);
+					sumj_dd = cplx_add(sumj_dd, temp_fi);
+				}
+				//up_down
+				if (inspin > 0.0 && outspin < 1.0){
+					comp_sld = cassign(b_sld.re_ud, b_sld.im_ud);
+					temp_fi = cplx_mult(comp_sld, ephase);
+					sumj_ud = cplx_add(sumj_ud, temp_fi);
+				}
+				//down_up
+				if (inspin < 1.0 && outspin > 0.0){
+					comp_sld = cassign(b_sld.re_du, b_sld.im_du);
+					temp_fi = cplx_mult(comp_sld, ephase);
+					sumj_du = cplx_add(sumj_du, temp_fi);
+				}
+
+
+				if (i == 0){
+					count += vol_pix[j];
+				}
+			}
+		}
+		//printf("aa%d=%g %g %d\n", i, (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im), (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im), count);
+
+		I_out[i] = (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im);
+		I_out[i] += (sumj_ud.re*sumj_ud.re + sumj_ud.im*sumj_ud.im);
+		I_out[i] += (sumj_du.re*sumj_du.re + sumj_du.im*sumj_du.im);
+		I_out[i] += (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im);
+
+		I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
+	}
+	//printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
+}
+/**
+ * Compute 1D isotropic
+ * Isotropic: Assumes all slds are real (no magnetic)
+ * Also assumes there is no polarization: No dependency on spin
+ */
+void GenI :: genicom(int npoints, double *q, double *I_out){
+	//npoints is given negative for angular averaging 
+	// Assumes that q doesn't have qz component and sld_n is all real
+	//double Pi = 4.0*atan(1.0);
+	int is_sym = 0;
+	double qr = 0.0;
+	double sumj;
+	double sld_j = 0.0;
+	double count = 0.0;
+	if (n_pix < 0 ){
+		is_sym = 1;
+		n_pix = n_pix * -1;
+	}
+	//Assume that pixel volumes are given in vol_pix in A^3 unit
+	// Loop over q-values and multiply apply matrix
+	for(size_t i=0; i<npoints; i++){
+		sumj =0.0;		
+		for(size_t j=0; j<n_pix; j++){
+			//Isotropic: Assumes all slds are real (no magnetic)
+			//Also assumes there is no polarization: No dependency on spin
+			if (is_sym == 1){
+				// approximation for a spherical symmetric particle
+				qr = sqrt(x_val[j]*x_val[j]+y_val[j]*y_val[j]+z_val[j]*z_val[j])*q[i];
+				if (qr > 0.0){
+					qr = sin(qr) / qr;
+					sumj += sldn_val[j] * vol_pix[j] * qr;
+				}
+				else{
+					sumj += sldn_val[j] * vol_pix[j];
+				}
+			}
+			else{
+				//full calculation
+				//pragma omp parallel for
+				for(size_t k=0; k<n_pix; k++){
+					sld_j =  sldn_val[j] * sldn_val[k] * vol_pix[j] * vol_pix[k];
+					qr = (x_val[j]-x_val[k])*(x_val[j]-x_val[k])+
+						      (y_val[j]-y_val[k])*(y_val[j]-y_val[k])+ 
+						      (z_val[j]-z_val[k])*(z_val[j]-z_val[k]);
+					qr = sqrt(qr) * q[i];
+					if (qr > 0.0){
+						sumj += sld_j*sin(qr)/qr;
+					}
+					else{
+						sumj += sld_j;
+					}
+				}
+			}
+			if (i == 0){
+				count += vol_pix[j];
+			}
+		}
+		I_out[i] = sumj;
+		if (is_sym == 1){
+			I_out[i] *= sumj;
+		}
+		I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
+	}
+	//printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
+}
Index: src/sas/sascalc/calculator/c_extensions/sld2i.hh
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i.hh	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/sld2i.hh	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,45 @@
+/**
+Computes the (magnetic) scattering form sld (n and m) profile
+ */
+#ifndef SLD2I_CLASS_H
+#define SLD2I_CLASS_H
+
+#include <vector>
+
+using namespace std;
+
+/**
+ * Base class
+ */
+class GenI {
+protected:
+	vector<double>* I_out;
+	// vectors
+	int n_pix;
+	double* qx;
+	double* qy;
+	double* x_val;
+	double* y_val;
+	double* z_val;
+	double* sldn_val;
+	double* mx_val;
+	double* my_val;
+	double* mz_val;
+	double* vol_pix;
+    // spin ratios
+    double inspin;
+    double outspin;
+    double stheta;
+
+public:
+    // Constructor
+	GenI(int npix, double* x, double* y, double* z, double* sldn, 
+			double* mx, double* my, double* mz, double* voli,
+			double in_spin, double out_spin,
+			double s_theta);
+	// compute function
+	virtual void genicomXY(int npoints, double* qx, double* qy, double *I_out);
+	virtual void genicom(int npoints, double* q, double *I_out);
+};
+
+#endif
Index: src/sas/sascalc/calculator/c_extensions/sld2i_module.cpp
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i_module.cpp	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/sld2i_module.cpp	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,159 @@
+/**
+  SLD2I module to perform point and I calculations
+ */
+#include <Python.h>
+#include <stdio.h>
+#include <sld2i.hh>
+
+// Utilities
+#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)
+
+
+/**
+ * Delete a GenI object
+ */
+void del_sld2i(void *ptr){
+	GenI* sld2i = static_cast<GenI *>(ptr);
+	delete sld2i;
+	return;
+}
+
+/**
+ * Create a GenI as a python object by supplying arrays
+ */
+PyObject * new_GenI(PyObject *, PyObject *args) {
+	PyObject *x_val_obj;
+	PyObject *y_val_obj;
+	PyObject *z_val_obj;
+	PyObject *sldn_val_obj;
+	PyObject *mx_val_obj;
+	PyObject *my_val_obj;
+	PyObject *mz_val_obj;
+	PyObject *vol_pix_obj;
+	Py_ssize_t n_x;
+	//PyObject rlimit_obj;
+	//PyObject npoints_obj;
+	//PyObject nrbins_obj;
+	//PyObject nphibins_obj;
+	int n_pix;
+	double* x_val;
+	double* y_val;
+	double* z_val;
+	double* sldn_val;
+	double* mx_val;
+	double* my_val;
+	double* mz_val;
+	double* vol_pix;
+	double inspin;
+	double outspin;
+	double stheta;
+
+	if (!PyArg_ParseTuple(args, "iOOOOOOOOddd", &n_pix, &x_val_obj, &y_val_obj, &z_val_obj, &sldn_val_obj, &mx_val_obj, &my_val_obj, &mz_val_obj, &vol_pix_obj, &inspin, &outspin, &stheta)) return NULL;
+	OUTVECTOR(x_val_obj, x_val, n_x);
+	OUTVECTOR(y_val_obj, y_val, n_x);
+	OUTVECTOR(z_val_obj, z_val, n_x);
+	OUTVECTOR(sldn_val_obj, sldn_val, n_x);
+	OUTVECTOR(mx_val_obj, mx_val, n_x);
+	OUTVECTOR(my_val_obj, my_val, n_x);
+	OUTVECTOR(mz_val_obj, mz_val, n_x);
+	OUTVECTOR(vol_pix_obj, vol_pix, n_x);
+	GenI* sld2i = new GenI(n_pix,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta);
+	return PyCObject_FromVoidPtr(sld2i, del_sld2i);
+}
+
+/**
+ * GenI the given input (2D) according to a given object
+ */
+PyObject * genicom_inputXY(PyObject *, PyObject *args) {
+	int npoints;
+	PyObject *qx_obj;
+	double *qx;
+	PyObject *qy_obj;
+	double *qy;
+	PyObject *I_out_obj;
+	Py_ssize_t n_out;
+	double *I_out;
+	PyObject *gen_obj;
+
+	if (!PyArg_ParseTuple(args, "OiOOO",  &gen_obj, &npoints, &qx_obj, &qy_obj, &I_out_obj)) return NULL;
+	OUTVECTOR(qx_obj, qx, n_out);
+	OUTVECTOR(qy_obj, qy, n_out);
+	OUTVECTOR(I_out_obj, I_out, n_out);
+
+	// Sanity check
+	//if(n_in!=n_out) return Py_BuildValue("i",-1);
+
+	// Set the array pointers
+	void *temp = PyCObject_AsVoidPtr(gen_obj);
+	GenI* s = static_cast<GenI *>(temp);
+
+	s->genicomXY(npoints, qx, qy, I_out);
+	//return PyCObject_FromVoidPtr(s, del_genicom);
+	return Py_BuildValue("i",1);
+}
+
+/**
+ * GenI the given 1D input according to a given object
+ */
+PyObject * genicom_input(PyObject *, PyObject *args) {
+	int npoints;
+	PyObject *q_obj;
+	double *q;
+	PyObject *I_out_obj;
+	Py_ssize_t n_out;
+	double *I_out;
+	PyObject *gen_obj;
+
+	if (!PyArg_ParseTuple(args, "OiOO",  &gen_obj, &npoints, &q_obj, &I_out_obj)) return NULL;
+	OUTVECTOR(q_obj, q, n_out);
+	OUTVECTOR(I_out_obj, I_out, n_out);
+
+	// Sanity check
+	//if(n_in!=n_out) return Py_BuildValue("i",-1);
+
+	// Set the array pointers
+	void *temp = PyCObject_AsVoidPtr(gen_obj);
+	GenI* s = static_cast<GenI *>(temp);
+
+	s->genicom(npoints, q, I_out);
+	//return PyCObject_FromVoidPtr(s, del_genicom);
+	return Py_BuildValue("i",1);
+}
+
+/**
+ * Define module methods
+ */
+static PyMethodDef module_methods[] = {
+	{"new_GenI", (PyCFunction)new_GenI, METH_VARARGS,
+		  "Create a new GenI object"},
+	{"genicom",(PyCFunction)genicom_input, METH_VARARGS,
+		  "genicom the given 1d input arrays"},
+	{"genicomXY",(PyCFunction)genicom_inputXY, METH_VARARGS,
+		  "genicomXY the given 2d input arrays"},
+    {NULL}
+};
+
+
+#ifndef PyMODINIT_FUNC	/* declarations for DLL import/export */
+#define PyMODINIT_FUNC void
+#endif
+PyMODINIT_FUNC
+initsld2i(void)
+{
+    PyObject* m;
+
+    m = Py_InitModule3("sld2i", module_methods,
+                       "Sld2i module");
+}
Index: src/sas/sascalc/calculator/c_extensions/winFuncs.c
===================================================================
--- src/sas/sascalc/calculator/c_extensions/winFuncs.c	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/winFuncs.c	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,146 @@
+/* winFuncs.c
+   Adding functions missing from windows math library
+   Andrew Jackson, October 2007
+   */
+
+#include <stdio.h>
+#include <math.h>
+
+#include "winFuncs.h"
+
+
+double fmax(double x, double y){
+	//Doesn't exactly match BSD as if one or other value is NaN behaviour is undefined
+	//BSD returns the other value if one is NaN and NaN if both are NaN
+	//probably wouldn't want to rely on that in any case.
+	if (x > y) {
+		//x is greater than y
+		return x;
+	} else if (y > x) {
+		//y is greater than x
+		return y;
+	} else {
+		//equal
+		return x;
+	}
+}
+
+double fmin(double x, double y){
+	//Doesn't exactly match BSD as if one or other value is NaN behaviour is undefined
+	//BSD returns the other value if one is NaN and NaN if both are NaN
+	//probably wouldn't want to rely on that in any case.
+	if (x < y) {
+		//x is less than y
+		return x;
+	} else if (y < x) {
+		//y is less than x
+		return y;
+	} else {
+		//equal
+		return x;
+	}
+}
+
+double trunc(double x){
+	//This is probably slow as hell
+	if (x > 0){
+		//positive - need floor
+		return floor(x);
+	} else if (x < 0) {
+		//negative - need ceiling
+		return ceil(x);
+	} else {
+		//x is zero or infinity, return x
+		return x;
+	}
+
+}
+
+
+/* erf.c  - public domain implementation of error function erf(3m)
+
+reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
+            (New Algorithm handbook in C language) (Gijyutsu hyouron
+            sha, Tokyo, 1991) p.227 [in Japanese]                 */
+
+
+#ifdef _WIN32
+# include <float.h>
+# if !defined __MINGW32__ || defined __NO_ISOCEXT
+#  ifndef isnan
+#   define isnan(x) _isnan(x)
+#  endif
+#  ifndef isinf
+#   define isinf(x) (!_finite(x) && !_isnan(x))
+#  endif
+#  ifndef finite
+#   define finite(x) _finite(x)
+#  endif
+# endif
+#endif
+
+static double q_gamma(double, double, double);
+
+/* Incomplete gamma function
+   1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt  */
+static double p_gamma(double a, double x, double loggamma_a)
+{
+    int k;
+    double result, term, previous;
+
+    if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
+    if (x == 0)     return 0;
+    result = term = exp(a * log(x) - x - loggamma_a) / a;
+    for (k = 1; k < 1000; k++) {
+        term *= x / (a + k);
+        previous = result;  result += term;
+        if (result == previous) return result;
+    }
+    fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__);
+    return result;
+}
+
+/* Incomplete gamma function
+   1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt  */
+static double q_gamma(double a, double x, double loggamma_a)
+{
+    int k;
+    double result, w, temp, previous;
+    double la = 1, lb = 1 + x - a;  /* Laguerre polynomial */
+
+    if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
+    w = exp(a * log(x) - x - loggamma_a);
+    result = w / lb;
+    for (k = 2; k < 1000; k++) {
+        temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
+        la = lb;  lb = temp;
+        w *= (k - 1 - a) / k;
+        temp = w / (la * lb);
+        previous = result;  result += temp;
+        if (result == previous) return result;
+    }
+    fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__);
+    return result;
+}
+
+#define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */
+
+double erf(double x)
+{
+    if (!finite(x)) {
+        if (isnan(x)) return x;      /* erf(NaN)   = NaN   */
+        return (x>0 ? 1.0 : -1.0);   /* erf(+-inf) = +-1.0 */
+    }
+    if (x >= 0) return   p_gamma(0.5, x * x, LOG_PI_OVER_2);
+    else        return - p_gamma(0.5, x * x, LOG_PI_OVER_2);
+}
+
+double erfc(double x)
+{
+    if (!finite(x)) {
+        if (isnan(x)) return x;      /* erfc(NaN)   = NaN      */
+        return (x>0 ? 0.0 : 2.0);    /* erfc(+-inf) = 0.0, 2.0 */
+    }
+    if (x >= 0) return  q_gamma(0.5, x * x, LOG_PI_OVER_2);
+    else        return  1 + p_gamma(0.5, x * x, LOG_PI_OVER_2);
+}
Index: src/sas/sascalc/calculator/c_extensions/winFuncs.h
===================================================================
--- src/sas/sascalc/calculator/c_extensions/winFuncs.h	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
+++ src/sas/sascalc/calculator/c_extensions/winFuncs.h	(revision 9e531f2d45ba748a75200d6849494152c0ad9383)
@@ -0,0 +1,50 @@
+#if !defined(WINFUNCS_H)
+#define WINFUNCS_H
+/* winFuncs.h
+   Definitions for missing math lib functions
+   Andrew Jackson, October 2007
+   */
+
+double fmax(double x, double y);
+//long double fmaxl (long double x, long double y);
+//float fmaxf (float x, float y);
+
+double fmin(double x, double y);
+
+double trunc(double x);
+//long double truncl(long double x);
+//float truncf(float x);
+
+double erf(double x);
+//long double erfl(long double x);
+//float erff (float x);
+//double erfc(double x);
+//long double erfcl(long double x);
+//float erfcf(float x;
+
+
+// Define INFINITY and NAN
+typedef union { unsigned char __c[4]; float __f; } __huge_valf_t;
+
+# if __BYTE_ORDER == __BIG_ENDIAN
+#  define __HUGE_VALF_bytes { 0x7f, 0x80, 0, 0 }
+# endif
+# if __BYTE_ORDER == __LITTLE_ENDIAN
+#  define __HUGE_VALF_bytes { 0, 0, 0x80, 0x7f }
+# endif
+
+static __huge_valf_t __huge_valf = { __HUGE_VALF_bytes };
+# define INFINITY  (__huge_valf.__f)
+
+
+# if __BYTE_ORDER == __BIG_ENDIAN
+#  define __nan_bytes   { 0x7f, 0xc0, 0, 0 }
+# endif
+# if __BYTE_ORDER == __LITTLE_ENDIAN
+#  define __nan_bytes   { 0, 0, 0xc0, 0x7f }
+# endif
+
+static union { unsigned char __c[4]; float __d; } __nan_union;
+# define NAN  (__nan_union.__d)
+
+#endif
