onion.cpp @ 6d48919

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

Last change on this file since 6d48919 was e096270, checked in by Jae Cho <jhjcho@…>, 14 years ago
added a new model
Property mode set to `100644`
File size: 14.0 KB

Line
1	/**
2	This software was developed by the University of Tennessee as part of the
3	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
4	project funded by the US National Science Foundation.
5
6	If you use DANSE applications to do scientific research that leads to
7	publication, we ask that you acknowledge the use of the software with the
8	following sentence:
9
10	"This work benefited from DANSE software developed under NSF award DMR-0520547."
11
12	copyright 2008, University of Tennessee
13	*/
14
15	/**
16	* Scattering model classes
17	* The classes use the IGOR library found in
18	* sansmodels/src/libigor
19	*
20	*/
21
22	#include <math.h>
23	#include "models.hh"
24	#include "parameters.hh"
25	#include <stdio.h>
26	using namespace std;
27
28	extern "C" {
29	#include "onion.h"
30	}
31
32	OnionModel :: OnionModel() {
33	n_shells = Parameter(1.0);
34	scale = Parameter(1.0);
35	rad_core = Parameter(200.0);
36	sld_core = Parameter(1e-06);
37	sld_solv = Parameter(6.4e-06);
38	background = Parameter(0.0);
39
40
41	sld_out_shell1 = Parameter(1.0e-06);
42	sld_out_shell2 = Parameter(1.0e-06);
43	sld_out_shell3 = Parameter(1.0e-06);
44	sld_out_shell4 = Parameter(1.0e-06);
45	sld_out_shell5 = Parameter(1.0e-06);
46	sld_out_shell6 = Parameter(1.0e-06);
47	sld_out_shell7 = Parameter(1.0e-06);
48	sld_out_shell8 = Parameter(1.0e-06);
49	sld_out_shell9 = Parameter(1.0e-06);
50	sld_out_shell10 = Parameter(1.0e-06);
51
52
53	sld_in_shell1 = Parameter(2.3e-06);
54	sld_in_shell2 = Parameter(2.6e-06);
55	sld_in_shell3 = Parameter(2.9e-06);
56	sld_in_shell4 = Parameter(3.2e-06);
57	sld_in_shell5 = Parameter(3.5e-06);
58	sld_in_shell6 = Parameter(3.8e-06);
59	sld_in_shell7 = Parameter(4.1e-06);
60	sld_in_shell8 = Parameter(4.4e-06);
61	sld_in_shell9 = Parameter(4.7e-06);
62	sld_in_shell10 = Parameter(5.0e-06);
63
64
65	A_shell1 = Parameter(1.0);
66	A_shell2 = Parameter(1.0);
67	A_shell3 = Parameter(1.0);
68	A_shell4 = Parameter(1.0);
69	A_shell5 = Parameter(1.0);
70	A_shell6 = Parameter(1.0);
71	A_shell7 = Parameter(1.0);
72	A_shell8 = Parameter(1.0);
73	A_shell9 = Parameter(1.0);
74	A_shell10 = Parameter(1.0);
75
76
77	thick_shell1 = Parameter(50.0);
78	thick_shell2 = Parameter(50.0);
79	thick_shell3 = Parameter(50.0);
80	thick_shell4 = Parameter(50.0);
81	thick_shell5 = Parameter(50.0);
82	thick_shell6 = Parameter(50.0);
83	thick_shell7 = Parameter(50.0);
84	thick_shell8 = Parameter(50.0);
85	thick_shell9 = Parameter(50.0);
86	thick_shell10 = Parameter(50.0);
87
88
89	func_shell1 = Parameter(2);
90	func_shell2 = Parameter(2);
91	func_shell3 = Parameter(2);
92	func_shell4 = Parameter(2);
93	func_shell5 = Parameter(2);
94	func_shell6 = Parameter(2);
95	func_shell7 = Parameter(2);
96	func_shell8 = Parameter(2);
97	func_shell9 = Parameter(2);
98	func_shell10 = Parameter(2);
99
100	}
101
102	/**
103	* Function to evaluate 1D scattering function
104	* The NIST IGOR library is used for the actual calculation.
105	* @param q: q-value
106	* @return: function value
107	*/
108	double OnionModel :: operator()(double q) {
109	double dp[56];
110	// Fill parameter array for IGOR library
111	// Add the background after averaging
112	dp[0] = n_shells();
113	dp[1] = scale();
114	dp[2] = rad_core();
115	dp[3] = sld_core();
116	dp[4] = sld_solv();
117	dp[5] = 0.0;
118
119	dp[6] = sld_out_shell1();
120	dp[7] = sld_out_shell2();
121	dp[8] = sld_out_shell3();
122	dp[9] = sld_out_shell4();
123	dp[10] = sld_out_shell5();
124	dp[11] = sld_out_shell6();
125	dp[12] = sld_out_shell7();
126	dp[13] = sld_out_shell8();
127	dp[14] = sld_out_shell9();
128	dp[15] = sld_out_shell10();
129
130	dp[16] = sld_in_shell1();
131	dp[17] = sld_in_shell2();
132	dp[18] = sld_in_shell3();
133	dp[19] = sld_in_shell4();
134	dp[20] = sld_in_shell5();
135	dp[21] = sld_in_shell6();
136	dp[22] = sld_in_shell7();
137	dp[23] = sld_in_shell8();
138	dp[24] = sld_in_shell9();
139	dp[25] = sld_in_shell10();
140
141	dp[26] = A_shell1();
142	dp[27] = A_shell2();
143	dp[28] = A_shell3();
144	dp[29] = A_shell4();
145	dp[30] = A_shell5();
146	dp[31] = A_shell6();
147	dp[32] = A_shell7();
148	dp[33] = A_shell8();
149	dp[34] = A_shell9();
150	dp[35] = A_shell10();
151
152	dp[36] = thick_shell1();
153	dp[37] = thick_shell2();
154	dp[38] = thick_shell3();
155	dp[39] = thick_shell4();
156	dp[40] = thick_shell5();
157	dp[41] = thick_shell6();
158	dp[42] = thick_shell7();
159	dp[43] = thick_shell8();
160	dp[44] = thick_shell9();
161	dp[45] = thick_shell10();
162
163	dp[46] = func_shell1();
164	dp[47] = func_shell2();
165	dp[48] = func_shell3();
166	dp[49] = func_shell4();
167	dp[50] = func_shell5();
168	dp[51] = func_shell6();
169	dp[52] = func_shell7();
170	dp[53] = func_shell8();
171	dp[54] = func_shell9();
172	dp[55] = func_shell10();
173
174
175	// Get the dispersion points for the radius
176	vector<WeightPoint> weights_rad;
177	rad_core.get_weights(weights_rad);
178
179	// Get the dispersion points for the thick 1
180	vector<WeightPoint> weights_s1;
181	thick_shell1.get_weights(weights_s1);
182
183	// Get the dispersion points for the thick 2
184	vector<WeightPoint> weights_s2;
185	thick_shell2.get_weights(weights_s2);
186
187	// Get the dispersion points for the thick 3
188	vector<WeightPoint> weights_s3;
189	thick_shell3.get_weights(weights_s3);
190
191	// Get the dispersion points for the thick 4
192	vector<WeightPoint> weights_s4;
193	thick_shell4.get_weights(weights_s4);
194
195	// Get the dispersion points for the thick 5
196	vector<WeightPoint> weights_s5;
197	thick_shell5.get_weights(weights_s5);
198
199	// Get the dispersion points for the thick 6
200	vector<WeightPoint> weights_s6;
201	thick_shell6.get_weights(weights_s6);
202
203	// Get the dispersion points for the thick 7
204	vector<WeightPoint> weights_s7;
205	thick_shell7.get_weights(weights_s7);
206
207	// Get the dispersion points for the thick 8
208	vector<WeightPoint> weights_s8;
209	thick_shell8.get_weights(weights_s8);
210	// Get the dispersion points for the thick 9
211	vector<WeightPoint> weights_s9;
212	thick_shell9.get_weights(weights_s9);
213
214	// Get the dispersion points for the thick 10
215	vector<WeightPoint> weights_s10;
216	thick_shell10.get_weights(weights_s10);
217
218
219	// Perform the computation, with all weight points
220	double sum = 0.0;
221	double norm = 0.0;
222	double vol = 0.0;
223
224	// Loop over radius weight points
225	for(int i=0; i<weights_rad.size(); i++) {
226	dp[2] = weights_rad[i].value;
227	// Loop over radius weight points
228	for(int j=0; j<weights_s1.size(); j++) {
229	dp[36] = weights_s1[j].value;
230	// Loop over radius weight points
231	for(int k=0; k<weights_s2.size(); k++) {
232	dp[37] = weights_s2[k].value;
233	// Loop over radius weight points
234	for(int l=0; l<weights_s3.size(); l++) {
235	dp[38] = weights_s3[l].value;
236	// Loop over radius weight points
237	for(int m=0; m<weights_s4.size(); m++) {
238	dp[39] = weights_s4[m].value;
239	for(int n=0; n<weights_s5.size(); n++) {
240	dp[40] = weights_s5[n].value;
241	for(int o=0; o<weights_s6.size(); o++) {
242	dp[41] = weights_s6[o].value;
243	for(int p=0; p<weights_s7.size(); p++) {
244	dp[42] = weights_s7[p].value;
245	for(int t=0; t<weights_s8.size(); t++) {
246	dp[43] = weights_s8[t].value;
247	for(int r=0; r<weights_s9.size(); r++) {
248	dp[44] = weights_s9[r].value;
249	for(int s=0; s<weights_s10.size(); s++) {
250	dp[45] = weights_s10[s].value;
251	//Un-normalize Shells by volume
252	sum += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weightweights_s3[l].weightweights_s4[m].weight
253	weights_s5[n].weightweights_s6[o].weightweights_s7[p].weightweights_s8[t].weight
254	weights_s9[r].weightweights_s10[s].weight
255	* so_kernel(dp,q) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
256	+weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
257	+weights_s9[r].value+weights_s10[s].value),3.0);
258	//Find average volume
259	vol += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weightweights_s3[l].weightweights_s4[m].weight
260	weights_s5[n].weightweights_s6[o].weightweights_s7[p].weightweights_s8[t].weight
261	weights_s9[r].weightweights_s10[s].weight
262	* pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
263	+weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
264	+weights_s9[r].value+weights_s10[s].value),3.0);
265	norm += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weightweights_s3[l].weightweights_s4[m].weight
266	weights_s5[n].weightweights_s6[o].weightweights_s7[p].weightweights_s8[t].weight
267	weights_s9[r].weightweights_s10[s].weight;
268	}
269	}
270	}
271	}
272	}
273	}
274	}
275	}
276	}
277	}
278	}
279
280	if (vol != 0.0 && norm != 0.0) {
281	//Re-normalize by avg volume
282	sum = sum/(vol/norm);}
283
284	return sum/norm + background();
285	}
286
287	/**
288	* Function to evaluate 2D scattering function
289	* @param q_x: value of Q along x
290	* @param q_y: value of Q along y
291	* @return: function value
292	*/
293	double OnionModel :: operator()(double qx, double qy) {
294	double q = sqrt(qxqx + qyqy);
295	return (*this).operator()(q);
296	}
297
298	/**
299	* Function to evaluate 2D scattering function
300	* @param pars: parameters of the sphere
301	* @param q: q-value
302	* @param phi: angle phi
303	* @return: function value
304	*/
305	double OnionModel :: evaluate_rphi(double q, double phi) {
306	return (*this).operator()(q);
307	}
308
309	/**
310	* Function to calculate effective radius
311	* @return: effective radius value
312	*/
313	double OnionModel :: calculate_ER() {
314	OnionParameters dp;
315	dp.rad_core = rad_core();
316	dp.thick_shell1 = thick_shell1();
317	dp.thick_shell2 = thick_shell2();
318	dp.thick_shell3 = thick_shell3();
319	dp.thick_shell4 = thick_shell4();
320	dp.thick_shell5 = thick_shell5();
321	dp.thick_shell6 = thick_shell6();
322	dp.thick_shell7 = thick_shell7();
323	dp.thick_shell8 = thick_shell8();
324	dp.thick_shell9 = thick_shell9();
325	dp.thick_shell10 = thick_shell10();
326
327
328	double rad_out = 0.0;
329	// Perform the computation, with all weight points
330	double sum = 0.0;
331	double norm = 0.0;
332
333	// Get the dispersion points for the radius
334	vector<WeightPoint> weights_rad;
335	rad_core.get_weights(weights_rad);
336
337	// Get the dispersion points for the thick 1
338	vector<WeightPoint> weights_s1;
339	thick_shell1.get_weights(weights_s1);
340
341	// Get the dispersion points for the thick 2
342	vector<WeightPoint> weights_s2;
343	thick_shell2.get_weights(weights_s2);
344
345	// Get the dispersion points for the thick 3
346	vector<WeightPoint> weights_s3;
347	thick_shell3.get_weights(weights_s3);
348
349	// Get the dispersion points for the thick 4
350	vector<WeightPoint> weights_s4;
351	thick_shell4.get_weights(weights_s4);
352	// Get the dispersion points for the thick 5
353	vector<WeightPoint> weights_s5;
354	thick_shell5.get_weights(weights_s5);
355
356	// Get the dispersion points for the thick 6
357	vector<WeightPoint> weights_s6;
358	thick_shell6.get_weights(weights_s6);
359
360	// Get the dispersion points for the thick 7
361	vector<WeightPoint> weights_s7;
362	thick_shell7.get_weights(weights_s7);
363
364	// Get the dispersion points for the thick 8
365	vector<WeightPoint> weights_s8;
366	thick_shell8.get_weights(weights_s8);
367	// Get the dispersion points for the thick 9
368	vector<WeightPoint> weights_s9;
369	thick_shell9.get_weights(weights_s9);
370
371	// Get the dispersion points for the thick 10
372	vector<WeightPoint> weights_s10;
373	thick_shell10.get_weights(weights_s10);
374
375
376	// Loop over radius weight points
377	for(int i=0; i<weights_rad.size(); i++) {
378	dp.rad_core = weights_rad[i].value;
379	// Loop over radius weight points
380	for(int j=0; j<weights_s1.size(); j++) {
381	dp.thick_shell1 = weights_s1[j].value;
382	// Loop over radius weight points
383	for(int k=0; k<weights_s2.size(); k++) {
384	dp.thick_shell2 = weights_s2[k].value;
385	// Loop over radius weight points
386	for(int l=0; l<weights_s3.size(); l++) {
387	dp.thick_shell3 = weights_s3[l].value;
388	// Loop over radius weight points
389	for(int m=0; m<weights_s4.size(); m++) {
390	dp.thick_shell4 = weights_s4[m].value;
391	// Loop over radius weight points
392	for(int n=0; j<weights_s5.size(); n++) {
393	dp.thick_shell5 = weights_s5[n].value;
394	// Loop over radius weight points
395	for(int o=0; k<weights_s6.size(); o++) {
396	dp.thick_shell6 = weights_s6[o].value;
397	// Loop over radius weight points
398	for(int p=0; l<weights_s7.size(); p++) {
399	dp.thick_shell7 = weights_s7[p].value;
400	// Loop over radius weight points
401	for(int t=0; m<weights_s8.size(); t++) {
402	dp.thick_shell8 = weights_s8[t].value;
403	// Loop over radius weight points
404	for(int r=0; l<weights_s9.size(); r++) {
405	dp.thick_shell8 = weights_s9[r].value;
406	// Loop over radius weight points
407	for(int s=0; m<weights_s10.size(); s++) {
408	dp.thick_shell10 = weights_s10[s].value;
409	//Un-normalize FourShell by volume
410	sum += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weightweights_s3[l].weightweights_s4[m].weight
411	weights_s5[n].weightweights_s6[o].weightweights_s7[p].weightweights_s8[t].weight
412	weights_s9[r].weightweights_s10[s].weight
413	* (dp.rad_core+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4+dp.thick_shell5
414	+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10);
415	norm += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weight*weights_s3[l].weight
416	weights_s4[m].weightweights_s5[n].weightweights_s6[o].weightweights_s7[p].weight
417	weights_s8[t].weightweights_s9[r].weight*weights_s10[s].weight;
418	}
419	}
420	}
421	}
422	}
423	}
424	}
425	}
426	}
427	}
428	}
429
430	if (norm != 0){
431	//return the averaged value
432	rad_out = sum/norm;}
433	else{
434	//return normal value
435	rad_out = dp.rad_core+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4
436	+dp.thick_shell5+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10;}
437	return rad_out;
438	}

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

SasView

source: sasview/sansmodels/src/sans/models/c_models/onion.cpp @ 6d48919

Download in other formats: