Index: sansmodels/src/libigor/libCylinder.c
===================================================================
--- sansmodels/src/libigor/libCylinder.c	(revision 96b593847bbd8eb6072269e750a93261a4f743b3)
+++ sansmodels/src/libigor/libCylinder.c	(revision 975ec8e8211549489759b5a5ffe3b60607de065b)
@@ -80,5 +80,5 @@
 	double va,vb;		//upper and lower integration limits
 	double summ,zi,yyy,answer,vell;			//running tally of integration
-	double summj,vaj,vbj,zij,arg;			//for the inner integration
+	double summj,vaj,vbj,zij,arg, si;			//for the inner integration
 
 	Pi = 4.0*atan(1.0);
@@ -114,5 +114,10 @@
 		//now calculate outer integral
 		arg = q*length*zi/2;
-		yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg;
+		if (arg == 0){
+			si = 1;
+		}else{
+			si = sin(arg) * sin(arg) / arg / arg;
+		}
+		yyy = Gauss76Wt[i] * answer * si;
 		summ += yyy;
 	}
@@ -155,5 +160,5 @@
 	double va,vb;		//upper and lower integration limits
 	double summ,zi,yyy,answer,vell;			//running tally of integration
-	double summj,vaj,vbj,zij,arg;			//for the inner integration
+	double summj,vaj,vbj,zij,arg,si;			//for the inner integration
 
 	Pi = 4.0*atan(1.0);
@@ -189,5 +194,10 @@
 		//now calculate outer integral
 		arg = q*length*zi/2;
-		yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg;
+		if (arg == 0){
+			si = 1;
+		}else{
+			si = sin(arg) * sin(arg) / arg / arg;
+		}
+		yyy = Gauss76Wt[i] * answer * si;
 		summ += yyy;
 	}
@@ -893,5 +903,32 @@
 	return(inten+bkg);
 }
-
+/*	LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included
+						-NO polydispersion included
+*/
+double lamellar_kernel(double dp[], double q){
+	double scale,del,sld_bi,sld_sol,contr,bkg;		//local variables of coefficient wave
+	double inten, qval,Pq;
+	double Pi;
+
+
+	Pi = 4.0*atan(1.0);
+	scale = dp[0];
+	del = dp[1];
+	sld_bi = dp[2];
+	sld_sol = dp[3];
+	bkg = dp[4];
+	qval = q;
+	contr = sld_bi -sld_sol;
+
+	Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
+
+	inten = 2.0*Pi*scale*Pq/(qval*qval);		//this is now dimensionless...
+
+	inten /= del;			//normalize by the thickness (in A)
+
+	inten *= 1.0e8;		// 1/A to 1/cm
+
+	return(inten+bkg);
+}
 /*	LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
 -------
@@ -1682,5 +1719,5 @@
 {
 	// local variables
-	double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,gfn2,pi43,gfn,Pi;
+	double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi;
 
 	Pi = 4.0*atan(1.0);
@@ -1695,6 +1732,16 @@
 	vc = pi43*aa*bb*bb;
 	vt = pi43*trmaj*trmin*trmin;
-	gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc;
-	gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps;
+   	if (uq == 0){
+   		siq = 1/3;
+   	}else{
+   		siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
+   	}
+   	if (ut == 0){
+   		sit = 1/3;
+   	}else{
+   		sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
+   	}
+	gfnc = 3.0*siq*vc*delpc;
+	gfnt = 3.0*sit*vt*delps;
 	gfn = gfnc+gfnt;
 	gfn2 = gfn*gfn;
@@ -1714,5 +1761,5 @@
 {
 	// local variables
-	double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
+	double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
 
 	Pi = 4.0*atan(1.0);
@@ -1726,6 +1773,16 @@
 	vc = pi43*aa*aa*bb;
    	vt = pi43*trmaj*trmaj*trmin;
-   	gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc;
-  	gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps;
+   	if (uq == 0){
+   		siq = 1/3;
+   	}else{
+   		siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
+   	}
+   	if (ut == 0){
+   		sit = 1/3;
+   	}else{
+   		sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
+   	}
+   	gfnc = 3.0*siq*vc*delpc;
+  	gfnt = 3.0*sit*vt*delps;
   	tgfn = gfnc+gfnt;
   	gfn4 = tgfn*tgfn;
@@ -1744,6 +1801,7 @@
 
     Pq = Sk_WR(q,zi,lb);		//does not have cross section term
-    Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
-
+    if (qr !=0){
+    	Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
+    } //else Pk *=1;
     vcyl=Pi*radius*radius*zi;
     Pq *= vcyl*vcyl;
@@ -1764,5 +1822,7 @@
 
     Pq = Sk_WR(q,Lc,Lb);		//does not have cross section term
-    Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
+    if (qr !=0){
+    	Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
+    }
 
     vcyl=Pi*zi*zi*Lc;
@@ -1811,5 +1871,5 @@
 	// dum is the dummy variable for the integration (theta)
 
-	double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval;
+	double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
 	double Pi;
 
@@ -1825,7 +1885,27 @@
 	sinarg1 = qq*length*cos(dum);
 	sinarg2 = qq*(length+facthick)*cos(dum);
-
-	t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
-	t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
+	if (besarg1 == 0){
+		be1 = 0.5;
+	}else{
+		be1 = NR_BessJ1(besarg1)/besarg1;
+	}
+	if (besarg2 == 0){
+		be2 = 0.5;
+	}else{
+		be2 = NR_BessJ1(besarg2)/besarg2;
+	}
+	if (sinarg1 == 0){
+		si1 = 1;
+	}else{
+		si1 = sin(sinarg1)/sinarg1;
+	}
+	if (besarg2 == 0){
+		si2 = 1;
+	}else{
+		si2 = sin(sinarg2)/sinarg2;
+	}
+
+	t1 = 2.0*vol1*dr1*si1*be1;
+	t2 = 2.0*vol2*dr2*si2*be2;
 
 	retval = ((t1+t2)*(t1+t2))*sin(dum);
@@ -1839,5 +1919,5 @@
 {
 
-    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval;
+    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
     double Pi;
 
@@ -1854,6 +1934,27 @@
     sinarg2 = qq*(length+thick)*cos(dum);
 
-    t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
-    t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
+	if (besarg1 == 0){
+		be1 = 0.5;
+	}else{
+		be1 = NR_BessJ1(besarg1)/besarg1;
+	}
+	if (besarg2 == 0){
+		be2 = 0.5;
+	}else{
+		be2 = NR_BessJ1(besarg2)/besarg2;
+	}
+	if (sinarg1 == 0){
+		si1 = 1;
+	}else{
+		si1 = sin(sinarg1)/sinarg1;
+	}
+	if (besarg2 == 0){
+		si2 = 1;
+	}else{
+		si2 = sin(sinarg2)/sinarg2;
+	}
+
+    t1 = 2.0*vol1*dr1*si1*be1;
+    t2 = 2.0*vol2*dr2*si2*be2;
 
     retval = ((t1+t2)*(t1+t2))*sin(dum);
@@ -1906,5 +2007,5 @@
 
 	//Local variables
-	double totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
+	double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
 	double Pi;
 	int kk;
@@ -1923,6 +2024,27 @@
 	sinarg2 = qq*(length+thick)*cos(dum);
 
-	t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg1)/besarg1);
-	t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg2)/besarg2);
+	if (besarg1 == 0){
+		be1 = 0.5;
+	}else{
+		be1 = NR_BessJ1(besarg1)/besarg1;
+	}
+	if (besarg2 == 0){
+		be2 = 0.5;
+	}else{
+		be2 = NR_BessJ1(besarg2)/besarg2;
+	}
+	if (sinarg1 == 0){
+		si1 = 1;
+	}else{
+		si1 = sin(sinarg1)/sinarg1;
+	}
+	if (besarg2 == 0){
+		si2 = 1;
+	}else{
+		si2 = sin(sinarg2)/sinarg2;
+	}
+
+	t1 = 2*area*(2*length)*dr1*(si1)*(be1);
+	t2 = 2*area*dr2*(totald*si2-2*length*si1)*(be2);
 
 	retval =((t1+t2)*(t1+t2))*sin(dum);
@@ -1937,4 +2059,5 @@
 	// end of loop for S(q)
 	sqq=1.0+2.0*sqq/N;
+
 	retval *= sqq;
 
@@ -2016,6 +2139,9 @@
     nu = nua/a;
     arg = qq*a*sqrt(1+dum*dum*(nu*nu-1));
-
-    retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
+    if (arg == 0){
+    	retval =1/3;
+    }else{
+    	retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
+    }
     retval *= retval;
     retval *= 9;
@@ -2032,10 +2158,22 @@
     arg1 = qq*rshell*sqrt(1-dum*dum);		//1=shell (outer radius)
     arg2 = qq*rcore*sqrt(1-dum*dum);			//2=core (inner radius)
-    lam1 = 2*NR_BessJ1(arg1)/arg1;
-    lam2 = 2*NR_BessJ1(arg2)/arg2;
+    if (arg1 == 0){
+    	lam1 = 1;
+    }else{
+    	lam1 = 2*NR_BessJ1(arg1)/arg1;
+    }
+    if (arg2 == 0){
+    	lam2 = 1;
+    }else{
+    	lam2 = 2*NR_BessJ1(arg2)/arg2;
+    }
+    //Todo: Need to check psi behavior as gamma goes to 1.
     psi = 1/(1-gamma*gamma)*(lam1 -  gamma*gamma*lam2);		//SRK 10/19/00
-
     sinarg = qq*length*dum/2;
-    t2 = sin(sinarg)/sinarg;
+    if (sinarg == 0){
+    	t2 = 1;
+    }else{
+    	t2 = sin(sinarg)/sinarg;
+    }
 
     retval = psi*psi*t2*t2;
@@ -2084,6 +2222,9 @@
     arg += cc*cc*dy*dy;
     arg = q*sqrt(arg);
-    val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
-
+    if (arg == 0){
+    	val = 1; // as arg --> 0, val should go to 1.0
+    }else{
+    	val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
+    }
     return (val);
 
@@ -2099,13 +2240,13 @@
 	// h is the HALF-LENGTH of the cylinder = L/2 (A)
 
-    double besarg,bj,retval,d1,t1,b1,t2,b2;		 //Local variables
+    double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si;		 //Local variables
 
 
     besarg = qq*rr*sin(theta);
-
+    siarg = qq * h * cos(theta);
     bj =NR_BessJ1(besarg);
 
 	//* Computing 2nd power */
-    d1 = sin(qq * h * cos(theta));
+    d1 = sin(siarg);
     t1 = d1 * d1;
 	//* Computing 2nd power */
@@ -2113,10 +2254,20 @@
     t2 = d1 * d1 * 4.0 * sin(theta);
 	//* Computing 2nd power */
-    d1 = qq * h * cos(theta);
+    d1 = siarg;
     b1 = d1 * d1;
 	//* Computing 2nd power */
     d1 = qq * rr * sin(theta);
     b2 = d1 * d1;
-    retval = t1 * t2 / b1 / b2;
+    if (besarg == 0){
+    	be = sin(theta);
+    }else{
+    	be = t2 / b2;
+    }
+    if (siarg == 0){
+    	si = 1;
+    }else{
+    	si = t1 / b1;
+    }
+    retval = be * si;
 
     return (retval);
@@ -2136,6 +2287,9 @@
 
 	arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2);
-
-	retval = 2*NR_BessJ1(arg)/arg;
+	if (arg == 0){
+		retval = 1;
+	}else{
+		retval = 2*NR_BessJ1(arg)/arg;
+	}
 
 	//square it
Index: sansmodels/src/libigor/libCylinder.h
===================================================================
--- sansmodels/src/libigor/libCylinder.h	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/libigor/libCylinder.h	(revision 975ec8e8211549489759b5a5ffe3b60607de065b)
@@ -34,4 +34,5 @@
 double EllipCylKernel(double qq, double ra,double nu, double theta);
 double TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy);
+double lamellar_kernel(double dp[], double q);
 double PPKernel(double aa, double mu, double uu);
 double HolCylKernel(double qq, double rcore, double rshell, double length, double dum);
Index: sansmodels/src/libigor/libSphere.c
===================================================================
--- sansmodels/src/libigor/libSphere.c	(revision ae3ce4eaa593202ebf49122e196897c0d45a25b6)
+++ sansmodels/src/libigor/libSphere.c	(revision 975ec8e8211549489759b5a5ffe3b60607de065b)
@@ -14,6 +14,6 @@
 {
 	double scale,radius,delrho,bkg;		//my local names
-	double bes,f,vol,f2,pi;
-	
+	double bes,f,vol,f2,pi,qr;
+
 	pi = 4.0*atan(1.0);
 	scale = dp[0];
@@ -21,16 +21,17 @@
 	delrho = dp[2];
 	bkg = dp[3];
-	//handle q==0 separately
-	if(q==0){
-		f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*scale*1.0e8 + bkg;
-		return(f);
+
+	qr = q*radius;
+	if(qr == 0){
+		bes = 1.0;
+	}else{
+		bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
 	}
-	
-	bes = 3.0*(sin(q*radius)-q*radius*cos(q*radius))/(q*q*q)/(radius*radius*radius);
+
 	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
 }
@@ -43,8 +44,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];
@@ -57,5 +58,9 @@
 	qr=x*rcore;
 	contr = rhocore-rhoshel;
-	bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	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;
@@ -63,16 +68,20 @@
 	qr=x*(rcore+thick);
 	contr = rhoshel-rhosolv;
-	bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	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 [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);
 }
@@ -88,5 +97,5 @@
 	pi = 4.0*atan(1.0);
 	x= q;
-	
+
 	scale = dp[0];
 	rcore = dp[1];
@@ -99,5 +108,9 @@
 	qr=x*rcore;
 	contr = rhocore-rhoshel;
-	bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	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;
@@ -105,18 +118,22 @@
 	qr=x*(rcore+thick);
 	contr = rhoshel-rhosolv;
-	bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+	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 the particle volume and rescale from [A-1] to [cm-1]
 	//note that for the vesicle model, the volume is ONLY the shell volume
 	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);
 }
@@ -135,5 +152,5 @@
 	double zp3,vpoly;
 	double s2y, arg1, arg2, arg3, drh1, drh2;
-	
+
 	pi = 4.0*atan(1.0);
 	qq= q;
@@ -145,9 +162,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;
@@ -157,6 +174,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;
@@ -179,10 +196,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));
@@ -199,5 +216,5 @@
 	// then add in the background
 	form += bkg;
-	
+
 	return(form);
 }
@@ -218,8 +235,8 @@
 	double ka2,r1,heff;
 	double hh,k;
-	
+
 	pi = 4.0*atan(1.0);
 	k= q;
-	
+
 	rad = dp[0];		// radius (A)
 	z2 = dp[1];		//polydispersity (0<z2<1)
@@ -228,9 +245,9 @@
 	bkg = dp[4];
 	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;
@@ -264,30 +281,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);
 	//
@@ -297,10 +314,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
 }
@@ -315,7 +332,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];
@@ -328,10 +345,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);
 }
@@ -344,7 +361,7 @@
 	double x,pq;
 	double scale,ravg,pd,bkg,rho,rhos;		//my local names
-	
+
 	x= q;
-	
+
 	scale = dp[0];
 	ravg = dp[1];
@@ -356,5 +373,5 @@
 	// add in the background
 	pq += bkg;
-	
+
 	return(pq);
 }
@@ -367,9 +384,9 @@
 	double aa,at1,at2,rt1,rt2,rt3,t1,t2,t3;
 	double v1,v2,v3,g1,pq,pi,delrho,zz;
-	
+
 	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;
@@ -381,5 +398,5 @@
 	//
 	aa = (zz+1)/x/ravg;
-	
+
 	at1 = atan(1.0/aa);
 	at2 = atan(2.0/aa);
@@ -387,5 +404,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);
@@ -399,23 +416,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;
-	
+
 	//
-	// 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;
@@ -423,5 +440,5 @@
 	//scale, convert to cm^-1
 	pq *= scale * 1.0e8;
-    
+
     return(pq);
 }
@@ -435,8 +452,8 @@
 	double scale,rad,pd,cont,bkg;		//my local names
 	double inten,h1,qw,qr,width,sig,averad3;
-	
+
 	pi = 4.0*atan(1.0);
 	x= q;
-	
+
 	scale = dp[0];
 	rad = dp[1];		// radius (A)
@@ -444,12 +461,12 @@
 	cont = dp[3];		// contrast (A^-2)
 	bkg = dp[4];
-	
+
 	// 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;
@@ -463,16 +480,16 @@
 	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;
@@ -483,5 +500,5 @@
 	// then add in the background
 	inten += bkg;
-	
+
 	return(inten);
 }
@@ -497,8 +514,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];
@@ -509,5 +526,5 @@
 	delrho=rho-rhos;
 	bkg=dp[5];
-	
+
 	va = -4.0*sig + rad;
 	if (va<0) {
@@ -515,5 +532,5 @@
 	}
 	vb = 4.0*sig +rad;
-	
+
 	summ = 0.0;		// initialize integral
 	for(ii=0;ii<nord;ii+=1) {
@@ -525,17 +542,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
 }
@@ -550,8 +567,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
@@ -562,5 +579,5 @@
 	delrho=rho-rhos;
 	bkg=dp[5];
-	
+
 	va = -3.5*sig + mu;
 	va = exp(va);
@@ -570,5 +587,5 @@
 	vb = 3.5*sig*(1.0+sig) +mu;
 	vb = exp(vb);
-	
+
 	summ = 0.0;		// initialize integral
 	for(ii=0;ii<nord;ii+=1) {
@@ -580,18 +597,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);
 }
@@ -599,7 +616,7 @@
 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 );
@@ -609,7 +626,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);
@@ -627,8 +644,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];
@@ -639,6 +656,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
@@ -648,10 +665,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;
@@ -659,14 +676,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
@@ -674,8 +691,8 @@
 	form *= 1.0e8;
 	form *= scale;
-	
+
 	//add in the background
 	form += bkg;
-	
+
 	return(form);
 }
@@ -689,5 +706,5 @@
 {
 	double z1,z2,z3,u,ww,term1,term2,term3,ans;
-	
+
 	z1=zz+1.0;
 	z2=zz+2.0;
@@ -699,5 +716,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);
 }
@@ -709,7 +726,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;
@@ -729,5 +746,5 @@
 	ans=4.5/z1/pow(yy,6);
 	ans *=(z1*(term1-term2)+yy*yy*pp*z2*(term3+term4)+z1*(term5-term6));
-	
+
 	return(ans);
 }
@@ -752,7 +769,7 @@
 	double psf11,psf12,psf22;
 	double phi1,phi2,phr,a3;
-	double v1,v2,n1,n2,qr1,qr2,b1,b2;
+	double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2;
 	int err;
-	
+
 	pi = 4.0*atan(1.0);
 	x= q;
@@ -765,6 +782,6 @@
 	rhos = dp[6];
 	bgd = dp[7];
-	
-	
+
+
 	phi = phi1 + phi2;
 	aa = r1/r2;
@@ -775,18 +792,28 @@
 	// 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;
-	
-	b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*(sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
-	b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*(sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
+
+	if (qr1 == 0){
+		sc1 = 1.0/3.0;
+	}else{
+		sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
+	}
+	if (qr2 == 0){
+		sc2 = 1.0/3.0;
+	}else{
+		sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
+	}
+	b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*sc1;
+	b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*sc2;
 	inten = n1*b1*b1*psf11;
 	inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
@@ -794,7 +821,7 @@
 	///* convert I(1/A) to (1/cm)  */
 	inten *= 1.0e8;
-	
+
 	inten += bgd;
-	
+
 	return(inten);
 }
@@ -808,5 +835,5 @@
 	double phi1,phi2,phr,a3;
 	int err;
-	
+
 	pi = 4.0*atan(1.0);
 	x= q;
@@ -827,5 +854,5 @@
 	// calculate the PSF's here
 	err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-	
+
     return(psf11);	//scale, and add in the background
 }
@@ -839,5 +866,5 @@
 	double phi1,phi2,phr,a3;
 	int err;
-	
+
 	pi = 4.0*atan(1.0);
 	x= q;
@@ -858,5 +885,5 @@
 	// calculate the PSF's here
 	err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-	
+
     return(psf12);	//scale, and add in the background
 }
@@ -870,8 +897,8 @@
 	double phi1,phi2,phr,a3;
 	int err;
-	
+
 	pi = 4.0*atan(1.0);
 	x= q;
-	
+
 	r2 = dp[0];
 	r1 = dp[1];
@@ -890,5 +917,5 @@
 	// calculate the PSF's here
 	err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-	
+
     return(psf22);	//scale, and add in the background
 }
@@ -898,5 +925,5 @@
 {
 	//	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;
@@ -905,5 +932,5 @@
 	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;
@@ -926,9 +953,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
@@ -941,5 +968,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;
@@ -949,5 +976,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);
@@ -969,12 +996,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);
 }
@@ -989,5 +1016,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
@@ -1004,7 +1031,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];
@@ -1016,10 +1043,10 @@
 	num = dp[6];
 	bkg = dp[7];
-	
+
 	//calculate with a loop, two shells at a time
-	
+
 	ii=0;
 	fval=0;
-	
+
 	do {
 		ri = rcore + (double)ii*ts + (double)ii*tw;
@@ -1033,10 +1060,10 @@
 		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 += bkg;
-	
+
 	return(fval);
 }
@@ -1068,8 +1095,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
@@ -1081,15 +1108,15 @@
 	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;
@@ -1104,5 +1131,5 @@
 		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);
@@ -1116,10 +1143,10 @@
 		first = 0;
 	} while(ii<=maxPairs);
-	
+
 	fval /= 4*pi/3*avg*avg*avg;		//normalize by the overall volume
 	fval /= distr;
 	fval *= scale;
 	fval += bkg;
-	
+
 	return(fval);
 }
@@ -1127,12 +1154,12 @@
 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;
-    
+
     do {
         ri = rcore + (double)ii*ts + (double)ii*tw;
@@ -1146,9 +1173,9 @@
         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;
-    
+
     return(fval);	// this result still needs to be multiplied by scale and have background added
 }
@@ -1156,7 +1183,7 @@
 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));
     return (exp(dr));
@@ -1165,5 +1192,5 @@
 static double
 gammln(double xx) {
-	
+
     double x,y,tmp,ser;
     static double cof[6]={76.18009172947146,-86.50532032941677,
@@ -1171,5 +1198,5 @@
 		0.1208650973866179e-2,-0.5395239384953e-5};
     int j;
-	
+
     y=x=xx;
     tmp=x+5.5;
@@ -1182,5 +1209,11 @@
 double
 F_func(double qr) {
-	return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
-}
-
+	double sc;
+	if (qr == 0){
+		sc = 1.0;
+	}else{
+		sc=(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
+	}
+	return sc;
+}
+
