Index: sansmodels/src/libigor/libCylinder.c
===================================================================
--- sansmodels/src/libigor/libCylinder.c	(revision 34c30203f24453aa528815f48e84a90b0bbd95be)
+++ sansmodels/src/libigor/libCylinder.c	(revision 8e91f01e64a5994e9308c38bccb548b77ea00bed)
@@ -25,11 +25,11 @@
 	double uplim,lolim;		//upper and lower integration limits
 	double summ,zi,yyy,answer,vcyl;			//running tally of integration
-	
+
 	Pi = 4.0*atan(1.0);
 	lolim = 0;
 	uplim = Pi/2.0;
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	radius = dp[1];
@@ -43,5 +43,5 @@
 		summ += yyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
 	// Multiply by contrast^2
@@ -57,5 +57,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -81,5 +81,5 @@
 	double summ,zi,yyy,answer,vell;			//running tally of integration
 	double summj,vaj,vbj,zij,arg;			//for the inner integration
-	
+
 	Pi = 4.0*atan(1.0);
 	va = 0;
@@ -87,7 +87,7 @@
 	vaj=0;
 	vbj=Pi;		//endpoints of inner integral
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	ra = dp[1];
@@ -111,5 +111,5 @@
 		//divide integral by Pi
 		answer /=Pi;
-		
+
 		//now calculate outer integral
 		arg = q*length*zi/2;
@@ -130,5 +130,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -156,5 +156,5 @@
 	double summ,zi,yyy,answer,vell;			//running tally of integration
 	double summj,vaj,vbj,zij,arg;			//for the inner integration
-	
+
 	Pi = 4.0*atan(1.0);
 	va = 0;
@@ -162,7 +162,7 @@
 	vaj=0;
 	vbj=Pi;		//endpoints of inner integral
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	ra = dp[1];
@@ -186,5 +186,5 @@
 		//divide integral by Pi
 		answer /=Pi;
-		
+
 		//now calculate outer integral
 		arg = q*length*zi/2;
@@ -192,5 +192,5 @@
 		summ += yyy;
 	}
-	
+
 	answer = (vb-va)/2.0*summ;
 	// Multiply by contrast^2
@@ -205,6 +205,6 @@
 	answer *= scale;
 	// add in the background
-	answer += bkg;	
-	
+	answer += bkg;
+
 	return answer;
 }
@@ -231,5 +231,5 @@
 	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;
@@ -237,7 +237,7 @@
 	vaj = 0;
 	vbj = 1;		//endpoints of inner integral
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	aa = dp[1];
@@ -258,10 +258,10 @@
 		//now calculate the value of the inner integral
 		answer = (vbj-vaj)/2.0*summj;
-		
+
 		//now calculate outer integral
 		yyy = Gauss76Wt[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
@@ -275,5 +275,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -301,6 +301,6 @@
 	double summj,vaj,vbj;			//for the inner integration
 	double mu,mudum,arg,sigma,uu,vol;
-	
-	
+
+
 	//	Pi = 4.0*atan(1.0);
 	va = 0;
@@ -308,7 +308,7 @@
 	vaj = 0;
 	vbj = 1;		//endpoints of inner integral
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	aa = dp[1];
@@ -317,5 +317,5 @@
 	delrho = dp[4];
 	bkg = dp[5];
-	
+
 	mu = q*bb;
 	vol = aa*bb*cc;
@@ -323,10 +323,10 @@
 	aa = aa/bb;
 	cc = cc/bb;
-	
+
 	for(i=0;i<nordi;i++) {
 		//setup inner integral over the ellipsoidal cross-section
 		summj=0;
 		sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy
-		
+
 		for(j=0;j<nordj;j++) {
 			//76 gauss points for the inner integral
@@ -338,5 +338,5 @@
 		//now calculate the value of the inner integral
 		answer = (vbj-vaj)/2.0*summj;
-		
+
 		arg = mu*cc*sigma/2;
 		if ( arg == 0 ) {
@@ -345,10 +345,10 @@
 			answer *= sin(arg)*sin(arg)/arg/arg;
 		}
-		
+
 		//now sum up the outer integral
 		yyy = Gauss76Wt[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
@@ -362,5 +362,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -385,11 +385,11 @@
 	double va,vb,zi;		//upper and lower integration limits
 	double summ,answer,pi;			//running tally of integration
-	
+
 	pi = 4.0*atan(1.0);
 	va = 0;
 	vb = 1;		//limits of numerical integral
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];		//make local copies in case memory moves
 	rcore = dp[1];
@@ -398,10 +398,10 @@
 	delrho = dp[4];
 	bkg = dp[5];
-	
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
 		summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi);
 	}
- 	
+
 	answer = (vb-va)/2.0*summ;
 	// Multiply by contrast^2
@@ -415,5 +415,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -438,11 +438,11 @@
 	double va,vb,zi;		//upper and lower integration limits
 	double summ,answer,pi;			//running tally of integration
-	
+
 	pi = 4.0*atan(1.0);
 	va = 0;
 	vb = 1;		//limits of numerical integral
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	nua = dp[1];
@@ -450,10 +450,10 @@
 	delrho = dp[3];
 	bkg = dp[4];
-	
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
 		summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi);
 	}
-	
+
 	answer = (vb-va)/2.0*summ;
 	// Multiply by contrast^2
@@ -467,5 +467,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -485,10 +485,10 @@
 	double summ,zi,yyy,answer,Vpoly;			//running tally of integration
 	double range,zz,Pi;
-	
+
 	Pi = 4.0*atan(1.0);
 	range = 3.4;
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	radius = dp[1];
@@ -497,7 +497,7 @@
 	delrho = dp[4];
 	bkg = dp[5];
-	
+
 	zz = (1.0/pd)*(1.0/pd) - 1.0;
-	
+
 	lolim = radius*(1.0-range*pd);		//set the upper/lower limits to cover the distribution
 	if(lolim<0) {
@@ -508,5 +508,5 @@
 	}
 	uplim = radius*(1.0+range*pd);
-	
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -514,5 +514,5 @@
 		summ += yyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
 	//normalize by average cylinder volume
@@ -526,5 +526,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -543,11 +543,11 @@
 	double summ,zi,yyy,answer,Vpoly;			//running tally of integration
 	double range,zz,Pi;
-	
-	
+
+
 	Pi = 4.0*atan(1.0);
 	range = 3.4;
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	radius = dp[1];
@@ -556,7 +556,7 @@
 	delrho = dp[4];
 	bkg = dp[5];
-	
+
 	zz = (1.0/pd)*(1.0/pd) - 1.0;
-	
+
 	lolim = length*(1.0-range*pd);		//set the upper/lower limits to cover the distribution
 	if(lolim<0) {
@@ -567,5 +567,5 @@
 	}
 	uplim = length*(1.0+range*pd);
-	
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -573,5 +573,5 @@
 		summ += yyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
 	//normalize by average cylinder volume (first moment)
@@ -585,5 +585,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -603,12 +603,12 @@
 	double summ,zi,yyy,answer,Vcyl;			//running tally of integration
 	double Pi;
-	
-	Pi = 4.0*atan(1.0);
-	
+
+	Pi = 4.0*atan(1.0);
+
 	lolim = 0.0;
 	uplim = Pi/2.0;
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];		//make local copies in case memory moves
 	rcore = dp[1];
@@ -619,7 +619,7 @@
 	rhosolv = dp[6];
 	bkg = dp[7];
-	
+
 	halfheight = length/2.0;
-	
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -627,7 +627,7 @@
 		summ += yyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
-	// length is the total core length 
+	// length is the total core length
 	Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick);
 	answer /= Vcyl;
@@ -638,5 +638,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -657,10 +657,10 @@
 	double summ,yyy,answer,Vpoly;			//running tally of integration
 	double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr;
-	
-	Pi = 4.0*atan(1.0);
-	
+
+	Pi = 4.0*atan(1.0);
+
 	summ = 0.0;			//initialize intergral
 	Rsqrsumm = 0.0;
-	
+
 	scale = dp[0];
 	radius = dp[1];
@@ -673,5 +673,5 @@
 	rhosolv = dp[8];
 	bkg = dp[9];
-	
+
 	lolim = exp(log(radius)-(4.*sigma));
 	if (lolim<0) {
@@ -679,5 +679,5 @@
 	}
 	uplim = exp(log(radius)+(4.*sigma));
-	
+
 	for(i=0;i<nord;i++) {
 		rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -688,5 +688,5 @@
 		Rsqrsumm += Rsqryyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
 	Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
@@ -700,5 +700,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -717,13 +717,13 @@
 	double summ,zi,yyy,answer,oblatevol;			//running tally of integration
 	double Pi;
-	
-	Pi = 4.0*atan(1.0);
-	
+
+	Pi = 4.0*atan(1.0);
+
 	lolim = 0.0;
 	uplim = 1.0;
-	
+
 	summ = 0.0;			//initialize intergral
-	
-	
+
+
 	scale = dp[0];		//make local copies in case memory moves
 	crmaj = dp[1];
@@ -733,6 +733,6 @@
 	delpc = dp[5];
 	delps = dp[6];
-	bkg = dp[7];                        
-	
+	bkg = dp[7];
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -740,10 +740,10 @@
 		summ += yyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
 	// normalize by particle volume
 	oblatevol = 4*Pi/3*trmaj*trmaj*trmin;
 	answer /= oblatevol;
-	
+
 	//convert to [cm-1]
 	answer *= 1.0e8;
@@ -752,5 +752,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -769,12 +769,12 @@
 	double summ,zi,yyy,answer,prolatevol;			//running tally of integration
 	double Pi;
-	
-	Pi = 4.0*atan(1.0);
-	
+
+	Pi = 4.0*atan(1.0);
+
 	lolim = 0.0;
 	uplim = 1.0;
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];		//make local copies in case memory moves
 	crmaj = dp[1];
@@ -784,6 +784,6 @@
 	delpc = dp[5];
 	delps = dp[6];
-	bkg = dp[7];                        
-	
+	bkg = dp[7];
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -791,10 +791,10 @@
 		summ += yyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
 	// normalize by particle volume
 	prolatevol = 4*Pi/3*trmaj*trmin*trmin;
 	answer /= prolatevol;
-	
+
 	//convert to [cm-1]
 	answer *= 1.0e8;
@@ -803,5 +803,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -820,13 +820,13 @@
 	int nord=76;			//order of integration
 	double Pi;
-	
-	
-	Pi = 4.0*atan(1.0);
-	
+
+
+	Pi = 4.0*atan(1.0);
+
 	va = 0.0;
 	vb = Pi/2.0;
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];
 	rcore = dp[1];
@@ -839,8 +839,8 @@
 	gsd = dp[8];
 	bkg = dp[9];
-	
+
 	d=2.0*thick+length;
 	halfheight = length/2.0;
-	
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0;
@@ -848,7 +848,7 @@
 		summ += yyy;
 	}
-	
+
 	answer = (vb-va)/2.0*summ;
-	// length is the total core length 
+	// length is the total core length
 	vcyl=Pi*rcore*rcore*(2.0*thick+length)*N;
 	answer /= vcyl;
@@ -859,5 +859,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -873,6 +873,6 @@
 	double inten, qval,Pq;
 	double Pi;
-	
-	
+
+
 	Pi = 4.0*atan(1.0);
 	scale = dp[0];
@@ -882,13 +882,13 @@
 	bkg = dp[4];
 	qval = q;
-	
+
 	Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
-	
+
 	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);
 }
@@ -906,12 +906,12 @@
 	double Pi,Euler,dQDefault,fii;
 	int ii,NNint;
-	
+
 	Euler = 0.5772156649;		// Euler's constant
-	dQDefault = 0.0025;		//[=] 1/A, q-resolution, default value
+	dQDefault = 0.0;//0.0025;		//[=] 1/A, q-resolution, default value
 	dQ = dQDefault;
-	
+
 	Pi = 4.0*atan(1.0);
 	qval = q;
-	
+
 	scale = dp[0];
 	dd = dp[1];
@@ -922,14 +922,14 @@
 	Cp = dp[6];
 	bkg = dp[7];
-	
+
 	Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
-	
+
 	NNint = (int)NN;		//cast to an integer for the loop
 	ii=0;
 	Sq = 0.0;
 	for(ii=1;ii<(NNint-1);ii+=1) {
-        
+
 		fii = (double)ii;		//do I really need to do this?
-		
+
 		temp = 0.0;
 		alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
@@ -937,20 +937,20 @@
 		t2 = 2.0*qval*qval*dd*dd*alpha;
 		t3 = dQ*dQ*dd*dd*ii*ii;
-		
+
 		temp = 1.0-ii/NN;
 		temp *= cos(dd*qval*ii/(1.0+t1));
 		temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
 		temp /= sqrt(1.0+t1);
-		
+
 		Sq += temp;
 	}
-	
+
 	Sq *= 2.0;
 	Sq += 1.0;
-	
+
 	inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
-	
+
 	inten *= 1.0e8;		// 1/A to 1/cm
-	
+
     return(inten+bkg);
 }
@@ -969,13 +969,13 @@
 	double Pi,Euler,dQDefault,fii;
 	int ii,NNint;
-	
-	
+
+
 	Euler = 0.5772156649;		// Euler's constant
 	dQDefault = 0.0025;		//[=] 1/A, q-resolution, default value
 	dQ = dQDefault;
-	
+
 	Pi = 4.0*atan(1.0);
 	qval= q;
-	
+
 	scale = dp[0];
 	dd = dp[1];
@@ -988,20 +988,20 @@
 	Cp = dp[8];
 	bkg = dp[9];
-	
-	
+
+
 	drh = SLD_H - SLD_S;
 	drt = SLD_T - SLD_S;	//correction 13FEB06 by L.Porcar
-	
+
 	Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
 	Pq *= Pq;
 	Pq *= 4.0/(qval*qval);
-	
+
 	NNint = (int)NN;		//cast to an integer for the loop
 	ii=0;
 	Sq = 0.0;
 	for(ii=1;ii<(NNint-1);ii+=1) {
-        
+
 		fii = (double)ii;		//do I really need to do this?
-		
+
 		temp = 0.0;
 		alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
@@ -1009,20 +1009,20 @@
 		t2 = 2.0*qval*qval*dd*dd*alpha;
 		t3 = dQ*dQ*dd*dd*ii*ii;
-		
+
 		temp = 1.0-ii/NN;
 		temp *= cos(dd*qval*ii/(1.0+t1));
 		temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
 		temp /= sqrt(1.0+t1);
-		
+
 		Sq += temp;
 	}
-	
+
 	Sq *= 2.0;
 	Sq += 1.0;
-	
+
 	inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
-	
+
 	inten *= 1.0e8;		// 1/A to 1/cm
-	
+
 	return(inten+bkg);
 }
@@ -1038,6 +1038,6 @@
 	double inten, qval,Pq,drh,drt;
 	double Pi;
-	
-	
+
+
 	Pi = 4.0*atan(1.0);
 	qval= q;
@@ -1049,19 +1049,19 @@
 	slds = dp[5];
 	bkg = dp[6];
-	
-	
+
+
 	drh = sldh - slds;
 	drt = sldt - slds;		//correction 13FEB06 by L.Porcar
-	
+
 	Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
 	Pq *= Pq;
 	Pq *= 4.0/(qval*qval);
-	
+
 	inten = 2.0*Pi*scale*Pq/(qval*qval);		//dimensionless...
-	
+
 	inten /= 2.0*(delT+delH);			//normalize by the bilayer thickness
-	
+
 	inten *= 1.0e8;		// 1/A to 1/cm
-	
+
 	return(inten+bkg);
 }
@@ -1076,8 +1076,8 @@
 	double scale,L,B,bkg,rad,qr,cont;
 	double Pi,flex,crossSect,answer;
-	
-	
-	Pi = 4.0*atan(1.0);
-	
+
+
+	Pi = 4.0*atan(1.0);
+
 	scale = dp[0];		//make local copies in case memory moves
 	L = dp[1];
@@ -1086,10 +1086,10 @@
 	cont = dp[4];
 	bkg = dp[5];
-	
-    
+
+
 	qr = q*rad;
-	
+
 	flex = Sk_WR(q,L,B);
-	
+
 	crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
 	flex *= crossSect;
@@ -1098,5 +1098,5 @@
 	flex *= 1.0e8;
 	answer = scale*flex + bkg;
-	
+
 	return answer;
 }
@@ -1111,6 +1111,6 @@
 	double scale,L,B,bkg,rad,qr,cont,ellRatio;
 	double Pi,flex,crossSect,answer;
-	
-	
+
+
 	Pi = 4.0*atan(1.0);
 	scale = dp[0];		//make local copies in case memory moves
@@ -1121,9 +1121,9 @@
 	cont = dp[5];
 	bkg = dp[6];
-	
+
 	qr = q*rad;
-	
+
 	flex = Sk_WR(q,L,B);
-	
+
 	crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio));
 	flex *= crossSect;
@@ -1132,5 +1132,5 @@
 	flex *= 1.0e8;
 	answer = scale*flex + bkg;
-	
+
 	return answer;
 }
@@ -1141,10 +1141,10 @@
     double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
     int i,nord=76;
-    
+
     Pi = 4.0*atan(1.0);
     lolim=0.0;
     uplim=Pi/2.0;
     summ=0.0;
-    
+
     for(i=0;i<nord;i++) {
 		zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -1157,5 +1157,5 @@
     answer *= 2.0/Pi;
     return(answer);
-	
+
 }
 /*	FlexCyl_PolyLenX  :  calculates the form factor of a flecible cylinder at the given x-value p->x
@@ -1172,8 +1172,8 @@
 	double summ,zi,yyy,answer,Vpoly;			//running tally of integration
 	double range,zz,Pi;
-	
+
 	Pi = 4.0*atan(1.0);
 	range = 3.4;
-	
+
 	summ = 0.0;			//initialize intergral
 	scale = dp[0];			//make local copies in case memory moves
@@ -1184,7 +1184,7 @@
 	delrho = dp[5];
 	bkg = dp[6];
-	
+
 	zz = (1.0/pd)*(1.0/pd) - 1.0;
-	
+
 	lolim = length*(1.0-range*pd);		//set the upper/lower limits to cover the distribution
 	if(lolim<0) {
@@ -1195,5 +1195,5 @@
 	}
 	uplim = length*(1.0+range*pd);
-	
+
 	for(i=0;i<nord;i++) {
 		zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -1201,12 +1201,12 @@
 		summ += yyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
 	//normalize by average cylinder volume (first moment), using the average length
 	Vpoly=Pi*radius*radius*length;
 	answer /= Vpoly;
-	
+
 	answer *=delrho*delrho;
-	
+
 	//convert to [cm-1]
 	answer *= 1.0e8;
@@ -1215,5 +1215,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -1232,11 +1232,11 @@
 	double summ,zi,yyy,answer,Vpoly;			//running tally of integration
 	double range,zz,Pi;
-	
-	
+
+
 	Pi = 4.0*atan(1.0);
 	range = 3.4;
-	
+
 	summ = 0.0;			//initialize intergral
-	
+
 	scale = dp[0];			//make local copies in case memory moves
 	length = dp[1];			//radius
@@ -1246,7 +1246,7 @@
 	delrho = dp[5];
 	bkg = dp[6];
-	
+
 	zz = (1.0/pd)*(1.0/pd) - 1.0;
-	
+
 	lolim = radius*(1.0-range*pd);		//set the upper/lower limits to cover the distribution
 	if(lolim<0) {
@@ -1257,5 +1257,5 @@
 	}
 	uplim = radius*(1.0+range*pd);
-	
+
 	for(i=0;i<nord;i++) {
 		//zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -1265,12 +1265,12 @@
 		summ += yyy;
 	}
-	
+
 	answer = (uplim-lolim)/2.0*summ;
 	//normalize by average cylinder volume (second moment), using the average radius
 	Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
 	answer /= Vpoly;
-	
+
 	answer *=delrho*delrho;
-	
+
 	//convert to [cm-1]
 	answer *= 1.0e8;
@@ -1279,5 +1279,5 @@
 	// add in the background
 	answer += bkg;
-	
+
 	return answer;
 }
@@ -1290,7 +1290,7 @@
 	double p1,p2,p1short,p2short,q0,qconnect;
 	double C,epsilon,ans,q0short,Sexvmodify,pi;
-	
+
 	pi = 4.0*atan(1.0);
-	
+
 	p1 = 4.12;
 	p2 = 4.42;
@@ -1299,7 +1299,7 @@
 	q0 = 3.1;
 	qconnect = q0/b;
-	//	
+	//
 	q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
-	
+
 	//
 	if(L/b > 10.0) {
@@ -1311,5 +1311,5 @@
 	}
 	//
-	
+
 	if( L > 4*b ) { // Longer Chains
 		if (q*b <= 3.1) {		//Modified by Yun on Oct. 15,
@@ -1318,5 +1318,5 @@
 		} else { //q(i)*b > 3.1
 			ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L);
-		} 
+		}
 	} else { //L <= 4*b Shorter Chains
 		if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
@@ -1330,5 +1330,5 @@
 		}
 	}
-	
+
 	return(ans);
 	//return(a2long(q, L, b, p1, p2, q0));
@@ -1341,5 +1341,5 @@
     double yy;
     yy = 0.5*(1 + tanh((x - 1.523)/0.1477));
-	
+
     return (yy);
 }
@@ -1350,7 +1350,7 @@
 {
     double yy;
-	
+
     yy = Rgsquareshort(q,L,b)*q*q;
-    
+
     return (yy);
 }
@@ -1370,8 +1370,8 @@
 static double
 Rgsquarezero(double q, double L, double b)
-{    
+{
     double yy;
     yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b))));
-    
+
     return (yy);
 }
@@ -1380,8 +1380,8 @@
 static double
 Rgsquareshort(double q, double L, double b)
-{    
+{
     double yy;
     yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);
-    
+
     return (yy);
 }
@@ -1393,5 +1393,5 @@
     double yy;
     yy = AlphaSquare(L/b)*L*b/6.0;
-    
+
     return (yy);
 }
@@ -1400,8 +1400,8 @@
 static double
 AlphaSquare(double x)
-{    
+{
     double yy;
     yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
-	
+
     return (yy);
 }
@@ -1410,8 +1410,8 @@
 static double
 miu(double x)
-{    
+{
     double yy;
     yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x)));
-    
+
     return (yy);
 }
@@ -1420,8 +1420,8 @@
 static double
 Sdebye(double q, double L, double b)
-{    
+{
     double yy;
     yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
-	
+
     return (yy);
 }
@@ -1430,8 +1430,8 @@
 static double
 Sdebye1(double q, double L, double b)
-{    
+{
     double yy;
     yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
-    
+
     return (yy);
 }
@@ -1440,5 +1440,5 @@
 static double
 Sexv(double q, double L, double b)
-{    
+{
     double yy,C1,C2,C3,miu,Rg2;
     C1=1.22;
@@ -1446,9 +1446,9 @@
     C3=-1.651;
     miu = 0.585;
-	
+
     Rg2 = Rgsquare(q,L,b);
-    
+
     yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +  C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +    C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
-    
+
     return (yy);
 }
@@ -1457,16 +1457,16 @@
 static double
 Sexvnew(double q, double L, double b)
-{    
+{
     double yy,C1,C2,C3,miu;
     double del=1.05,C_star2,Rg2;
-    
+
     C1=1.22;
     C2=0.4288;
     C3=-1.651;
     miu = 0.585;
-	
+
     //calculating the derivative to decide on the corection (cutoff) term?
     // I have modified this from WRs original code
-    
+
     if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
         C_star2 = 0.0;
@@ -1474,9 +1474,9 @@
         C_star2 = 1.0;
     }
-	
+
     Rg2 = Rgsquare(q,L,b);
-    
+
 	yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
-	
+
     return (yy);
 }
@@ -1485,9 +1485,9 @@
 static double
 a2short(double q, double L, double b, double p1short, double p2short, double q0)
-{    
+{
     double yy,Rg2_sh;
     double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
     double pi;
-	
+
     E = 2.718281828459045091;
 	pi = 4.0*atan(1.0);
@@ -1496,11 +1496,11 @@
     t1 = ((q0*q0*Rg2_sh)/(b*b));
     Et1 = pow(E,t1);
-    Emt1 =pow(E,-t1); 
+    Emt1 =pow(E,-t1);
     q02 = q0*q0;
     q0p = pow(q0,(-4.0 + p2short) );
-    
+
     //E is the number e
     yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2)))))));
-	
+
     return (yy);
 }
@@ -1509,9 +1509,9 @@
 static double
 a1short(double q, double L, double b, double p1short, double p2short, double q0)
-{    
+{
     double yy,Rg2_sh;
     double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
 	double pi;
-    
+
     E = 2.718281828459045091;
 	pi = 4.0*atan(1.0);
@@ -1521,10 +1521,10 @@
     t1 = ((q0*q0*Rg2_sh)/(b*b));
     Et1 = pow(E,t1);
-    Emt1 =pow(E,-t1); 
+    Emt1 =pow(E,-t1);
     q02 = q0*q0;
     q0p = pow(q0,(-4.0 + p1short) );
-    
-    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2)))))); 
-	
+
+    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2))))));
+
     return(yy);
 }
@@ -1537,5 +1537,5 @@
     double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
     double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
-    
+
     pi = 4.0*atan(1.0);
     E = 2.718281828459045091;
@@ -1545,5 +1545,5 @@
         C = 1.0;
     }
-	
+
     C1 = 1.22;
     C2 = 0.4288;
@@ -1552,5 +1552,5 @@
     C5 = 0.1477;
     miu = 0.585;
-	
+
     Rg2 = Rgsquare(q,L,b);
     Rg22 = Rg2*Rg2;
@@ -1562,28 +1562,28 @@
     q04 = q03*q0;
     q05 = q04*q0;
-    
+
     t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));
-    
+
     t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14*b3*pow(E,(-((q02*Rg2)/b2))))/(15*q03*Rg2) + (2*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15*q02*Rg2)))*Rg2)/b)))/L;
-    
+
     t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(2*C5);
-    
+
     t4 = (b4*sqrt(Rg2)*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);
-    
+
     t5 = (2*b4*(((2*q0*Rg2)/b - (2*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
-    
+
     t6 = (8*b4*b*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);
-    
+
     t7 = (((-((3*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 3/miu)))/miu)) - (2*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 2/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 1/miu)))/miu));
-    
+
     t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-    
+
     t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7*b2)/(15*q02*Rg2))) + (7*b2)/(15*q02*Rg2))))/L;
-    
+
     t10 = (2*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
-	
-    
+
+
     yy = ((-1*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1)/miu))))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
-	
+
     return (yy);
 }
@@ -1592,5 +1592,5 @@
 static double
 sech_WR(double x)
-{	
+{
 	return(1/cosh(x));
 }
@@ -1599,13 +1599,13 @@
 static double
 a1long(double q, double L, double b, double p1, double p2, double q0)
-{    
+{
     double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
     double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
     double E,pi;
     double b2,b3,b4,q02,q03,q04,q05,Rg22;
-	
+
     pi = 4.0*atan(1.0);
     E = 2.718281828459045091;
-    
+
     if( L/b > 10.0) {
         C = 3.06/pow((L/b),0.44);
@@ -1613,5 +1613,5 @@
         C = 1.0;
     }
-	
+
     C1 = 1.22;
     C2 = 0.4288;
@@ -1620,5 +1620,5 @@
     C5 = 0.1477;
     miu = 0.585;
-	
+
     Rg2 = Rgsquare(q,L,b);
     Rg22 = Rg2*Rg2;
@@ -1630,38 +1630,38 @@
     q04 = q03*q0;
     q05 = q04*q0;
-    
+
     t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
-    
+
     t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
-    
+
     t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
-    
+
     t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-	
+
     t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));
-    
+
     t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b)));
-    
+
     t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
-    
+
     t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
-    
+
     t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
-    
+
     t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
-    
+
     t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));
-    
+
     t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-    
+
     t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
-    
+
     t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
-    
+
     t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
-	
-    
+
+
     yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))));
-	
+
     return (yy);
 }
@@ -1673,5 +1673,5 @@
 //
 //     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
-//                       BY (53) AND (56,57) IN CHEN AND 
+//                       BY (53) AND (56,57) IN CHEN AND
 //                       KOTLARCHYK REFERENCE
 //
@@ -1683,7 +1683,7 @@
 	// local variables
 	double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,gfn2,pi43,gfn,Pi;
-	
-	Pi = 4.0*atan(1.0);
-	
+
+	Pi = 4.0*atan(1.0);
+
 	pi43=4.0/3.0*Pi;
 	aa = crmaj;
@@ -1699,5 +1699,5 @@
 	gfn = gfnc+gfnt;
 	gfn2 = gfn*gfn;
-	
+
 	return (gfn2);
 }
@@ -1709,5 +1709,5 @@
 //
 //       <OBLATE ELLIPSOID>
-// function gfn4 for oblate ellipsoids 
+// function gfn4 for oblate ellipsoids
 double
 gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
@@ -1715,5 +1715,5 @@
 	// local variables
 	double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
-	
+
 	Pi = 4.0*atan(1.0);
 	pi43=4.0/3.0*Pi;
@@ -1730,5 +1730,5 @@
   	tgfn = gfnc+gfnt;
   	gfn4 = tgfn*tgfn;
-  	
+
   	return (gfn4);
 }
@@ -1739,17 +1739,17 @@
     double Pq,vcyl,dl;
     double Pi,qr;
-    
+
     Pi = 4.0*atan(1.0);
     qr = q*radius;
-    
+
     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);
-    
+
     vcyl=Pi*radius*radius*zi;
     Pq *= vcyl*vcyl;
-    
+
     dl = SchulzPoint_cpr(zi,length,zz);
-    return (Pq*dl);	
-	
+    return (Pq*dl);
+
 }
 
@@ -1759,17 +1759,17 @@
     double Pq,vcyl,dr;
     double Pi,qr;
-    
+
     Pi = 4.0*atan(1.0);
     qr = q*zi;
-    
+
     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);
-    
+
     vcyl=Pi*zi*zi*Lc;
     Pq *= vcyl*vcyl;
-    
+
     dr = SchulzPoint_cpr(zi,ravg,zz);
-    return (Pq*dr);	
-	
+    return (Pq*dr);
+
 }
 
@@ -1780,6 +1780,6 @@
 	double lolim,uplim,summ,yyy,zi;
 	int nord,i;
-	
-	// set up the integration end points 
+
+	// set up the integration end points
 	Pi = 4.0*atan(1.0);
 	nord = 76;
@@ -1787,5 +1787,5 @@
 	uplim = Pi/2;
 	halfheight = length/2.0;
-	
+
 	summ = 0.0;				// initialize integral
 	i=0;
@@ -1795,5 +1795,5 @@
 		summ += yyy;
 	}
-	
+
 	// calculate value of integral to return
 	answer = (uplim-lolim)/2.0*summ;
@@ -1803,33 +1803,33 @@
 double
 CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum)
-{	
+{
 	// qq is the q-value for the calculation (1/A)
 	// radius is the core radius of the cylinder (A)
 	//  radthick and facthick are the radial and face layer thicknesses
 	// rho(n) are the respective SLD's
-	// length is the *Half* CORE-LENGTH of the cylinder 
+	// length is the *Half* CORE-LENGTH of the cylinder
 	// dum is the dummy variable for the integration (theta)
-	
+
 	double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval;
 	double Pi;
-	
-	Pi = 4.0*atan(1.0); 
-	
+
+	Pi = 4.0*atan(1.0);
+
 	dr1 = rhoc-rhos;
 	dr2 = rhos-rhosolv;
 	vol1 = Pi*rad*rad*(2.0*length);
 	vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
-	
+
 	besarg1 = qq*rad*sin(dum);
 	besarg2 = qq*(rad+radthick)*sin(dum);
 	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;
-	
+
 	retval = ((t1+t2)*(t1+t2))*sin(dum);
 	return (retval);
-    
+
 }
 
@@ -1838,25 +1838,25 @@
 CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum)
 {
-	
+
     double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval;
     double Pi;
-    
+
     Pi = 4.0*atan(1.0);
-    
+
     dr1 = rhoc-rhos;
     dr2 = rhos-rhosolv;
     vol1 = Pi*rcore*rcore*(2.0*length);
     vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick);
-    
+
     besarg1 = qq*rcore*sin(dum);
     besarg2 = qq*(rcore+thick)*sin(dum);
     sinarg1 = qq*length*cos(dum);
     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;
-    
+
     retval = ((t1+t2)*(t1+t2))*sin(dum);
-	
+
     return (retval);
 }
@@ -1865,9 +1865,9 @@
 Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen)
 {
-	
+
     double halfheight,uplim,lolim,zi,summ,yyy,Pi;
     double answer,dr,Vcyl;
     int i,nord;
-    
+
     Pi = 4.0*atan(1.0);
     lolim = 0;
@@ -1876,5 +1876,5 @@
     nord=20;
     summ = 0.0;
-    
+
     //do the cylinder orientational average
     for(i=0;i<nord;i++) {
@@ -1890,5 +1890,5 @@
     Vcyl = Pi*radius*radius*dumLen;
     answer *= Vcyl*Vcyl;
-    
+
     dr = SchulzPoint_cpr(dumLen,len_avg,zz);
     return(dr*answer);
@@ -1898,5 +1898,5 @@
 double
 Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N)
-{         	
+{
 	// qq is the q-value for the calculation (1/A)
 	// rcore is the core radius of the cylinder (A)
@@ -1904,28 +1904,28 @@
 	// length is the *Half* CORE-LENGTH of the cylinder = L (A)
 	// dum is the dummy variable for the integration (x in Feigin's notation)
-	
-	//Local variables 
+
+	//Local variables
 	double totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
 	double Pi;
 	int kk;
-	
-	Pi = 4.0*atan(1.0);
-	
+
+	Pi = 4.0*atan(1.0);
+
 	dr1 = rhoc-rhosolv;
 	dr2 = rhol-rhosolv;
 	area = Pi*rcore*rcore;
 	totald=2.0*(thick+length);
-	
+
 	besarg1 = qq*rcore*sin(dum);
 	besarg2 = qq*rcore*sin(dum);
-	
+
 	sinarg1 = qq*length*cos(dum);
 	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);
-	
+
 	retval =((t1+t2)*(t1+t2))*sin(dum);
-	
+
 	// loop for the structure facture S(q)
 	sqq=0.0;
@@ -1933,10 +1933,10 @@
 		dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0;
 		sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt);
-	}			
-	
+	}
+
 	// end of loop for S(q)
 	sqq=1.0+2.0*sqq/N;
 	retval *= sqq;
-    
+
 	return(retval);
 }
@@ -1946,9 +1946,9 @@
 Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad)
 {
-	
+
     double halfheight,uplim,lolim,zi,summ,yyy,Pi;
     double answer,dr,Vcyl;
     int i,nord;
-    
+
     Pi = 4.0*atan(1.0);
     lolim = 0;
@@ -1958,5 +1958,5 @@
     nord=76;
     summ = 0.0;
-    
+
     //do the cylinder orientational average
 	//    for(i=0;i<nord;i++) {
@@ -1977,5 +1977,5 @@
     Vcyl = Pi*dumRad*dumRad*length;
     answer *= Vcyl*Vcyl;
-    
+
     dr = SchulzPoint_cpr(dumRad,radius,zz);
     return(dr*answer);
@@ -1986,5 +1986,5 @@
 {
     double dr;
-    
+
     dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0));
     return(exp(dr));
@@ -1999,5 +1999,5 @@
 		0.1208650973866179e-2,-0.5395239384953e-5};
 	int j;
-	
+
 	y=x=xx;
 	tmp=x+5.5;
@@ -2013,12 +2013,12 @@
 {
     double arg,nu,retval;		//local variables
-	
+
     nu = nua/a;
     arg = qq*a*sqrt(1+dum*dum*(nu*nu-1));
-    
+
     retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
     retval *= retval;
     retval *= 9;
-    
+
     return(retval);
 }//Function EllipsoidKernel()
@@ -2028,5 +2028,5 @@
 {
     double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;		//local variables
-    
+
     gamma = rcore/rshell;
     arg1 = qq*rshell*sqrt(1-dum*dum);		//1=shell (outer radius)
@@ -2035,10 +2035,10 @@
     lam2 = 2*NR_BessJ1(arg2)/arg2;
     psi = 1/(1-gamma*gamma)*(lam1 -  gamma*gamma*lam2);		//SRK 10/19/00
-    
+
     sinarg = qq*length*dum/2;
     t2 = sin(sinarg)/sinarg;
-    
+
     retval = psi*psi*t2*t2;
-    
+
     return(retval);
 }//Function HolCylKernel()
@@ -2049,7 +2049,7 @@
     // mu passed in is really mu*sqrt(1-sig^2)
     double arg1,arg2,Pi,tmp1,tmp2;			//local variables
-	
+
     Pi = 4.0*atan(1.0);
-    
+
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
     arg1 = (mu/2)*cos(Pi*uu/2);
@@ -2060,5 +2060,5 @@
 		tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
     }
-    
+
     if (arg2==0) {
 		tmp2 = 1;
@@ -2066,7 +2066,7 @@
 		tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
     }
-	
+
     return (tmp1*tmp2);
-	
+
 }//Function PPKernel()
 
@@ -2075,9 +2075,9 @@
 TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy)
 {
-	
+
     double arg,val,pi;			//local variables
-	
+
     pi = 4.0*atan(1.0);
-    
+
     arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2);
     arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy);
@@ -2085,7 +2085,7 @@
     arg = q*sqrt(arg);
     val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
-    
+
     return (val);
-	
+
 }//Function TriaxialKernel()
 
@@ -2094,16 +2094,16 @@
 CylKernel(double qq, double rr,double h, double theta)
 {
-	
+
 	// qq is the q-value for the calculation (1/A)
 	// rr is the radius of the cylinder (A)
 	// 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;		 //Local variables
+
+
     besarg = qq*rr*sin(theta);
-	
+
     bj =NR_BessJ1(besarg);
-	
+
 	//* Computing 2nd power */
     d1 = sin(qq * h * cos(theta));
@@ -2119,7 +2119,7 @@
     b2 = d1 * d1;
     retval = t1 * t2 / b1 / b2;
-	
+
     return (retval);
-	
+
 }//Function CylKernel()
 
@@ -2127,21 +2127,21 @@
 EllipCylKernel(double qq, double ra,double nu, double theta)
 {
-	//this is the function LAMBDA1^2 in Feigin's notation	
+	//this is the function LAMBDA1^2 in Feigin's notation
 	// qq is the q-value for the calculation (1/A)
 	// ra is the transformed radius"a" in Feigin's notation
 	// nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
 	// theta is the dummy variable of the integration
-	
-	double retval,arg;		 //Local variables 
-	
+
+	double retval,arg;		 //Local variables
+
 	arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2);
-	
+
 	retval = 2*NR_BessJ1(arg)/arg;
-	
+
 	//square it
 	retval *= retval;
-	
+
     return(retval);
-	
+
 }//Function EllipCylKernel()
 
@@ -2150,5 +2150,5 @@
 	double ax,z;
 	double xx,y,ans,ans1,ans2;
-	
+
 	if ((ax=fabs(x)) < 8.0) {
 		y=x*x;
@@ -2170,5 +2170,5 @@
 		if (x < 0.0) ans = -ans;
 	}
-	
+
 	return(ans);
 }
