Context Navigation

source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/TwoPhase.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.
File size: 19.3 KB

Line
1	/* TwoPhaseFit.c
2
3	*/
4
5	#pragma XOP_SET_STRUCT_PACKING // All structures are 2-byte-aligned.
6
7	#include "XOPStandardHeaders.h" // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
8	#include "SANSAnalysis.h"
9	#include "TwoPhase.h"
10
11	// scattering from the Teubner-Strey model for microemulsions - hardly needs to be an XOP...
12	int
13	TeubnerStreyModelX(FitParamsPtr p)
14	{
15	DOUBLE *dp; // Pointer to double precision wave data.
16	float *fp; // Pointer to single precision wave data.
17	DOUBLE q;
18	DOUBLE inten,q2,q4; //my local names
19
20	if (p->waveHandle == NIL) {
21	SetNaN64(&p->result);
22	return NON_EXISTENT_WAVE;
23	}
24
25	q= p->x;
26	q2 = q*q;
27	q4 = q2*q2;
28
29	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
30	case NT_FP32:
31	fp= WaveData(p->waveHandle);
32	// scale = fp[0];
33	// radius = fp[1];
34	// delrho = fp[2];
35	// bkg = fp[3];
36
37	inten = 1.0/(fp[0]+fp[1]q2+fp[2]q4);
38	inten += fp[3];
39
40	break;
41	case NT_FP64:
42	dp= WaveData(p->waveHandle);
43
44	inten = 1.0/(dp[0]+dp[1]q2+dp[2]q4);
45	inten += dp[3];
46
47	break;
48	default: // We can't handle this wave data type.
49	SetNaN64(&p->result);
50	return REQUIRES_SP_OR_DP_WAVE;
51	}
52
53	p->result= (inten); //scale, and add in the background
54	return 0;
55	}
56
57	int
58	Power_Law_ModelX(FitParamsPtr p)
59	{
60	DOUBLE *dp; // Pointer to double precision wave data.
61	float *fp; // Pointer to single precision wave data.
62	DOUBLE qval;
63	DOUBLE inten,A,m,bgd; //my local names
64
65	if (p->waveHandle == NIL) {
66	SetNaN64(&p->result);
67	return NON_EXISTENT_WAVE;
68	}
69
70	qval= p->x;
71
72	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
73	case NT_FP32:
74	fp= WaveData(p->waveHandle);
75	A = fp[0];
76	m = fp[1];
77	bgd = fp[2];
78
79	break;
80	case NT_FP64:
81	dp= WaveData(p->waveHandle);
82	A = dp[0];
83	m = dp[1];
84	bgd = dp[2];
85
86	break;
87	default: // We can't handle this wave data type.
88	SetNaN64(&p->result);
89	return REQUIRES_SP_OR_DP_WAVE;
90	}
91
92	inten = A*pow(qval,-m) + bgd;
93	p->result= (inten); //scale, and add in the background
94	return 0;
95	}
96
97
98	int
99	Peak_Lorentz_ModelX(FitParamsPtr p)
100	{
101	DOUBLE *dp; // Pointer to double precision wave data.
102	float *fp; // Pointer to single precision wave data.
103	DOUBLE qval;
104	DOUBLE inten,I0, qpk, dq,bgd; //my local names
105
106	if (p->waveHandle == NIL) {
107	SetNaN64(&p->result);
108	return NON_EXISTENT_WAVE;
109	}
110
111	qval= p->x;
112
113	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
114	case NT_FP32:
115	fp= WaveData(p->waveHandle);
116	I0 = fp[0];
117	qpk = fp[1];
118	dq = fp[2];
119	bgd = fp[3];
120
121	break;
122	case NT_FP64:
123	dp= WaveData(p->waveHandle);
124	I0 = dp[0];
125	qpk = dp[1];
126	dq = dp[2];
127	bgd = dp[3];
128
129	break;
130	default: // We can't handle this wave data type.
131	SetNaN64(&p->result);
132	return REQUIRES_SP_OR_DP_WAVE;
133	}
134
135	inten = I0/(1.0 + pow( (qval-qpk)/dq,2) ) + bgd;
136
137	p->result= (inten); //scale, and add in the background
138	return 0;
139	}
140
141	int
142	Peak_Gauss_ModelX(FitParamsPtr p)
143	{
144	DOUBLE *dp; // Pointer to double precision wave data.
145	float *fp; // Pointer to single precision wave data.
146	DOUBLE qval;
147	DOUBLE inten,I0, qpk, dq,bgd; //my local names
148
149	if (p->waveHandle == NIL) {
150	SetNaN64(&p->result);
151	return NON_EXISTENT_WAVE;
152	}
153
154	qval= p->x;
155
156	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
157	case NT_FP32:
158	fp= WaveData(p->waveHandle);
159	I0 = fp[0];
160	qpk = fp[1];
161	dq = fp[2];
162	bgd = fp[3];
163
164	break;
165	case NT_FP64:
166	dp= WaveData(p->waveHandle);
167	I0 = dp[0];
168	qpk = dp[1];
169	dq = dp[2];
170	bgd = dp[3];
171
172	break;
173	default: // We can't handle this wave data type.
174	SetNaN64(&p->result);
175	return REQUIRES_SP_OR_DP_WAVE;
176	}
177
178	inten = I0exp(-0.5pow((qval-qpk)/dq,2))+ bgd;
179
180	p->result= (inten); //scale, and add in the background
181	return 0;
182	}
183
184	int
185	Lorentz_ModelX(FitParamsPtr p)
186	{
187	DOUBLE *dp; // Pointer to double precision wave data.
188	float *fp; // Pointer to single precision wave data.
189	DOUBLE qval;
190	DOUBLE inten,I0, L,bgd; //my local names
191
192	if (p->waveHandle == NIL) {
193	SetNaN64(&p->result);
194	return NON_EXISTENT_WAVE;
195	}
196
197	qval= p->x;
198
199	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
200	case NT_FP32:
201	fp= WaveData(p->waveHandle);
202	I0 = fp[0];
203	L = fp[1];
204	bgd = fp[2];
205
206	break;
207	case NT_FP64:
208	dp= WaveData(p->waveHandle);
209	I0 = dp[0];
210	L = dp[1];
211	bgd = dp[2];
212
213	break;
214	default: // We can't handle this wave data type.
215	SetNaN64(&p->result);
216	return REQUIRES_SP_OR_DP_WAVE;
217	}
218
219	inten = I0/(1.0 + (qvalL)(qval*L)) + bgd;
220
221	p->result= (inten); //scale, and add in the background
222	return 0;
223	}
224
225	int
226	FractalX(FitParamsPtr p)
227	{
228	DOUBLE *dp; // Pointer to double precision wave data.
229	float *fp; // Pointer to single precision wave data.
230	DOUBLE x,pi;
231	DOUBLE r0,Df,corr,phi,sldp,sldm,bkg;
232	DOUBLE pq,sq,ans;
233
234	if (p->waveHandle == NIL) {
235	SetNaN64(&p->result);
236	return NON_EXISTENT_WAVE;
237	}
238
239	pi = 4.0*atan(1.0);
240	x= p->x;
241
242	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
243	case NT_FP32:
244	fp= WaveData(p->waveHandle);
245	phi = fp[0]; // volume fraction of building block spheres...
246	r0 = fp[1]; // radius of building block
247	Df = fp[2]; // fractal dimension
248	corr = fp[3]; // correlation length of fractal-like aggregates
249	sldp = fp[4]; // SLD of building block
250	sldm = fp[5]; // SLD of matrix or solution
251	bkg = fp[6]; // flat background
252
253	break;
254	case NT_FP64:
255	dp= WaveData(p->waveHandle);
256	phi = dp[0]; // volume fraction of building block spheres...
257	r0 = dp[1]; // radius of building block
258	Df = dp[2]; // fractal dimension
259	corr = dp[3]; // correlation length of fractal-like aggregates
260	sldp = dp[4]; // SLD of building block
261	sldm = dp[5]; // SLD of matrix or solution
262	bkg = dp[6]; // flat background
263
264	break;
265	default: // We can't handle this wave data type.
266	SetNaN64(&p->result);
267	return REQUIRES_SP_OR_DP_WAVE;
268	}
269
270	//calculate P(q) for the spherical subunits, units cm-1 sr-1
271	pq = 1.0e8phi4.0/3.0pir0r0r0(sldp-sldm)(sldp-sldm)pow((3(sin(xr0) - xr0cos(xr0))/pow((x*r0),3)),2);
272
273	//calculate S(q)
274	sq = Dfexp(gammln(Df-1.0))sin((Df-1.0)atan(xcorr));
275	sq /= pow((xr0),Df) pow((1.0 + 1.0/(xcorr)/(xcorr)),((Df-1.0)/2.0));
276	sq += 1.0;
277	//combine and return
278	ans = pq*sq + bkg;
279
280	p->result= (ans); //scale, and add in the background
281	return 0;
282	}
283
284	int
285	DAB_ModelX(FitParamsPtr p)
286	{
287	DOUBLE *dp; // Pointer to double precision wave data.
288	float *fp; // Pointer to single precision wave data.
289	DOUBLE qval,inten;
290	DOUBLE Izero, range, incoh;
291
292	if (p->waveHandle == NIL) {
293	SetNaN64(&p->result);
294	return NON_EXISTENT_WAVE;
295	}
296
297	qval= p->x;
298
299
300	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
301	case NT_FP32:
302	fp= WaveData(p->waveHandle);
303	Izero = fp[0];
304	range = fp[1];
305	incoh = fp[2];
306
307	break;
308	case NT_FP64:
309	dp= WaveData(p->waveHandle);
310	Izero = dp[0];
311	range = dp[1];
312	incoh = dp[2];
313
314	break;
315	default: // We can't handle this wave data type.
316	SetNaN64(&p->result);
317	return REQUIRES_SP_OR_DP_WAVE;
318	}
319
320	inten = Izero/pow((1.0 + (qvalrange)(qval*range)),2) + incoh;
321
322	p->result= (inten); //scale, and add in the background
323	return 0;
324	}
325
326	// G. Beaucage's Unified Model (1-4 levels)
327	//
328	int
329	OneLevelX(FitParamsPtr p)
330	{
331	DOUBLE *dp; // Pointer to double precision wave data.
332	float *fp; // Pointer to single precision wave data.
333	DOUBLE x,ans,erf1;
334	DOUBLE G1,Rg1,B1,Pow1,bkg,scale;
335
336	if (p->waveHandle == NIL) {
337	SetNaN64(&p->result);
338	return NON_EXISTENT_WAVE;
339	}
340
341	x= p->x;
342
343	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
344	case NT_FP32:
345	fp= WaveData(p->waveHandle);
346	scale = fp[0];
347	G1 = fp[1];
348	Rg1 = fp[2];
349	B1 = fp[3];
350	Pow1 = fp[4];
351	bkg = fp[5];
352
353
354	break;
355	case NT_FP64:
356	dp= WaveData(p->waveHandle);
357	scale = dp[0];
358	G1 = dp[1];
359	Rg1 = dp[2];
360	B1 = dp[3];
361	Pow1 = dp[4];
362	bkg = dp[5];
363
364
365	break;
366	default: // We can't handle this wave data type.
367	SetNaN64(&p->result);
368	return REQUIRES_SP_OR_DP_WAVE;
369	}
370
371	erf1 = erf( (x*Rg1/sqrt(6.0)));
372
373	ans = G1exp(-xxRg1Rg1/3.0);
374	ans += B1pow((erf1erf1*erf1/x),Pow1);
375
376	ans *= scale;
377	ans += bkg;
378
379	p->result= ans; //scale, and add in the background
380	return 0;
381	}
382
383	// G. Beaucage's Unified Model (1-4 levels)
384	//
385	int
386	TwoLevelX(FitParamsPtr p)
387	{
388	DOUBLE *dp; // Pointer to double precision wave data.
389	float *fp; // Pointer to single precision wave data.
390	DOUBLE x;
391	DOUBLE ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
392	DOUBLE erf1,erf2,scale;
393
394	if (p->waveHandle == NIL) {
395	SetNaN64(&p->result);
396	return NON_EXISTENT_WAVE;
397	}
398
399	x= p->x;
400
401	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
402	case NT_FP32:
403	fp= WaveData(p->waveHandle);
404	scale = fp[0];
405	G1 = fp[1]; //equivalent to I(0)
406	Rg1 = fp[2];
407	B1 = fp[3];
408	Pow1 = fp[4];
409	G2 = fp[5];
410	Rg2 = fp[6];
411	B2 = fp[7];
412	Pow2 = fp[8];
413	bkg = fp[9];
414
415	break;
416	case NT_FP64:
417	dp= WaveData(p->waveHandle);
418	scale = dp[0];
419	G1 = dp[1]; //equivalent to I(0)
420	Rg1 = dp[2];
421	B1 = dp[3];
422	Pow1 = dp[4];
423	G2 = dp[5];
424	Rg2 = dp[6];
425	B2 = dp[7];
426	Pow2 = dp[8];
427	bkg = dp[9];
428
429	break;
430	default: // We can't handle this wave data type.
431	SetNaN64(&p->result);
432	return REQUIRES_SP_OR_DP_WAVE;
433	}
434
435	erf1 = erf( (x*Rg1/sqrt(6.0)) );
436	erf2 = erf( (x*Rg2/sqrt(6.0)) );
437	//Print erf1
438
439	ans = G1exp(-xxRg1Rg1/3.0);
440	ans += B1exp(-xxRg2Rg2/3.0)pow((erf1erf1*erf1/x),Pow1);
441	ans += G2exp(-xxRg2Rg2/3.0);
442	ans += B2pow((erf2erf2*erf2/x),Pow2);
443
444	ans *= scale;
445	ans += bkg;
446
447	p->result= ans; //scale, and add in the background
448	return 0;
449	}
450
451	// G. Beaucage's Unified Model (1-4 levels)
452	//
453	int
454	ThreeLevelX(FitParamsPtr p)
455	{
456	DOUBLE *dp; // Pointer to double precision wave data.
457	float *fp; // Pointer to single precision wave data.
458	DOUBLE x;
459	DOUBLE ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
460	DOUBLE G3,Rg3,B3,Pow3,erf3;
461	DOUBLE erf1,erf2,scale;
462
463	if (p->waveHandle == NIL) {
464	SetNaN64(&p->result);
465	return NON_EXISTENT_WAVE;
466	}
467
468	x= p->x;
469
470	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
471	case NT_FP32:
472	fp= WaveData(p->waveHandle);
473	scale = fp[0];
474	G1 = fp[1]; //equivalent to I(0)
475	Rg1 = fp[2];
476	B1 = fp[3];
477	Pow1 = fp[4];
478	G2 = fp[5];
479	Rg2 = fp[6];
480	B2 = fp[7];
481	Pow2 = fp[8];
482	G3 = fp[9];
483	Rg3 = fp[10];
484	B3 = fp[11];
485	Pow3 = fp[12];
486	bkg = fp[13];
487
488	break;
489	case NT_FP64:
490	dp= WaveData(p->waveHandle);
491	scale = dp[0];
492	G1 = dp[1]; //equivalent to I(0)
493	Rg1 = dp[2];
494	B1 = dp[3];
495	Pow1 = dp[4];
496	G2 = dp[5];
497	Rg2 = dp[6];
498	B2 = dp[7];
499	Pow2 = dp[8];
500	G3 = dp[9];
501	Rg3 = dp[10];
502	B3 = dp[11];
503	Pow3 = dp[12];
504	bkg = dp[13];
505
506	break;
507	default: // We can't handle this wave data type.
508	SetNaN64(&p->result);
509	return REQUIRES_SP_OR_DP_WAVE;
510	}
511
512	erf1 = erf( (x*Rg1/sqrt(6.0)) );
513	erf2 = erf( (x*Rg2/sqrt(6.0)) );
514	erf3 = erf( (x*Rg3/sqrt(6.0)) );
515	//Print erf1
516
517	ans = G1exp(-xxRg1Rg1/3.0) + B1exp(-xxRg2Rg2/3.0)pow((erf1erf1*erf1/x),Pow1);
518	ans += G2exp(-xxRg2Rg2/3.0) + B2exp(-xxRg3Rg3/3.0)pow((erf2erf2*erf2/x),Pow2);
519	ans += G3exp(-xxRg3Rg3/3.0) + B3pow((erf3erf3*erf3/x),Pow3);
520
521	ans *= scale;
522	ans += bkg;
523
524	p->result= ans; //scale, and add in the background
525	return 0;
526	}
527
528	// G. Beaucage's Unified Model (1-4 levels)
529	//
530	int
531	FourLevelX(FitParamsPtr p)
532	{
533	DOUBLE *dp; // Pointer to double precision wave data.
534	float *fp; // Pointer to single precision wave data.
535	DOUBLE x;
536	DOUBLE ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
537	DOUBLE G3,Rg3,B3,Pow3,erf3;
538	DOUBLE G4,Rg4,B4,Pow4,erf4;
539	DOUBLE erf1,erf2,scale;
540
541	if (p->waveHandle == NIL) {
542	SetNaN64(&p->result);
543	return NON_EXISTENT_WAVE;
544	}
545
546	x= p->x;
547
548	switch(WaveType(p->waveHandle)){ // We can handle single and double precision coefficient waves.
549	case NT_FP32:
550	fp= WaveData(p->waveHandle);
551	scale = fp[0];
552	G1 = fp[1]; //equivalent to I(0)
553	Rg1 = fp[2];
554	B1 = fp[3];
555	Pow1 = fp[4];
556	G2 = fp[5];
557	Rg2 = fp[6];
558	B2 = fp[7];
559	Pow2 = fp[8];
560	G3 = fp[9];
561	Rg3 = fp[10];
562	B3 = fp[11];
563	Pow3 = fp[12];
564	G4 = fp[13];
565	Rg4 = fp[14];
566	B4 = fp[15];
567	Pow4 = fp[16];
568	bkg = fp[17];
569
570	break;
571	case NT_FP64:
572	dp= WaveData(p->waveHandle);
573	scale = dp[0];
574	G1 = dp[1]; //equivalent to I(0)
575	Rg1 = dp[2];
576	B1 = dp[3];
577	Pow1 = dp[4];
578	G2 = dp[5];
579	Rg2 = dp[6];
580	B2 = dp[7];
581	Pow2 = dp[8];
582	G3 = dp[9];
583	Rg3 = dp[10];
584	B3 = dp[11];
585	Pow3 = dp[12];
586	G4 = dp[13];
587	Rg4 = dp[14];
588	B4 = dp[15];
589	Pow4 = dp[16];
590	bkg = dp[17];
591
592	break;
593	default: // We can't handle this wave data type.
594	SetNaN64(&p->result);
595	return REQUIRES_SP_OR_DP_WAVE;
596	}
597
598	erf1 = erf( (x*Rg1/sqrt(6.0)) );
599	erf2 = erf( (x*Rg2/sqrt(6.0)) );
600	erf3 = erf( (x*Rg3/sqrt(6.0)) );
601	erf4 = erf( (x*Rg4/sqrt(6.0)) );
602
603	ans = G1exp(-xxRg1Rg1/3.0) + B1exp(-xxRg2Rg2/3.0)pow((erf1erf1*erf1/x),Pow1);
604	ans += G2exp(-xxRg2Rg2/3.0) + B2exp(-xxRg3Rg3/3.0)pow((erf2erf2*erf2/x),Pow2);
605	ans += G3exp(-xxRg3Rg3/3.0) + B3exp(-xxRg4Rg4/3.0)pow((erf3erf3*erf3/x),Pow3);
606	ans += G4exp(-xxRg4Rg4/3.0) + B4pow((erf4erf4*erf4/x),Pow4);
607
608	ans *= scale;
609	ans += bkg;
610
611	p->result= ans; //scale, and add in the background
612	return 0;
613	}
614
615
616	static DOUBLE
617	gammln(double xx) {
618
619	double x,y,tmp,ser;
620	static double cof[6]={76.18009172947146,-86.50532032941677,
621	24.01409824083091,-1.231739572450155,
622	0.1208650973866179e-2,-0.5395239384953e-5};
623	int j;
624
625	y=x=xx;
626	tmp=x+5.5;
627	tmp -= (x+0.5)*log(tmp);
628	ser=1.000000000190015;
629	for (j=0;j<=5;j++) ser += cof[j]/++y;
630	return -tmp+log(2.5066282746310005*ser/x);
631	}
632
633
634
635	#pragma XOP_RESET_STRUCT_PACKING // All structures are 2-byte-aligned.
636
637	///////////end of XOP
638
639

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

Download in other formats: