source: sasview/sansmodels/src/sans/models/c_models/spheresld.cpp @ bdc25e2

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 bdc25e2 was 7a8faf8, checked in by Jae Cho <jhjcho@…>, 14 years ago

Added polydispersion in a couple of params in spheresldmodel

  • Property mode set to 100644
File size: 9.5 KB
RevLine 
[0164899a]1
2#include <math.h>
3#include "models.hh"
4#include "parameters.hh"
5#include <stdio.h>
6using namespace std;
7
8extern "C" {
9        #include "spheresld.h"
10}
11
12SphereSLDModel :: SphereSLDModel() {
13        n_shells = Parameter(1.0);
14        scale = Parameter(1.0);
[7a8faf8]15        thick_inter0 = Parameter(1.0, true);
16        thick_inter0.set_min(0.0);
[0164899a]17        func_inter0 = Parameter(0);
18        sld_core0 = Parameter(2.07e-06);
19        sld_solv = Parameter(1.0e-06);
20    background = Parameter(0.0);
21
22
23    sld_flat1 = Parameter(2.7e-06);
24    sld_flat2 = Parameter(3.5e-06);
25    sld_flat3 = Parameter(4.0e-06);
26    sld_flat4 = Parameter(3.5e-06);
27    sld_flat5 = Parameter(4.0e-06);
28    sld_flat6 = Parameter(3.5e-06);
29    sld_flat7 = Parameter(4.0e-06);
30    sld_flat8 = Parameter(3.5e-06);
31    sld_flat9 = Parameter(4.0e-06);
32    sld_flat10 = Parameter(3.5e-06);
33
34
35    thick_inter1 = Parameter(1.0);
36    thick_inter2 = Parameter(1.0);
37    thick_inter3 = Parameter(1.0);
38    thick_inter4 = Parameter(1.0);
39    thick_inter5 = Parameter(1.0);
40    thick_inter6 = Parameter(1.0);
41    thick_inter7 = Parameter(1.0);
42    thick_inter8 = Parameter(1.0);
43    thick_inter9 = Parameter(1.0);
44    thick_inter10 = Parameter(1.0);
45
46
47    thick_flat1 = Parameter(100.0);
48    thick_flat2 = Parameter(100.0);
49    thick_flat3 = Parameter(100.0);
50    thick_flat4 = Parameter(100.0);
51    thick_flat5 = Parameter(100.0);
52    thick_flat6 = Parameter(100.0);
53    thick_flat7 = Parameter(100.0);
54    thick_flat8 = Parameter(100.0);
55    thick_flat9 = Parameter(100.0);
56    thick_flat10 = Parameter(100.0);
57
58
59    func_inter1 = Parameter(0);
60    func_inter2 = Parameter(0);
61    func_inter3 = Parameter(0);
62    func_inter4 = Parameter(0);
63    func_inter5 = Parameter(0);
64    func_inter6 = Parameter(0);
65    func_inter7 = Parameter(0);
66    func_inter8 = Parameter(0);
67    func_inter9 = Parameter(0);
68    func_inter10 = Parameter(0);
69
70    nu_inter1 = Parameter(2.5);
71    nu_inter2 = Parameter(2.5);
72    nu_inter3 = Parameter(2.5);
73    nu_inter4 = Parameter(2.5);
74    nu_inter5 = Parameter(2.5);
75    nu_inter6 = Parameter(2.5);
76    nu_inter7 = Parameter(2.5);
77    nu_inter8 = Parameter(2.5);
78    nu_inter9 = Parameter(2.5);
79    nu_inter10 = Parameter(2.5);
80
81    npts_inter = Parameter(35.0);
82    nu_inter0 = Parameter(2.5);
[7a8faf8]83    rad_core0 = Parameter(60.0, true);
84    rad_core0.set_min(0.0);
[0164899a]85}
86
87/**
88 * Function to evaluate 1D SphereSLD function
89 * @param q: q-value
90 * @return: function value
91 */
92double SphereSLDModel :: operator()(double q) {
93        double dp[60];
94        // Fill parameter array for IGOR library
95        // Add the background after averaging
96        dp[0] = n_shells();
97        dp[1] = scale();
98        dp[2] = thick_inter0();
99        dp[3] = func_inter0();
100        dp[4] = sld_core0();
101        dp[5] = sld_solv();
[7a8faf8]102        dp[6] = 0.0;
[0164899a]103
104        dp[7] = sld_flat1();
105        dp[8] = sld_flat2();
106        dp[9] = sld_flat3();
107        dp[10] = sld_flat4();
108        dp[11] = sld_flat5();
109        dp[12] = sld_flat6();
110        dp[13] = sld_flat7();
111        dp[14] = sld_flat8();
112        dp[15] = sld_flat9();
113        dp[16] = sld_flat10();
114
115        dp[17] = thick_inter1();
116        dp[18] = thick_inter2();
117        dp[19] = thick_inter3();
118        dp[20] = thick_inter4();
119        dp[21] = thick_inter5();
120        dp[22] = thick_inter6();
121        dp[23] = thick_inter7();
122        dp[24] = thick_inter8();
123        dp[25] = thick_inter9();
124        dp[26] = thick_inter10();
125
126        dp[27] = thick_flat1();
127        dp[28] = thick_flat2();
128        dp[29] = thick_flat3();
129        dp[30] = thick_flat4();
130        dp[31] = thick_flat5();
131        dp[32] = thick_flat6();
132        dp[33] = thick_flat7();
133        dp[34] = thick_flat8();
134        dp[35] = thick_flat9();
135        dp[36] = thick_flat10();
136
137        dp[37] = func_inter1();
138        dp[38] = func_inter2();
139        dp[39] = func_inter3();
140        dp[40] = func_inter4();
141        dp[41] = func_inter5();
142        dp[42] = func_inter6();
143        dp[43] = func_inter7();
144        dp[44] = func_inter8();
145        dp[45] = func_inter9();
146        dp[46] = func_inter10();
147
148        dp[47] = nu_inter1();
149        dp[48] = nu_inter2();
150        dp[49] = nu_inter3();
151        dp[50] = nu_inter4();
152        dp[51] = nu_inter5();
153        dp[52] = nu_inter6();
154        dp[53] = nu_inter7();
155        dp[54] = nu_inter8();
156        dp[55] = nu_inter9();
157        dp[56] = nu_inter10();
158
159
160        dp[57] = npts_inter();
161        dp[58] = nu_inter0();
162        dp[59] = rad_core0();
163
164        // No polydispersion supported in this model.
[7a8faf8]165        // Get the dispersion points for the radius
166        vector<WeightPoint> weights_rad_core0;
167        rad_core0.get_weights(weights_rad_core0);
168        vector<WeightPoint> weights_thick_inter0;
169        thick_inter0.get_weights(weights_thick_inter0);
170        // Perform the computation, with all weight points
171        double sum = 0.0;
172        double norm = 0.0;
173        double vol = 0.0;
174
175        // Loop over core weight points
176        for(int i=0; i<weights_rad_core0.size(); i++) {
177                dp[59] = weights_rad_core0[i].value;
178                // Loop over thick_inter0 weight points
179                for(int j=0; j<weights_thick_inter0.size(); j++) {
180                        dp[2] = weights_thick_inter0[j].value;
181
182                        //Un-normalize Sphere by volume
183                        sum += weights_rad_core0[i].weight * weights_thick_inter0[j].weight
184                                * sphere_sld_kernel(dp,q) * pow((weights_rad_core0[i].value +
185                                                weights_thick_inter0[j].value),3.0);
186                        //Find average volume
187                        vol += weights_rad_core0[i].weight * weights_thick_inter0[j].weight
188                                * pow((weights_rad_core0[i].value+weights_thick_inter0[j].value),3.0);
189
190                        norm += weights_rad_core0[i].weight * weights_thick_inter0[j].weight;
191                }
192        }
193
194        if (vol != 0.0 && norm != 0.0) {
195                //Re-normalize by avg volume
196                sum = sum/(vol/norm);}
197
198        return sum/norm + background();
[0164899a]199}
200
201/**
202 * Function to evaluate 2D SphereSLD function
203 * @param q_x: value of Q along x
204 * @param q_y: value of Q along y
205 * @return: function value
206 */
207double SphereSLDModel :: operator()(double qx, double qy) {
208        double q = sqrt(qx*qx + qy*qy);
209        return (*this).operator()(q);
210}
211
212/**
213 * Function to evaluate SphereSLD function
214 * @param pars: parameters of the SphereSLD
215 * @param q: q-value
216 * @param phi: angle phi
217 * @return: function value
218 */
219double SphereSLDModel :: evaluate_rphi(double q, double phi) {
220        return (*this).operator()(q);
221}
222
223/**
224 * Function to calculate TOTAL radius
225 * ToDo: Find What is the effective radius for this model.
226 * @return: effective radius value
227 */
228// No polydispersion supported in this model.
229// Calculate max radius assumming max_radius = effective radius
230// Note that this max radius is not affected by sld of layer, sld of interface, or
231// sld of solvent.
232double SphereSLDModel :: calculate_ER() {
233        SphereSLDParameters dp;
234
235        dp.n_shells = n_shells();
236
237        dp.rad_core0 = rad_core0();
238        dp.thick_flat1 = thick_flat1();
239        dp.thick_flat2 = thick_flat2();
240        dp.thick_flat3 = thick_flat3();
241        dp.thick_flat4 = thick_flat4();
242        dp.thick_flat5 = thick_flat5();
243        dp.thick_flat6 = thick_flat6();
244        dp.thick_flat7 = thick_flat7();
245        dp.thick_flat8 = thick_flat8();
246        dp.thick_flat9 = thick_flat9();
247        dp.thick_flat10 = thick_flat10();
248
249        dp.thick_inter0 = thick_inter0();
250        dp.thick_inter1 = thick_inter1();
251        dp.thick_inter2 = thick_inter2();
252        dp.thick_inter3 = thick_inter3();
253        dp.thick_inter4 = thick_inter4();
254        dp.thick_inter5 = thick_inter5();
255        dp.thick_inter6 = thick_inter6();
256        dp.thick_inter7 = thick_inter7();
257        dp.thick_inter8 = thick_inter8();
258        dp.thick_inter9 = thick_inter9();
259        dp.thick_inter10 = thick_inter10();
260
[7a8faf8]261        double rad_out = 0.0;
262        double out = 0.0;
263        // Perform the computation, with all weight points
264        double sum = 0.0;
265        double norm = 0.0;
266
267        // Get the dispersion points for the radius
268        vector<WeightPoint> weights_rad_core0;
269        rad_core0.get_weights(weights_rad_core0);
270
271        // Get the dispersion points for the thick 1
272        vector<WeightPoint> weights_thick_inter0;
273        thick_inter0.get_weights(weights_thick_inter0);
274        // Loop over radius weight points
275        for(int i=0; i<weights_rad_core0.size(); i++) {
276                dp.rad_core0 = weights_rad_core0[i].value;
277                // Loop over radius weight points
278                for(int j=0; j<weights_thick_inter0.size(); j++) {
279                        dp.thick_inter0 = weights_thick_inter0[j].value;
280                        rad_out = dp.rad_core0 + dp.thick_inter0;
281                        if (dp.n_shells > 0)
282                                rad_out += dp.thick_flat1 + dp.thick_inter1;
283                        if (dp.n_shells > 1)
284                                rad_out += dp.thick_flat2 + dp.thick_inter2;
285                        if (dp.n_shells > 2)
286                                rad_out += dp.thick_flat3 + dp.thick_inter3;
287                        if (dp.n_shells > 3)
288                                rad_out += dp.thick_flat4 + dp.thick_inter4;
289                        if (dp.n_shells > 4)
290                                rad_out += dp.thick_flat5 + dp.thick_inter5;
291                        if (dp.n_shells > 5)
292                                rad_out += dp.thick_flat6 + dp.thick_inter6;
293                        if (dp.n_shells > 6)
294                                rad_out += dp.thick_flat7 + dp.thick_inter7;
295                        if (dp.n_shells > 7)
296                                rad_out += dp.thick_flat8 + dp.thick_inter8;
297                        if (dp.n_shells > 8)
298                                rad_out += dp.thick_flat9 + dp.thick_inter9;
299                        if (dp.n_shells > 9)
300                                rad_out += dp.thick_flat10 + dp.thick_inter10;
301                        sum += weights_rad_core0[i].weight*weights_thick_inter0[j].weight
302                                * (rad_out);
303                        norm += weights_rad_core0[i].weight*weights_thick_inter0[j].weight;
304                }
305        }
306        if (norm != 0){
307                //return the averaged value
308                out =  sum/norm;}
309        else{
310                //return normal value
311                out = dp.rad_core0 + dp.thick_inter0;
312                if (dp.n_shells > 0)
313                        out += dp.thick_flat1 + dp.thick_inter1;
314                if (dp.n_shells > 1)
315                        out += dp.thick_flat2 + dp.thick_inter2;
316                if (dp.n_shells > 2)
317                        out += dp.thick_flat3 + dp.thick_inter3;
318                if (dp.n_shells > 3)
319                        out += dp.thick_flat4 + dp.thick_inter4;
320                if (dp.n_shells > 4)
321                        out += dp.thick_flat5 + dp.thick_inter5;
322                if (dp.n_shells > 5)
323                        out += dp.thick_flat6 + dp.thick_inter6;
324                if (dp.n_shells > 6)
325                        out += dp.thick_flat7 + dp.thick_inter7;
326                if (dp.n_shells > 7)
327                        out += dp.thick_flat8 + dp.thick_inter8;
328                if (dp.n_shells > 8)
329                        out += dp.thick_flat9 + dp.thick_inter9;
330                if (dp.n_shells > 9)
331                        out += dp.thick_flat10 + dp.thick_inter10;
332                }
333
334        return out;
[0164899a]335
336}
Note: See TracBrowser for help on using the repository browser.