libCylinder.c @ cd2f75d

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since cd2f75d was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago
Reorganizing models in preparation of cpp cleanup
Property mode set to `100644`
File size: 80.4 KB

Line
1	/* CylinderFit.c
2
3	A simplified project designed to act as a template for your curve fitting function.
4	The fitting function is a Cylinder form factor. No resolution effects are included (yet)
5	*/
6
7	#include "StandardHeaders.h" // Include ANSI headers, Mac headers
8	#include "GaussWeights.h"
9	#include "libCylinder.h"
10
11	/////////functions for WRC implementation of flexible cylinders
12	static double
13	gammaln(double xx)
14	{
15	double x,y,tmp,ser;
16	static double cof[6]={76.18009172947146,-86.50532032941677,
17	24.01409824083091,-1.231739572450155,
18	0.1208650973866179e-2,-0.5395239384953e-5};
19	int j;
20
21	y=x=xx;
22	tmp=x+5.5;
23	tmp -= (x+0.5)*log(tmp);
24	ser=1.000000000190015;
25	for (j=0;j<=5;j++) ser += cof[j]/++y;
26	return -tmp+log(2.5066282746310005*ser/x);
27	}
28	/////////functions for WRC implementation of flexible cylinders
29
30	//
31	static double
32	AlphaSquare(double x)
33	{
34	double yy;
35	yy = pow( (1.0 + (x/3.12)(x/3.12) + (x/8.67)(x/8.67)*(x/8.67)),(0.176/3.0) );
36
37	return (yy);
38	}
39
40	//
41	static double
42	Rgsquarezero(double q, double L, double b)
43	{
44	double yy;
45	yy = (Lb/6.0) (1.0 - 1.5(b/L) + 1.5pow((b/L),2) - 0.75pow((b/L),3)(1.0 - exp(-2.0*(L/b))));
46
47	return (yy);
48	}
49
50	//
51	static double
52	Rgsquareshort(double q, double L, double b)
53	{
54	double yy;
55	yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);
56
57	return (yy);
58	}
59
60	//
61	static double
62	Rgsquare(double q, double L, double b)
63	{
64	double yy;
65	yy = AlphaSquare(L/b)Lb/6.0;
66
67	return (yy);
68	}
69
70	// ?? funciton is not used - but should the log actually be log10???
71	/*
72	static double
73	miu(double x)
74	{
75	double yy;
76	yy = (1.0/8.0)(9.0x - 2.0 + 2.0log(1.0 + x)/x)exp(1.0/2.565(1.0/x + (1.0 - 1.0/(xx))*log(1.0 + x)));
77
78	return (yy);
79	}
80	*/
81
82	//WR named this w (too generic)
83	static double
84	w_WR(double x)
85	{
86	double yy;
87	yy = 0.5*(1 + tanh((x - 1.523)/0.1477));
88
89	return (yy);
90	}
91
92	//
93	static double
94	u1(double q, double L, double b)
95	{
96	double yy;
97
98	yy = Rgsquareshort(q,L,b)qq;
99
100	return (yy);
101	}
102
103	// was named u
104	static double
105	u_WR(double q, double L, double b)
106	{
107	double yy;
108	yy = Rgsquare(q,L,b)qq;
109	return (yy);
110	}
111
112
113	//
114	static double
115	Sdebye1(double q, double L, double b)
116	{
117	double yy;
118	yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
119
120	return (yy);
121	}
122
123
124	//
125	static double
126	Sdebye(double q, double L, double b)
127	{
128	double yy;
129	yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
130
131	return (yy);
132	}
133
134	//
135	static double
136	Sexv(double q, double L, double b)
137	{
138	double yy,C1,C2,C3,miu,Rg2;
139	C1=1.22;
140	C2=0.4288;
141	C3=-1.651;
142	miu = 0.585;
143
144	Rg2 = Rgsquare(q,L,b);
145
146	yy = (1.0 - w_WR(qsqrt(Rg2)))Sdebye(q,L,b) + w_WR(qsqrt(Rg2))(C1pow((qsqrt(Rg2)),(-1.0/miu)) + C2pow((qsqrt(Rg2)),(-2.0/miu)) + C3pow((qsqrt(Rg2)),(-3.0/miu)));
147
148	return (yy);
149	}
150
151	// this must be WR modified version
152	static double
153	Sexvnew(double q, double L, double b)
154	{
155	double yy,C1,C2,C3,miu;
156	double del=1.05,C_star2,Rg2;
157
158	C1=1.22;
159	C2=0.4288;
160	C3=-1.651;
161	miu = 0.585;
162
163	//calculating the derivative to decide on the corection (cutoff) term?
164	// I have modified this from WRs original code
165
166	if( (Sexv(qdel,L,b)-Sexv(q,L,b))/(qdel - q) >= 0.0 ) {
167	C_star2 = 0.0;
168	} else {
169	C_star2 = 1.0;
170	}
171
172	Rg2 = Rgsquare(q,L,b);
173
174	yy = (1.0 - w_WR(qsqrt(Rg2)))Sdebye(q,L,b) + C_star2w_WR(qsqrt(Rg2))(C1pow((qsqrt(Rg2)),(-1.0/miu)) + C2pow((qsqrt(Rg2)),(-2.0/miu)) + C3pow((q*sqrt(Rg2)),(-3.0/miu)));
175
176	return (yy);
177	}
178
179	// these are the messy ones
180	static double
181	a2short(double q, double L, double b, double p1short, double p2short, double q0)
182	{
183	double yy,Rg2_sh;
184	double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
185	double pi;
186
187	E = 2.718281828459045091;
188	pi = 4.0*atan(1.0);
189	Rg2_sh = Rgsquareshort(q,L,b);
190	Rg2_sh2 = Rg2_sh*Rg2_sh;
191	t1 = ((q0q0Rg2_sh)/(b*b));
192	Et1 = pow(E,t1);
193	Emt1 =pow(E,-t1);
194	q02 = q0*q0;
195	q0p = pow(q0,(-4.0 + p2short) );
196
197	//E is the number e
198	yy = ((-(1.0/(L((p1short - p2short))Rg2_sh2)((bEmt1q0p((8.0bbbL - 8.0bbbEt1L - 2.0bbbLp1short + 2.0bbbEt1Lp1short + 4.0bLq02Rg2_sh + 4.0bEt1Lq02Rg2_sh - 2.0bEt1Lp1shortq02Rg2_sh - Et1piq02q0Rg2_sh2 + Et1p1shortpiq02q0Rg2_sh2)))))));
199
200	return (yy);
201	}
202
203	//
204	static double
205	a1short(double q, double L, double b, double p1short, double p2short, double q0)
206	{
207	double yy,Rg2_sh;
208	double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
209	double pi;
210
211	E = 2.718281828459045091;
212	pi = 4.0*atan(1.0);
213	Rg2_sh = Rgsquareshort(q,L,b);
214	Rg2_sh2 = Rg2_sh*Rg2_sh;
215	b3 = bbb;
216	t1 = ((q0q0Rg2_sh)/(b*b));
217	Et1 = pow(E,t1);
218	Emt1 =pow(E,-t1);
219	q02 = q0*q0;
220	q0p = pow(q0,(-4.0 + p1short) );
221
222	yy = ((1.0/(L((p1short - p2short))Rg2_sh2)((bEmt1q0p((8.0b3L - 8.0b3Et1L - 2.0b3Lp2short + 2.0b3Et1Lp2short + 4.0bLq02Rg2_sh + 4.0bEt1Lq02Rg2_sh - 2.0bEt1Lp2shortq02Rg2_sh - Et1piq02q0Rg2_sh2 + Et1p2shortpiq02q0Rg2_sh2))))));
223
224	return(yy);
225	}
226
227
228	//need to define this on my own
229	static double
230	sech_WR(double x)
231	{
232	return(1/cosh(x));
233	}
234
235	// this one will be lots of trouble
236	static double
237	a2long(double q, double L, double b, double p1, double p2, double q0)
238	{
239	double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
240	double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
241	double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
242
243	pi = 4.0*atan(1.0);
244	E = 2.718281828459045091;
245	if( L/b > 10.0) {
246	C = 3.06/pow((L/b),0.44);
247	} else {
248	C = 1.0;
249	}
250
251	C1 = 1.22;
252	C2 = 0.4288;
253	C3 = -1.651;
254	C4 = 1.523;
255	C5 = 0.1477;
256	miu = 0.585;
257
258	Rg2 = Rgsquare(q,L,b);
259	Rg22 = Rg2*Rg2;
260	b2 = b*b;
261	b3 = bbb;
262	b4 = b3*b;
263	q02 = q0*q0;
264	q03 = q0q0q0;
265	q04 = q03*q0;
266	q05 = q04*q0;
267
268	t1 = (1.0/(b* p1pow(q0,((-1.0) - p1 - p2)) - bp2*pow(q0,((-1.0) - p1 - p2)) ));
269
270	t2 = (bC(((-1.0((14.0b3)/(15.0q03Rg2))) + (14.0b3pow(E,(-((q02Rg2)/b2))))/(15.0q03Rg2) + (2.0pow(E,(-((q02Rg2)/b2)))q0((11.0/15.0 + (7b2)/(15.0q02Rg2)))*Rg2)/b)))/L;
271
272	t3 = (sqrt(Rg2)((C3pow((((sqrt(Rg2)q0)/b)),((-3.0)/miu)) + C2pow((((sqrt(Rg2)q0)/b)),((-2.0)/miu)) + C1pow((((sqrt(Rg2)q0)/b)),((-1.0)/miu))))pow(sech_WR(((-C4) + (sqrt(Rg2)q0)/b)/C5),2.0))/(2.0C5);
273
274	t4 = (b4sqrt(Rg2)(((-1.0) + pow(E,(-((q02Rg2)/b2))) + (q02Rg2)/b2))pow(sech_WR(((-C4) + (sqrt(Rg2)q0)/b)/C5),2))/(C5q04Rg22);
275
276	t5 = (2.0b4(((2.0q0Rg2)/b - (2.0pow(E,(-((q02Rg2)/b2)))q0Rg2)/b))((1.0 + 1.0/2.0(((-1.0) - tanh(((-C4) + (sqrt(Rg2)q0)/b)/C5))))))/(q04Rg22);
277
278	t6 = (8.0b4b(((-1.0) + pow(E,(-((q02Rg2)/b2))) + (q02Rg2)/b2))((1.0 + 1.0/2.0(((-1) - tanh(((-C4) + (sqrt(Rg2)q0)/b)/C5))))))/(q05*Rg22);
279
280	t7 = (((-((3.0C3sqrt(Rg2)pow((((sqrt(Rg2)q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0C2sqrt(Rg2)pow((((sqrt(Rg2)q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1sqrt(Rg2)pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));
281
282	t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
283
284	t9 = (bC((4.0/15.0 - pow(E,(-((q02Rg2)/b2)))((11.0/15.0 + (7.0b2)/(15q02Rg2))) + (7.0b2)/(15.0q02Rg2))))/L;
285
286	t10 = (2.0b4(((-1) + pow(E,(-((q02Rg2)/b2))) + (q02Rg2)/b2))((1.0 + 1.0/2.0(((-1) - tanh(((-C4) + (sqrt(Rg2)q0)/b)/C5))))))/(q04Rg22);
287
288
289	yy = ((-1.0(t1 ((-pow(q0,-p1)(((b2pi)/(Lq02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0t7t8)) - bp1pow(q0,((-1.0) - p1))(((-((bpi)/(Lq0))) + t9 + t10 + 1.0/2.0((C3pow((((sqrt(Rg2)q0)/b)),((-3.0)/miu)) + C2pow((((sqrt(Rg2)q0)/b)),((-2.0)/miu)) + C1pow((((sqrt(Rg2)q0)/b)),((-1.0)/miu))))((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
290
291	return (yy);
292	}
293	//
294	static double
295	a1long(double q, double L, double b, double p1, double p2, double q0)
296	{
297	double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
298	double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
299	double E,pi;
300	double b2,b3,b4,q02,q03,q04,q05,Rg22;
301
302	pi = 4.0*atan(1.0);
303	E = 2.718281828459045091;
304
305	if( L/b > 10.0) {
306	C = 3.06/pow((L/b),0.44);
307	} else {
308	C = 1.0;
309	}
310
311	C1 = 1.22;
312	C2 = 0.4288;
313	C3 = -1.651;
314	C4 = 1.523;
315	C5 = 0.1477;
316	miu = 0.585;
317
318	Rg2 = Rgsquare(q,L,b);
319	Rg22 = Rg2*Rg2;
320	b2 = b*b;
321	b3 = bbb;
322	b4 = b3*b;
323	q02 = q0*q0;
324	q03 = q0q0q0;
325	q04 = q03*q0;
326	q05 = q04*q0;
327
328	t1 = (bC((4.0/15.0 - pow(E,(-((q02Rg2)/b2)))((11.0/15.0 + (7.0b2)/(15.0q02Rg2))) + (7.0b2)/(15.0q02Rg2))));
329
330	t2 = (2.0b4(((-1.0) + pow(E,(-((q02Rg2)/b2))) + (q02Rg2)/b2))((1.0 + 1.0/2.0(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
331
332	t3 = ((C3pow((((sqrt(Rg2)q0)/b)),((-3.0)/miu)) + C2pow((((sqrt(Rg2)q0)/b)),((-2.0)/miu)) + C1pow((((sqrt(Rg2)q0)/b)),((-1.0)/miu))));
333
334	t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
335
336	t5 = (1.0/(bp1pow(q0,((-1.0) - p1 - p2)) - bp2pow(q0,((-1.0) - p1 - p2))));
337
338	t6 = (bC(((-((14.0b3)/(15.0q03Rg2))) + (14.0b3pow(E,(-((q02Rg2)/b2))))/(15.0q03Rg2) + (2.0pow(E,(-((q02Rg2)/b2)))q0((11.0/15.0 + (7.0b2)/(15.0q02Rg2)))Rg2)/b)));
339
340	t7 = (sqrt(Rg2)((C3pow((((sqrt(Rg2)q0)/b)),((-3.0)/miu)) + C2pow((((sqrt(Rg2)q0)/b)),((-2.0)/miu)) + C1pow((((sqrt(Rg2)q0)/b)),((-1.0)/miu))))pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
341
342	t8 = (b4sqrt(Rg2)(((-1.0) + pow(E,(-((q02Rg2)/b2))) + (q02Rg2)/b2))pow(sech_WR(((-C4) + (sqrt(Rg2)q0)/b)/C5),2));
343
344	t9 = (2.0b4(((2.0q0Rg2)/b - (2.0pow(E,(-((q02Rg2)/b2)))q0Rg2)/b))((1.0 + 1.0/2.0(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
345
346	t10 = (8.0b4b(((-1.0) + pow(E,(-((q02Rg2)/b2))) + (q02Rg2)/b2))((1.0 + 1.0/2.0(((-1.0) - tanh(((-C4) + (sqrt(Rg2)q0)/b)/C5))))));
347
348	t11 = (((-((3.0C3sqrt(Rg2)pow((((sqrt(Rg2)q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0C2sqrt(Rg2)pow((((sqrt(Rg2)q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1sqrt(Rg2)pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));
349
350	t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
351
352	t13 = (bC((4.0/15.0 - pow(E,(-((q02Rg2)/b2)))((11.0/15.0 + (7.0b2)/(15.0q02* Rg2))) + (7.0b2)/(15.0q02*Rg2))));
353
354	t14 = (2.0b4(((-1.0) + pow(E,(-((q02Rg2)/b2))) + (q02Rg2)/b2))((1.0 + 1.0/2.0(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
355
356	t15 = ((C3pow((((sqrt(Rg2)q0)/b)),((-3.0)/miu)) + C2pow((((sqrt(Rg2)q0)/b)),((-2.0)/miu)) + C1pow((((sqrt(Rg2)q0)/b)),((-1.0)/miu))));
357
358
359	yy = (pow(q0,p1)(((-((bpi)/(Lq0))) +t1/L +t2/(q04Rg22) + 1.0/2.0t3t4)) + (t5((pow(q0,(p1 - p2))(((-pow(q0,(-p1)))(((b2pi)/(Lq02) +t6/L +t7/(2.0C5) -t8/(C5q04Rg22) +t9/(q04Rg22) -t10/(q05Rg22) + 1.0/2.0t11t12)) - bp1pow(q0,((-1.0) - p1))(((-((bpi)/(Lq0))) +t13/L +t14/(q04Rg22) + 1.0/2.0t15((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))));
360
361	return (yy);
362	}
363
364
365	static double
366	Sk_WR(double q, double L, double b)
367	{
368	//
369	double p1,p2,p1short,p2short,q0;
370	double C,ans,q0short,Sexvmodify,pi;
371
372	pi = 4.0*atan(1.0);
373
374	p1 = 4.12;
375	p2 = 4.42;
376	p1short = 5.36;
377	p2short = 5.62;
378	q0 = 3.1;
379	//
380	q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
381
382	//
383	if(L/b > 10.0) {
384	C = 3.06/pow((L/b),0.44);
385	} else {
386	C = 1.0;
387	}
388	//
389
390	if( L > 4*b ) { // Longer Chains
391	if (q*b <= 3.1) { //Modified by Yun on Oct. 15,
392	Sexvmodify = Sexvnew(q, L, b);
393	ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0u_WR(q,L,b)))exp(-u_WR(q,L,b)))(b/L);
394	} else { //q(i)*b > 3.1
395	ans = a1long(q, L, b, p1, p2, q0)/(pow((qb),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((qb),p2)) + pi/(q*L);
396	}
397	} else { //L <= 4*b Shorter Chains
398	if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
399	if (q*b<=0.01) {
400	ans = 1.0 - Rgsquareshort(q,L,b)(qq)/3.0;
401	} else {
402	ans = Sdebye1(q,L,b);
403	}
404	} else { //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
405	ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((qb),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((qb),p2short)) + pi/(q*L);
406	}
407	}
408
409	return(ans);
410	//return(a2long(q, L, b, p1, p2, q0));
411	}
412
413	/* CylinderForm : calculates the form factor of a cylinder at the give x-value p->x
414
415	Warning:
416	The call to WaveData() below returns a pointer to the middle
417	of an unlocked Macintosh handle. In the unlikely event that your
418	calculations could cause memory to move, you should copy the coefficient
419	values to local variables or an array before such operations.
420	*/
421	double
422	CylinderForm(double dp[], double q)
423	{
424
425	int i;
426	double Pi;
427	double scale,radius,length,delrho,bkg,halfheight,sldCyl,sldSolv; //local variables of coefficient wave
428	int nord=76; //order of integration
429	double uplim,lolim; //upper and lower integration limits
430	double summ,zi,yyy,answer,vcyl; //running tally of integration
431
432	Pi = 4.0*atan(1.0);
433	lolim = 0.0;
434	uplim = Pi/2.0;
435
436	summ = 0.0; //initialize intergral
437
438	scale = dp[0]; //make local copies in case memory moves
439	radius = dp[1];
440	length = dp[2];
441	sldCyl = dp[3];
442	sldSolv = dp[4];
443	bkg = dp[5];
444
445	delrho = sldCyl-sldSolv;
446	halfheight = length/2.0;
447	for(i=0;i<nord;i++) {
448	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
449	yyy = Gauss76Wt[i] * CylKernel(q, radius, halfheight, zi);
450	summ += yyy;
451	}
452
453	answer = (uplim-lolim)/2.0*summ;
454	// Multiply by contrast^2
455	answer = delrhodelrho;
456	//normalize by cylinder volume
457	//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
458	vcyl=Piradiusradius*length;
459	answer *= vcyl;
460	//convert to [cm-1]
461	answer *= 1.0e8;
462	//Scale
463	answer *= scale;
464	// add in the background
465	answer += bkg;
466
467	return answer;
468	}
469
470	/* EllipCyl76X : calculates the form factor of a elliptical cylinder at the given x-value p->x
471
472	Uses 76 pt Gaussian quadrature for both integrals
473
474	Warning:
475	The call to WaveData() below returns a pointer to the middle
476	of an unlocked Macintosh handle. In the unlikely event that your
477	calculations could cause memory to move, you should copy the coefficient
478	values to local variables or an array before such operations.
479	*/
480	double
481	EllipCyl76(double dp[], double q)
482	{
483	int i,j;
484	double Pi,slde,sld;
485	double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave
486	int nord=76; //order of integration
487	double va,vb; //upper and lower integration limits
488	double summ,zi,yyy,answer,vell; //running tally of integration
489	double summj,vaj,vbj,zij,arg, si; //for the inner integration
490
491	Pi = 4.0*atan(1.0);
492	va = 0.0;
493	vb = 1.0; //orintational average, outer integral
494	vaj=0.0;
495	vbj=Pi; //endpoints of inner integral
496
497	summ = 0.0; //initialize intergral
498
499	scale = dp[0]; //make local copies in case memory moves
500	ra = dp[1];
501	nu = dp[2];
502	length = dp[3];
503	slde = dp[4];
504	sld = dp[5];
505	delrho = slde - sld;
506	bkg = dp[6];
507
508	for(i=0;i<nord;i++) {
509	//setup inner integral over the ellipsoidal cross-section
510	summj=0;
511	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy
512	arg = rasqrt(1.0-zizi);
513	for(j=0;j<nord;j++) {
514	//76 gauss points for the inner integral as well
515	zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy
516	yyy = Gauss76Wt[j] * EllipCylKernel(q,arg,nu,zij);
517	summj += yyy;
518	}
519	//now calculate the value of the inner integral
520	answer = (vbj-vaj)/2.0*summj;
521	//divide integral by Pi
522	answer /=Pi;
523
524	//now calculate outer integral
525	arg = qlengthzi/2.0;
526	if (arg == 0.0){
527	si = 1.0;
528	}else{
529	si = sin(arg) * sin(arg) / arg / arg;
530	}
531	yyy = Gauss76Wt[i] * answer * si;
532	summ += yyy;
533	}
534	answer = (vb-va)/2.0*summ;
535	// Multiply by contrast^2
536	answer = delrhodelrho;
537	//normalize by cylinder volume
538	//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
539	vell = Pira(nura)length;
540	answer *= vell;
541	//convert to [cm-1]
542	answer *= 1.0e8;
543	//Scale
544	answer *= scale;
545	// add in the background
546	answer += bkg;
547
548	return answer;
549	}
550
551	/* EllipCyl20X : calculates the form factor of a elliptical cylinder at the given x-value p->x
552
553	Uses 76 pt Gaussian quadrature for orientational integral
554	Uses 20 pt quadrature for the inner integral over the elliptical cross-section
555
556	Warning:
557	The call to WaveData() below returns a pointer to the middle
558	of an unlocked Macintosh handle. In the unlikely event that your
559	calculations could cause memory to move, you should copy the coefficient
560	values to local variables or an array before such operations.
561	*/
562	double
563	EllipCyl20(double dp[], double q)
564	{
565	int i,j;
566	double Pi,slde,sld;
567	double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave
568	int nordi=76; //order of integration
569	int nordj=20;
570	double va,vb; //upper and lower integration limits
571	double summ,zi,yyy,answer,vell; //running tally of integration
572	double summj,vaj,vbj,zij,arg,si; //for the inner integration
573
574	Pi = 4.0*atan(1.0);
575	va = 0.0;
576	vb = 1.0; //orintational average, outer integral
577	vaj=0.0;
578	vbj=Pi; //endpoints of inner integral
579
580	summ = 0.0; //initialize intergral
581
582	scale = dp[0]; //make local copies in case memory moves
583	ra = dp[1];
584	nu = dp[2];
585	length = dp[3];
586	slde = dp[4];
587	sld = dp[5];
588	delrho = slde - sld;
589	bkg = dp[6];
590
591	for(i=0;i<nordi;i++) {
592	//setup inner integral over the ellipsoidal cross-section
593	summj=0;
594	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy
595	arg = rasqrt(1.0-zizi);
596	for(j=0;j<nordj;j++) {
597	//20 gauss points for the inner integral
598	zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy
599	yyy = Gauss20Wt[j] * EllipCylKernel(q,arg,nu,zij);
600	summj += yyy;
601	}
602	//now calculate the value of the inner integral
603	answer = (vbj-vaj)/2.0*summj;
604	//divide integral by Pi
605	answer /=Pi;
606
607	//now calculate outer integral
608	arg = qlengthzi/2.0;
609	if (arg == 0.0){
610	si = 1.0;
611	}else{
612	si = sin(arg) * sin(arg) / arg / arg;
613	}
614	yyy = Gauss76Wt[i] * answer * si;
615	summ += yyy;
616	}
617
618	answer = (vb-va)/2.0*summ;
619	// Multiply by contrast^2
620	answer = delrhodelrho;
621	//normalize by cylinder volume
622	//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
623	vell = Pira(nura)length;
624	answer *= vell;
625	//convert to [cm-1]
626	answer *= 1.0e8;
627	//Scale
628	answer *= scale;
629	// add in the background
630	answer += bkg;
631
632	return answer;
633	}
634
635	/* TriaxialEllipsoidX : calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
636
637	Uses 76 pt Gaussian quadrature for both integrals
638
639	Warning:
640	The call to WaveData() below returns a pointer to the middle
641	of an unlocked Macintosh handle. In the unlikely event that your
642	calculations could cause memory to move, you should copy the coefficient
643	values to local variables or an array before such operations.
644	*/
645	double
646	TriaxialEllipsoid(double dp[], double q)
647	{
648	int i,j;
649	double Pi;
650	double scale,aa,bb,cc,delrho,bkg; //local variables of coefficient wave
651	int nordi=76; //order of integration
652	int nordj=76;
653	double va,vb; //upper and lower integration limits
654	double summ,zi,yyy,answer; //running tally of integration
655	double summj,vaj,vbj,zij,slde,sld; //for the inner integration
656
657	Pi = 4.0*atan(1.0);
658	va = 0.0;
659	vb = 1.0; //orintational average, outer integral
660	vaj = 0.0;
661	vbj = 1.0; //endpoints of inner integral
662
663	summ = 0.0; //initialize intergral
664
665	scale = dp[0]; //make local copies in case memory moves
666	aa = dp[1];
667	bb = dp[2];
668	cc = dp[3];
669	slde = dp[4];
670	sld = dp[5];
671	delrho = slde - sld;
672	bkg = dp[6];
673	for(i=0;i<nordi;i++) {
674	//setup inner integral over the ellipsoidal cross-section
675	summj=0.0;
676	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy
677	for(j=0;j<nordj;j++) {
678	//20 gauss points for the inner integral
679	zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy
680	yyy = Gauss76Wt[j] * TriaxialKernel(q,aa,bb,cc,zi,zij);
681	summj += yyy;
682	}
683	//now calculate the value of the inner integral
684	answer = (vbj-vaj)/2.0*summj;
685
686	//now calculate outer integral
687	yyy = Gauss76Wt[i] * answer;
688	summ += yyy;
689	} //final scaling is done at the end of the function, after the NT_FP64 case
690
691	answer = (vb-va)/2.0*summ;
692	// Multiply by contrast^2
693	answer = delrhodelrho;
694	//normalize by ellipsoid volume
695	answer = 4.0Pi/3.0aabb*cc;
696	//convert to [cm-1]
697	answer *= 1.0e8;
698	//Scale
699	answer *= scale;
700	// add in the background
701	answer += bkg;
702
703	return answer;
704	}
705
706	/* ParallelepipedX : calculates the form factor of a Parallelepiped (a rectangular solid)
707	at the given x-value p->x
708
709	Uses 76 pt Gaussian quadrature for both integrals
710
711	Warning:
712	The call to WaveData() below returns a pointer to the middle
713	of an unlocked Macintosh handle. In the unlikely event that your
714	calculations could cause memory to move, you should copy the coefficient
715	values to local variables or an array before such operations.
716	*/
717	double
718	Parallelepiped(double dp[], double q)
719	{
720	int i,j;
721	double scale,aa,bb,cc,delrho,bkg; //local variables of coefficient wave
722	int nordi=76; //order of integration
723	int nordj=76;
724	double va,vb; //upper and lower integration limits
725	double summ,yyy,answer; //running tally of integration
726	double summj,vaj,vbj; //for the inner integration
727	double mu,mudum,arg,sigma,uu,vol,sldp,sld;
728
729
730	// Pi = 4.0*atan(1.0);
731	va = 0.0;
732	vb = 1.0; //orintational average, outer integral
733	vaj = 0.0;
734	vbj = 1.0; //endpoints of inner integral
735
736	summ = 0.0; //initialize intergral
737
738	scale = dp[0]; //make local copies in case memory moves
739	aa = dp[1];
740	bb = dp[2];
741	cc = dp[3];
742	sldp = dp[4];
743	sld = dp[5];
744	delrho = sldp - sld;
745	bkg = dp[6];
746
747	mu = q*bb;
748	vol = aabbcc;
749	// normalize all WRT bb
750	aa = aa/bb;
751	cc = cc/bb;
752
753	for(i=0;i<nordi;i++) {
754	//setup inner integral over the ellipsoidal cross-section
755	summj=0.0;
756	sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy
757
758	for(j=0;j<nordj;j++) {
759	//76 gauss points for the inner integral
760	uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy
761	mudum = musqrt(1.0-sigmasigma);
762	yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu);
763	summj += yyy;
764	}
765	//now calculate the value of the inner integral
766	answer = (vbj-vaj)/2.0*summj;
767
768	arg = muccsigma/2.0;
769	if ( arg == 0.0 ) {
770	answer *= 1.0;
771	} else {
772	answer = sin(arg)sin(arg)/arg/arg;
773	}
774
775	//now sum up the outer integral
776	yyy = Gauss76Wt[i] * answer;
777	summ += yyy;
778	} //final scaling is done at the end of the function, after the NT_FP64 case
779
780	answer = (vb-va)/2.0*summ;
781	// Multiply by contrast^2
782	answer = delrhodelrho;
783	//normalize by volume
784	answer *= vol;
785	//convert to [cm-1]
786	answer *= 1.0e8;
787	//Scale
788	answer *= scale;
789	// add in the background
790	answer += bkg;
791
792	return answer;
793	}
794
795	/* HollowCylinderX : calculates the form factor of a Hollow Cylinder
796	at the given x-value p->x
797
798	Uses 76 pt Gaussian quadrature for the single integral
799
800	Warning:
801	The call to WaveData() below returns a pointer to the middle
802	of an unlocked Macintosh handle. In the unlikely event that your
803	calculations could cause memory to move, you should copy the coefficient
804	values to local variables or an array before such operations.
805	*/
806	double
807	HollowCylinder(double dp[], double q)
808	{
809	int i;
810	double scale,rcore,rshell,length,delrho,bkg; //local variables of coefficient wave
811	int nord=76; //order of integration
812	double va,vb,zi; //upper and lower integration limits
813	double summ,answer,pi,sldc,sld; //running tally of integration
814
815	pi = 4.0*atan(1.0);
816	va = 0.0;
817	vb = 1.0; //limits of numerical integral
818
819	summ = 0.0; //initialize intergral
820
821	scale = dp[0]; //make local copies in case memory moves
822	rcore = dp[1];
823	rshell = dp[2];
824	length = dp[3];
825	sldc = dp[4];
826	sld = dp[5];
827	delrho = sldc - sld;
828	bkg = dp[6];
829
830	for(i=0;i<nord;i++) {
831	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
832	summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi);
833	}
834
835	answer = (vb-va)/2.0*summ;
836	// Multiply by contrast^2
837	answer = delrhodelrho;
838	//normalize by volume
839	answer = pi(rshellrshell-rcorercore)*length;
840	//convert to [cm-1]
841	answer *= 1.0e8;
842	//Scale
843	answer *= scale;
844	// add in the background
845	answer += bkg;
846
847	return answer;
848	}
849
850	/* EllipsoidFormX : calculates the form factor of an ellipsoid of revolution with semiaxes a:a:nua
851	at the given x-value p->x
852
853	Uses 76 pt Gaussian quadrature for the single integral
854
855	Warning:
856	The call to WaveData() below returns a pointer to the middle
857	of an unlocked Macintosh handle. In the unlikely event that your
858	calculations could cause memory to move, you should copy the coefficient
859	values to local variables or an array before such operations.
860	*/
861	double
862	EllipsoidForm(double dp[], double q)
863	{
864	int i;
865	double scale,a,nua,delrho,bkg; //local variables of coefficient wave
866	int nord=76; //order of integration
867	double va,vb,zi; //upper and lower integration limits
868	double summ,answer,pi,slde,sld; //running tally of integration
869
870	pi = 4.0*atan(1.0);
871	va = 0.0;
872	vb = 1.0; //limits of numerical integral
873
874	summ = 0.0; //initialize intergral
875
876	scale = dp[0]; //make local copies in case memory moves
877	nua = dp[1];
878	a = dp[2];
879	slde = dp[3];
880	sld = dp[4];
881	delrho = slde - sld;
882	bkg = dp[5];
883
884	for(i=0;i<nord;i++) {
885	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
886	summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi);
887	}
888
889	answer = (vb-va)/2.0*summ;
890	// Multiply by contrast^2
891	answer = delrhodelrho;
892	//normalize by volume
893	answer = 4.0pi/3.0aa*nua;
894	//convert to [cm-1]
895	answer *= 1.0e8;
896	//Scale
897	answer *= scale;
898	// add in the background
899	answer += bkg;
900
901	return answer;
902	}
903
904
905	/* Cyl_PolyRadiusX : calculates the form factor of a cylinder at the given x-value p->x
906	the cylinder has a polydisperse cross section
907
908	*/
909	double
910	Cyl_PolyRadius(double dp[], double q)
911	{
912	int i;
913	double scale,radius,length,pd,delrho,bkg; //local variables of coefficient wave
914	int nord=20; //order of integration
915	double uplim,lolim; //upper and lower integration limits
916	double summ,zi,yyy,answer,Vpoly; //running tally of integration
917	double range,zz,Pi,sldc,sld;
918
919	Pi = 4.0*atan(1.0);
920	range = 3.4;
921
922	summ = 0.0; //initialize intergral
923
924	scale = dp[0]; //make local copies in case memory moves
925	radius = dp[1];
926	length = dp[2];
927	pd = dp[3];
928	sldc = dp[4];
929	sld = dp[5];
930	delrho = sldc - sld;
931	bkg = dp[6];
932
933	zz = (1.0/pd)*(1.0/pd) - 1.0;
934
935	lolim = radius(1.0-rangepd); //set the upper/lower limits to cover the distribution
936	if(lolim<0.0) {
937	lolim = 0.0;
938	}
939	if(pd>0.3) {
940	range = 3.4 + (pd-0.3)*18.0;
941	}
942	uplim = radius(1.0+rangepd);
943
944	for(i=0;i<nord;i++) {
945	zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
946	yyy = Gauss20Wt[i] * Cyl_PolyRadKernel(q, radius, length, zz, delrho, zi);
947	summ += yyy;
948	}
949
950	answer = (uplim-lolim)/2.0*summ;
951	//normalize by average cylinder volume
952	//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
953	Vpoly=Piradiusradiuslength(zz+2.0)/(zz+1.0);
954	answer /= Vpoly;
955	//convert to [cm-1]
956	answer *= 1.0e8;
957	//Scale
958	answer *= scale;
959	// add in the background
960	answer += bkg;
961
962	return answer;
963	}
964
965	/* Cyl_PolyLengthX : calculates the form factor of a cylinder at the given x-value p->x
966	the cylinder has a polydisperse Length
967
968	*/
969	double
970	Cyl_PolyLength(double dp[], double q)
971	{
972	int i;
973	double scale,radius,length,pd,delrho,bkg; //local variables of coefficient wave
974	int nord=20; //order of integration
975	double uplim,lolim; //upper and lower integration limits
976	double summ,zi,yyy,answer,Vpoly; //running tally of integration
977	double range,zz,Pi,sldc,sld;
978
979
980	Pi = 4.0*atan(1.0);
981	range = 3.4;
982
983	summ = 0.0; //initialize intergral
984
985	scale = dp[0]; //make local copies in case memory moves
986	radius = dp[1];
987	length = dp[2];
988	pd = dp[3];
989	sldc = dp[4];
990	sld = dp[5];
991	delrho = sldc - sld;
992	bkg = dp[6];
993
994	zz = (1.0/pd)*(1.0/pd) - 1.0;
995
996	lolim = length(1.0-rangepd); //set the upper/lower limits to cover the distribution
997	if(lolim<0.0) {
998	lolim = 0.0;
999	}
1000	if(pd>0.3) {
1001	range = 3.4 + (pd-0.3)*18.0;
1002	}
1003	uplim = length(1.0+rangepd);
1004
1005	for(i=0;i<nord;i++) {
1006	zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1007	yyy = Gauss20Wt[i] * Cyl_PolyLenKernel(q, radius, length, zz, delrho, zi);
1008	summ += yyy;
1009	}
1010
1011	answer = (uplim-lolim)/2.0*summ;
1012	//normalize by average cylinder volume (first moment)
1013	//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
1014	Vpoly=Piradiusradius*length;
1015	answer /= Vpoly;
1016	//convert to [cm-1]
1017	answer *= 1.0e8;
1018	//Scale
1019	answer *= scale;
1020	// add in the background
1021	answer += bkg;
1022
1023	return answer;
1024	}
1025
1026	/* CoreShellCylinderX : calculates the form factor of a cylinder at the given x-value p->x
1027	the cylinder has a core-shell structure
1028
1029	*/
1030	double
1031	CoreShellCylinder(double dp[], double q)
1032	{
1033	int i;
1034	double scale,rcore,length,bkg; //local variables of coefficient wave
1035	double thick,rhoc,rhos,rhosolv;
1036	int nord=76; //order of integration
1037	double uplim,lolim,halfheight; //upper and lower integration limits
1038	double summ,zi,yyy,answer,Vcyl; //running tally of integration
1039	double Pi;
1040
1041	Pi = 4.0*atan(1.0);
1042
1043	lolim = 0.0;
1044	uplim = Pi/2.0;
1045
1046	summ = 0.0; //initialize intergral
1047
1048	scale = dp[0]; //make local copies in case memory moves
1049	rcore = dp[1];
1050	thick = dp[2];
1051	length = dp[3];
1052	rhoc = dp[4];
1053	rhos = dp[5];
1054	rhosolv = dp[6];
1055	bkg = dp[7];
1056
1057	halfheight = length/2.0;
1058
1059	for(i=0;i<nord;i++) {
1060	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1061	yyy = Gauss76Wt[i] * CoreShellCylKernel(q, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi);
1062	summ += yyy;
1063	}
1064
1065	answer = (uplim-lolim)/2.0*summ;
1066	// length is the total core length
1067	Vcyl=Pi(rcore+thick)(rcore+thick)(length+2.0thick);
1068	answer /= Vcyl;
1069	//convert to [cm-1]
1070	answer *= 1.0e8;
1071	//Scale
1072	answer *= scale;
1073	// add in the background
1074	answer += bkg;
1075
1076	return answer;
1077	}
1078
1079
1080	/* PolyCoShCylinderX : calculates the form factor of a core-shell cylinder at the given x-value p->x
1081	the cylinder has a polydisperse CORE radius
1082
1083	*/
1084	double
1085	PolyCoShCylinder(double dp[], double q)
1086	{
1087	int i;
1088	double scale,radius,length,sigma,bkg; //local variables of coefficient wave
1089	double rad,radthick,facthick,rhoc,rhos,rhosolv;
1090	int nord=20; //order of integration
1091	double uplim,lolim; //upper and lower integration limits
1092	double summ,yyy,answer,Vpoly; //running tally of integration
1093	double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr;
1094
1095	Pi = 4.0*atan(1.0);
1096
1097	summ = 0.0; //initialize intergral
1098	Rsqrsumm = 0.0;
1099
1100	scale = dp[0];
1101	radius = dp[1];
1102	sigma = dp[2]; //sigma is the standard mean deviation
1103	length = dp[3];
1104	radthick = dp[4];
1105	facthick= dp[5];
1106	rhoc = dp[6];
1107	rhos = dp[7];
1108	rhosolv = dp[8];
1109	bkg = dp[9];
1110
1111	lolim = exp(log(radius)-(4.*sigma));
1112	if (lolim<0.0) {
1113	lolim=0.0; //to avoid numerical error when va<0 (-ve r value)
1114	}
1115	uplim = exp(log(radius)+(4.*sigma));
1116
1117	for(i=0;i<nord;i++) {
1118	rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1119	AR=(1.0/(radsigmasqrt(2.0Pi)))exp(-(0.5((log(radius/rad))/sigma)((log(radius/rad))/sigma)));
1120	yyy = AR* Gauss20Wt[i] * CSCylIntegration(q,rad,radthick,facthick,rhoc,rhos,rhosolv,length);
1121	Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick); //SRK normalize to total dimensions
1122	summ += yyy;
1123	Rsqrsumm += Rsqryyy;
1124	}
1125
1126	answer = (uplim-lolim)/2.0*summ;
1127	Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
1128	//normalize by average cylinder volume
1129	Vpoly = PiRsqr(length+2*facthick);
1130	answer /= Vpoly;
1131	//convert to [cm-1]
1132	answer *= 1.0e8;
1133	//Scale
1134	answer *= scale;
1135	// add in the background
1136	answer += bkg;
1137
1138	return answer;
1139	}
1140
1141	/* OblateFormX : calculates the form factor of a core-shell Oblate ellipsoid at the given x-value p->x
1142	the ellipsoid has a core-shell structure
1143
1144	*/
1145	double
1146	OblateForm(double dp[], double q)
1147	{
1148	int i;
1149	double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
1150	int nord=76; //order of integration
1151	double uplim,lolim; //upper and lower integration limits
1152	double summ,zi,yyy,answer,oblatevol; //running tally of integration
1153	double Pi,sldc,slds,sld;
1154
1155	Pi = 4.0*atan(1.0);
1156
1157	lolim = 0.0;
1158	uplim = 1.0;
1159
1160	summ = 0.0; //initialize intergral
1161
1162
1163	scale = dp[0]; //make local copies in case memory moves
1164	crmaj = dp[1];
1165	crmin = dp[2];
1166	trmaj = dp[3];
1167	trmin = dp[4];
1168	sldc = dp[5];
1169	slds = dp[6];
1170	sld = dp[7];
1171	delpc = sldc - slds; //core - shell
1172	delps = slds - sld; //shell - solvent
1173	bkg = dp[8];
1174
1175	for(i=0;i<nord;i++) {
1176	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1177	yyy = Gauss76Wt[i] * gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
1178	summ += yyy;
1179	}
1180
1181	answer = (uplim-lolim)/2.0*summ;
1182	// normalize by particle volume
1183	oblatevol = 4.0Pi/3.0trmajtrmajtrmin;
1184	answer /= oblatevol;
1185
1186	//convert to [cm-1]
1187	answer *= 1.0e8;
1188	//Scale
1189	answer *= scale;
1190	// add in the background
1191	answer += bkg;
1192
1193	return answer;
1194	}
1195
1196	/* ProlateFormX : calculates the form factor of a core-shell Prolate ellipsoid at the given x-value p->x
1197	the ellipsoid has a core-shell structure
1198
1199	*/
1200	double
1201	ProlateForm(double dp[], double q)
1202	{
1203	int i;
1204	double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
1205	int nord=76; //order of integration
1206	double uplim,lolim; //upper and lower integration limits
1207	double summ,zi,yyy,answer,prolatevol; //running tally of integration
1208	double Pi,sldc,slds,sld;
1209
1210	Pi = 4.0*atan(1.0);
1211
1212	lolim = 0.0;
1213	uplim = 1.0;
1214
1215	summ = 0.0; //initialize intergral
1216
1217	scale = dp[0]; //make local copies in case memory moves
1218	crmaj = dp[1];
1219	crmin = dp[2];
1220	trmaj = dp[3];
1221	trmin = dp[4];
1222	sldc = dp[5];
1223	slds = dp[6];
1224	sld = dp[7];
1225	delpc = sldc - slds; //core - shell
1226	delps = slds - sld; //shell - sovent
1227	bkg = dp[8];
1228
1229	for(i=0;i<nord;i++) {
1230	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1231	yyy = Gauss76Wt[i] * gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
1232	summ += yyy;
1233	}
1234
1235	answer = (uplim-lolim)/2.0*summ;
1236	// normalize by particle volume
1237	prolatevol = 4.0Pi/3.0trmajtrmintrmin;
1238	answer /= prolatevol;
1239
1240	//convert to [cm-1]
1241	answer *= 1.0e8;
1242	//Scale
1243	answer *= scale;
1244	// add in the background
1245	answer += bkg;
1246
1247	return answer;
1248	}
1249
1250
1251	/* StackedDiscsX : calculates the form factor of a stacked "tactoid" of core shell disks
1252	like clay platelets that are not exfoliated
1253
1254	*/
1255	double
1256	StackedDiscs(double dp[], double q)
1257	{
1258	int i;
1259	double scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd; //local variables of coefficient wave
1260	double va,vb,vcyl,summ,yyy,zi,halfheight,d,answer;
1261	int nord=76; //order of integration
1262	double Pi;
1263
1264
1265	Pi = 4.0*atan(1.0);
1266
1267	va = 0.0;
1268	vb = Pi/2.0;
1269
1270	summ = 0.0; //initialize intergral
1271
1272	scale = dp[0];
1273	rcore = dp[1];
1274	length = dp[2];
1275	thick = dp[3];
1276	rhoc = dp[4];
1277	rhol = dp[5];
1278	rhosolv = dp[6];
1279	N = dp[7];
1280	gsd = dp[8];
1281	bkg = dp[9];
1282
1283	d=2.0*thick+length;
1284	halfheight = length/2.0;
1285
1286	for(i=0;i<nord;i++) {
1287	zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0;
1288	yyy = Gauss76Wt[i] * Stackdisc_kern(q, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N);
1289	summ += yyy;
1290	}
1291
1292	answer = (vb-va)/2.0*summ;
1293	// length is the total core length
1294	vcyl=Pircorercore(2.0thick+length)*N;
1295	answer /= vcyl;
1296	//Convert to [cm-1]
1297	answer *= 1.0e8;
1298	//Scale
1299	answer *= scale;
1300	// add in the background
1301	answer += bkg;
1302
1303	return answer;
1304	}
1305
1306
1307	/* LamellarFFX : calculates the form factor of a lamellar structure - no S(q) effects included
1308
1309	*/
1310	double
1311	LamellarFF(double dp[], double q)
1312	{
1313	double scale,del,sig,contr,bkg; //local variables of coefficient wave
1314	double inten, qval,Pq;
1315	double Pi,sldb,sld;
1316
1317
1318	Pi = 4.0*atan(1.0);
1319	scale = dp[0];
1320	del = dp[1];
1321	sig = dp[2]*del;
1322	sldb = dp[3];
1323	sld = dp[4];
1324	contr = sldb - sld;
1325	bkg = dp[5];
1326	qval=q;
1327
1328	Pq = 2.0contrcontr/qval/qval(1.0-cos(qvaldel)exp(-0.5qvalqvalsig*sig));
1329
1330	inten = 2.0PiscalePq/(qvalqval); //this is now dimensionless...
1331
1332	inten /= del; //normalize by the thickness (in A)
1333
1334	inten *= 1.0e8; // 1/A to 1/cm
1335
1336	return(inten+bkg);
1337	}
1338
1339	/* LamellarPSX : calculates the form factor of a lamellar structure - with S(q) effects included
1340	--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
1341	model is smeared just like any other function
1342	*/
1343	double
1344	LamellarPS(double dp[], double q)
1345	{
1346	double scale,dd,del,sig,contr,NN,Cp,bkg; //local variables of coefficient wave
1347	double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
1348	double Pi,Euler,dQDefault,fii,sldb,sld;
1349	int ii,NNint;
1350	// char buf[256];
1351
1352
1353	Euler = 0.5772156649; // Euler's constant
1354	// dQDefault = 0.0025; //[=] 1/A, q-resolution, default value
1355	dQDefault = 0.0;
1356	dQ = dQDefault;
1357
1358	Pi = 4.0*atan(1.0);
1359	qval = q;
1360
1361	scale = dp[0];
1362	dd = dp[1];
1363	del = dp[2];
1364	sig = dp[3]*del;
1365	sldb = dp[4];
1366	sld = dp[5];
1367	contr = sldb - sld;
1368	NN = trunc(dp[6]); //be sure that NN is an integer
1369	Cp = dp[7];
1370	bkg = dp[8];
1371
1372	Pq = 2.0contrcontr/qval/qval(1.0-cos(qvaldel)exp(-0.5qvalqvalsig*sig));
1373
1374	NNint = (int)NN; //cast to an integer for the loop
1375
1376	// sprintf(buf, "qval = %g\r", qval);
1377	// XOPNotice(buf);
1378
1379	ii=0;
1380	Sq = 0.0;
1381	for(ii=1;ii<(NNint-1);ii+=1) {
1382
1383	fii = (double)ii; //do I really need to do this?
1384
1385	temp = 0.0;
1386	alpha = Cp/4.0/Pi/Pi(log(Pifii) + Euler);
1387	t1 = 2.0dQdQdddd*alpha;
1388	t2 = 2.0qvalqvaldddd*alpha;
1389	t3 = dQdQddddfii*fii;
1390
1391	temp = 1.0-fii/NN;
1392	temp = cos(ddqval*fii/(1.0+t1));
1393	temp = exp(-1.0(t2 + t3)/(2.0*(1.0+t1)) );
1394	temp /= sqrt(1.0+t1);
1395
1396	Sq += temp;
1397	}
1398
1399	Sq *= 2.0;
1400	Sq += 1.0;
1401
1402	inten = 2.0PiscalePqSq/(ddqvalqval);
1403
1404	inten *= 1.0e8; // 1/A to 1/cm
1405
1406	return(inten+bkg);
1407	}
1408
1409
1410	/* LamellarPS_HGX : calculates the form factor of a lamellar structure - with S(q) effects included
1411	--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
1412	model is smeared just like any other function
1413	*/
1414	double
1415	LamellarPS_HG(double dp[], double q)
1416	{
1417	double scale,dd,delT,delH,SLD_T,SLD_H,SLD_S,NN,Cp,bkg; //local variables of coefficient wave
1418	double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ,drh,drt;
1419	double Pi,Euler,dQDefault,fii;
1420	int ii,NNint;
1421
1422
1423	Euler = 0.5772156649; // Euler's constant
1424	// dQDefault = 0.0025; //[=] 1/A, q-resolution, default value
1425	dQDefault = 0.0;
1426	dQ = dQDefault;
1427
1428	Pi = 4.0*atan(1.0);
1429	qval= q;
1430
1431	scale = dp[0];
1432	dd = dp[1];
1433	delT = dp[2];
1434	delH = dp[3];
1435	SLD_T = dp[4];
1436	SLD_H = dp[5];
1437	SLD_S = dp[6];
1438	NN = trunc(dp[7]); //be sure that NN is an integer
1439	Cp = dp[8];
1440	bkg = dp[9];
1441
1442
1443	drh = SLD_H - SLD_S;
1444	drt = SLD_T - SLD_S; //correction 13FEB06 by L.Porcar
1445
1446	Pq = drh(sin(qval(delH+delT))-sin(qvaldelT)) + drtsin(qval*delT);
1447	Pq *= Pq;
1448	Pq = 4.0/(qvalqval);
1449
1450	NNint = (int)NN; //cast to an integer for the loop
1451	ii=0;
1452	Sq = 0.0;
1453	for(ii=1;ii<(NNint-1);ii+=1) {
1454
1455	fii = (double)ii; //do I really need to do this?
1456
1457	temp = 0.0;
1458	alpha = Cp/4.0/Pi/Pi(log(Piii) + Euler);
1459	t1 = 2.0dQdQdddd*alpha;
1460	t2 = 2.0qvalqvaldddd*alpha;
1461	t3 = dQdQddddii*ii;
1462
1463	temp = 1.0-ii/NN;
1464	temp = cos(ddqval*ii/(1.0+t1));
1465	temp = exp(-1.0(t2 + t3)/(2.0*(1.0+t1)) );
1466	temp /= sqrt(1.0+t1);
1467
1468	Sq += temp;
1469	}
1470
1471	Sq *= 2.0;
1472	Sq += 1.0;
1473
1474	inten = 2.0PiscalePqSq/(ddqvalqval);
1475
1476	inten *= 1.0e8; // 1/A to 1/cm
1477
1478	return(inten+bkg);
1479	}
1480
1481	/* LamellarFF_HGX : calculates the form factor of a lamellar structure - no S(q) effects included
1482	but extra SLD for head groups is included
1483
1484	*/
1485	double
1486	LamellarFF_HG(double dp[], double q)
1487	{
1488	double scale,delT,delH,slds,sldh,sldt,bkg; //local variables of coefficient wave
1489	double inten, qval,Pq,drh,drt;
1490	double Pi;
1491
1492
1493	Pi = 4.0*atan(1.0);
1494	qval= q;
1495	scale = dp[0];
1496	delT = dp[1];
1497	delH = dp[2];
1498	sldt = dp[3];
1499	sldh = dp[4];
1500	slds = dp[5];
1501	bkg = dp[6];
1502
1503
1504	drh = sldh - slds;
1505	drt = sldt - slds; //correction 13FEB06 by L.Porcar
1506
1507	Pq = drh(sin(qval(delH+delT))-sin(qvaldelT)) + drtsin(qval*delT);
1508	Pq *= Pq;
1509	Pq = 4.0/(qvalqval);
1510
1511	inten = 2.0PiscalePq/(qvalqval); //dimensionless...
1512
1513	inten /= 2.0*(delT+delH); //normalize by the bilayer thickness
1514
1515	inten *= 1.0e8; // 1/A to 1/cm
1516
1517	return(inten+bkg);
1518	}
1519
1520	/* FlexExclVolCylX : calculates the form factor of a flexible cylinder with a circular cross section
1521	-- incorporates Wei-Ren Chen's fixes - 2006
1522
1523	*/
1524	double
1525	FlexExclVolCyl(double dp[], double q)
1526	{
1527	double scale,L,B,bkg,rad,qr,cont,sldc,slds;
1528	double Pi,flex,crossSect,answer;
1529
1530
1531	Pi = 4.0*atan(1.0);
1532
1533	scale = dp[0]; //make local copies in case memory moves
1534	L = dp[1];
1535	B = dp[2];
1536	rad = dp[3];
1537	sldc = dp[4];
1538	slds = dp[5];
1539	cont = sldc-slds;
1540	bkg = dp[6];
1541
1542
1543	qr = q*rad;
1544
1545	flex = Sk_WR(q,L,B);
1546
1547	crossSect = (2.0NR_BessJ1(qr)/qr)(2.0*NR_BessJ1(qr)/qr);
1548	flex *= crossSect;
1549	flex = PiradradL;
1550	flex = contcont;
1551	flex *= 1.0e8;
1552	answer = scale*flex + bkg;
1553
1554	return answer;
1555	}
1556
1557	/* FlexCyl_EllipX : calculates the form factor of a flexible cylinder with an elliptical cross section
1558	-- incorporates Wei-Ren Chen's fixes - 2006
1559
1560	*/
1561	double
1562	FlexCyl_Ellip(double dp[], double q)
1563	{
1564	double scale,L,B,bkg,rad,qr,cont,ellRatio,slds,sldc;
1565	double Pi,flex,crossSect,answer;
1566
1567
1568	Pi = 4.0*atan(1.0);
1569	scale = dp[0]; //make local copies in case memory moves
1570	L = dp[1];
1571	B = dp[2];
1572	rad = dp[3];
1573	ellRatio = dp[4];
1574	sldc = dp[5];
1575	slds = dp[6];
1576	bkg = dp[7];
1577
1578	cont = sldc - slds;
1579	qr = q*rad;
1580
1581	flex = Sk_WR(q,L,B);
1582
1583	crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio));
1584	flex *= crossSect;
1585	flex = PiradradellRatio*L;
1586	flex = contcont;
1587	flex *= 1.0e8;
1588	answer = scale*flex + bkg;
1589
1590	return answer;
1591	}
1592
1593	double
1594	EllipticalCross_fn(double qq, double a, double b)
1595	{
1596	double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
1597	int i,nord=76;
1598
1599	Pi = 4.0*atan(1.0);
1600	lolim=0.0;
1601	uplim=Pi/2.0;
1602	summ=0.0;
1603
1604	for(i=0;i<nord;i++) {
1605	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1606	arg = qqsqrt(aasin(zi)sin(zi)+bbcos(zi)*cos(zi));
1607	yyy = pow((2.0 * NR_BessJ1(arg) / arg),2);
1608	yyy *= Gauss76Wt[i];
1609	summ += yyy;
1610	}
1611	answer = (uplim-lolim)/2.0*summ;
1612	answer *= 2.0/Pi;
1613	return(answer);
1614
1615	}
1616	/* FlexCyl_PolyLenX : calculates the form factor of a flecible cylinder at the given x-value p->x
1617	the cylinder has a polydisperse Length
1618
1619	*/
1620	double
1621	FlexCyl_PolyLen(double dp[], double q)
1622	{
1623	int i;
1624	double scale,radius,length,pd,bkg,lb,delrho,sldc,slds; //local variables of coefficient wave
1625	int nord=20; //order of integration
1626	double uplim,lolim; //upper and lower integration limits
1627	double summ,zi,yyy,answer,Vpoly; //running tally of integration
1628	double range,zz,Pi;
1629
1630	Pi = 4.0*atan(1.0);
1631	range = 3.4;
1632
1633	summ = 0.0; //initialize intergral
1634	scale = dp[0]; //make local copies in case memory moves
1635	length = dp[1]; //radius
1636	pd = dp[2]; // average length
1637	lb = dp[3];
1638	radius = dp[4];
1639	sldc = dp[5];
1640	slds = dp[6];
1641	bkg = dp[7];
1642
1643	delrho = sldc - slds;
1644	zz = (1.0/pd)*(1.0/pd) - 1.0;
1645
1646	lolim = length(1.0-rangepd); //set the upper/lower limits to cover the distribution
1647	if(lolim<0.0) {
1648	lolim = 0.0;
1649	}
1650	if(pd>0.3) {
1651	range = 3.4 + (pd-0.3)*18.0;
1652	}
1653	uplim = length(1.0+rangepd);
1654
1655	for(i=0;i<nord;i++) {
1656	zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1657	yyy = Gauss20Wt[i] * FlePolyLen_kernel(q,radius,length,lb,zz,delrho,zi);
1658	summ += yyy;
1659	}
1660
1661	answer = (uplim-lolim)/2.0*summ;
1662	//normalize by average cylinder volume (first moment), using the average length
1663	Vpoly=Piradiusradius*length;
1664	answer /= Vpoly;
1665
1666	answer =delrhodelrho;
1667
1668	//convert to [cm-1]
1669	answer *= 1.0e8;
1670	//Scale
1671	answer *= scale;
1672	// add in the background
1673	answer += bkg;
1674
1675	return answer;
1676	}
1677
1678	/* FlexCyl_PolyLenX : calculates the form factor of a flexible cylinder at the given x-value p->x
1679	the cylinder has a polydisperse cross sectional radius
1680
1681	*/
1682	double
1683	FlexCyl_PolyRad(double dp[], double q)
1684	{
1685	int i;
1686	double scale,radius,length,pd,delrho,bkg,lb,sldc,slds; //local variables of coefficient wave
1687	int nord=76; //order of integration
1688	double uplim,lolim; //upper and lower integration limits
1689	double summ,zi,yyy,answer,Vpoly; //running tally of integration
1690	double range,zz,Pi;
1691
1692
1693	Pi = 4.0*atan(1.0);
1694	range = 3.4;
1695
1696	summ = 0.0; //initialize intergral
1697
1698	scale = dp[0]; //make local copies in case memory moves
1699	length = dp[1]; //radius
1700	lb = dp[2]; // average length
1701	radius = dp[3];
1702	pd = dp[4];
1703	sldc = dp[5];
1704	slds = dp[6];
1705	bkg = dp[7];
1706
1707	delrho = sldc-slds;
1708	zz = (1.0/pd)*(1.0/pd) - 1.0;
1709
1710	lolim = radius(1.0-rangepd); //set the upper/lower limits to cover the distribution
1711	if(lolim<0.0) {
1712	lolim = 0.0;
1713	}
1714	if(pd>0.3) {
1715	range = 3.4 + (pd-0.3)*18.0;
1716	}
1717	uplim = radius(1.0+rangepd);
1718
1719	for(i=0;i<nord;i++) {
1720	//zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1721	//yyy = Gauss20Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
1722	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1723	yyy = Gauss76Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
1724	summ += yyy;
1725	}
1726
1727	answer = (uplim-lolim)/2.0*summ;
1728	//normalize by average cylinder volume (second moment), using the average radius
1729	Vpoly = Piradiusradiuslength(zz+2.0)/(zz+1.0);
1730	answer /= Vpoly;
1731
1732	answer =delrhodelrho;
1733
1734	//convert to [cm-1]
1735	answer *= 1.0e8;
1736	//Scale
1737	answer *= scale;
1738	// add in the background
1739	answer += bkg;
1740
1741	return answer;
1742	}
1743
1744
1745
1746	///////////////
1747
1748	//
1749	// FUNCTION gfn2: CONTAINS F(Q,A,B,mu)**2 AS GIVEN
1750	// BY (53) AND (56,57) IN CHEN AND
1751	// KOTLARCHYK REFERENCE
1752	//
1753	// <PROLATE ELLIPSOIDS>
1754	//
1755	double
1756	gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
1757	{
1758	// local variables
1759	double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi;
1760
1761	Pi = 4.0*atan(1.0);
1762
1763	pi43=4.0/3.0*Pi;
1764	aa = crmaj;
1765	bb = crmin;
1766	u2 = (aaaaxxxx + bbbb(1.0-xxxx));
1767	ut2 = (trmajtrmajxxxx + trmintrmin(1.0-xxxx));
1768	uq = sqrt(u2)*qq;
1769	ut= sqrt(ut2)*qq;
1770	vc = pi43aabb*bb;
1771	vt = pi43trmajtrmin*trmin;
1772	if (uq == 0.0){
1773	siq = 1.0/3.0;
1774	}else{
1775	siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
1776	}
1777	if (ut == 0.0){
1778	sit = 1.0/3.0;
1779	}else{
1780	sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
1781	}
1782	gfnc = 3.0siqvc*delpc;
1783	gfnt = 3.0sitvt*delps;
1784	gfn = gfnc+gfnt;
1785	gfn2 = gfn*gfn;
1786
1787	return (gfn2);
1788	}
1789
1790	//
1791	// FUNCTION gfn4: CONTAINS F(Q,A,B,MU)**2 AS GIVEN
1792	// BY (53) & (58-59) IN CHEN AND
1793	// KOTLARCHYK REFERENCE
1794	//
1795	// <OBLATE ELLIPSOID>
1796	// function gfn4 for oblate ellipsoids
1797	double
1798	gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
1799	{
1800	// local variables
1801	double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
1802
1803	Pi = 4.0*atan(1.0);
1804	pi43=4.0/3.0*Pi;
1805	aa = crmaj;
1806	bb = crmin;
1807	u2 = (bbbbxxxx + aaaa(1.0-xxxx));
1808	ut2 = (trmintrminxxxx + trmajtrmaj(1.0-xxxx));
1809	uq = sqrt(u2)*qq;
1810	ut= sqrt(ut2)*qq;
1811	vc = pi43aaaa*bb;
1812	vt = pi43trmajtrmaj*trmin;
1813	if (uq == 0.0){
1814	siq = 1.0/3.0;
1815	}else{
1816	siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
1817	}
1818	if (ut == 0.0){
1819	sit = 1.0/3.0;
1820	}else{
1821	sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
1822	}
1823	gfnc = 3.0siqvc*delpc;
1824	gfnt = 3.0sitvt*delps;
1825	tgfn = gfnc+gfnt;
1826	gfn4 = tgfn*tgfn;
1827
1828	return (gfn4);
1829	}
1830
1831	double
1832	FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi)
1833	{
1834	double Pq,vcyl,dl;
1835	double Pi,qr;
1836
1837	Pi = 4.0*atan(1.0);
1838	qr = q*radius;
1839
1840	Pq = Sk_WR(q,zi,lb); //does not have cross section term
1841	if (qr !=0){
1842	Pq = (2.0NR_BessJ1(qr)/qr)(2.0NR_BessJ1(qr)/qr);
1843	}
1844	vcyl=Piradiusradius*zi;
1845	Pq = vcylvcyl;
1846
1847	dl = SchulzPoint_cpr(zi,length,zz);
1848	return (Pq*dl);
1849
1850	}
1851
1852	double
1853	FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi)
1854	{
1855	double Pq,vcyl,dr;
1856	double Pi,qr;
1857
1858	Pi = 4.0*atan(1.0);
1859	qr = q*zi;
1860
1861	Pq = Sk_WR(q,Lc,Lb); //does not have cross section term
1862	if (qr !=0){
1863	Pq = (2.0NR_BessJ1(qr)/qr)(2.0NR_BessJ1(qr)/qr);
1864	}
1865
1866	vcyl=Pizizi*Lc;
1867	Pq = vcylvcyl;
1868
1869	dr = SchulzPoint_cpr(zi,ravg,zz);
1870	return (Pq*dr);
1871
1872	}
1873
1874	double
1875	CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length)
1876	{
1877	double answer,halfheight,Pi;
1878	double lolim,uplim,summ,yyy,zi;
1879	int nord,i;
1880
1881	// set up the integration end points
1882	Pi = 4.0*atan(1.0);
1883	nord = 76;
1884	lolim = 0.0;
1885	uplim = Pi/2.0;
1886	halfheight = length/2.0;
1887
1888	summ = 0.0; // initialize integral
1889	i=0;
1890	for(i=0;i<nord;i++) {
1891	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1892	yyy = Gauss76Wt[i] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi);
1893	summ += yyy;
1894	}
1895
1896	// calculate value of integral to return
1897	answer = (uplim-lolim)/2.0*summ;
1898	return (answer);
1899	}
1900
1901	double
1902	CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum)
1903	{
1904	// qq is the q-value for the calculation (1/A)
1905	// radius is the core radius of the cylinder (A)
1906	// radthick and facthick are the radial and face layer thicknesses
1907	// rho(n) are the respective SLD's
1908	// length is the Half CORE-LENGTH of the cylinder
1909	// dum is the dummy variable for the integration (theta)
1910
1911	double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
1912	double Pi;
1913
1914	Pi = 4.0*atan(1.0);
1915
1916	dr1 = rhoc-rhos;
1917	dr2 = rhos-rhosolv;
1918	vol1 = Piradrad(2.0length);
1919	vol2 = Pi(rad+radthick)(rad+radthick)(2.0length+2.0*facthick);
1920
1921	besarg1 = qqradsin(dum);
1922	besarg2 = qq(rad+radthick)sin(dum);
1923	sinarg1 = qqlengthcos(dum);
1924	sinarg2 = qq(length+facthick)cos(dum);
1925	if (besarg1 == 0.0){
1926	be1 = 0.5;
1927	}else{
1928	be1 = NR_BessJ1(besarg1)/besarg1;
1929	}
1930	if (besarg2 == 0.0){
1931	be2 = 0.5;
1932	}else{
1933	be2 = NR_BessJ1(besarg2)/besarg2;
1934	}
1935	if (sinarg1 == 0.0){
1936	si1 = 1.0;
1937	}else{
1938	si1 = sin(sinarg1)/sinarg1;
1939	}
1940	if (sinarg2 == 0.0){
1941	si2 = 1.0;
1942	}else{
1943	si2 = sin(sinarg2)/sinarg2;
1944	}
1945
1946	t1 = 2.0vol1dr1si1be1;
1947	t2 = 2.0vol2dr2si2be2;
1948
1949	retval = ((t1+t2)(t1+t2))sin(dum);
1950	return (retval);
1951
1952	}
1953
1954
1955	double
1956	CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum)
1957	{
1958
1959	double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
1960	double Pi;
1961
1962	Pi = 4.0*atan(1.0);
1963
1964	dr1 = rhoc-rhos;
1965	dr2 = rhos-rhosolv;
1966	vol1 = Pircorercore(2.0length);
1967	vol2 = Pi(rcore+thick)(rcore+thick)(2.0length+2.0*thick);
1968
1969	besarg1 = qqrcoresin(dum);
1970	besarg2 = qq(rcore+thick)sin(dum);
1971	sinarg1 = qqlengthcos(dum);
1972	sinarg2 = qq(length+thick)cos(dum);
1973
1974	if (besarg1 == 0.0){
1975	be1 = 0.5;
1976	}else{
1977	be1 = NR_BessJ1(besarg1)/besarg1;
1978	}
1979	if (besarg2 == 0.0){
1980	be2 = 0.5;
1981	}else{
1982	be2 = NR_BessJ1(besarg2)/besarg2;
1983	}
1984	if (sinarg1 == 0.0){
1985	si1 = 1.0;
1986	}else{
1987	si1 = sin(sinarg1)/sinarg1;
1988	}
1989	if (sinarg2 == 0.0){
1990	si2 = 1.0;
1991	}else{
1992	si2 = sin(sinarg2)/sinarg2;
1993	}
1994
1995	t1 = 2.0vol1dr1si1be1;
1996	t2 = 2.0vol2dr2si2be2;
1997
1998	retval = ((t1+t2)(t1+t2))sin(dum);
1999
2000	return (retval);
2001	}
2002
2003	double
2004	Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen)
2005	{
2006
2007	double halfheight,uplim,lolim,zi,summ,yyy,Pi;
2008	double answer,dr,Vcyl;
2009	int i,nord;
2010
2011	Pi = 4.0*atan(1.0);
2012	lolim = 0.0;
2013	uplim = Pi/2.0;
2014	halfheight = dumLen/2.0;
2015	nord=20;
2016	summ = 0.0;
2017
2018	//do the cylinder orientational average
2019	for(i=0;i<nord;i++) {
2020	zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2021	yyy = Gauss20Wt[i] * CylKernel(q, radius, halfheight, zi);
2022	summ += yyy;
2023	}
2024	answer = (uplim-lolim)/2.0*summ;
2025	// Multiply by contrast^2
2026	answer = delrhodelrho;
2027	// don't do the normal scaling to volume here
2028	// instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
2029	Vcyl = Piradiusradius*dumLen;
2030	answer = VcylVcyl;
2031
2032	dr = SchulzPoint_cpr(dumLen,len_avg,zz);
2033	return(dr*answer);
2034	}
2035
2036
2037	double
2038	Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N)
2039	{
2040	// qq is the q-value for the calculation (1/A)
2041	// rcore is the core radius of the cylinder (A)
2042	// rho(n) are the respective SLD's
2043	// length is the Half CORE-LENGTH of the cylinder = L (A)
2044	// dum is the dummy variable for the integration (x in Feigin's notation)
2045
2046	//Local variables
2047	double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
2048	double Pi;
2049	int kk;
2050
2051	Pi = 4.0*atan(1.0);
2052
2053	dr1 = rhoc-rhosolv;
2054	dr2 = rhol-rhosolv;
2055	area = Pircorercore;
2056	totald=2.0*(thick+length);
2057
2058	besarg1 = qqrcoresin(dum);
2059	besarg2 = qqrcoresin(dum);
2060
2061	sinarg1 = qqlengthcos(dum);
2062	sinarg2 = qq(length+thick)cos(dum);
2063
2064	if (besarg1 == 0.0){
2065	be1 = 0.5;
2066	}else{
2067	be1 = NR_BessJ1(besarg1)/besarg1;
2068	}
2069	if (besarg2 == 0.0){
2070	be2 = 0.5;
2071	}else{
2072	be2 = NR_BessJ1(besarg2)/besarg2;
2073	}
2074	if (sinarg1 == 0.0){
2075	si1 = 1.0;
2076	}else{
2077	si1 = sin(sinarg1)/sinarg1;
2078	}
2079	if (sinarg2 == 0.0){
2080	si2 = 1.0;
2081	}else{
2082	si2 = sin(sinarg2)/sinarg2;
2083	}
2084
2085	t1 = 2.0area(2.0length)dr1(si1)(be1);
2086	t2 = 2.0areadr2(totaldsi2-2.0lengthsi1)*(be2);
2087
2088	retval =((t1+t2)(t1+t2))sin(dum);
2089
2090	// loop for the structure facture S(q)
2091	sqq=0.0;
2092	for(kk=1;kk<N;kk+=1) {
2093	dexpt=qqcos(dum)qqcos(dum)ddgsdgsdkk/2.0;
2094	sqq=sqq+(N-kk)cos(qqcos(dum)dkk)exp(-1.dexpt);
2095	}
2096
2097	// end of loop for S(q)
2098	sqq=1.0+2.0*sqq/N;
2099	retval *= sqq;
2100
2101	return(retval);
2102	}
2103
2104
2105	double
2106	Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad)
2107	{
2108
2109	double halfheight,uplim,lolim,zi,summ,yyy,Pi;
2110	double answer,dr,Vcyl;
2111	int i,nord;
2112
2113	Pi = 4.0*atan(1.0);
2114	lolim = 0.0;
2115	uplim = Pi/2.0;
2116	halfheight = length/2.0;
2117	// nord=20;
2118	nord=76;
2119	summ = 0.0;
2120
2121	//do the cylinder orientational average
2122	// for(i=0;i<nord;i++) {
2123	// zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2124	// yyy = Gauss20Wt[i] * CylKernel(q, dumRad, halfheight, zi);
2125	// summ += yyy;
2126	// }
2127	for(i=0;i<nord;i++) {
2128	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2129	yyy = Gauss76Wt[i] * CylKernel(q, dumRad, halfheight, zi);
2130	summ += yyy;
2131	}
2132	answer = (uplim-lolim)/2.0*summ;
2133	// Multiply by contrast^2
2134	answer = delrhodelrho;
2135	// don't do the normal scaling to volume here
2136	// instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
2137	Vcyl = PidumRaddumRad*length;
2138	answer = VcylVcyl;
2139
2140	dr = SchulzPoint_cpr(dumRad,radius,zz);
2141	return(dr*answer);
2142	}
2143
2144	double
2145	SchulzPoint_cpr(double dumRad, double radius, double zz)
2146	{
2147	double dr;
2148
2149	dr = zzlog(dumRad) - gammaln(zz+1.0) + (zz+1.0)log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0));
2150	return(exp(dr));
2151	}
2152
2153
2154	double
2155	EllipsoidKernel(double qq, double a, double nua, double dum)
2156	{
2157	double arg,nu,retval; //local variables
2158
2159	nu = nua/a;
2160	arg = qqasqrt(1.0+dumdum(nu*nu-1.0));
2161	if (arg == 0.0){
2162	retval =1.0/3.0;
2163	}else{
2164	retval = (sin(arg)-argcos(arg))/(argarg*arg);
2165	}
2166	retval *= retval;
2167	retval *= 9.0;
2168
2169	return(retval);
2170	}//Function EllipsoidKernel()
2171
2172	double
2173	HolCylKernel(double qq, double rcore, double rshell, double length, double dum)
2174	{
2175	double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval; //local variables
2176
2177	gamma = rcore/rshell;
2178	arg1 = qqrshellsqrt(1.0-dum*dum); //1=shell (outer radius)
2179	arg2 = qqrcoresqrt(1.0-dum*dum); //2=core (inner radius)
2180	if (arg1 == 0.0){
2181	lam1 = 1.0;
2182	}else{
2183	lam1 = 2.0*NR_BessJ1(arg1)/arg1;
2184	}
2185	if (arg2 == 0.0){
2186	lam2 = 1.0;
2187	}else{
2188	lam2 = 2.0*NR_BessJ1(arg2)/arg2;
2189	}
2190	//Todo: Need to check psi behavior as gamma goes to 1.
2191	psi = 1.0/(1.0-gammagamma)(lam1 - gammagammalam2); //SRK 10/19/00
2192	sinarg = qqlengthdum/2.0;
2193	if (sinarg == 0.0){
2194	t2 = 1.0;
2195	}else{
2196	t2 = sin(sinarg)/sinarg;
2197	}
2198
2199	retval = psipsit2*t2;
2200
2201	return(retval);
2202	}//Function HolCylKernel()
2203
2204	double
2205	PPKernel(double aa, double mu, double uu)
2206	{
2207	// mu passed in is really mu*sqrt(1-sig^2)
2208	double arg1,arg2,Pi,tmp1,tmp2; //local variables
2209
2210	Pi = 4.0*atan(1.0);
2211
2212	//handle arg=0 separately, as sin(t)/t -> 1 as t->0
2213	arg1 = (mu/2.0)cos(Piuu/2.0);
2214	arg2 = (muaa/2.0)sin(Pi*uu/2.0);
2215	if(arg1==0.0) {
2216	tmp1 = 1.0;
2217	} else {
2218	tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
2219	}
2220
2221	if (arg2==0.0) {
2222	tmp2 = 1.0;
2223	} else {
2224	tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
2225	}
2226
2227	return (tmp1*tmp2);
2228
2229	}//Function PPKernel()
2230
2231
2232	double
2233	TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy)
2234	{
2235
2236	double arg,val,pi; //local variables
2237
2238	pi = 4.0*atan(1.0);
2239
2240	arg = aaaacos(pidx/2)cos(pi*dx/2);
2241	arg += bbbbsin(pidx/2)sin(pidx/2)(1-dy*dy);
2242	arg += ccccdy*dy;
2243	arg = q*sqrt(arg);
2244	if (arg == 0.0){
2245	val = 1.0; // as arg --> 0, val should go to 1.0
2246	}else{
2247	val = 9.0 * ( (sin(arg) - argcos(arg))/(argargarg) ) ( (sin(arg) - argcos(arg))/(argarg*arg) );
2248	}
2249	return (val);
2250
2251	}//Function TriaxialKernel()
2252
2253
2254	double
2255	CylKernel(double qq, double rr,double h, double theta)
2256	{
2257
2258	// qq is the q-value for the calculation (1/A)
2259	// rr is the radius of the cylinder (A)
2260	// h is the HALF-LENGTH of the cylinder = L/2 (A)
2261
2262	double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si; //Local variables
2263
2264
2265	besarg = qqrrsin(theta);
2266	siarg = qq * h * cos(theta);
2267	bj =NR_BessJ1(besarg);
2268
2269	//* Computing 2nd power */
2270	d1 = sin(siarg);
2271	t1 = d1 * d1;
2272	//* Computing 2nd power */
2273	d1 = bj;
2274	t2 = d1 * d1 * 4.0 * sin(theta);
2275	//* Computing 2nd power */
2276	d1 = siarg;
2277	b1 = d1 * d1;
2278	//* Computing 2nd power */
2279	d1 = qq * rr * sin(theta);
2280	b2 = d1 * d1;
2281	if (besarg == 0.0){
2282	be = sin(theta);
2283	}else{
2284	be = t2 / b2;
2285	}
2286	if (siarg == 0.0){
2287	si = 1.0;
2288	}else{
2289	si = t1 / b1;
2290	}
2291	retval = be * si;
2292
2293	return (retval);
2294
2295	}//Function CylKernel()
2296
2297	double
2298	EllipCylKernel(double qq, double ra,double nu, double theta)
2299	{
2300	//this is the function LAMBDA1^2 in Feigin's notation
2301	// qq is the q-value for the calculation (1/A)
2302	// ra is the transformed radius"a" in Feigin's notation
2303	// nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
2304	// theta is the dummy variable of the integration
2305
2306	double retval,arg; //Local variables
2307
2308	arg = qqrasqrt((1.0+nunu)/2+(1.0-nunu)*cos(theta)/2);
2309	if (arg == 0.0){
2310	retval = 1.0;
2311	}else{
2312	retval = 2.0*NR_BessJ1(arg)/arg;
2313	}
2314
2315	//square it
2316	retval *= retval;
2317
2318	return(retval);
2319
2320	}//Function EllipCylKernel()
2321
2322	double NR_BessJ1(double x)
2323	{
2324	double ax,z;
2325	double xx,y,ans,ans1,ans2;
2326
2327	if ((ax=fabs(x)) < 8.0) {
2328	y=x*x;
2329	ans1=x(72362614232.0+y(-7895059235.0+y*(242396853.1
2330	+y(-2972611.439+y(15704.48260+y*(-30.16036606))))));
2331	ans2=144725228442.0+y(2300535178.0+y(18583304.74
2332	+y(99447.43394+y(376.9991397+y*1.0))));
2333	ans=ans1/ans2;
2334	} else {
2335	z=8.0/ax;
2336	y=z*z;
2337	xx=ax-2.356194491;
2338	ans1=1.0+y(0.183105e-2+y(-0.3516396496e-4
2339	+y(0.2457520174e-5+y(-0.240337019e-6))));
2340	ans2=0.04687499995+y*(-0.2002690873e-3
2341	+y(0.8449199096e-5+y(-0.88228987e-6
2342	+y*0.105787412e-6)));
2343	ans=sqrt(0.636619772/ax)(cos(xx)ans1-zsin(xx)ans2);
2344	if (x < 0.0) ans = -ans;
2345	}
2346
2347	return(ans);
2348	}
2349
2350	/* Lamellar_ParaCrystal - Pedersen's model
2351
2352	*/
2353	double
2354	Lamellar_ParaCrystal(double w[], double q)
2355	{
2356	// Input (fitting) variables are:
2357	//[0] scale factor
2358	//[1] thickness
2359	//[2] number of layers
2360	//[3] spacing between layers
2361	//[4] polydispersity of spacing
2362	//[5] SLD lamellar
2363	//[6] SLD solvent
2364	//[7] incoherent background
2365	// give them nice names
2366	double inten,qval,scale,th,nl,davg,pd,contr,bkg,xn;
2367	double xi,ww,Pbil,Znq,Snq,an,sldLayer,sldSolvent,pi;
2368	long n1,n2;
2369
2370	pi = 4.0*atan(1.0);
2371	scale = w[0];
2372	th = w[1];
2373	nl = w[2];
2374	davg = w[3];
2375	pd = w[4];
2376	sldLayer = w[5];
2377	sldSolvent = w[6];
2378	bkg = w[7];
2379
2380	contr = w[5] - w[6];
2381	qval = q;
2382
2383	//get the fractional part of nl, to determine the "mixing" of N's
2384
2385	n1 = trunc(nl); //rounds towards zero
2386	n2 = n1 + 1;
2387	xn = (double)n2 - nl; //fractional contribution of n1
2388
2389	ww = exp(-qvalqvalpdpddavg*davg/2.0);
2390
2391	//calculate the n1 contribution
2392	an = paraCryst_an(ww,qval,davg,n1);
2393	Snq = paraCryst_sn(ww,qval,davg,n1,an);
2394
2395	Znq = xn*Snq;
2396
2397	//calculate the n2 contribution
2398	an = paraCryst_an(ww,qval,davg,n2);
2399	Snq = paraCryst_sn(ww,qval,davg,n2,an);
2400
2401	Znq += (1.0-xn)*Snq;
2402
2403	//and the independent contribution
2404	Znq += (1.0-wwww)/(1.0+wwww-2.0wwcos(qval*davg));
2405
2406	//the limit when NL approaches infinity
2407	// Zq = (1-ww^2)/(1+ww^2-2wwcos(qval*davg))
2408
2409	xi = th/2.0; //use 1/2 the bilayer thickness
2410	Pbil = (sin(qvalxi)/(qvalxi))(sin(qvalxi)/(qval*xi));
2411
2412	inten = 2.0picontrcontrPbilZnq/(qvalqval);
2413	inten *= 1.0e8;
2414
2415	return(scale*inten+bkg);
2416	}
2417
2418	// functions for the lamellar paracrystal model
2419	double
2420	paraCryst_sn(double ww, double qval, double davg, long nl, double an) {
2421
2422	double Snq;
2423
2424	Snq = an/( (double)nlpow((1.0+wwww-2.0wwcos(qval*davg)),2) );
2425
2426	return(Snq);
2427	}
2428
2429
2430	double
2431	paraCryst_an(double ww, double qval, double davg, long nl) {
2432
2433	double an;
2434
2435	an = 4.0wwww - 2.0(wwwwww+ww)cos(qval*davg);
2436	an -= 4.0pow(ww,(nl+2))cos((double)nlqvaldavg);
2437	an += 2.0pow(ww,(nl+3))cos((double)(nl-1)qvaldavg);
2438	an += 2.0pow(ww,(nl+1))cos((double)(nl+1)qvaldavg);
2439
2440	return(an);
2441	}
2442
2443
2444	/* Spherocylinder :
2445
2446	Uses 76 pt Gaussian quadrature for both integrals
2447	*/
2448	double
2449	Spherocylinder(double w[], double x)
2450	{
2451	int i,j;
2452	double Pi;
2453	double scale,contr,bkg,sldc,slds;
2454	double len,rad,hDist,endRad;
2455	int nordi=76; //order of integration
2456	int nordj=76;
2457	double va,vb; //upper and lower integration limits
2458	double summ,zi,yyy,answer; //running tally of integration
2459	double summj,vaj,vbj,zij; //for the inner integration
2460	double SphCyl_tmp[7],arg1,arg2,inner,be;
2461
2462
2463	scale = w[0];
2464	rad = w[1];
2465	len = w[2];
2466	sldc = w[3];
2467	slds = w[4];
2468	bkg = w[5];
2469
2470	SphCyl_tmp[0] = w[0];
2471	SphCyl_tmp[1] = w[1];
2472	SphCyl_tmp[2] = w[2];
2473	SphCyl_tmp[3] = w[1]; //end radius is same as cylinder radius
2474	SphCyl_tmp[4] = w[3];
2475	SphCyl_tmp[5] = w[4];
2476	SphCyl_tmp[6] = w[5];
2477
2478	hDist = 0; //by definition for this model
2479	endRad = rad;
2480
2481	contr = sldc-slds;
2482
2483	Pi = 4.0*atan(1.0);
2484	va = 0.0;
2485	vb = Pi/2.0; //orintational average, outer integral
2486	vaj = -1.0*hDist/endRad;
2487	vbj = 1.0; //endpoints of inner integral
2488
2489	summ = 0.0; //initialize intergral
2490
2491	for(i=0;i<nordi;i++) {
2492	//setup inner integral over the ellipsoidal cross-section
2493	summj=0.0;
2494	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy
2495
2496	for(j=0;j<nordj;j++) {
2497	//20 gauss points for the inner integral
2498	zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy
2499	yyy = Gauss76Wt[j] * SphCyl_kernel(SphCyl_tmp,x,zij,zi);
2500	summj += yyy;
2501	}
2502	//now calculate the value of the inner integral
2503	inner = (vbj-vaj)/2.0*summj;
2504	inner = 4.0PiendRadendRad*endRad;
2505
2506	//now calculate outer integrand
2507	arg1 = xlen/2.0cos(zi);
2508	arg2 = xradsin(zi);
2509	yyy = inner;
2510
2511	if(arg2 == 0) {
2512	be = 0.5;
2513	} else {
2514	be = NR_BessJ1(arg2)/arg2;
2515	}
2516
2517	if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h
2518	yyy += Piradradlen2.0*be;
2519	} else {
2520	yyy += Piradradlensin(arg1)/arg12.0be;
2521	}
2522	yyy *= yyy;
2523	yyy = sin(zi); // = \|A(q)\|^2sin(theta)
2524	yyy *= Gauss76Wt[i];
2525	summ += yyy;
2526	} //final scaling is done at the end of the function, after the NT_FP64 case
2527
2528	answer = (vb-va)/2.0*summ;
2529
2530	answer /= Piradradlen + Pi4.0endRadendRad*endRad/3.0; //divide by volume
2531	answer *= 1.0e8; //convert to cm^-1
2532	answer = contrcontr;
2533	answer *= scale;
2534	answer += bkg;
2535
2536	return answer;
2537	}
2538
2539
2540	// inner integral of the sphereocylinder model, special case of hDist = 0
2541	//
2542	double
2543	SphCyl_kernel(double w[], double x, double tt, double theta) {
2544
2545	double val,arg1,arg2;
2546	double scale,bkg,sldc,slds;
2547	double len,rad,hDist,endRad,be;
2548	scale = w[0];
2549	rad = w[1];
2550	len = w[2];
2551	endRad = w[3];
2552	sldc = w[4];
2553	slds = w[5];
2554	bkg = w[6];
2555
2556	hDist = 0.0;
2557
2558	arg1 = xcos(theta)(endRad*tt+hDist+len/2.0);
2559	arg2 = xendRadsin(theta)sqrt(1.0-tttt);
2560
2561	if(arg2 == 0) {
2562	be = 0.5;
2563	} else {
2564	be = NR_BessJ1(arg2)/arg2;
2565	}
2566	val = cos(arg1)(1.0-tttt)*be;
2567
2568	return(val);
2569	}
2570
2571
2572	/* Convex Lens : special case where L ~ 0 and hDist < 0
2573
2574	Uses 76 pt Gaussian quadrature for both integrals
2575	*/
2576	double
2577	ConvexLens(double w[], double x)
2578	{
2579	int i,j;
2580	double Pi;
2581	double scale,contr,bkg,sldc,slds;
2582	double len,rad,hDist,endRad;
2583	int nordi=76; //order of integration
2584	int nordj=76;
2585	double va,vb; //upper and lower integration limits
2586	double summ,zi,yyy,answer; //running tally of integration
2587	double summj,vaj,vbj,zij; //for the inner integration
2588	double CLens_tmp[7],arg1,arg2,inner,hh,be;
2589
2590
2591	scale = w[0];
2592	rad = w[1];
2593	// len = w[2]
2594	endRad = w[2];
2595	sldc = w[3];
2596	slds = w[4];
2597	bkg = w[5];
2598
2599	len = 0.01;
2600
2601	CLens_tmp[0] = w[0];
2602	CLens_tmp[1] = w[1];
2603	CLens_tmp[2] = len; //length is some small number, essentially zero
2604	CLens_tmp[3] = w[2];
2605	CLens_tmp[4] = w[3];
2606	CLens_tmp[5] = w[4];
2607	CLens_tmp[6] = w[5];
2608
2609	hDist = -1.0sqrt(fabs(endRadendRad-rad*rad)); //by definition for this model
2610
2611	contr = sldc-slds;
2612
2613	Pi = 4.0*atan(1.0);
2614	va = 0.0;
2615	vb = Pi/2.0; //orintational average, outer integral
2616	vaj = -1.0*hDist/endRad;
2617	vbj = 1.0; //endpoints of inner integral
2618
2619	summ = 0.0; //initialize intergral
2620
2621	for(i=0;i<nordi;i++) {
2622	//setup inner integral over the ellipsoidal cross-section
2623	summj=0.0;
2624	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy
2625
2626	for(j=0;j<nordj;j++) {
2627	//20 gauss points for the inner integral
2628	zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy
2629	yyy = Gauss76Wt[j] * ConvLens_kernel(CLens_tmp,x,zij,zi);
2630	summj += yyy;
2631	}
2632	//now calculate the value of the inner integral
2633	inner = (vbj-vaj)/2.0*summj;
2634	inner = 4.0PiendRadendRad*endRad;
2635
2636	//now calculate outer integrand
2637	arg1 = xlen/2.0cos(zi);
2638	arg2 = xradsin(zi);
2639	yyy = inner;
2640
2641	if(arg2 == 0) {
2642	be = 0.5;
2643	} else {
2644	be = NR_BessJ1(arg2)/arg2;
2645	}
2646
2647	if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h
2648	yyy += Piradradlen2.0*be;
2649	} else {
2650	yyy += Piradradlensin(arg1)/arg12.0be;
2651	}
2652	yyy *= yyy;
2653	yyy = sin(zi); // = \|A(q)\|^2sin(theta)
2654	yyy *= Gauss76Wt[i];
2655	summ += yyy;
2656	} //final scaling is done at the end of the function, after the NT_FP64 case
2657
2658	answer = (vb-va)/2.0*summ;
2659
2660	hh = fabs(hDist); //need positive value for spherical cap volume
2661	answer /= 2.0(1.0/3.0Pi(endRad-hh)(endRad-hh)(2.0endRad+hh)); //divide by volume
2662	answer *= 1.0e8; //convert to cm^-1
2663	answer = contrcontr;
2664	answer *= scale;
2665	answer += bkg;
2666
2667	return answer;
2668	}
2669
2670	/* Capped Cylinder : special case where L is nonzero and hDist < 0
2671
2672	-- uses the same Kernel as the Convex Lens
2673
2674	Uses 76 pt Gaussian quadrature for both integrals
2675	*/
2676	double
2677	CappedCylinder(double w[], double x)
2678	{
2679	int i,j;
2680	double Pi;
2681	double scale,contr,bkg,sldc,slds;
2682	double len,rad,hDist,endRad;
2683	int nordi=76; //order of integration
2684	int nordj=76;
2685	double va,vb; //upper and lower integration limits
2686	double summ,zi,yyy,answer; //running tally of integration
2687	double summj,vaj,vbj,zij; //for the inner integration
2688	double arg1,arg2,inner,hh,be;
2689
2690
2691	scale = w[0];
2692	rad = w[1];
2693	len = w[2];
2694	endRad = w[3];
2695	sldc = w[4];
2696	slds = w[5];
2697	bkg = w[6];
2698
2699	hDist = -1.0sqrt(fabs(endRadendRad-rad*rad)); //by definition for this model
2700
2701	contr = sldc-slds;
2702
2703	Pi = 4.0*atan(1.0);
2704	va = 0.0;
2705	vb = Pi/2.0; //orintational average, outer integral
2706	vaj = -1.0*hDist/endRad;
2707	vbj = 1.0; //endpoints of inner integral
2708
2709	summ = 0.0; //initialize intergral
2710
2711	for(i=0;i<nordi;i++) {
2712	//setup inner integral over the ellipsoidal cross-section
2713	summj=0.0;
2714	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy
2715
2716	for(j=0;j<nordj;j++) {
2717	//20 gauss points for the inner integral
2718	zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy
2719	yyy = Gauss76Wt[j] * ConvLens_kernel(w,x,zij,zi); //uses the same kernel as ConvexLens, except here L != 0
2720	summj += yyy;
2721	}
2722	//now calculate the value of the inner integral
2723	inner = (vbj-vaj)/2.0*summj;
2724	inner = 4.0PiendRadendRad*endRad;
2725
2726	//now calculate outer integrand
2727	arg1 = xlen/2.0cos(zi);
2728	arg2 = xradsin(zi);
2729	yyy = inner;
2730
2731	if(arg2 == 0) {
2732	be = 0.5;
2733	} else {
2734	be = NR_BessJ1(arg2)/arg2;
2735	}
2736
2737	if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h
2738	yyy += Piradradlen2.0*be;
2739	} else {
2740	yyy += Piradradlensin(arg1)/arg12.0be;
2741	}
2742
2743
2744
2745	yyy *= yyy;
2746	yyy = sin(zi); // = \|A(q)\|^2sin(theta)
2747	yyy *= Gauss76Wt[i];
2748	summ += yyy;
2749	} //final scaling is done at the end of the function, after the NT_FP64 case
2750
2751	answer = (vb-va)/2.0*summ;
2752
2753	hh = fabs(hDist); //need positive value for spherical cap volume
2754	answer /= Piradradlen + 2.0(1.0/3.0Pi(endRad-hh)(endRad-hh)(2.0*endRad+hh)); //divide by volume
2755	answer *= 1.0e8; //convert to cm^-1
2756	answer = contrcontr;
2757	answer *= scale;
2758	answer += bkg;
2759
2760	return answer;
2761	}
2762
2763
2764
2765	// inner integral of the ConvexLens model, special case where L ~ 0 and hDist < 0
2766	//
2767	double
2768	ConvLens_kernel(double w[], double x, double tt, double theta) {
2769
2770	double val,arg1,arg2;
2771	double scale,bkg,sldc,slds;
2772	double len,rad,hDist,endRad,be;
2773	scale = w[0];
2774	rad = w[1];
2775	len = w[2];
2776	endRad = w[3];
2777	sldc = w[4];
2778	slds = w[5];
2779	bkg = w[6];
2780
2781	hDist = -1.0sqrt(fabs(endRadendRad-rad*rad));
2782
2783	arg1 = xcos(theta)(endRad*tt+hDist+len/2.0);
2784	arg2 = xendRadsin(theta)sqrt(1.0-tttt);
2785
2786	if(arg2 == 0) {
2787	be = 0.5;
2788	} else {
2789	be = NR_BessJ1(arg2)/arg2;
2790	}
2791
2792	val = cos(arg1)(1.0-tttt)*be;
2793
2794	return(val);
2795	}
2796
2797
2798	/* Dumbbell : special case where L ~ 0 and hDist > 0
2799
2800	Uses 76 pt Gaussian quadrature for both integrals
2801	*/
2802	double
2803	Dumbbell(double w[], double x)
2804	{
2805	int i,j;
2806	double Pi;
2807	double scale,contr,bkg,sldc,slds;
2808	double len,rad,hDist,endRad;
2809	int nordi=76; //order of integration
2810	int nordj=76;
2811	double va,vb; //upper and lower integration limits
2812	double summ,zi,yyy,answer; //running tally of integration
2813	double summj,vaj,vbj,zij; //for the inner integration
2814	double Dumb_tmp[7],arg1,arg2,inner,be;
2815
2816
2817	scale = w[0];
2818	rad = w[1];
2819	// len = w[2]
2820	endRad = w[2];
2821	sldc = w[3];
2822	slds = w[4];
2823	bkg = w[5];
2824
2825	len = 0.01;
2826
2827	Dumb_tmp[0] = w[0];
2828	Dumb_tmp[1] = w[1];
2829	Dumb_tmp[2] = len; //length is some small number, essentially zero
2830	Dumb_tmp[3] = w[2];
2831	Dumb_tmp[4] = w[3];
2832	Dumb_tmp[5] = w[4];
2833	Dumb_tmp[6] = w[5];
2834
2835	hDist = sqrt(fabs(endRadendRad-radrad)); //by definition for this model
2836
2837	contr = sldc-slds;
2838
2839	Pi = 4.0*atan(1.0);
2840	va = 0.0;
2841	vb = Pi/2.0; //orintational average, outer integral
2842	vaj = -1.0*hDist/endRad;
2843	vbj = 1.0; //endpoints of inner integral
2844
2845	summ = 0.0; //initialize intergral
2846
2847	for(i=0;i<nordi;i++) {
2848	//setup inner integral over the ellipsoidal cross-section
2849	summj=0.0;
2850	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy
2851
2852	for(j=0;j<nordj;j++) {
2853	//20 gauss points for the inner integral
2854	zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy
2855	yyy = Gauss76Wt[j] * Dumb_kernel(Dumb_tmp,x,zij,zi);
2856	summj += yyy;
2857	}
2858	//now calculate the value of the inner integral
2859	inner = (vbj-vaj)/2.0*summj;
2860	inner = 4.0PiendRadendRad*endRad;
2861
2862	//now calculate outer integrand
2863	arg1 = xlen/2.0cos(zi);
2864	arg2 = xradsin(zi);
2865	yyy = inner;
2866
2867	if(arg2 == 0) {
2868	be = 0.5;
2869	} else {
2870	be = NR_BessJ1(arg2)/arg2;
2871	}
2872
2873	if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h
2874	yyy += Piradradlen2.0*be;
2875	} else {
2876	yyy += Piradradlensin(arg1)/arg12.0be;
2877	}
2878	yyy *= yyy;
2879	yyy = sin(zi); // = \|A(q)\|^2sin(theta)
2880	yyy *= Gauss76Wt[i];
2881	summ += yyy;
2882	} //final scaling is done at the end of the function, after the NT_FP64 case
2883
2884	answer = (vb-va)/2.0*summ;
2885
2886	answer /= 2.0Pi(2.0endRadendRadendRad/3.0+endRadendRadhDist-hDisthDist*hDist/3.0); //divide by volume
2887	answer *= 1.0e8; //convert to cm^-1
2888	answer = contrcontr;
2889	answer *= scale;
2890	answer += bkg;
2891
2892	return answer;
2893	}
2894
2895
2896	/* Barbell : "normal" case where L is nonzero 0 and hDist > 0
2897
2898	-- uses the same kernel as the Dumbbell case
2899
2900	Uses 76 pt Gaussian quadrature for both integrals
2901	*/
2902	double
2903	Barbell(double w[], double x)
2904	{
2905	int i,j;
2906	double Pi;
2907	double scale,contr,bkg,sldc,slds;
2908	double len,rad,hDist,endRad;
2909	int nordi=76; //order of integration
2910	int nordj=76;
2911	double va,vb; //upper and lower integration limits
2912	double summ,zi,yyy,answer; //running tally of integration
2913	double summj,vaj,vbj,zij; //for the inner integration
2914	double arg1,arg2,inner,be;
2915
2916
2917	scale = w[0];
2918	rad = w[1];
2919	len = w[2];
2920	endRad = w[3];
2921	sldc = w[4];
2922	slds = w[5];
2923	bkg = w[6];
2924
2925	hDist = sqrt(fabs(endRadendRad-radrad)); //by definition for this model
2926
2927	contr = sldc-slds;
2928
2929	Pi = 4.0*atan(1.0);
2930	va = 0.0;
2931	vb = Pi/2.0; //orintational average, outer integral
2932	vaj = -1.0*hDist/endRad;
2933	vbj = 1.0; //endpoints of inner integral
2934
2935	summ = 0.0; //initialize intergral
2936
2937	for(i=0;i<nordi;i++) {
2938	//setup inner integral over the ellipsoidal cross-section
2939	summj=0.0;
2940	zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy
2941
2942	for(j=0;j<nordj;j++) {
2943	//20 gauss points for the inner integral
2944	zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy
2945	yyy = Gauss76Wt[j] * Dumb_kernel(w,x,zij,zi); //uses the same Kernel as the Dumbbell, here L>0
2946	summj += yyy;
2947	}
2948	//now calculate the value of the inner integral
2949	inner = (vbj-vaj)/2.0*summj;
2950	inner = 4.0PiendRadendRad*endRad;
2951
2952	//now calculate outer integrand
2953	arg1 = xlen/2.0cos(zi);
2954	arg2 = xradsin(zi);
2955	yyy = inner;
2956
2957	if(arg2 == 0) {
2958	be = 0.5;
2959	} else {
2960	be = NR_BessJ1(arg2)/arg2;
2961	}
2962
2963	if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h
2964	yyy += Piradradlen2.0*be;
2965	} else {
2966	yyy += Piradradlensin(arg1)/arg12.0be;
2967	}
2968	yyy *= yyy;
2969	yyy = sin(zi); // = \|A(q)\|^2sin(theta)
2970	yyy *= Gauss76Wt[i];
2971	summ += yyy;
2972	} //final scaling is done at the end of the function, after the NT_FP64 case
2973
2974	answer = (vb-va)/2.0*summ;
2975
2976	answer /= Piradradlen + 2.0Pi(2.0endRadendRadendRad/3.0+endRadendRadhDist-hDisthDisthDist/3.0); //divide by volume
2977	answer *= 1.0e8; //convert to cm^-1
2978	answer = contrcontr;
2979	answer *= scale;
2980	answer += bkg;
2981
2982	return answer;
2983	}
2984
2985
2986
2987	// inner integral of the Dumbbell model, special case where L ~ 0 and hDist > 0
2988	//
2989	// inner integral of the Barbell model if L is nonzero
2990	//
2991	double
2992	Dumb_kernel(double w[], double x, double tt, double theta) {
2993
2994	double val,arg1,arg2;
2995	double scale,bkg,sldc,slds;
2996	double len,rad,hDist,endRad,be;
2997	scale = w[0];
2998	rad = w[1];
2999	len = w[2];
3000	endRad = w[3];
3001	sldc = w[4];
3002	slds = w[5];
3003	bkg = w[6];
3004
3005	hDist = sqrt(fabs(endRadendRad-radrad));
3006
3007	arg1 = xcos(theta)(endRad*tt+hDist+len/2.0);
3008	arg2 = xendRadsin(theta)sqrt(1.0-tttt);
3009
3010	if(arg2 == 0) {
3011	be = 0.5;
3012	} else {
3013	be = NR_BessJ1(arg2)/arg2;
3014	}
3015	val = cos(arg1)(1.0-tttt)*be;
3016
3017	return(val);
3018	}
3019
3020	double PolyCoreBicelle(double dp[], double q)
3021	{
3022	int i;
3023	int nord = 20;
3024	double scale, length, sigma, bkg, radius, radthick, facthick;
3025	double rhoc, rhoh, rhor, rhosolv;
3026	double answer, Vpoly;
3027	double Pi,lolim,uplim,summ,yyy,rad,AR,Rsqr,Rsqrsumm,Rsqryyy;
3028
3029	scale = dp[0];
3030	radius = dp[1];
3031	sigma = dp[2]; //sigma is the standard mean deviation
3032	length = dp[3];
3033	radthick = dp[4];
3034	facthick= dp[5];
3035	rhoc = dp[6];
3036	rhoh = dp[7];
3037	rhor=dp[8];
3038	rhosolv = dp[9];
3039	bkg = dp[10];
3040
3041	Pi = 4.0*atan(1.0);
3042
3043	lolim = exp(log(radius)-(4.*sigma));
3044	if (lolim<0.0) {
3045	lolim=0.0; //to avoid numerical error when va<0 (-ve r value)
3046	}
3047	uplim = exp(log(radius)+(4.*sigma));
3048
3049	summ = 0.0;
3050	Rsqrsumm = 0.0;
3051
3052	for(i=0;i<nord;i++) {
3053	rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
3054	AR=(1.0/(radsigmasqrt(2.0Pi)))exp(-(0.5((log(radius/rad))/sigma)((log(radius/rad))/sigma)));
3055	yyy = AR* Gauss20Wt[i] * BicelleIntegration(q,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length);
3056	Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick); //SRK normalize to total dimensions
3057	summ += yyy;
3058	Rsqrsumm += Rsqryyy;
3059	}
3060
3061	answer = (uplim-lolim)/2.0*summ;
3062	Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
3063	//normalize by average cylinder volume
3064	Vpoly = PiRsqr(length+2*facthick);
3065	answer /= Vpoly;
3066	//convert to [cm-1]
3067	answer *= 1.0e8;
3068	//Scale
3069	answer *= scale;
3070	// add in the background
3071	answer += bkg;
3072
3073	return answer;
3074
3075	}
3076
3077	double
3078	BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length){
3079
3080	double answer,halfheight,Pi;
3081	double lolim,uplim,summ,yyy,zi;
3082	int nord,i;
3083
3084	// set up the integration end points
3085	Pi = 4.0*atan(1.0);
3086	nord = 76;
3087	lolim = 0.0;
3088	uplim = Pi/2;
3089	halfheight = length/2.0;
3090
3091	summ = 0.0; // initialize integral
3092	i=0;
3093	for(i=0;i<nord;i++) {
3094	zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
3095	yyy = Gauss76Wt[i] * BicelleKernel(qq, rad, radthick, facthick, rhoc, rhoh, rhor,rhosolv, halfheight, zi);
3096	summ += yyy;
3097	}
3098
3099	// calculate value of integral to return
3100	answer = (uplim-lolim)/2.0*summ;
3101	return(answer);
3102	}
3103
3104	double
3105	BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum)
3106	{
3107	double dr1,dr2,dr3;
3108	double besarg1,besarg2;
3109	double vol1,vol2,vol3;
3110	double sinarg1,sinarg2;
3111	double t1,t2,t3;
3112	double retval,si1,si2,be1,be2;
3113
3114	double Pi = 4.0*atan(1.0);
3115
3116	dr1 = rhoc-rhoh;
3117	dr2 = rhor-rhosolv;
3118	dr3= rhoh-rhor;
3119	vol1 = Piradrad(2.0length);
3120	vol2 = Pi(rad+radthick)(rad+radthick)(2.0length+2.0*facthick);
3121	vol3= Pi(rad)(rad)(2.0length+2.0*facthick);
3122	besarg1 = qqradsin(dum);
3123	besarg2 = qq(rad+radthick)sin(dum);
3124	sinarg1 = qqlengthcos(dum);
3125	sinarg2 = qq(length+facthick)cos(dum);
3126
3127	if(besarg1 == 0) {
3128	be1 = 0.5;
3129	} else {
3130	be1 = NR_BessJ1(besarg1)/besarg1;
3131	}
3132	if(besarg2 == 0) {
3133	be2 = 0.5;
3134	} else {
3135	be2 = NR_BessJ1(besarg2)/besarg2;
3136	}
3137	if(sinarg1 == 0) {
3138	si1 = 1.0;
3139	} else {
3140	si1 = sin(sinarg1)/sinarg1;
3141	}
3142	if(sinarg2 == 0) {
3143	si2 = 1.0;
3144	} else {
3145	si2 = sin(sinarg2)/sinarg2;
3146	}
3147	t1 = 2.0vol1dr1si1be1;
3148	t2 = 2.0vol2dr2si2be2;
3149	t3 = 2.0vol3dr3si2be1;
3150
3151	retval = ((t1+t2+t3)(t1+t2+t3))sin(dum);
3152	return(retval);
3153
3154	}
3155
3156
3157	double
3158	CSPPKernel(double dp[], double mu, double uu)
3159	{
3160	double aa,bb,cc, ta,tb,tc;
3161	double Vin,Vot,V1,V2;
3162	double rhoA,rhoB,rhoC, rhoP, rhosolv;
3163	double dr0, drA,drB, drC;
3164	double arg1,arg2,arg3,arg4,t1,t2, t3, t4;
3165	double Pi,retVal;
3166
3167	aa = dp[1];
3168	bb = dp[2];
3169	cc = dp[3];
3170	ta = dp[4];
3171	tb = dp[5];
3172	tc = dp[6];
3173	rhoA=dp[7];
3174	rhoB=dp[8];
3175	rhoC=dp[9];
3176	rhoP=dp[10];
3177	rhosolv=dp[11];
3178	dr0=rhoP-rhosolv;
3179	drA=rhoA-rhosolv;
3180	drB=rhoB-rhosolv;
3181	drC=rhoC-rhosolv;
3182	Vin=(aabbcc);
3183	Vot=(aabbcc+2.0tabbcc+2.0aatbcc+2.0aabb*tc);
3184	V1=(2.0tabbcc); // incorrect V1 (aabbcc+2tabbcc)
3185	V2=(2.0aatbcc); // incorrect V2(aabbcc+2aatbcc)
3186	aa = aa/bb;
3187	ta=(aa+2.0*ta)/bb;
3188	tb=(aa+2.0*tb)/bb;
3189
3190	Pi = 4.0*atan(1.0);
3191
3192	arg1 = (muaa/2.0)sin(Pi*uu/2.0);
3193	arg2 = (mu/2.0)cos(Piuu/2.0);
3194	arg3= (muta/2.0)sin(Pi*uu/2.0);
3195	arg4= (mutb/2.0)cos(Pi*uu/2.0);
3196
3197	if(arg1==0.0){
3198	t1 = 1.0;
3199	} else {
3200	t1 = (sin(arg1)/arg1); //defn for CSPP model sin(arg1)/arg1 test: (sin(arg1)/arg1)*(sin(arg1)/arg1)
3201	}
3202	if(arg2==0.0){
3203	t2 = 1.0;
3204	} else {
3205	t2 = (sin(arg2)/arg2); //defn for CSPP model sin(arg2)/arg2 test: (sin(arg2)/arg2)*(sin(arg2)/arg2)
3206	}
3207	if(arg3==0.0){
3208	t3 = 1.0;
3209	} else {
3210	t3 = sin(arg3)/arg3;
3211	}
3212	if(arg4==0.0){
3213	t4 = 1.0;
3214	} else {
3215	t4 = sin(arg4)/arg4;
3216	}
3217	retVal =( dr0t1t2Vin + drA(t3-t1)t2V1+ drBt1(t4-t2)V2 )( dr0t1t2Vin + drA(t3-t1)t2V1+ drBt1(t4-t2)*V2 ); // correct FF : square of sum of phase factors
3218	return(retVal);
3219
3220	}
3221
3222	/* CSParallelepiped : calculates the form factor of a Parallelepiped with a core-shell structure
3223	-- different SLDs can be used for the face and rim
3224
3225	Uses 76 pt Gaussian quadrature for both integrals
3226	*/
3227	double
3228	CSParallelepiped(double dp[], double q)
3229	{
3230	int i,j;
3231	double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg; //local variables of coefficient wave
3232	int nordi=76; //order of integration
3233	int nordj=76;
3234	double va,vb; //upper and lower integration limits
3235	double summ,yyy,answer; //running tally of integration
3236	double summj,vaj,vbj; //for the inner integration
3237	double mu,mudum,arg,sigma,uu,vol;
3238
3239
3240	// Pi = 4.0*atan(1.0);
3241	va = 0.0;
3242	vb = 1.0; //orintational average, outer integral
3243	vaj = 0.0;
3244	vbj = 1.0; //endpoints of inner integral
3245
3246	summ = 0.0; //initialize intergral
3247
3248	scale = dp[0];
3249	aa = dp[1];
3250	bb = dp[2];
3251	cc = dp[3];
3252	ta = dp[4];
3253	tb = dp[5];
3254	tc = dp[6]; // is 0 at the moment
3255	rhoA=dp[7]; //rim A SLD
3256	rhoB=dp[8]; //rim B SLD
3257	rhoC=dp[9]; //rim C SLD
3258	rhoP = dp[10]; //Parallelpiped core SLD
3259	rhosolv=dp[11]; // Solvent SLD
3260	bkg = dp[12];
3261
3262	mu = q*bb;
3263	vol = aabbcc+2.0tabbcc+2.0aatbcc+2.0aabb*tc; //calculate volume before rescaling
3264
3265	// do the rescaling here, not in the kernel
3266	// normalize all WRT bb
3267	aa = aa/bb;
3268	cc = cc/bb;
3269
3270	for(i=0;i<nordi;i++) {
3271	//setup inner integral over the ellipsoidal cross-section
3272	summj=0.0;
3273	sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy
3274
3275	for(j=0;j<nordj;j++) {
3276	//76 gauss points for the inner integral
3277	uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy
3278	mudum = musqrt(1.0-sigmasigma);
3279	yyy = Gauss76Wt[j] * CSPPKernel(dp,mudum,uu);
3280	summj += yyy;
3281	}
3282	//now calculate the value of the inner integral
3283	answer = (vbj-vaj)/2.0*summj;
3284
3285	//finish the outer integral cc already scaled
3286	arg = muccsigma/2.0;
3287	if ( arg == 0.0 ) {
3288	answer *= 1.0;
3289	} else {
3290	answer = sin(arg)sin(arg)/arg/arg;
3291	}
3292
3293	//now sum up the outer integral
3294	yyy = Gauss76Wt[i] * answer;
3295	summ += yyy;
3296	} //final scaling is done at the end of the function, after the NT_FP64 case
3297
3298	answer = (vb-va)/2.0*summ;
3299
3300	//normalize by volume
3301	answer /= vol;
3302	//convert to [cm-1]
3303	answer *= 1.0e8;
3304	//Scale
3305	answer *= scale;
3306	// add in the background
3307	answer += bkg;
3308
3309	return answer;
3310	}
3311

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/sansmodels/src/libigor/libCylinder.c @ cd2f75d

Download in other formats: