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 a1daf86c0d7d0aba84e29d778ffd6148e933d93a)
+++ src/sas/sascalc/calculator/c_extensions/librefl.c	(revision 144e032af2559d8208e907563f404ab479862dbf)
@@ -112,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;
 
@@ -193,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);
@@ -205,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)));
@@ -213,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 a1daf86c0d7d0aba84e29d778ffd6148e933d93a)
+++ src/sas/sascalc/calculator/c_extensions/sld2i.c	(revision 144e032af2559d8208e907563f404ab479862dbf)
@@ -54,16 +54,20 @@
 	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;
+	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
@@ -77,13 +81,12 @@
 	for(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);
+		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++){
-			//printf("j: %d\n", j);
 			if (this->sldn_val[j]!=0.0
 				||this->mx_val[j]!=0.0
@@ -91,40 +94,40 @@
 				||this->mz_val[j]!=0.0)
 			{
+			    // printf("i,j: %d,%d\n", i,j);
 				//anisotropic
-				temp_fi = cassign(0.0, 0.0);
-				b_sld = cal_msld(0, qx[i], qy[i], this->sldn_val[j],
+				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]);
-				iqr = cassign(0.0, qr);
-				ephase = cplx_exp(iqr);
+				cassign(&iqr, 0.0, qr);
+				cplx_exp(&ephase, iqr);
 
 				//Let's multiply pixel(atomic) volume here
-				ephase = rcmult(this->vol_pix[j], ephase);
+				rcmult(&ephase, this->vol_pix[j], ephase);
 				//up_up
 				if (this->inspin > 0.0 && this->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);
+					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){
-					comp_sld = cassign(b_sld.dd, 0.0);
-					temp_fi = cplx_mult(comp_sld, ephase);
-					sumj_dd = cplx_add(sumj_dd, temp_fi);
+					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){
-					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);
+					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){
-					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);
-				}
-
+					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){
Index: src/sas/sascalc/calculator/sas_gen.py
===================================================================
--- src/sas/sascalc/calculator/sas_gen.py	(revision a1daf86c0d7d0aba84e29d778ffd6148e933d93a)
+++ src/sas/sascalc/calculator/sas_gen.py	(revision 144e032af2559d8208e907563f404ab479862dbf)
@@ -133,11 +133,17 @@
         sldn = copy.deepcopy(self.data_sldn)
         sldn -= self.params['solvent_SLD']
-        model = mod.new_GenI((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'])
+        # **** 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)
