Index: sansmodels/src/libigor/libSphere.c
===================================================================
--- sansmodels/src/libigor/libSphere.c	(revision 975ec8e8211549489759b5a5ffe3b60607de065b)
+++ sansmodels/src/libigor/libSphere.c	(revision 6e93a029d775ebbe413e7c484e2b69f68be665ab)
@@ -13,5 +13,5 @@
 SphereForm(double dp[], double q)
 {
-	double scale,radius,delrho,bkg;		//my local names
+	double scale,radius,delrho,bkg,sldSph,sldSolv;		//my local names
 	double bes,f,vol,f2,pi,qr;
 
@@ -19,19 +19,22 @@
 	scale = dp[0];
 	radius = dp[1];
-	delrho = dp[2];
-	bkg = dp[3];
-
+	sldSph = dp[2];
+	sldSolv = dp[3];
+	bkg = dp[4];
+	
+	delrho = sldSph - sldSolv;
+	//handle qr==0 separately
 	qr = q*radius;
-	if(qr == 0){
+	if(qr == 0.0){
 		bes = 1.0;
 	}else{
 		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
 	}
-
+	
 	vol = 4.0*pi/3.0*radius*radius*radius;
 	f = vol*bes*delrho;		// [=] A-1
 							// normalize to single particle volume, convert to 1/cm
 	f2 = f * f / vol * 1.0e8;		// [=] 1/cm
-
+	
 	return(scale*f2+bkg);	//scale, and add in the background
 }
@@ -44,8 +47,8 @@
 	double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;		//my local names
 	double bes,f,vol,qr,contr,f2;
-
+	
 	pi = 4.0*atan(1.0);
 	x=q;
-
+	
 	scale = dp[0];
 	rcore = dp[1];
@@ -58,5 +61,5 @@
 	qr=x*rcore;
 	contr = rhocore-rhoshel;
-	if(qr == 0){
+	if(qr == 0.0){
 		bes = 1.0;
 	}else{
@@ -68,5 +71,5 @@
 	qr=x*(rcore+thick);
 	contr = rhoshel-rhosolv;
-	if(qr == 0){
+	if(qr == 0.0){
 		bes = 1.0;
 	}else{
@@ -78,10 +81,10 @@
 	// normalize to particle volume and rescale from [A-1] to [cm-1]
 	f2 = f*f/vol*1.0e8;
-
+	
 	//scale if desired
 	f2 *= scale;
 	// then add in the background
 	f2 += bkg;
-
+	
 	return(f2);
 }
@@ -97,5 +100,5 @@
 	pi = 4.0*atan(1.0);
 	x= q;
-
+	
 	scale = dp[0];
 	rcore = dp[1];
@@ -118,5 +121,5 @@
 	qr=x*(rcore+thick);
 	contr = rhoshel-rhosolv;
-	if(qr == 0){
+	if(qr == 0.0){
 		bes = 1.0;
 	}else{
@@ -130,10 +133,10 @@
 	vol = 4.0*pi/3.0*(pow((rcore+thick),3)-pow(rcore,3));
 	f2 = f*f/vol*1.0e8;
-
+	
 	//scale if desired
 	f2 *= scale;
 	// then add in the background
 	f2 += bkg;
-
+	
 	return(f2);
 }
@@ -152,5 +155,5 @@
 	double zp3,vpoly;
 	double s2y, arg1, arg2, arg3, drh1, drh2;
-
+	
 	pi = 4.0*atan(1.0);
 	qq= q;
@@ -162,9 +165,9 @@
 	drho2 = dp[5]-dp[6];		//shell-solvent
 	bkg = dp[7];
-
-	zz = (1.0/sig)*(1.0/sig) - 1.0;
-
+	
+	zz = (1.0/sig)*(1.0/sig) - 1.0;   
+	
 	h=qq;
-
+	
 	drh1 = drho1;
 	drh2 = drho2;
@@ -174,6 +177,6 @@
 	zp3 = zz + 3.;
 	vpoly = 4*pi/3*zp3*zp2/zp1/zp1*pow((corrad+del),3);
-
-
+	
+	
 	// remember that h is the passed in value of q for the calculation
 	y = h *del;
@@ -196,10 +199,10 @@
 	c8 = c4 - .5 + g * cy - g * .5 * g * (y * y + 1.) * c2y;
 	c9 = g * sy * (1. - g * cy);
-
+	
 	tb = log(zp1 * zp1 / (zp1 * zp1 + x * 4. * x));
 	tb1 = exp(zp1 * .5 * tb);
 	tb2 = exp(zp2 * .5 * tb);
 	tb3 = exp(zp3 * .5 * tb);
-
+	
 	t1 = c1 + c2 * x + c3 * x * x * zp2 / zp1;
 	t2 = tb1 * (c4 * cos(arg1) + c7 * sin(arg1));
@@ -216,5 +219,5 @@
 	// then add in the background
 	form += bkg;
-
+	
 	return(form);
 }
@@ -234,20 +237,22 @@
 	double capl,capl1,capmu,capmu1,r3,pq;
 	double ka2,r1,heff;
-	double hh,k;
-
+	double hh,k,slds,sld;
+	
 	pi = 4.0*atan(1.0);
 	k= q;
-
+	
 	rad = dp[0];		// radius (A)
 	z2 = dp[1];		//polydispersity (0<z2<1)
 	phi = dp[2];		// volume fraction (0<phi<1)
-	cont = dp[3]*1.0e4;		// contrast (odd units)
-	bkg = dp[4];
+	slds = dp[3];
+	sld = dp[4];
+	cont = (slds - sld)*1.0e4;		// contrast (odd units)
+	bkg = dp[5];
 	sigma = 2*rad;
-
+	
 	zz=1.0/(z2*z2)-1.0;
 	bb = sigma/(zz+1.0);
 	cc = zz+1.0;
-
+	
 	//*c   Compute the number density by <r-cubed>, not <r> cubed*/
 	r1 = sigma/2.0;
@@ -281,30 +286,30 @@
 	d5=pow((pi/capd),2)*((rho/k)*(chi-1.0)+0.5*k*t2);
 	d6=pow((pi/capd),2)*(rho/k)*(chi2-s2);
-
+	
 	e1=pow((pi/capd),2)*pow((rho/k/k),2)*((chi-1.0)*(chi2-s2)-(chi1-s1)*(chi1-s1)-(k*s1-pp)*(k*s3-p2)+pow((k*s2-p1),2));
 	ee=1.0-(2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(k*s1-pp)-(2.0*pi/capd)*rho/k/k*((chi1-s1)+(0.25*pi*t2/capd)*(chi2-s2))-e1;
 	y1=pow((pi/capd),2)*pow((rho/k/k),2)*((k*s1-pp)*(chi2-s2)-2.0*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1.0));
-	yy = (2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1.0)-(2.0*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1;
-
+	yy = (2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1.0)-(2.0*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1;    
+	
 	capl=2.0*pi*cont*rho/k/k/k*(pp-0.5*k*(s1+chi1));
 	capl1=2.0*pi*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2));
 	capmu=2.0*pi*cont*rho/k/k/k*(1.0-chi-0.5*k*p1);
 	capmu1=2.0*pi*cont*rho/k/k/k*(s1-chi1-0.5*k*p2);
-
+	
 	h1=capl*(capl*(yy*d1-ee*d6)+capl1*(yy*d2-ee*d4)+capmu*(ee*d1+yy*d6)+capmu1*(ee*d2+yy*d4));
 	h2=capl1*(capl*(yy*d2-ee*d4)+capl1*(yy*d3-ee*d5)+capmu*(ee*d2+yy*d4)+capmu1*(ee*d3+yy*d5));
 	h3=capmu*(capl*(ee*d1+yy*d6)+capl1*(ee*d2+yy*d4)+capmu*(ee*d6-yy*d1)+capmu1*(ee*d4-yy*d2));
 	h4=capmu1*(capl*(ee*d2+yy*d4)+capl1*(ee*d3+yy*d5)+capmu*(ee*d4-yy*d2)+capmu1*(ee*d5-yy*d3));
-
+	
 	//*  This part computes the second integral in equation (1) of the paper.*/
-
+	
 	hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(ee*ee+yy*yy));
-
+	
 	//*  This part computes the first integral in equation (1).  It also
 	// generates the KC approximated effective structure factor.*/
-
+	
 	pq=4.0*pi*cont*(sin(k*sigma/2.0)-0.5*k*sigma*cos(k*sigma/2.0));
 	hint2=8.0*pi*pi*rho*cont*cont/(k*k*k*k*k*k)*(1.0-chi-k*p1+0.25*k*k*(s2+chi2));
-
+	
 	ka=k*(sigma/2.0);
 	//
@@ -314,10 +319,10 @@
 	ka2=ka*ka;
 	//*
-	//  heff is PY analytical solution for intensity divided by the
+	//  heff is PY analytical solution for intensity divided by the 
 	//   form factor.  happ is the KC approximated effective S(q)
-
+	 
 	 //*******************
 	 //  add in the background then return the intensity value
-
+	 
 	 return(hh+bkg);	//scale, and add in the background
 }
@@ -332,7 +337,7 @@
 	double scale,ravg,pd,bkg,rho,rhos;		//my local names
 	double scale2,ravg2,pd2,rho2;		//my local names
-
+	
 	x= q;
-
+	
 	scale = dp[0];
 	ravg = dp[1];
@@ -345,10 +350,10 @@
 	rhos = dp[8];
 	bkg = dp[9];
-
+	
 	pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x);
 	pq += SchulzSphere_Fn( scale2,  ravg2,  pd2,  rho2,  rhos, x);
 	// add in the background
 	pq += bkg;
-
+	
 	return (pq);
 }
@@ -361,7 +366,7 @@
 	double x,pq;
 	double scale,ravg,pd,bkg,rho,rhos;		//my local names
-
+	
 	x= q;
-
+	
 	scale = dp[0];
 	ravg = dp[1];
@@ -373,5 +378,5 @@
 	// add in the background
 	pq += bkg;
-
+	
 	return(pq);
 }
@@ -384,9 +389,10 @@
 	double aa,at1,at2,rt1,rt2,rt3,t1,t2,t3;
 	double v1,v2,v3,g1,pq,pi,delrho,zz;
-
+	double i_zero,Rg2,zp8;
+	
 	pi = 4.0*atan(1.0);
 	delrho = rho-rhos;
-	zz = (1.0/pd)*(1.0/pd) - 1.0;
-
+	zz = (1.0/pd)*(1.0/pd) - 1.0;   
+	
 	zp1 = zz + 1.0;
 	zp2 = zz + 2.0;
@@ -397,6 +403,19 @@
 	zp7 = zz + 7.0;
 	//
-	aa = (zz+1)/x/ravg;
-
+	
+	//small QR limit - use Guinier approx
+	zp8 = zz+8.0;
+	if(x*ravg < 0.1) {
+		i_zero = scale*delrho*delrho*1.e8*4.*pi/3.*pow(ravg,3);
+		i_zero *= zp6*zp5*zp4/zp1/zp1/zp1;		//6th moment / 3rd moment
+		Rg2 = 3.*zp8*zp7/5./(zp1*zp1)*ravg*ravg;
+		pq = i_zero*exp(-x*x*Rg2/3.);
+		//pq += bkg;		//unlike the Igor code, the backgorund is added in the wrapper (above)
+		return(pq);
+	}
+	//
+
+	aa = (zz+1.0)/x/ravg;
+	
 	at1 = atan(1.0/aa);
 	at2 = atan(2.0/aa);
@@ -404,5 +423,5 @@
 	//  calculations are performed to avoid  large # errors
 	// - trick is to propogate the a^(z+7) term through the g1
-	//
+	// 
 	t1 = zp7*log10(aa) - zp1/2.0*log10(aa*aa+4.0);
 	t2 = zp7*log10(aa) - zp3/2.0*log10(aa*aa+4.0);
@@ -416,23 +435,23 @@
 	v3 = -2.0*zp1*rt3*sin(zp2*at2);
 	g1 = (v1+v2+v3);
-
+	
 	pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg);
-	pq = pow(10,pq)*8*pi*pi*delrho*delrho;
-
+	pq = pow(10,pq)*8.0*pi*pi*delrho*delrho;
+	
 	//
-	// beta factor is not used here, but could be for the
+	// beta factor is not used here, but could be for the 
 	// decoupling approximation
-	//
+	// 
 	//	g11 = g1
 	//	gd = -zp7*log(aa)
 	//	g1 = log(g11) + gd
-	//
+	//                       
 	//	t1 = zp1*at1
 	//	t2 = zp2*at1
 	//	g2 = sin( t1 ) - zp1/sqrt(aa*aa+1)*cos( t2 )
 	//	g22 = g2*g2
-	//	beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22)
+	//	beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22) 
 	//	beta = 2*alog(beta)
-
+	
 	//re-normalize by the average volume
 	vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*ravg*ravg*ravg;
@@ -440,5 +459,5 @@
 	//scale, convert to cm^-1
 	pq *= scale * 1.0e8;
-
+    
     return(pq);
 }
@@ -451,45 +470,65 @@
 	double pi,x;
 	double scale,rad,pd,cont,bkg;		//my local names
-	double inten,h1,qw,qr,width,sig,averad3;
-
+	double inten,h1,qw,qr,width,sig,averad3,Rg2,slds,sld;
+	
 	pi = 4.0*atan(1.0);
 	x= q;
-
+	
 	scale = dp[0];
 	rad = dp[1];		// radius (A)
 	pd = dp[2];		//polydispersity of rectangular distribution
-	cont = dp[3];		// contrast (A^-2)
-	bkg = dp[4];
-
+	slds = dp[3];
+	sld = dp[4];
+	cont = slds - sld;		// contrast (A^-2)
+	bkg = dp[5];
+	
 	// as usual, poly = sig/ravg
 	// for the rectangular distribution, sig = width/sqrt(3)
 	// width is the HALF- WIDTH of the rectangular distrubution
-
+	
 	sig = pd*rad;
 	width = sqrt(3.0)*sig;
-
+	
 	//x is the q-value
 	qw = x*width;
 	qr = x*rad;
+	
+	// as for the numerical inatabilities at low QR, the function is calculating the sines and cosines
+	// just fine - the problem seems to be that the 
+	// leading terms nearly cancel with the last term (the -6*qr... term), to within machine
+	// precision - the difference is on the order of 10^-20
+	// so just use the limiting Guiner value
+	if(qr<0.1) {
+		h1 = scale*cont*cont*1.e8*4.*pi/3.0*pow(rad,3);
+		h1 *= 	(1. + 15.*pow(pd,2) + 27.*pow(pd,4) +27./7.*pow(pd,6) );				//6th moment
+		h1 /= (1.+3.*pd*pd);	//3rd moment
+		Rg2 = 3.0/5.0*rad*rad*( 1.+28.*pow(pd,2)+126.*pow(pd,4)+108.*pow(pd,6)+27.*pow(pd,8) );
+		Rg2 /= (1.+15.*pow(pd,2)+27.*pow(pd,4)+27./7.*pow(pd,6));
+		h1 *= exp(-1./3.*Rg2*x*x);
+		h1 += bkg;
+		return(h1);
+	}
+	
+	// normal calculation
 	h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0;
-	h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw);
-	h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw);
-	h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw);
-	h1 += qw*qr*sin(2*qr)*cos(2*qw);
+	h1 -= 5.0/2.0*cos(2.0*qr)*sin(qw)*cos(qw);
+	h1 += 0.5*qr*qr*cos(2.0*qr)*sin(2.0*qw);
+	h1 += 0.5*qw*qw*cos(2.0*qr)*sin(2.0*qw);
+	h1 += qw*qr*sin(2.0*qr)*cos(2.0*qw);
 	h1 += 3.0*qw*(cos(qr)*cos(qw))*(cos(qr)*cos(qw));
 	h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw));
 	h1 -= 6.0*qr*cos(qr)*sin(qr)*cos(qw)*sin(qw);
-
+	
 	// calculate P(q) = <f^2>
 	inten = 8.0*pi*pi*cont*cont/width/pow(x,7)*h1;
-
+	
 	// beta(q) would be calculated as 2/width/x/h1*h2*h2
-	// with
+	// with 
 	// h2 = 2*sin(x*rad)*sin(x*width)-x*rad*cos(x*rad)*sin(x*width)-x*width*sin(x*rad)*cos(x*width)
-
+	
 	// normalize to the average volume
 	// <R^3> = ravg^3*(1+3*pd^2)
 	// or... "zf"  = (1 + 3*p^2), which will be greater than one
-
+	
 	averad3 =  rad*rad*rad*(1.0+3.0*pd*pd);
 	inten /= 4.0*pi/3.0*averad3;
@@ -500,5 +539,5 @@
 	// then add in the background
 	inten += bkg;
-
+	
 	return(inten);
 }
@@ -514,8 +553,8 @@
 	double va,vb,zi,yy,summ,inten;
 	int nord=20,ii;
-
+	
 	pi = 4.0*atan(1.0);
 	x= q;
-
+	
 	scale=dp[0];
 	rad=dp[1];
@@ -526,11 +565,11 @@
 	delrho=rho-rhos;
 	bkg=dp[5];
-
+	
 	va = -4.0*sig + rad;
-	if (va<0) {
-		va=0;		//to avoid numerical error when  va<0 (-ve q-value)
+	if (va<0.0) {
+		va=0.0;		//to avoid numerical error when  va<0 (-ve q-value)
 	}
 	vb = 4.0*sig +rad;
-
+	
 	summ = 0.0;		// initialize integral
 	for(ii=0;ii<nord;ii+=1) {
@@ -542,17 +581,17 @@
 		yy *= yy;
 		yy *= Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi);
-
+		
 		summ += yy;		//add to the running total of the quadrature
    	}
 	// calculate value of integral to return
 	inten = (vb-va)/2.0*summ;
-
+	
 	//re-normalize by polydisperse sphere volume
 	inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
-
+	
 	inten *= 1.0e8;
 	inten *= scale;
 	inten += bkg;
-
+	
     return(inten);	//scale, and add in the background
 }
@@ -567,8 +606,8 @@
 	double va,vb,zi,yy,summ,inten;
 	int nord=76,ii;
-
+	
 	pi = 4.0*atan(1.0);
 	x= q;
-
+	
 	scale=dp[0];
 	rad=dp[1];		//rad is the median radius
@@ -579,13 +618,13 @@
 	delrho=rho-rhos;
 	bkg=dp[5];
-
+	
 	va = -3.5*sig + mu;
 	va = exp(va);
-	if (va<0) {
-		va=0;		//to avoid numerical error when  va<0 (-ve q-value)
+	if (va<0.0) {
+		va=0.0;		//to avoid numerical error when  va<0 (-ve q-value)
 	}
 	vb = 3.5*sig*(1.0+sig) +mu;
 	vb = exp(vb);
-
+	
 	summ = 0.0;		// initialize integral
 	for(ii=0;ii<nord;ii+=1) {
@@ -597,18 +636,18 @@
 		yy *= yy;
 		yy *= Gauss76Wt[ii] *  LogNormal_distr(sig,mu,zi);
-
+		
 		summ += yy;		//add to the running total of the quadrature
    	}
 	// calculate value of integral to return
 	inten = (vb-va)/2.0*summ;
-
+	
 	//re-normalize by polydisperse sphere volume
 	r3 = exp(3.0*mu + 9.0/2.0*sig*sig);		// <R^3> directly
 	inten /= (4.0*pi/3.0*r3);		//polydisperse volume
-
+	
 	inten *= 1.0e8;
 	inten *= scale;
 	inten += bkg;
-
+	
 	return(inten);
 }
@@ -616,9 +655,9 @@
 static double
 LogNormal_distr(double sig, double mu, double pt)
-{
+{	
 	double retval,pi;
-
-	pi = 4.0*atan(1.0);
-	retval = (1/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );
+	
+	pi = 4.0*atan(1.0);
+	retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );
 	return(retval);
 }
@@ -626,7 +665,7 @@
 static double
 Gauss_distr(double sig, double avg, double pt)
-{
+{	
 	double retval,Pi;
-
+	
 	Pi = 4.0*atan(1.0);
 	retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0);
@@ -644,8 +683,8 @@
 	double sld1,sld2,sld3,zp1,zp2,zp3,vpoly;
 	double pi43,c1,c2,form,volume,arg1,arg2;
-
+	
 	pi = 4.0*atan(1.0);
 	x= q;
-
+	
 	scale = dp[0];
 	corrad = dp[1];
@@ -656,6 +695,6 @@
 	sld3 = dp[6];
 	bkg = dp[7];
-
-	zz = (1.0/sig)*(1.0/sig) - 1.0;
+	
+	zz = (1.0/sig)*(1.0/sig) - 1.0;   
 	shlrad = corrad + thick;
 	drho1 = sld1-sld2;		//core-shell
@@ -665,10 +704,10 @@
 	zp3 = zz + 3.;
 	vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((corrad+thick),3);
-
+	
 	// the beta factor is not calculated
  	// the calculated form factor <f^2> has units [length^2]
  	// and must be multiplied by number density [l^-3] and the correct unit
  	// conversion to get to absolute scale
-
+	
 	pi43=4.0/3.0*pi;
  	pp=corrad/shlrad;
@@ -676,14 +715,14 @@
  	c1=drho1*volume;
  	c2=drho2*volume;
-
+	
 	arg1 = x*shlrad*pp;
 	arg2 = x*shlrad;
-
+	
 	form=pow(pp,6)*c1*c1*fnt2(arg1,zz);
 	form += c2*c2*fnt2(arg2,zz);
 	form += 2.0*c1*c2*fnt3(arg2,pp,zz);
-
+	
 	//convert the result to [cm^-1]
-
+	
 	//scale the result
 	// - divide by the polydisperse volume, mult by 10^8
@@ -691,8 +730,8 @@
 	form *= 1.0e8;
 	form *= scale;
-
+	
 	//add in the background
 	form += bkg;
-
+	
 	return(form);
 }
@@ -706,5 +745,5 @@
 {
 	double z1,z2,z3,u,ww,term1,term2,term3,ans;
-
+	
 	z1=zz+1.0;
 	z2=zz+2.0;
@@ -716,5 +755,5 @@
 	term3=1.0+cos(z3*ww)/pow((1.0+4.0*u*u),(z3/2.0));
 	ans=(4.50/z1/pow(yy,6))*(z1*(1.0-term1-term2)+yy*yy*z2*term3);
-
+	
 	return(ans);
 }
@@ -726,7 +765,7 @@
 double
 fnt3(double yy, double pp, double zz)
-{
+{     	
 	double z1,z2,z3,yp,yn,up,un,vp,vn,term1,term2,term3,term4,term5,term6,ans;
-
+	
 	z1=zz+1.0;
 	z2=zz+2.0;
@@ -746,5 +785,5 @@
 	ans=4.5/z1/pow(yy,6);
 	ans *=(z1*(term1-term2)+yy*yy*pp*z2*(term3+term4)+z1*(term5-term6));
-
+	
 	return(ans);
 }
@@ -771,5 +810,5 @@
 	double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2;
 	int err;
-
+	
 	pi = 4.0*atan(1.0);
 	x= q;
@@ -782,6 +821,6 @@
 	rhos = dp[6];
 	bgd = dp[7];
-
-
+	
+	
 	phi = phi1 + phi2;
 	aa = r1/r2;
@@ -792,13 +831,13 @@
 	// calculate the PSF's here
 	err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-
+	
 	// /* do form factor calculations  */
-
+	
 	v1 = 4.0*pi/3.0*r1*r1*r1;
 	v2 = 4.0*pi/3.0*r2*r2*r2;
-
+	
 	n1 = phi1/v1;
 	n2 = phi2/v2;
-
+	
 	qr1 = r1*x;
 	qr2 = r2*x;
@@ -821,7 +860,7 @@
 	///* convert I(1/A) to (1/cm)  */
 	inten *= 1.0e8;
-
+	
 	inten += bgd;
-
+	
 	return(inten);
 }
@@ -835,5 +874,5 @@
 	double phi1,phi2,phr,a3;
 	int err;
-
+	
 	pi = 4.0*atan(1.0);
 	x= q;
@@ -854,5 +893,5 @@
 	// calculate the PSF's here
 	err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-
+	
     return(psf11);	//scale, and add in the background
 }
@@ -866,5 +905,5 @@
 	double phi1,phi2,phr,a3;
 	int err;
-
+	
 	pi = 4.0*atan(1.0);
 	x= q;
@@ -885,5 +924,5 @@
 	// calculate the PSF's here
 	err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-
+	
     return(psf12);	//scale, and add in the background
 }
@@ -897,8 +936,8 @@
 	double phi1,phi2,phr,a3;
 	int err;
-
+	
 	pi = 4.0*atan(1.0);
 	x= q;
-
+	
 	r2 = dp[0];
 	r1 = dp[1];
@@ -917,5 +956,5 @@
 	// calculate the PSF's here
 	err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-
+	
     return(psf22);	//scale, and add in the background
 }
@@ -925,12 +964,12 @@
 {
 	//	variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
-
+	
 	//   calculate constant terms
 	double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
 	double a1,a2i,a2,b1,b2,b12,gm1,gm12;
-	double err=0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
+	double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
 	double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
 	double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
-
+	
 	s2 = 2.0*r2;
 	s1 = aa*s2;
@@ -953,9 +992,9 @@
 	gm1=(v1*a1+a3*v2*a2)*.5;
 	gm12=2.*gm1*(1.-aa)/aa;
-	//c
+	//c  
 	//c   calculate the direct correlation functions and print results
 	//c
 	//	do 20 j=1,npts
-
+	
 	yy=qval*s2;
 	//c   calculate direct correlation functions
@@ -968,5 +1007,5 @@
 	t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
 	f11=24.*v1*(t1+t2+t3)/ay3;
-
+	
 	//c ------c22
 	y2=yy*yy;
@@ -976,5 +1015,5 @@
 	tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
 	f22=24.*v2*(tt1+tt2+tt3)/y3;
-
+	
 	//c   -----c12
 	yl=.5*yy*(1.-aa);
@@ -996,12 +1035,12 @@
 	ttt4=a1*(t41+t42)/y1;
 	f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
-
+	
 	c11=f11;
 	c22=f22;
 	c12=f12;
 	*s11=1./(1.+c11-(c12)*c12/(1.+c22));
-	*s22=1./(1.+c22-(c12)*c12/(1.+c11));
-	*s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));
-
+	*s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
+	*s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));   
+	
 	return(err);
 }
@@ -1016,5 +1055,5 @@
  // bragg peaks arise naturally from the periodicity of the sample
  // resolution smeared version gives he most appropriate view of the model
-
+ 
 	Warning:
  The call to WaveData() below returns a pointer to the middle
@@ -1031,7 +1070,7 @@
 	double fval,voli,ri,sldi;
 	double pi;
-
-	pi = 4.0*atan(1.0);
-
+	
+	pi = 4.0*atan(1.0);
+	
 	x= q;
 	scale = dp[0];
@@ -1043,27 +1082,27 @@
 	num = dp[6];
 	bkg = dp[7];
-
+	
 	//calculate with a loop, two shells at a time
-
+	
 	ii=0;
-	fval=0;
-
+	fval=0.0;
+	
 	do {
 		ri = rcore + (double)ii*ts + (double)ii*tw;
-		voli = 4*pi/3*ri*ri*ri;
+		voli = 4.0*pi/3.0*ri*ri*ri;
 		sldi = rhocore-rhoshel;
 		fval += voli*sldi*F_func(ri*x);
 		ri += ts;
-		voli = 4*pi/3*ri*ri*ri;
+		voli = 4.0*pi/3.0*ri*ri*ri;
 		sldi = rhoshel-rhocore;
 		fval += voli*sldi*F_func(ri*x);
 		ii+=1;		//do 2 layers at a time
 	} while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
-
+	
 	fval *= fval;		//square it
 	fval /= voli;		//normalize by the overall volume
-	fval *= scale*1e8;
+	fval *= scale*1.0e8;
 	fval += bkg;
-
+	
 	return(fval);
 }
@@ -1095,8 +1134,8 @@
 	double fval,ri,pi;
 	double avg,pd,zz,lo,hi,r1,r2,d1,d2,distr;
-
-	pi = 4.0*atan(1.0);
+	
+	pi = 4.0*atan(1.0);	
 	x= q;
-
+	
 	scale = dp[0];
 	avg = dp[1];		// average (total) outer radius
@@ -1108,28 +1147,28 @@
 	rhoshel = dp[7];
 	bkg = dp[8];
-
+	
 	zz = (1.0/pd)*(1.0/pd)-1.0;
-
+	
 	//max radius set to be 5 std deviations past mean
 	hi = avg + pd*avg*5.0;
 	lo = avg - pd*avg*5.0;
-
+	
 	maxPairs = trunc( (hi-rcore+tw)/(ts+tw) );
 	minPairs = trunc( (lo-rcore+tw)/(ts+tw) );
 	minPairs = (minPairs < 1) ? 1 : minPairs;	// need a minimum of one
-
+	
 	ii=minPairs;
-	fval=0;
-	d1 = 0;
-	d2 = 0;
-	r1 = 0;
-	r2 = 0;
-	distr = 0;
-	first = 1;
+	fval=0.0;
+	d1 = 0.0;
+	d2 = 0.0;
+	r1 = 0.0;
+	r2 = 0.0;
+	distr = 0.0;
+	first = 1.0;
 	do {
 		//make the current values old
 		r1 = r2;
 		d1 = d2;
-
+		
 		ri = (double)ii*(ts+tw) - tw + rcore;
 		fval += SchulzPoint(ri,avg,zz) * MultiShellGuts(x, rcore, ts, tw, rhocore, rhoshel, ii) * (4*pi/3*ri*ri*ri);
@@ -1143,10 +1182,10 @@
 		first = 0;
 	} while(ii<=maxPairs);
-
-	fval /= 4*pi/3*avg*avg*avg;		//normalize by the overall volume
+	
+	fval /= 4.0*pi/3.0*avg*avg*avg;		//normalize by the overall volume
 	fval /= distr;
 	fval *= scale;
 	fval += bkg;
-
+	
 	return(fval);
 }
@@ -1154,28 +1193,28 @@
 double
 MultiShellGuts(double x,double rcore,double ts,double tw,double rhocore,double rhoshel,int num) {
-
+	
     double ri,sldi,fval,voli,pi;
     int ii;
-
+    
 	pi = 4.0*atan(1.0);
     ii=0;
-    fval=0;
-
+    fval=0.0;
+    
     do {
         ri = rcore + (double)ii*ts + (double)ii*tw;
-        voli = 4*pi/3*ri*ri*ri;
+        voli = 4.0*pi/3.0*ri*ri*ri;
         sldi = rhocore-rhoshel;
         fval += voli*sldi*F_func(ri*x);
         ri += ts;
-        voli = 4*pi/3*ri*ri*ri;
+        voli = 4.0*pi/3.0*ri*ri*ri;
         sldi = rhoshel-rhocore;
         fval += voli*sldi*F_func(ri*x);
         ii+=1;		//do 2 layers at a time
     } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
-
+    
     fval *= fval;
     fval /= voli;
-    fval *= 1e8;
-
+    fval *= 1.0e8;
+    
     return(fval);	// this result still needs to be multiplied by scale and have background added
 }
@@ -1183,8 +1222,8 @@
 static double
 SchulzPoint(double x, double avg, double zz) {
-
+	
     double dr;
-
-    dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1));
+    
+    dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0));
     return (exp(dr));
 }
@@ -1192,5 +1231,5 @@
 static double
 gammln(double xx) {
-
+	
     double x,y,tmp,ser;
     static double cof[6]={76.18009172947146,-86.50532032941677,
@@ -1198,5 +1237,5 @@
 		0.1208650973866179e-2,-0.5395239384953e-5};
     int j;
-
+	
     y=x=xx;
     tmp=x+5.5;
@@ -1210,10 +1249,1036 @@
 F_func(double qr) {
 	double sc;
-	if (qr == 0){
+	if (qr == 0.0){
 		sc = 1.0;
 	}else{
-		sc=(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
+		sc=(3.0*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
 	}
 	return sc;
 }
 
+double
+OneShell(double dp[], double q)
+{
+	// variables are:
+	//[0] scale factor
+	//[1] radius of core []
+	//[2] SLD of the core	[-2]
+	//[3] thickness of the shell	[]
+	//[4] SLD of the shell
+	//[5] SLD of the solvent
+	//[6] background	[cm-1]
+	
+	double x,pi;
+	double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;		//my local names
+	double bes,f,vol,qr,contr,f2;
+	
+	pi = 4.0*atan(1.0);
+	x=q;
+	
+	scale = dp[0];
+	rcore = dp[1];
+	rhocore = dp[2];
+	thick = dp[3];
+	rhoshel = dp[4];
+	rhosolv = dp[5];
+	bkg = dp[6];
+	
+	// core first, then add in shell
+	qr=x*rcore;
+	contr = rhocore-rhoshel;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*rcore*rcore*rcore;
+	f = vol*bes*contr;
+	//now the shell
+	qr=x*(rcore+thick);
+	contr = rhoshel-rhosolv;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*pow((rcore+thick),3);
+	f += vol*bes*contr;
+	
+	// normalize to particle volume and rescale from [-1] to [cm-1]
+	f2 = f*f/vol*1.0e8;
+	
+	//scale if desired
+	f2 *= scale;
+	// then add in the background
+	f2 += bkg;
+	
+	return(f2);
+}
+
+double
+TwoShell(double dp[], double q)
+{
+	// variables are:
+	//[0] scale factor
+	//[1] radius of core []
+	//[2] SLD of the core	[-2]
+	//[3] thickness of shell 1 []
+	//[4] SLD of shell 1
+	//[5] thickness of shell 2 []
+	//[6] SLD of shell 2
+	//[7] SLD of the solvent
+	//[8] background	[cm-1]
+	
+	double x,pi;
+	double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;		//my local names
+	double bes,f,vol,qr,contr,f2;
+	double rhoshel2,thick2;
+	
+	pi = 4.0*atan(1.0);
+	x=q;
+	
+	scale = dp[0];
+	rcore = dp[1];
+	rhocore = dp[2];
+	thick1 = dp[3];
+	rhoshel1 = dp[4];
+	thick2 = dp[5];
+	rhoshel2 = dp[6];	
+	rhosolv = dp[7];
+	bkg = dp[8];
+		// core first, then add in shells
+	qr=x*rcore;
+	contr = rhocore-rhoshel1;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*rcore*rcore*rcore;
+	f = vol*bes*contr;
+	//now the shell (1)
+	qr=x*(rcore+thick1);
+	contr = rhoshel1-rhoshel2;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
+	f += vol*bes*contr;
+	//now the shell (2)
+	qr=x*(rcore+thick1+thick2);
+	contr = rhoshel2-rhosolv;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
+	f += vol*bes*contr;
+	
+		
+	// normalize to particle volume and rescale from [-1] to [cm-1]
+	f2 = f*f/vol*1.0e8;
+	
+	//scale if desired
+	f2 *= scale;
+	// then add in the background
+	f2 += bkg;
+	
+	return(f2);
+}
+
+double
+ThreeShell(double dp[], double q)
+{
+	// variables are:
+	//[0] scale factor
+	//[1] radius of core []
+	//[2] SLD of the core	[-2]
+	//[3] thickness of shell 1 []
+	//[4] SLD of shell 1
+	//[5] thickness of shell 2 []
+	//[6] SLD of shell 2
+	//[7] thickness of shell 3
+	//[8] SLD of shell 3
+	//[9] SLD of solvent
+	//[10] background	[cm-1]
+	
+	double x,pi;
+	double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;		//my local names
+	double bes,f,vol,qr,contr,f2;
+	double rhoshel2,thick2,rhoshel3,thick3;
+	
+	pi = 4.0*atan(1.0);
+	x=q;
+	
+	scale = dp[0];
+	rcore = dp[1];
+	rhocore = dp[2];
+	thick1 = dp[3];
+	rhoshel1 = dp[4];
+	thick2 = dp[5];
+	rhoshel2 = dp[6];	
+	thick3 = dp[7];
+	rhoshel3 = dp[8];	
+	rhosolv = dp[9];
+	bkg = dp[10];
+	
+		// core first, then add in shells
+	qr=x*rcore;
+	contr = rhocore-rhoshel1;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*rcore*rcore*rcore;
+	f = vol*bes*contr;
+	//now the shell (1)
+	qr=x*(rcore+thick1);
+	contr = rhoshel1-rhoshel2;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
+	f += vol*bes*contr;
+	//now the shell (2)
+	qr=x*(rcore+thick1+thick2);
+	contr = rhoshel2-rhoshel3;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
+	f += vol*bes*contr;
+	//now the shell (3)
+	qr=x*(rcore+thick1+thick2+thick3);
+	contr = rhoshel3-rhosolv;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
+	f += vol*bes*contr;
+		
+	// normalize to particle volume and rescale from [-1] to [cm-1]
+	f2 = f*f/vol*1.0e8;
+	
+	//scale if desired
+	f2 *= scale;
+	// then add in the background
+	f2 += bkg;
+	
+	return(f2);
+}
+
+double
+FourShell(double dp[], double q)
+{
+	// variables are:
+	//[0] scale factor
+	//[1] radius of core []
+	//[2] SLD of the core	[-2]
+	//[3] thickness of shell 1 []
+	//[4] SLD of shell 1
+	//[5] thickness of shell 2 []
+	//[6] SLD of shell 2
+	//[7] thickness of shell 3
+	//[8] SLD of shell 3
+	//[9] thickness of shell 3
+	//[10] SLD of shell 3
+	//[11] SLD of solvent
+	//[12] background	[cm-1]
+	
+	double x,pi;
+	double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;		//my local names
+	double bes,f,vol,qr,contr,f2;
+	double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4;
+	
+	pi = 4.0*atan(1.0);
+	x=q;
+	
+	scale = dp[0];
+	rcore = dp[1];
+	rhocore = dp[2];
+	thick1 = dp[3];
+	rhoshel1 = dp[4];
+	thick2 = dp[5];
+	rhoshel2 = dp[6];	
+	thick3 = dp[7];
+	rhoshel3 = dp[8];
+	thick4 = dp[9];
+	rhoshel4 = dp[10];	
+	rhosolv = dp[11];
+	bkg = dp[12];
+	
+		// core first, then add in shells
+	qr=x*rcore;
+	contr = rhocore-rhoshel1;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*rcore*rcore*rcore;
+	f = vol*bes*contr;
+	//now the shell (1)
+	qr=x*(rcore+thick1);
+	contr = rhoshel1-rhoshel2;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
+	f += vol*bes*contr;
+	//now the shell (2)
+	qr=x*(rcore+thick1+thick2);
+	contr = rhoshel2-rhoshel3;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
+	f += vol*bes*contr;
+	//now the shell (3)
+	qr=x*(rcore+thick1+thick2+thick3);
+	contr = rhoshel3-rhoshel4;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
+	f += vol*bes*contr;
+	//now the shell (4)
+	qr=x*(rcore+thick1+thick2+thick3+thick4);
+	contr = rhoshel4-rhosolv;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	}
+	vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4);
+	f += vol*bes*contr;
+	
+		
+	// normalize to particle volume and rescale from [-1] to [cm-1]
+	f2 = f*f/vol*1.0e8;
+	
+	//scale if desired
+	f2 *= scale;
+	// then add in the background
+	f2 += bkg;
+	
+	return(f2);
+}
+
+double
+PolyOneShell(double dp[], double x)
+{
+	double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg,pd,zz;		//my local names
+	double va,vb,summ,yyy,zi;
+	double answer,zp1,zp2,zp3,vpoly,range,temp_1sf[7],pi;
+	int nord=76,ii;
+	
+	pi = 4.0*atan(1.0);
+	
+	scale = dp[0];
+	rcore = dp[1];
+	pd = dp[2];
+	rhocore = dp[3];
+	thick = dp[4];
+	rhoshel = dp[5];
+	rhosolv = dp[6];
+	bkg = dp[7];
+		
+	zz = (1.0/pd)*(1.0/pd)-1.0;		//polydispersity of the core only
+	
+	range = 8.0;			//std deviations for the integration
+	va = rcore*(1.0-range*pd);
+	if (va<0.0) {
+		va=0.0;		//otherwise numerical error when pd >= 0.3, making a<0
+	}
+	if (pd>0.3) {
+		range = range + (pd-0.3)*18.0;		//stretch upper range to account for skewed tail
+	}
+	vb = rcore*(1.0+range*pd);		// is this far enough past avg radius?
+	
+	//temp set scale=1 and bkg=0 for quadrature calc
+	temp_1sf[0] = 1.0;
+	temp_1sf[1] = dp[1];	//the core radius will be changed in the loop
+	temp_1sf[2] = dp[3];
+	temp_1sf[3] = dp[4];
+	temp_1sf[4] = dp[5];
+	temp_1sf[5] = dp[6];
+	temp_1sf[6] = 0.0;
+	
+	summ = 0.0;		// initialize integral
+	for(ii=0;ii<nord;ii+=1) {
+		// calculate Gauss points on integration interval (r-value for evaluation)
+		zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
+		temp_1sf[1] = zi;
+		yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * OneShell(temp_1sf,x);
+		//un-normalize by volume
+		yyy *= 4.0*pi/3.0*pow((zi+thick),3);
+		summ += yyy;		//add to the running total of the quadrature
+   	}
+	// calculate value of integral to return
+	answer = (vb-va)/2.0*summ;
+   	
+	//re-normalize by the average volume
+	zp1 = zz + 1.0;
+	zp2 = zz + 2.0;
+	zp3 = zz + 3.0;
+	vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick),3);
+	answer /= vpoly;
+//scale
+	answer *= scale;
+// add in the background
+	answer += bkg;
+		
+    return(answer);
+}
+
+double
+PolyTwoShell(double dp[], double x)
+{
+	double scale,rcore,rhocore,rhosolv,bkg,pd,zz;		//my local names
+	double va,vb,summ,yyy,zi;
+	double answer,zp1,zp2,zp3,vpoly,range,temp_2sf[9],pi;
+	int nord=76,ii;
+	double thick1,thick2;
+	double rhoshel1,rhoshel2;
+	
+	scale = dp[0];
+	rcore = dp[1];
+	pd = dp[2];
+	rhocore = dp[3];
+	thick1 = dp[4];
+	rhoshel1 = dp[5];
+	thick2 = dp[6];
+	rhoshel2 = dp[7];
+	rhosolv = dp[8];
+	bkg = dp[9];
+	
+	pi = 4.0*atan(1.0);
+		
+	zz = (1.0/pd)*(1.0/pd)-1.0;		//polydispersity of the core only
+	
+	range = 8.0;			//std deviations for the integration
+	va = rcore*(1.0-range*pd);
+	if (va<0.0) {
+		va=0.0;		//otherwise numerical error when pd >= 0.3, making a<0
+	}
+	if (pd>0.3) {
+		range = range + (pd-0.3)*18.0;		//stretch upper range to account for skewed tail
+	}
+	vb = rcore*(1.0+range*pd);		// is this far enough past avg radius?
+	
+	//temp set scale=1 and bkg=0 for quadrature calc
+	temp_2sf[0] = 1.0;
+	temp_2sf[1] = dp[1];		//the core radius will be changed in the loop
+	temp_2sf[2] = dp[3];
+	temp_2sf[3] = dp[4];
+	temp_2sf[4] = dp[5];
+	temp_2sf[5] = dp[6];
+	temp_2sf[6] = dp[7];
+	temp_2sf[7] = dp[8];
+	temp_2sf[8] = 0.0;
+	
+	summ = 0.0;		// initialize integral
+	for(ii=0;ii<nord;ii+=1) {
+		// calculate Gauss points on integration interval (r-value for evaluation)
+		zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
+		temp_2sf[1] = zi;
+		yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * TwoShell(temp_2sf,x);
+		//un-normalize by volume
+		yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2),3);
+		summ += yyy;		//add to the running total of the quadrature
+   	}
+	// calculate value of integral to return
+	answer = (vb-va)/2.0*summ;
+   	
+	//re-normalize by the average volume
+	zp1 = zz + 1.0;
+	zp2 = zz + 2.0;
+	zp3 = zz + 3.0;
+	vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2),3);
+	answer /= vpoly;
+//scale
+	answer *= scale;
+// add in the background
+	answer += bkg;
+		
+    return(answer);
+}
+
+double
+PolyThreeShell(double dp[], double x)
+{
+	double scale,rcore,rhocore,rhosolv,bkg,pd,zz;		//my local names
+	double va,vb,summ,yyy,zi;
+	double answer,zp1,zp2,zp3,vpoly,range,temp_3sf[11],pi;
+	int nord=76,ii;
+	double thick1,thick2,thick3;
+	double rhoshel1,rhoshel2,rhoshel3;
+	
+	scale = dp[0];
+	rcore = dp[1];
+	pd = dp[2];
+	rhocore = dp[3];
+	thick1 = dp[4];
+	rhoshel1 = dp[5];
+	thick2 = dp[6];
+	rhoshel2 = dp[7];
+	thick3 = dp[8];
+	rhoshel3 = dp[9];
+	rhosolv = dp[10];
+	bkg = dp[11];
+	
+	pi = 4.0*atan(1.0);
+		
+	zz = (1.0/pd)*(1.0/pd)-1.0;		//polydispersity of the core only
+	
+	range = 8.0;			//std deviations for the integration
+	va = rcore*(1.0-range*pd);
+	if (va<0) {
+		va=0;		//otherwise numerical error when pd >= 0.3, making a<0
+	}
+	if (pd>0.3) {
+		range = range + (pd-0.3)*18.0;		//stretch upper range to account for skewed tail
+	}
+	vb = rcore*(1.0+range*pd);		// is this far enough past avg radius?
+	
+	//temp set scale=1 and bkg=0 for quadrature calc
+	temp_3sf[0] = 1.0;
+	temp_3sf[1] = dp[1];		//the core radius will be changed in the loop
+	temp_3sf[2] = dp[3];
+	temp_3sf[3] = dp[4];
+	temp_3sf[4] = dp[5];
+	temp_3sf[5] = dp[6];
+	temp_3sf[6] = dp[7];
+	temp_3sf[7] = dp[8];
+	temp_3sf[8] = dp[9];
+	temp_3sf[9] = dp[10];
+	temp_3sf[10] = 0.0;
+	
+	summ = 0.0;		// initialize integral
+	for(ii=0;ii<nord;ii+=1) {
+		// calculate Gauss points on integration interval (r-value for evaluation)
+		zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
+		temp_3sf[1] = zi;
+		yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * ThreeShell(temp_3sf,x);
+		//un-normalize by volume
+		yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3),3);
+		summ += yyy;		//add to the running total of the quadrature
+   	}
+	// calculate value of integral to return
+	answer = (vb-va)/2.0*summ;
+   	
+	//re-normalize by the average volume
+	zp1 = zz + 1.0;
+	zp2 = zz + 2.0;
+	zp3 = zz + 3.0;
+	vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3),3);
+	answer /= vpoly;
+//scale
+	answer *= scale;
+// add in the background
+	answer += bkg;
+		
+    return(answer);
+}
+
+double
+PolyFourShell(double dp[], double x)
+{
+	double scale,rcore,rhocore,rhosolv,bkg,pd,zz;		//my local names
+	double va,vb,summ,yyy,zi;
+	double answer,zp1,zp2,zp3,vpoly,range,temp_4sf[13],pi;
+	int nord=76,ii;
+	double thick1,thick2,thick3,thick4;
+	double rhoshel1,rhoshel2,rhoshel3,rhoshel4;
+	
+	scale = dp[0];
+	rcore = dp[1];
+	pd = dp[2];
+	rhocore = dp[3];
+	thick1 = dp[4];
+	rhoshel1 = dp[5];
+	thick2 = dp[6];
+	rhoshel2 = dp[7];
+	thick3 = dp[8];
+	rhoshel3 = dp[9];
+	thick4 = dp[10];
+	rhoshel4 = dp[11];
+	rhosolv = dp[12];
+	bkg = dp[13];
+	
+	pi = 4.0*atan(1.0);
+		
+	zz = (1.0/pd)*(1.0/pd)-1.0;		//polydispersity of the core only
+	
+	range = 8.0;			//std deviations for the integration
+	va = rcore*(1.0-range*pd);
+	if (va<0) {
+		va=0;		//otherwise numerical error when pd >= 0.3, making a<0
+	}
+	if (pd>0.3) {
+		range = range + (pd-0.3)*18.0;		//stretch upper range to account for skewed tail
+	}
+	vb = rcore*(1.0+range*pd);		// is this far enough past avg radius?
+	
+	//temp set scale=1 and bkg=0 for quadrature calc
+	temp_4sf[0] = 1.0;
+	temp_4sf[1] = dp[1];		//the core radius will be changed in the loop
+	temp_4sf[2] = dp[3];
+	temp_4sf[3] = dp[4];
+	temp_4sf[4] = dp[5];
+	temp_4sf[5] = dp[6];
+	temp_4sf[6] = dp[7];
+	temp_4sf[7] = dp[8];
+	temp_4sf[8] = dp[9];
+	temp_4sf[9] = dp[10];
+	temp_4sf[10] = dp[11];
+	temp_4sf[11] = dp[12];
+	temp_4sf[12] = 0.0;
+		
+	summ = 0.0;		// initialize integral
+	for(ii=0;ii<nord;ii+=1) {
+		// calculate Gauss points on integration interval (r-value for evaluation)
+		zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
+		temp_4sf[1] = zi;
+		yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * FourShell(temp_4sf,x);
+		//un-normalize by volume
+		yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3+thick4),3);
+		summ += yyy;		//add to the running total of the quadrature
+   	}
+	// calculate value of integral to return
+	answer = (vb-va)/2.0*summ;
+   	
+	//re-normalize by the average volume
+	zp1 = zz + 1.0;
+	zp2 = zz + 2.0;
+	zp3 = zz + 3.0;
+	vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3+thick4),3);
+	answer /= vpoly;
+//scale
+	answer *= scale;
+// add in the background
+	answer += bkg;
+		
+    return(answer);
+}
+
+
+/*	BCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
+
+Uses 150 pt Gaussian quadrature for both integrals
+
+*/
+double
+BCC_ParaCrystal(double w[], double x)
+{
+	int i,j;
+	double Pi;
+	double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;		//local variables of coefficient wave
+	int nordi=150;			//order of integration
+	int nordj=150;
+	double va,vb;		//upper and lower integration limits
+	double summ,zi,yyy,answer;			//running tally of integration
+	double summj,vaj,vbj,zij;			//for the inner integration
+	
+	Pi = 4.0*atan(1.0);
+	va = 0.0;
+	vb = 2.0*Pi;		//orintational average, outer integral
+	vaj = 0.0;
+	vbj = Pi;		//endpoints of inner integral
+	
+	summ = 0.0;			//initialize intergral
+	
+	scale = w[0];
+	Dnn = w[1];					//Nearest neighbor distance A
+	gg = w[2];					//Paracrystal distortion factor
+	Rad = w[3];					//Sphere radius
+	sld = w[4];
+	sldSolv = w[5];
+	background = w[6]; 
+	
+	contrast = sld - sldSolv;
+	
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	latticeScale = 2.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn/sqrt(3.0/4.0),3);
+	
+	for(i=0;i<nordi;i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		summj=0.0;
+		zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
+		for(j=0;j<nordj;j++) {
+			//20 gauss points for the inner integral
+			zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
+			yyy = Gauss150Wt[j] * BCC_Integrand(w,x,zi,zij);
+			summj += yyy;
+		}
+		//now calculate the value of the inner integral
+		answer = (vbj-vaj)/2.0*summj;
+		
+		//now calculate outer integral
+		yyy = Gauss150Wt[i] * answer;
+		summ += yyy;
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+	
+	answer = (vb-va)/2.0*summ;
+	// Multiply by contrast^2
+	answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
+	// add in the background
+	answer += background;
+	
+	return answer;
+}
+
+// xx is phi (outer)
+// yy is theta (inner)
+double
+BCC_Integrand(double w[], double qq, double xx, double yy) {
+	
+	double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi;
+	
+	Dnn = w[1]; //Nearest neighbor distance A
+	gg = w[2]; //Paracrystal distortion factor
+	aa = Dnn;
+	Da = gg*aa;
+	
+	Pi = 4.0*atan(1.0);
+	temp1 = qq*qq*Da*Da;
+	temp3 = qq*aa;	
+	
+	retVal = BCCeval(yy,xx,temp1,temp3);
+	retVal /=4.0*Pi;
+	
+	return(retVal);
+}
+
+double
+BCCeval(double Theta, double Phi, double temp1, double temp3) {
+
+	double temp6,temp7,temp8,temp9,temp10;
+	double result;
+	
+	temp6 = sin(Theta);
+	temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta);
+	temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta);
+	temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta);
+	temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
+	result = pow(1.0-(temp10*temp10),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
+	
+	return (result);
+}
+
+double
+SphereForm_Paracrystal(double radius, double delrho, double x) {					
+	
+	double bes,f,vol,f2,pi;
+	pi = 4.0*atan(1.0);
+	//
+	//handle q==0 separately
+	if(x==0) {
+		f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*1.0e8;
+		return(f);
+	}
+	
+	bes = 3.0*(sin(x*radius)-x*radius*cos(x*radius))/(x*x*x)/(radius*radius*radius);
+	vol = 4.0*pi/3.0*radius*radius*radius;
+	f = vol*bes*delrho	;	// [=] 
+	// normalize to single particle volume, convert to 1/cm
+	f2 = f * f / vol * 1.0e8;		// [=] 1/cm
+	
+	return (f2);
+}
+
+/*	FCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
+
+Uses 150 pt Gaussian quadrature for both integrals
+
+*/
+double
+FCC_ParaCrystal(double w[], double x)
+{
+	int i,j;
+	double Pi;
+	double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;		//local variables of coefficient wave
+	int nordi=150;			//order of integration
+	int nordj=150;
+	double va,vb;		//upper and lower integration limits
+	double summ,zi,yyy,answer;			//running tally of integration
+	double summj,vaj,vbj,zij;			//for the inner integration
+	
+	Pi = 4.0*atan(1.0);
+	va = 0.0;
+	vb = 2.0*Pi;		//orintational average, outer integral
+	vaj = 0.0;
+	vbj = Pi;		//endpoints of inner integral
+	
+	summ = 0.0;			//initialize intergral
+	
+	scale = w[0];
+	Dnn = w[1];					//Nearest neighbor distance A
+	gg = w[2];					//Paracrystal distortion factor
+	Rad = w[3];					//Sphere radius
+	sld = w[4];
+	sldSolv = w[5];
+	background = w[6]; 
+	
+	contrast = sld - sldSolv;
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	latticeScale = 4.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn*sqrt(2.0),3);
+	
+	for(i=0;i<nordi;i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		summj=0.0;
+		zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
+		for(j=0;j<nordj;j++) {
+			//20 gauss points for the inner integral
+			zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
+			yyy = Gauss150Wt[j] * FCC_Integrand(w,x,zi,zij);
+			summj += yyy;
+		}
+		//now calculate the value of the inner integral
+		answer = (vbj-vaj)/2.0*summj;
+		
+		//now calculate outer integral
+		yyy = Gauss150Wt[i] * answer;
+		summ += yyy;
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+	
+	answer = (vb-va)/2.0*summ;
+	// Multiply by contrast^2
+	answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
+	// add in the background
+	answer += background;
+	
+	return answer;
+}
+
+
+// xx is phi (outer)
+// yy is theta (inner)
+double
+FCC_Integrand(double w[], double qq, double xx, double yy) {
+	
+	double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi;
+	
+	Pi = 4.0*atan(1.0);
+	Dnn = w[1]; //Nearest neighbor distance A
+	gg = w[2]; //Paracrystal distortion factor
+	aa = Dnn;
+	Da = gg*aa;
+	
+	temp1 = qq*qq*Da*Da;
+	temp3 = qq*aa;
+	
+	retVal = FCCeval(yy,xx,temp1,temp3);
+	retVal /=4*Pi;
+	
+	return(retVal);
+}
+
+double
+FCCeval(double Theta, double Phi, double temp1, double temp3) {
+
+	double temp6,temp7,temp8,temp9,temp10;
+	double result;
+	
+	temp6 = sin(Theta);
+	temp7 = sin(Theta)*sin(Phi)+cos(Theta);
+	temp8 = -1.0*sin(Theta)*cos(Phi)+cos(Theta);
+	temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi);
+	temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
+	result = pow((1.0-(temp10*temp10)),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
+	
+	return (result);
+}
+
+
+/*	SC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
+
+Uses 150 pt Gaussian quadrature for both integrals
+
+*/
+double
+SC_ParaCrystal(double w[], double x)
+{
+	int i,j;
+	double Pi;
+	double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;		//local variables of coefficient wave
+	int nordi=150;			//order of integration
+	int nordj=150;
+	double va,vb;		//upper and lower integration limits
+	double summ,zi,yyy,answer;			//running tally of integration
+	double summj,vaj,vbj,zij;			//for the inner integration
+	
+	Pi = 4.0*atan(1.0);
+	va = 0.0;
+	vb = Pi/2.0;		//orintational average, outer integral
+	vaj = 0.0;
+	vbj = Pi/2.0;		//endpoints of inner integral
+	
+	summ = 0.0;			//initialize intergral
+	
+	scale = w[0];
+	Dnn = w[1];					//Nearest neighbor distance A
+	gg = w[2];					//Paracrystal distortion factor
+	Rad = w[3];					//Sphere radius
+	sld = w[4];
+	sldSolv = w[5];
+	background = w[6]; 
+	
+	contrast = sld - sldSolv;
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	latticeScale = (4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn,3);
+	
+	for(i=0;i<nordi;i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		summj=0.0;
+		zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
+		for(j=0;j<nordj;j++) {
+			//20 gauss points for the inner integral
+			zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
+			yyy = Gauss150Wt[j] * SC_Integrand(w,x,zi,zij);
+			summj += yyy;
+		}
+		//now calculate the value of the inner integral
+		answer = (vbj-vaj)/2.0*summj;
+		
+		//now calculate outer integral
+		yyy = Gauss150Wt[i] * answer;
+		summ += yyy;
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+	
+	answer = (vb-va)/2.0*summ;
+	// Multiply by contrast^2
+	answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
+	// add in the background
+	answer += background;
+	
+	return answer;
+}
+
+// xx is phi (outer)
+// yy is theta (inner)
+double
+SC_Integrand(double w[], double qq, double xx, double yy) {
+	
+	double retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg,Pi;
+	
+	Pi = 4.0*atan(1.0);
+	Dnn = w[1]; //Nearest neighbor distance A
+	gg = w[2]; //Paracrystal distortion factor
+	aa = Dnn;
+	Da = gg*aa;
+	
+	temp1 = qq*qq*Da*Da;
+	temp2 = pow( 1.0-exp(-1.0*temp1) ,3);
+	temp3 = qq*aa;
+	temp4 = 2.0*exp(-0.5*temp1);
+	temp5 = exp(-1.0*temp1);
+	
+	
+	retVal = temp2*SCeval(yy,xx,temp3,temp4,temp5);
+	retVal *= 2.0/Pi;
+	
+	return(retVal);
+}
+
+double
+SCeval(double Theta, double Phi, double temp3, double temp4, double temp5) { //Function to calculate integrand values for simple cubic structure
+
+	double temp6,temp7,temp8,temp9; //Theta and phi dependent parts of the equation
+	double result;
+	
+	temp6 = sin(Theta);
+	temp7 = -1.0*temp3*sin(Theta)*cos(Phi);
+	temp8 = temp3*sin(Theta)*sin(Phi);
+	temp9 = temp3*cos(Theta);
+	result = temp6/((1.0-temp4*cos((temp7))+temp5)*(1.0-temp4*cos((temp8))+temp5)*(1.0-temp4*cos((temp9))+temp5));
+	
+	return (result);
+}
+
+// scattering from a uniform sphere with a Gaussian size distribution
+//
+double
+FuzzySpheres(double dp[], double q)
+{
+	double pi,x;
+	double scale,rad,pd,sig,rho,rhos,bkg,delrho,sig_surf,f2,bes,vol,f;		//my local names
+	double va,vb,zi,yy,summ,inten;
+	int nord=20,ii;
+	
+	pi = 4.0*atan(1.0);
+	x= q;
+	
+	scale=dp[0];
+	rad=dp[1];
+	pd=dp[2];
+	sig=pd*rad;
+	sig_surf = dp[3];
+	rho=dp[4];
+	rhos=dp[5];
+	delrho=rho-rhos;
+	bkg=dp[6];
+	
+			
+	va = -4.0*sig + rad;
+	if (va<0) {
+		va=0;		//to avoid numerical error when  va<0 (-ve q-value)
+	}
+	vb = 4.0*sig +rad;
+	
+	summ = 0.0;		// initialize integral
+	for(ii=0;ii<nord;ii+=1) {
+		// calculate Gauss points on integration interval (r-value for evaluation)
+		zi = ( Gauss20Z[ii]*(vb-va) + vb + va )/2.0;
+		// calculate sphere scattering
+		//
+		//handle q==0 separately
+		if (x==0.0) {
+			f2 = 4.0/3.0*pi*zi*zi*zi*delrho*delrho*1.0e8;
+			f2 *= exp(-0.5*sig_surf*sig_surf*x*x);
+			f2 *= exp(-0.5*sig_surf*sig_surf*x*x);
+		} else {
+			bes = 3.0*(sin(x*zi)-x*zi*cos(x*zi))/(x*x*x)/(zi*zi*zi);
+			vol = 4.0*pi/3.0*zi*zi*zi;
+			f = vol*bes*delrho;		// [=] A
+			f *= exp(-0.5*sig_surf*sig_surf*x*x);
+			// normalize to single particle volume, convert to 1/cm
+			f2 = f * f / vol * 1.0e8;		// [=] 1/cm
+		}
+	
+		yy = Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi) * f2;
+		yy *= 4.0*pi/3.0*zi*zi*zi;		//un-normalize by current sphere volume
+		
+		summ += yy;		//add to the running total of the quadrature
+		
+		
+   	}
+	// calculate value of integral to return
+	inten = (vb-va)/2.0*summ;
+	
+	//re-normalize by polydisperse sphere volume
+	inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
+	
+	inten *= scale;
+	inten += bkg;
+	
+    return(inten);	//scale, and add in the background
+}
+
+
