Context Navigation

source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/Cylinder.c @ 94

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

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

Download in other formats: