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

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 db08737 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
RevLine 
[e096270]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>
26using namespace std;
27
28extern "C" {
29        #include "onion.h"
30}
31
32OnionModel :: 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 */
108double 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].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
253                                                                                                                *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
254                                                                                                                *weights_s9[r].weight*weights_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].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
260                                                                                                        *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
261                                                                                                        *weights_s9[r].weight*weights_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].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
266                                                                                                        *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
267                                                                                                        *weights_s9[r].weight*weights_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 */
293double OnionModel :: operator()(double qx, double qy) {
294        double q = sqrt(qx*qx + qy*qy);
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 */
305double 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 */
313double 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].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
411                                                                                                        *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
412                                                                                                        *weights_s9[r].weight*weights_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].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight
416                                                                                                        *weights_s4[m].weight*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight
417                                                                                                        *weights_s8[t].weight*weights_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.