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 a57ae0704487b63fbba07d89eff101d45ec4a470)
Index: src/sas/sascalc/calculator/c_extensions/libfunc.c
===================================================================
--- src/sas/sascalc/calculator/c_extensions/libfunc.c	(revision b6c8abea4ddb688eaf519b2c53bbd8aeebf55489)
+++ src/sas/sascalc/calculator/c_extensions/libfunc.c	(revision 144e032af2559d8208e907563f404ab479862dbf)
@@ -95,5 +95,5 @@
 // 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,
+void cal_msld(polar_sld *p_sld, int isangle, double qx, double qy, double bn,
 				double m01, double mtheta1, double mphi1,
 				double spinfraci, double spinfracf, double spintheta)
@@ -124,91 +124,91 @@
 	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;
+	double uu = sld;
+	double dd = sld;
+	double re_ud = 0.0;
+	double im_ud = 0.0;
+	double re_du = 0.0;
+	double im_du = 0.0;
 
 	//No mag means no further calculation
-	if (isangle>0){
+	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;
+			uu = sqrt(sqrt(in_spin * out_spin)) * uu;
+			dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd;
+		}
+	}
+	else if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){
+			uu = sqrt(sqrt(in_spin * out_spin)) * uu;
+			dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd;
+	} else {
+
+		//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
+		uu -= m_sigma_x;
+		dd += m_sigma_x;
+		re_ud = m_sigma_y;
+		re_du = m_sigma_y;
+		im_ud = m_sigma_z;
+		im_du = -m_sigma_z;
+
+		uu = sqrt(sqrt(in_spin * out_spin)) * uu;
+		dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd;
+
+		re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * re_ud;
+		im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * im_ud;
+		re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * re_du;
+		im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * im_du;
+	}
+	p_sld->uu = uu;
+	p_sld->dd = dd;
+	p_sld->re_ud = re_ud;
+	p_sld->im_ud = im_ud;
+	p_sld->re_du = re_du;
+	p_sld->im_du = im_du;
 }
 
Index: src/sas/sascalc/calculator/c_extensions/libfunc.h
===================================================================
--- src/sas/sascalc/calculator/c_extensions/libfunc.h	(revision b6c8abea4ddb688eaf519b2c53bbd8aeebf55489)
+++ src/sas/sascalc/calculator/c_extensions/libfunc.h	(revision 144e032af2559d8208e907563f404ab479862dbf)
@@ -18,5 +18,5 @@
 //double gamln(double x);
 
-polar_sld cal_msld(int isangle, double qx, double qy, double bn, double m01, double mtheta1, 
+void cal_msld(polar_sld*, int isangle, double qx, double qy, double bn, double m01, double mtheta1, 
 			double mphi1, double spinfraci, double spinfracf, double spintheta);
 
Index: src/sas/sascalc/calculator/c_extensions/librefl.c
===================================================================
--- src/sas/sascalc/calculator/c_extensions/librefl.c	(revision 4c29e4db880ce15b4400de7af0865fa99a32fe16)
+++ src/sas/sascalc/calculator/c_extensions/librefl.c	(revision 144e032af2559d8208e907563f404ab479862dbf)
@@ -7,5 +7,5 @@
 #include <stdio.h>
 #include <stdlib.h>
-#if defined(_MSC_VER)
+#if defined _MSC_VER  || defined __TINYCC__
 #define NEED_ERF
 #endif
@@ -21,5 +21,14 @@
 
 
-#ifdef _WIN32
+#ifdef __TINYCC__
+# ifdef isnan
+#   undef isnan
+# endif
+# ifdef isfinite
+#   undef isfinite
+# endif
+# define isnan(x) (x != x)
+# define isfinite(x) (x != INFINITY && x != -INFINITY)
+#elif defined _WIN32
 # include <float.h>
 # if !defined __MINGW32__ || defined __NO_ISOCEXT
@@ -30,6 +39,6 @@
 #   define isinf(x) (!_finite(x) && !_isnan(x))
 #  endif
-#  ifndef finite
-#   define finite(x) _finite(x)
+#  ifndef isfinite
+#   define isfinite(x) _finite(x)
 #  endif
 # endif
@@ -84,5 +93,5 @@
 double erf(double x)
 {
-    if (!finite(x)) {
+    if (!isfinite(x)) {
         if (isnan(x)) return x;      /* erf(NaN)   = NaN   */
         return (x>0 ? 1.0 : -1.0);   /* erf(+-inf) = +-1.0 */
@@ -94,5 +103,5 @@
 double erfc(double x)
 {
-    if (!finite(x)) {
+    if (!isfinite(x)) {
         if (isnan(x)) return x;      /* erfc(NaN)   = NaN      */
         return (x>0 ? 0.0 : 2.0);    /* erfc(+-inf) = 0.0, 2.0 */
@@ -103,78 +112,54 @@
 #endif // NEED_ERF
 
-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;
+void cassign(Cplx *x, double real, double imag)
+{
+	x->re = real;
+	x->im = imag;
+}
+
+
+void cplx_add(Cplx *z, Cplx x, Cplx y)
+{
+	z->re = x.re + y.re;
+	z->im = x.im + y.im;
+}
+
+void rcmult(Cplx *z, double x, Cplx y)
+{
+	z->re = x*y.re;
+	z->im = x*y.im;
+}
+
+void cplx_sub(Cplx *z, Cplx x, Cplx y)
+{
+	z->re = x.re - y.re;
+	z->im = x.im - y.im;
+}
+
+
+void cplx_mult(Cplx *z, Cplx x, Cplx y)
+{
+	z->re = x.re*y.re - x.im*y.im;
+	z->im = x.re*y.im + x.im*y.re;
+}
+
+void cplx_div(Cplx *z, Cplx x, Cplx y)
+{
+	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);
+}
+
+void cplx_exp(Cplx *z, Cplx b)
+{
 	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;
+	z->re = exp(br)*cos(bi);
+	z->im = exp(br)*sin(bi);
+}
+
+
+void cplx_sqrt(Cplx *c, Cplx z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
+{
 	double zr,zi,x,y,r,w;
 
@@ -184,10 +169,7 @@
 	if (zr==0.0 && zi==0.0)
 	{
-    c.re=0.0;
-	c.im=0.0;
-    return c;
-	}
-	else
-	{
+		c->re=0.0;
+		c->im=0.0;
+	} else {
 		x=fabs(zr);
 		y=fabs(zi);
@@ -196,7 +178,5 @@
 			r=y/x;
 			w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
-		}
-		else
-		{
+		} else {
 			r=x/y;
 			w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
@@ -204,27 +184,24 @@
 		if (zr >=0.0)
 		{
-			c.re=w;
-			c.im=zi/(2.0*w);
+			c->re=w;
+			c->im=zi/(2.0*w);
+		} else {
+			c->im=(zi >= 0) ? w : -w;
+			c->re=zi/(2.0*c->im);
 		}
-		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;
+	}
+}
+
+void cplx_cos(Cplx *z, Cplx b)
+{
+	// cos(b) = (e^bi + e^-bi)/2
+	//        = (e^b.im e^-i bi.re) + e^-b.im e^i b.re)/2
+	//        = (e^b.im cos(-b.re) + e^b.im sin(-b.re) i)/2 + (e^-b.im cos(b.re) + e^-b.im sin(b.re) i)/2
+	//        = e^b.im cos(b.re)/2 - e^b.im sin(b.re)/2 i + 1/e^b.im cos(b.re)/2 + 1/e^b.im sin(b.re)/2 i
+	//        = (e^b.im + 1/e^b.im)/2 cos(b.re) + (-e^b.im + 1/e^b.im)/2 sin(b.re) i
+	//        = cosh(b.im) cos(b.re) - sinh(b.im) sin(b.re) i
+	double exp_b_im = exp(b.im);
+	z->re = 0.5*(+exp_b_im + 1.0/exp_b_im) * cos(b.re);
+	z->im = -0.5*(exp_b_im - 1.0/exp_b_im) * sin(b.re);
 }
 
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 144e032af2559d8208e907563f404ab479862dbf)
@@ -5,30 +5,30 @@
 	double re;
 	double im;
-} complex;
+} Cplx;
 
 typedef struct {
-	complex a;
-	complex b;
-	complex c;
-	complex d;
+	Cplx a;
+	Cplx b;
+	Cplx c;
+	Cplx d;
 } matrix;
 
-complex cassign(double real, double imag);
+void cassign(Cplx*, double real, double imag);
 
-complex cplx_add(complex x,complex y);
+void cplx_add(Cplx*, Cplx x,Cplx y);
 
-complex rcmult(double x,complex y);
+void rcmult(Cplx*, double x,Cplx y);
 
-complex cplx_sub(complex x,complex y);
+void cplx_sub(Cplx*, Cplx x,Cplx y);
 
-complex cplx_mult(complex x,complex y);
+void cplx_mult(Cplx*, Cplx x,Cplx y);
 
-complex cplx_div(complex x,complex y);
+void cplx_div(Cplx*, Cplx x,Cplx y);
 
-complex cplx_exp(complex b);
+void cplx_exp(Cplx*, Cplx b);
 
-complex cplx_sqrt(complex z);
+void cplx_sqrt(Cplx*, Cplx z);
 
-complex cplx_cos(complex b);
+void cplx_cos(Cplx*, Cplx b);
 
 double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r);
Index: src/sas/sascalc/calculator/c_extensions/sld2i.c
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i.c	(revision 144e032af2559d8208e907563f404ab479862dbf)
+++ src/sas/sascalc/calculator/c_extensions/sld2i.c	(revision 144e032af2559d8208e907563f404ab479862dbf)
@@ -0,0 +1,210 @@
+/**
+Computes the (magnetic) scattering form sld (n and m) profile
+ */
+#include <stdio.h>
+#include <math.h>
+#include "sld2i.h"
+#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
+ */
+void initGenI(GenI* this, int is_avg, 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->is_avg = is_avg;
+	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 genicomXY(GenI* this, 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;
+	Cplx iqr;
+	Cplx ephase;
+	Cplx comp_sld;
+
+	Cplx sumj_uu;
+	Cplx sumj_ud;
+	Cplx sumj_du;
+	Cplx sumj_dd;
+	Cplx temp_fi;
+
+	double count = 0.0;
+	int i, j;
+
+	cassign(&iqr, 0.0, 0.0);
+	cassign(&ephase, 0.0, 0.0);
+	cassign(&comp_sld, 0.0, 0.0);
+
+	//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
+
+	//printf("npoints: %d, npix: %d\n", npoints, this->n_pix);
+	for(i=0; i<npoints; i++){
+		//I_out[i] = 0.0;
+		cassign(&sumj_uu, 0.0, 0.0);
+		cassign(&sumj_ud, 0.0, 0.0);
+		cassign(&sumj_du, 0.0, 0.0);
+		cassign(&sumj_dd, 0.0, 0.0);
+		//printf("i: %d\n", i);
+		//q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]);
+
+		for(j=0; j<this->n_pix; j++){
+			if (this->sldn_val[j]!=0.0
+				||this->mx_val[j]!=0.0
+				||this->my_val[j]!=0.0
+				||this->mz_val[j]!=0.0)
+			{
+			    // printf("i,j: %d,%d\n", i,j);
+				//anisotropic
+				cassign(&temp_fi, 0.0, 0.0);
+				cal_msld(&b_sld, 0, qx[i], qy[i], this->sldn_val[j],
+							 this->mx_val[j], this->my_val[j], this->mz_val[j],
+			 				 this->inspin, this->outspin, this->stheta);
+				qr = (qx[i]*this->x_val[j] + qy[i]*this->y_val[j]);
+				cassign(&iqr, 0.0, qr);
+				cplx_exp(&ephase, iqr);
+
+				//Let's multiply pixel(atomic) volume here
+				rcmult(&ephase, this->vol_pix[j], ephase);
+				//up_up
+				if (this->inspin > 0.0 && this->outspin > 0.0){
+					cassign(&comp_sld, b_sld.uu, 0.0);
+					cplx_mult(&temp_fi, comp_sld, ephase);
+					cplx_add(&sumj_uu, sumj_uu, temp_fi);
+				}
+				//down_down
+				if (this->inspin < 1.0 && this->outspin < 1.0){
+					cassign(&comp_sld, b_sld.dd, 0.0);
+					cplx_mult(&temp_fi, comp_sld, ephase);
+					cplx_add(&sumj_dd, sumj_dd, temp_fi);
+				}
+				//up_down
+				if (this->inspin > 0.0 && this->outspin < 1.0){
+					cassign(&comp_sld, b_sld.re_ud, b_sld.im_ud);
+					cplx_mult(&temp_fi, comp_sld, ephase);
+					cplx_add(&sumj_ud, sumj_ud, temp_fi);
+				}
+				//down_up
+				if (this->inspin < 1.0 && this->outspin > 0.0){
+					cassign(&comp_sld, b_sld.re_du, b_sld.im_du);
+					cplx_mult(&temp_fi, comp_sld, ephase);
+					cplx_add(&sumj_du, sumj_du, temp_fi);
+				}
+
+				if (i == 0){
+					count += this->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, this->sldn_val[0],this->mx_val[0], this->my_val[0], this->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 genicom(GenI* this, 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);
+	double qr = 0.0;
+	double sumj;
+	double sld_j = 0.0;
+	double count = 0.0;
+	int i, j, k;
+
+	//Assume that pixel volumes are given in vol_pix in A^3 unit
+	// Loop over q-values and multiply apply matrix
+	for(i=0; i<npoints; i++){
+		sumj =0.0;
+		for(j=0; j<this->n_pix; j++){
+			//Isotropic: Assumes all slds are real (no magnetic)
+			//Also assumes there is no polarization: No dependency on spin
+			if (this->is_avg == 1){
+				// approximation for a spherical symmetric particle
+				qr = sqrt(this->x_val[j]*this->x_val[j]+this->y_val[j]*this->y_val[j]+this->z_val[j]*this->z_val[j])*q[i];
+				if (qr > 0.0){
+					qr = sin(qr) / qr;
+					sumj += this->sldn_val[j] * this->vol_pix[j] * qr;
+				}
+				else{
+					sumj += this->sldn_val[j] * this->vol_pix[j];
+				}
+			}
+			else{
+				//full calculation
+				//pragma omp parallel for
+				for(k=0; k<this->n_pix; k++){
+					sld_j =  this->sldn_val[j] * this->sldn_val[k] * this->vol_pix[j] * this->vol_pix[k];
+					qr = (this->x_val[j]-this->x_val[k])*(this->x_val[j]-this->x_val[k])+
+						      (this->y_val[j]-this->y_val[k])*(this->y_val[j]-this->y_val[k])+
+						      (this->z_val[j]-this->z_val[k])*(this->z_val[j]-this->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 += this->vol_pix[j];
+			}
+		}
+		I_out[i] = sumj;
+		if (this->is_avg == 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.cpp
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i.cpp	(revision b523c0e4093c0c8e82cebd10672473d88adeaee9)
+++ 	(revision )
@@ -1,210 +1,0 @@
-/**
-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(int i=0; i<npoints; i++){
-		sumj =0.0;		
-		for(int 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(int 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.h
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i.h	(revision 54b0650fdadd5a49ef015c2834003a1ea44c39d2)
+++ src/sas/sascalc/calculator/c_extensions/sld2i.h	(revision 54b0650fdadd5a49ef015c2834003a1ea44c39d2)
@@ -0,0 +1,37 @@
+/**
+Computes the (magnetic) scattering form sld (n and m) profile
+ */
+#ifndef SLD2I_CLASS_H
+#define SLD2I_CLASS_H
+
+/**
+ * Base class
+ */
+typedef struct {
+	// vectors
+	int is_avg;
+	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;
+    // spin ratios
+    double inspin;
+    double outspin;
+    double stheta;
+} GenI;
+
+// Constructor
+void initGenI(GenI*, int is_avg, 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
+void genicomXY(GenI*, int npoints, double* qx, double* qy, double *I_out);
+void genicom(GenI*, int npoints, double* q, double *I_out);
+
+#endif
Index: src/sas/sascalc/calculator/c_extensions/sld2i.hh
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i.hh	(revision 4c29e4db880ce15b4400de7af0865fa99a32fe16)
+++ 	(revision )
@@ -1,45 +1,0 @@
-/**
-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
-	void genicomXY(int npoints, double* qx, double* qy, double *I_out);
-	void genicom(int npoints, double* q, double *I_out);
-};
-
-#endif
Index: src/sas/sascalc/calculator/c_extensions/sld2i_module.c
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i_module.c	(revision a1daf86c0d7d0aba84e29d778ffd6148e933d93a)
+++ src/sas/sascalc/calculator/c_extensions/sld2i_module.c	(revision a1daf86c0d7d0aba84e29d778ffd6148e933d93a)
@@ -0,0 +1,206 @@
+/**
+  SLD2I module to perform point and I calculations
+ */
+#include <Python.h>
+#include <stdio.h>
+#include "sld2i.h"
+
+#if PY_MAJOR_VERSION < 3
+typedef void (*PyCapsule_Destructor)(PyObject *);
+typedef void (*PyCObject_Destructor)(void *);
+#define PyCapsule_New(pointer, name, destructor) (PyCObject_FromVoidPtr(pointer, (PyCObject_Destructor)destructor))
+#define PyCapsule_GetPointer(capsule, name) (PyCObject_AsVoidPtr(capsule))
+#endif
+
+
+// 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(PyObject *obj){
+#if PY_MAJOR_VERSION < 3
+	GenI* sld2i = (GenI *)obj;
+#else
+	GenI* sld2i = (GenI *)(PyCapsule_GetPointer(obj, "GenI"));
+#endif
+	PyMem_Free((void *)sld2i);
+}
+
+/**
+ * Create a GenI as a python object by supplying arrays
+ */
+PyObject * new_GenI(PyObject *self, 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, n_y, n_z, n_sld, n_mx, n_my, n_mz, n_vol_pix;
+	int is_avg;
+	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;
+	PyObject *obj;
+	GenI* sld2i;
+
+	//printf("new GenI\n");
+	if (!PyArg_ParseTuple(args, "iOOOOOOOOddd", &is_avg, &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;
+	INVECTOR(x_val_obj, x_val, n_x);
+	INVECTOR(y_val_obj, y_val, n_y);
+	INVECTOR(z_val_obj, z_val, n_z);
+	INVECTOR(sldn_val_obj, sldn_val, n_sld);
+	INVECTOR(mx_val_obj, mx_val, n_mx);
+	INVECTOR(my_val_obj, my_val, n_my);
+	INVECTOR(mz_val_obj, mz_val, n_mz);
+	INVECTOR(vol_pix_obj, vol_pix, n_vol_pix);
+	sld2i = PyMem_Malloc(sizeof(GenI));
+	//printf("sldi:%p\n", sld2i);
+	if (sld2i != NULL) {
+		initGenI(sld2i,is_avg,(int)n_x,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta);
+	}
+	obj = PyCapsule_New(sld2i, "GenI", del_sld2i);
+	//printf("constructed %p\n", obj);
+	return obj;
+}
+
+/**
+ * GenI the given input (2D) according to a given object
+ */
+PyObject * genicom_inputXY(PyObject *self, PyObject *args) {
+	PyObject *gen_obj;
+	PyObject *qx_obj;
+	PyObject *qy_obj;
+	PyObject *I_out_obj;
+	Py_ssize_t n_qx, n_qy, n_out;
+	double *qx;
+	double *qy;
+	double *I_out;
+	GenI* sld2i;
+
+	//printf("in genicom_inputXY\n");
+	if (!PyArg_ParseTuple(args, "OOOO",  &gen_obj, &qx_obj, &qy_obj, &I_out_obj)) return NULL;
+	sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI");
+	INVECTOR(qx_obj, qx, n_qx);
+	INVECTOR(qy_obj, qy, n_qy);
+	OUTVECTOR(I_out_obj, I_out, n_out);
+	//printf("qx, qy, I_out: %d %d %d, %d %d %d\n", qx, qy, I_out, n_qx, n_qy, n_out);
+
+	// Sanity check
+	//if(n_q!=n_out) return Py_BuildValue("i",-1);
+
+	genicomXY(sld2i, (int)n_qx, qx, qy, I_out);
+	//printf("done calc\n");
+	//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 *self, PyObject *args) {
+	PyObject *gen_obj;
+	PyObject *q_obj;
+	PyObject *I_out_obj;
+	Py_ssize_t n_q, n_out;
+	double *q;
+	double *I_out;
+	GenI *sld2i;
+
+	if (!PyArg_ParseTuple(args, "OOO",  &gen_obj, &q_obj, &I_out_obj)) return NULL;
+	sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI");
+	INVECTOR(q_obj, q, n_q);
+	OUTVECTOR(I_out_obj, I_out, n_out);
+
+	// Sanity check
+	//if (n_q!=n_out) return Py_BuildValue("i",-1);
+
+	genicom(sld2i, (int)n_q, q, I_out);
+	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}
+};
+
+#define MODULE_DOC "Sld2i C Library"
+#define MODULE_NAME "sld2i"
+#define MODULE_INIT2 initsld2i
+#define MODULE_INIT3 PyInit_sld2i
+#define MODULE_METHODS module_methods
+
+/* ==== boilerplate python 2/3 interface bootstrap ==== */
+
+
+#if defined(WIN32) && !defined(__MINGW32__)
+    #define DLL_EXPORT __declspec(dllexport)
+#else
+    #define DLL_EXPORT
+#endif
+
+#if PY_MAJOR_VERSION >= 3
+
+  DLL_EXPORT PyMODINIT_FUNC MODULE_INIT3(void)
+  {
+    static struct PyModuleDef moduledef = {
+      PyModuleDef_HEAD_INIT,
+      MODULE_NAME,         /* m_name */
+      MODULE_DOC,          /* m_doc */
+      -1,                  /* m_size */
+      MODULE_METHODS,      /* m_methods */
+      NULL,                /* m_reload */
+      NULL,                /* m_traverse */
+      NULL,                /* m_clear */
+      NULL,                /* m_free */
+    };
+    return PyModule_Create(&moduledef);
+  }
+
+#else /* !PY_MAJOR_VERSION >= 3 */
+
+  DLL_EXPORT PyMODINIT_FUNC MODULE_INIT2(void)
+  {
+    Py_InitModule4(MODULE_NAME,
+		 MODULE_METHODS,
+		 MODULE_DOC,
+		 0,
+		 PYTHON_API_VERSION
+		 );
+  }
+
+#endif /* !PY_MAJOR_VERSION >= 3 */
Index: src/sas/sascalc/calculator/c_extensions/sld2i_module.cpp
===================================================================
--- src/sas/sascalc/calculator/c_extensions/sld2i_module.cpp	(revision 1d014cbb271bb612a92607b84c582f19f50b2b11)
+++ 	(revision )
@@ -1,202 +1,0 @@
-/**
-  SLD2I module to perform point and I calculations
- */
-#include <Python.h>
-#include <stdio.h>
-#include <sld2i.hh>
-
-#if PY_MAJOR_VERSION < 3
-typedef void (*PyCapsule_Destructor)(PyObject *);
-typedef void (*PyCObject_Destructor)(void *);
-#define PyCapsule_New(pointer, name, destructor) (PyCObject_FromVoidPtr(pointer, (PyCObject_Destructor)destructor))
-#define PyCapsule_GetPointer(capsule, name) (PyCObject_AsVoidPtr(capsule))
-#endif
-
-
-// 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(PyObject *obj){
-	GenI* sld2i = static_cast<GenI *>(PyCapsule_GetPointer(obj, "GenI"));
-	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 PyCapsule_New(sld2i, "GenI", 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 = PyCapsule_GetPointer(gen_obj, "GenI");
-	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 = PyCapsule_GetPointer(gen_obj, "GenI");
-	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}
-};
-
-#define MODULE_DOC "Sld2i C Library"
-#define MODULE_NAME "sld2i"
-#define MODULE_INIT2 initsld2i
-#define MODULE_INIT3 PyInit_sld2i
-#define MODULE_METHODS module_methods
-
-/* ==== boilerplate python 2/3 interface bootstrap ==== */
-
-
-#if defined(WIN32) && !defined(__MINGW32__)
-    #define DLL_EXPORT __declspec(dllexport)
-#else
-    #define DLL_EXPORT
-#endif
-
-#if PY_MAJOR_VERSION >= 3
-
-  DLL_EXPORT PyMODINIT_FUNC MODULE_INIT3(void)
-  {
-    static struct PyModuleDef moduledef = {
-      PyModuleDef_HEAD_INIT,
-      MODULE_NAME,         /* m_name */
-      MODULE_DOC,          /* m_doc */
-      -1,                  /* m_size */
-      MODULE_METHODS,      /* m_methods */
-      NULL,                /* m_reload */
-      NULL,                /* m_traverse */
-      NULL,                /* m_clear */
-      NULL,                /* m_free */
-    };
-    return PyModule_Create(&moduledef);
-  }
-
-#else /* !PY_MAJOR_VERSION >= 3 */
-
-  DLL_EXPORT PyMODINIT_FUNC MODULE_INIT2(void)
-  {
-    Py_InitModule4(MODULE_NAME,
-		 MODULE_METHODS,
-		 MODULE_DOC,
-		 0,
-		 PYTHON_API_VERSION
-		 );
-  }
-
-#endif /* !PY_MAJOR_VERSION >= 3 */
Index: src/sas/sascalc/calculator/sas_gen.py
===================================================================
--- src/sas/sascalc/calculator/sas_gen.py	(revision b58265c3c76f95c80a659bfb713f1605c80532a6)
+++ src/sas/sascalc/calculator/sas_gen.py	(revision 144e032af2559d8208e907563f404ab479862dbf)
@@ -118,5 +118,5 @@
         self.is_avg = is_avg
 
-    def _gen(self, x, y, i):
+    def _gen(self, qx, qy):
         """
         Evaluate the function
@@ -129,24 +129,35 @@
         pos_y = self.data_y
         pos_z = self.data_z
-        len_x = len(pos_x)
         if self.is_avg is None:
-            len_x *= -1
             pos_x, pos_y, pos_z = transform_center(pos_x, pos_y, pos_z)
-        len_q = len(x)
         sldn = copy.deepcopy(self.data_sldn)
         sldn -= self.params['solvent_SLD']
-        model = mod.new_GenI(len_x, pos_x, pos_y, pos_z,
-                             sldn, self.data_mx, self.data_my,
-                             self.data_mz, self.data_vol,
-                             self.params['Up_frac_in'],
-                             self.params['Up_frac_out'],
-                             self.params['Up_theta'])
-        if y == []:
-            mod.genicom(model, len_q, x, i)
-        else:
-            mod.genicomXY(model, len_q, x, y, i)
+        # **** WARNING **** new_GenI holds pointers to numpy vectors
+        # be sure that they are contiguous double precision arrays and make 
+        # sure the GC doesn't eat them before genicom is called.
+        # TODO: rewrite so that the parameters are passed directly to genicom
+        args = (
+            (1 if self.is_avg else 0),
+            pos_x, pos_y, pos_z,
+            sldn, self.data_mx, self.data_my,
+            self.data_mz, self.data_vol,
+            self.params['Up_frac_in'],
+            self.params['Up_frac_out'],
+            self.params['Up_theta'])
+        model = mod.new_GenI(*args)
+        if len(qy):
+            qx, qy = _vec(qx), _vec(qy)
+            I_out = np.empty_like(qx)
+            #print("npoints", qx.shape, "npixels", pos_x.shape)
+            mod.genicomXY(model, qx, qy, I_out)
+            #print("I_out after", I_out)
+        else:
+            qx = _vec(qx)
+            I_out = np.empty_like(qx)
+            mod.genicom(model, qx, I_out)
         vol_correction = self.data_total_volume / self.params['total_volume']
-        return  self.params['scale'] * vol_correction * i + \
-                        self.params['background']
+        result = (self.params['scale'] * vol_correction * I_out
+                  + self.params['background'])
+        return result
 
     def set_sld_data(self, sld_data=None):
@@ -156,12 +167,12 @@
         self.sld_data = sld_data
         self.data_pos_unit = sld_data.pos_unit
-        self.data_x = sld_data.pos_x
-        self.data_y = sld_data.pos_y
-        self.data_z = sld_data.pos_z
-        self.data_sldn = sld_data.sld_n
-        self.data_mx = sld_data.sld_mx
-        self.data_my = sld_data.sld_my
-        self.data_mz = sld_data.sld_mz
-        self.data_vol = sld_data.vol_pix
+        self.data_x = _vec(sld_data.pos_x)
+        self.data_y = _vec(sld_data.pos_y)
+        self.data_z = _vec(sld_data.pos_z)
+        self.data_sldn = _vec(sld_data.sld_n)
+        self.data_mx = _vec(sld_data.sld_mx)
+        self.data_my = _vec(sld_data.sld_my)
+        self.data_mz = _vec(sld_data.sld_mz)
+        self.data_vol = _vec(sld_data.vol_pix)
         self.data_total_volume = sum(sld_data.vol_pix)
         self.params['total_volume'] = sum(sld_data.vol_pix)
@@ -180,11 +191,10 @@
         :return: (I value)
         """
-        if x.__class__.__name__ == 'list':
+        if isinstance(x, list):
             if len(x[1]) > 0:
                 msg = "Not a 1D."
                 raise ValueError(msg)
-            i_out = np.zeros_like(x[0])
             # 1D I is found at y =0 in the 2D pattern
-            out = self._gen(x[0], [], i_out)
+            out = self._gen(x[0], [])
             return out
         else:
@@ -199,8 +209,6 @@
         :Use this runXY() for the computation
         """
-        if x.__class__.__name__ == 'list':
-            i_out = np.zeros_like(x[0])
-            out = self._gen(x[0], x[1], i_out)
-            return out
+        if isinstance(x, list):
+            return self._gen(x[0], x[1])
         else:
             msg = "Q must be given as list of qx's and qy's"
@@ -214,14 +222,13 @@
                       where qx,qy are 1D ndarrays (for 2D).
         """
-        if qdist.__class__.__name__ == 'list':
-            if len(qdist[1]) < 1:
-                out = self.run(qdist)
-            else:
-                out = self.runXY(qdist)
-            return out
+        if isinstance(qdist, list):
+            return self.run(qdist) if len(qdist[1]) < 1 else self.runXY(qdist)
         else:
             mesg = "evalDistribution is expecting an ndarray of "
             mesg += "a list [qx,qy] where qx,qy are arrays."
             raise RuntimeError(mesg)
+
+def _vec(v):
+    return np.ascontiguousarray(v, 'd')
 
 class OMF2SLD(object):
@@ -1041,4 +1048,10 @@
         self.line_z = line_z
 
+def _get_data_path(*path_parts):
+    from os.path import realpath, join as joinpath, dirname, abspath
+    # in sas/sascalc/calculator;  want sas/sasview/test
+    return joinpath(dirname(realpath(__file__)),
+                    '..', '..', 'sasview', 'test', *path_parts)
+
 def test_load():
     """
@@ -1046,18 +1059,12 @@
     """
     from mpl_toolkits.mplot3d import Axes3D
-    current_dir = os.path.abspath(os.path.curdir)
-    print(current_dir)
-    for i in range(6):
-        current_dir, _ = os.path.split(current_dir)
-        tfile = os.path.join(current_dir, "test", "CoreXY_ShellZ.txt")
-        ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf")
-        if os.path.isfile(tfile):
-            tfpath = tfile
-            ofpath = ofile
-            break
+    tfpath = _get_data_path("1d_data", "CoreXY_ShellZ.txt")
+    ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf")
+    if not os.path.isfile(tfpath) or not os.path.isfile(ofpath):
+        raise ValueError("file(s) not found: %r, %r"%(tfpath, ofpath))
     reader = SLDReader()
     oreader = OMFReader()
-    output = decode(reader.read(tfpath))
-    ooutput = decode(oreader.read(ofpath))
+    output = reader.read(tfpath)
+    ooutput = oreader.read(ofpath)
     foutput = OMF2SLD()
     foutput.set_data(ooutput)
@@ -1088,30 +1095,34 @@
     plt.show()
 
+def test_save():
+    ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf")
+    if not os.path.isfile(ofpath):
+        raise ValueError("file(s) not found: %r"%(ofpath,))
+    oreader = OMFReader()
+    omfdata = oreader.read(ofpath)
+    omf2sld = OMF2SLD()
+    omf2sld.set_data(omfdata)
+    writer = SLDReader()
+    writer.write("out.txt", omf2sld.output)
+
 def test():
     """
         Test code
     """
-    current_dir = os.path.abspath(os.path.curdir)
-    for i in range(3):
-        current_dir, _ = os.path.split(current_dir)
-        ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf")
-        if os.path.isfile(ofile):
-            ofpath = ofile
-            break
+    ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf")
+    if not os.path.isfile(ofpath):
+        raise ValueError("file(s) not found: %r"%(ofpath,))
     oreader = OMFReader()
-    ooutput = decode(oreader.read(ofpath))
-    foutput = OMF2SLD()
-    foutput.set_data(ooutput)
-    writer = SLDReader()
-    writer.write(os.path.join(os.path.dirname(ofpath), "out.txt"),
-                 foutput.output)
+    omfdata = oreader.read(ofpath)
+    omf2sld = OMF2SLD()
+    omf2sld.set_data(omfdata)
     model = GenSAS()
-    model.set_sld_data(foutput.output)
-    x = np.arange(1000)/10000. + 1e-5
-    y = np.arange(1000)/10000. + 1e-5
-    i = np.zeros(1000)
-    model.runXY([x, y, i])
+    model.set_sld_data(omf2sld.output)
+    x = np.linspace(0, 0.1, 11)[1:]
+    return model.runXY([x, x])
 
 if __name__ == "__main__":
+    #test_load()
+    #test_save()
+    #print(test())
     test()
-    test_load()
