source: sasview/sansmodels/src/c_models/corefourshell.cpp @ d555416

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 d555416 was 318b5bbb, checked in by Jae Cho <jhjcho@…>, 12 years ago

Added polarization and magnetic stuffs
  • Property mode set to 100644
File size: 17.2 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 "parameters.hh"
24#include <stdio.h>
25using namespace std;
26#include "corefourshell.h"
27
28extern "C" {
29        #include "libSphere.h"
30        #include "libmultifunc/libfunc.h"
31}
32
33typedef struct {
34  double scale;
35  double rad_core0;
36  double sld_core0;
37  double thick_shell1;
38  double sld_shell1;
39  double thick_shell2;
40  double sld_shell2;
41  double thick_shell3;
42  double sld_shell3;
43  double thick_shell4;
44  double sld_shell4;
45  double sld_solv;
46  double background;
47  double M0_sld_shell1;
48  double M_theta_shell1;
49  double M_phi_shell1;
50  double M0_sld_shell2;
51  double M_theta_shell2;
52  double M_phi_shell2;
53  double M0_sld_shell3;
54  double M_theta_shell3;
55  double M_phi_shell3;
56  double M0_sld_shell4;
57  double M_theta_shell4;
58  double M_phi_shell4;
59  double M0_sld_core0;
60  double M_theta_core0;
61  double M_phi_core0;
62  double M0_sld_solv;
63  double M_theta_solv;
64  double M_phi_solv;
65  double Up_frac_i;
66  double Up_frac_f;
67  double Up_theta;
68} CoreFourShellParameters;
69
70CoreFourShellModel :: CoreFourShellModel() {
71        scale      = Parameter(1.0);
72        rad_core0     = Parameter(60.0, true);
73        rad_core0.set_min(0.0);
74        sld_core0   = Parameter(6.4e-6);
75        thick_shell1     = Parameter(10.0, true);
76        thick_shell1.set_min(0.0);
77        sld_shell1   = Parameter(1.0e-6);
78        thick_shell2     = Parameter(10.0, true);
79        thick_shell2.set_min(0.0);
80        sld_shell2   = Parameter(2.0e-6);
81        thick_shell3     = Parameter(10.0, true);
82        thick_shell3.set_min(0.0);
83        sld_shell3   = Parameter(3.0e-6);
84        thick_shell4     = Parameter(10.0, true);
85        thick_shell4.set_min(0.0);
86        sld_shell4   = Parameter(4.0e-6);
87        sld_solv   = Parameter(6.4e-6);
88        background = Parameter(0.001);
89        M0_sld_shell1 = Parameter(0.0e-6);
90        M_theta_shell1 = Parameter(0.0);
91        M_phi_shell1 = Parameter(0.0);
92        M0_sld_shell2 = Parameter(0.0e-6);
93        M_theta_shell2 = Parameter(0.0);
94        M_phi_shell2 = Parameter(0.0);
95        M0_sld_shell3 = Parameter(0.0e-6);
96        M_theta_shell3 = Parameter(0.0);
97        M_phi_shell3 = Parameter(0.0);
98        M0_sld_shell4 = Parameter(0.0e-6);
99        M_theta_shell4 = Parameter(0.0);
100        M_phi_shell4 = Parameter(0.0);
101        M0_sld_core0 = Parameter(0.0e-6);
102        M_theta_core0 = Parameter(0.0);
103        M_phi_core0 = Parameter(0.0);
104        M0_sld_solv = Parameter(0.0e-6);
105        M_theta_solv = Parameter(0.0);
106        M_phi_solv = Parameter(0.0);
107        Up_frac_i = Parameter(0.5);
108        Up_frac_f = Parameter(0.5);
109        Up_theta = Parameter(0.0);
110
111}
112
113/**
114 * Function to evaluate 1D scattering function
115 * The NIST IGOR library is used for the actual calculation.
116 * @param q: q-value
117 * @return: function value
118 */
119double CoreFourShellModel :: operator()(double q) {
120        double dp[13];
121
122        // Fill parameter array for IGOR library
123        // Add the background after averaging
124
125        dp[0] = scale();
126        dp[1] = rad_core0();
127        dp[2] = sld_core0();
128        dp[3] = thick_shell1();
129        dp[4] = sld_shell1();
130        dp[5] = thick_shell2();
131        dp[6] = sld_shell2();
132        dp[7] = thick_shell3();
133        dp[8] = sld_shell3();
134        dp[9] = thick_shell4();
135        dp[10] = sld_shell4();
136        dp[11] = sld_solv();
137        dp[12] = 0.0;
138
139        // Get the dispersion points for the radius
140        vector<WeightPoint> weights_rad;
141        rad_core0.get_weights(weights_rad);
142
143        // Get the dispersion points for the thick 1
144        vector<WeightPoint> weights_s1;
145        thick_shell1.get_weights(weights_s1);
146
147        // Get the dispersion points for the thick 2
148        vector<WeightPoint> weights_s2;
149        thick_shell2.get_weights(weights_s2);
150
151        // Get the dispersion points for the thick 3
152        vector<WeightPoint> weights_s3;
153        thick_shell3.get_weights(weights_s3);
154
155        // Get the dispersion points for the thick 4
156        vector<WeightPoint> weights_s4;
157        thick_shell4.get_weights(weights_s4);
158
159        // Perform the computation, with all weight points
160        double sum = 0.0;
161        double norm = 0.0;
162        double vol = 0.0;
163
164        // Loop over radius weight points
165        for(size_t i=0; i<weights_rad.size(); i++) {
166                dp[1] = weights_rad[i].value;
167                // Loop over radius weight points
168                for(size_t j=0; j<weights_s1.size(); j++) {
169                        dp[3] = weights_s1[j].value;
170                        // Loop over radius weight points
171                        for(size_t k=0; k<weights_s2.size(); k++) {
172                                dp[5] = weights_s2[k].value;
173                                // Loop over radius weight points
174                                for(size_t l=0; l<weights_s3.size(); l++) {
175                                        dp[7] = weights_s3[l].value;
176                                        // Loop over radius weight points
177                                        for(size_t m=0; m<weights_s4.size(); m++) {
178                                                dp[9] = weights_s4[m].value;
179                                                //Un-normalize FourShell by volume
180                                                sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
181                                                        * FourShell(dp, q) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3);
182                                                //Find average volume
183                                                vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
184                                                        * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3);
185
186                                                norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight;
187                                        }
188                                }
189                        }
190                }
191        }
192
193        if (vol != 0.0 && norm != 0.0) {
194                //Re-normalize by avg volume
195                sum = sum/(vol/norm);}
196        return sum/norm + background();
197}
198
199/**
200 * Function to evaluate 2D scattering function
201 * @param q_x: value of Q along x
202 * @param q_y: value of Q along y
203 * @return: function value
204 */
205
206static double corefourshell_analytical_2D_scaled(CoreFourShellParameters *pars, double q, double q_x, double q_y) {
207        double dp[13];
208
209        // Fill parameter array for IGOR library
210        // Add the background after averaging
211
212        dp[0] = pars->scale;
213        dp[1] = pars->rad_core0;
214        dp[2] = 0.0; //sld_core0;
215        dp[3] = pars->thick_shell1;
216        dp[4] = 0.0; //sld_shell1;
217        dp[5] = pars->thick_shell2;
218        dp[6] = 0.0; //sld_shell2;
219        dp[7] = pars->thick_shell3;
220        dp[8] = 0.0; //sld_shell3;
221        dp[9] = pars->thick_shell4;
222        dp[10] = 0.0; //sld_shell4;
223        dp[11] = 0.0; //sld_solv;
224        dp[12] = 0.0;
225
226    double sld_core0 = pars->sld_core0;
227    double sld_shell1 = pars->sld_shell1;
228        double sld_shell2 = pars->sld_shell2;
229        double sld_shell3 = pars->sld_shell3;
230        double sld_shell4 = pars->sld_shell4;
231        double sld_solv = pars->sld_solv;
232        double answer = 0.0;
233        double m_max0 = pars->M0_sld_core0;
234        double m_max_shell1 = pars->M0_sld_shell1;
235        double m_max_shell2 = pars->M0_sld_shell2;
236        double m_max_shell3 = pars->M0_sld_shell3;
237        double m_max_shell4 = pars->M0_sld_shell4;
238        double m_max_solv = pars->M0_sld_solv;
239
240        if (m_max0 < 1.0e-32 && m_max_solv < 1.0e-32 && m_max_shell1 < 1.0e-32 && m_max_shell2 <
241                                                1.0e-32 && m_max_shell3 < 1.0e-32 && m_max_shell4 < 1.0e-32){
242                dp[2] = sld_core0;
243                dp[4] = sld_shell1;
244                dp[6] = sld_shell2;
245                dp[8] = sld_shell3;
246                dp[10] = sld_shell4;
247                dp[11] = sld_solv;
248                answer = FourShell(dp, q);
249        }
250        else{
251                double qx = q_x;
252                double qy = q_y;
253                double s_theta = pars->Up_theta;
254                double m_phi0 = pars->M_phi_core0;
255                double m_theta0 = pars->M_theta_core0;
256                double m_phi_shell1 = pars->M_phi_shell1;
257                double m_theta_shell1 = pars->M_theta_shell1;
258                double m_phi_shell2 = pars->M_phi_shell2;
259                double m_theta_shell2 = pars->M_theta_shell2;
260                double m_phi_shell3 = pars->M_phi_shell3;
261                double m_theta_shell3 = pars->M_theta_shell3;
262                double m_phi_shell4 = pars->M_phi_shell4;
263                double m_theta_shell4 = pars->M_theta_shell4;
264                double m_phi_solv = pars->M_phi_solv;
265                double m_theta_solv = pars->M_theta_solv;
266                double in_spin = pars->Up_frac_i;
267                double out_spin = pars->Up_frac_f;
268                polar_sld p_sld_core0;
269                polar_sld p_sld_shell1;
270                polar_sld p_sld_shell2;
271                polar_sld p_sld_shell3;
272                polar_sld p_sld_shell4;
273                polar_sld p_sld_solv;
274                //Find (b+m) slds
275                p_sld_core0 = cal_msld(1, qx, qy, sld_core0, m_max0, m_theta0, m_phi0,
276                                                                in_spin, out_spin, s_theta);
277                p_sld_shell1 = cal_msld(1, qx, qy, sld_shell1, m_max_shell1, m_theta_shell1, m_phi_shell1,
278                                                                in_spin, out_spin, s_theta);
279                p_sld_shell2 = cal_msld(1, qx, qy, sld_shell2, m_max_shell2, m_theta_shell2, m_phi_shell2,
280                                                                in_spin, out_spin, s_theta);
281                p_sld_shell3 = cal_msld(1, qx, qy, sld_shell3, m_max_shell3, m_theta_shell3, m_phi_shell3,
282                                                                in_spin, out_spin, s_theta);
283                p_sld_shell4 = cal_msld(1, qx, qy, sld_shell4, m_max_shell4, m_theta_shell4, m_phi_shell4,
284                                                                in_spin, out_spin, s_theta);
285                p_sld_solv = cal_msld(1, qx, qy, sld_solv, m_max_solv, m_theta_solv, m_phi_solv,
286                                                                in_spin, out_spin, s_theta);
287                //up_up
288                if (in_spin > 0.0 && out_spin > 0.0){
289                        dp[2] = p_sld_core0.uu;
290                        dp[4] = p_sld_shell1.uu;
291                        dp[6] = p_sld_shell2.uu;
292                        dp[8] = p_sld_shell3.uu;
293                        dp[10] = p_sld_shell4.uu;
294                        dp[11] = p_sld_solv.uu;
295                        answer += FourShell(dp, q);
296                        }
297                //down_down
298                if (in_spin < 1.0 && out_spin < 1.0){
299                        dp[2] = p_sld_core0.dd;
300                        dp[4] = p_sld_shell1.dd;
301                        dp[6] = p_sld_shell2.dd;
302                        dp[8] = p_sld_shell3.dd;
303                        dp[10] = p_sld_shell4.dd;
304                        dp[11] = p_sld_solv.dd;
305                        answer += FourShell(dp, q);
306                        }
307                //up_down
308                if (in_spin > 0.0 && out_spin < 1.0){
309                        dp[2] = p_sld_core0.re_ud;
310                        dp[4] = p_sld_shell1.re_ud;
311                        dp[6] = p_sld_shell2.re_ud;
312                        dp[8] = p_sld_shell3.re_ud;
313                        dp[10] = p_sld_shell4.re_ud;
314                        dp[11] = p_sld_solv.re_ud;
315                        answer += FourShell(dp, q);
316                        dp[2] = p_sld_core0.im_ud;
317                        dp[4] = p_sld_shell1.im_ud;
318                        dp[6] = p_sld_shell2.im_ud;
319                        dp[8] = p_sld_shell3.im_ud;
320                        dp[10] = p_sld_shell4.im_ud;
321                        dp[11] = p_sld_solv.im_ud;
322                        answer += FourShell(dp, q);
323                        }
324                //down_up
325                if (in_spin < 1.0 && out_spin > 0.0){
326                        dp[2] = p_sld_core0.re_du;
327                        dp[4] = p_sld_shell1.re_du;
328                        dp[6] = p_sld_shell2.re_du;
329                        dp[8] = p_sld_shell3.re_du;
330                        dp[10] = p_sld_shell4.re_du;
331                        dp[11] = p_sld_solv.re_du;
332                        answer += FourShell(dp, q);
333                        dp[2] = p_sld_core0.im_du;
334                        dp[4] = p_sld_shell1.im_du;
335                        dp[6] = p_sld_shell2.im_du;
336                        dp[8] = p_sld_shell3.im_du;
337                        dp[10] = p_sld_shell4.im_du;
338                        dp[11] = p_sld_solv.im_du;
339                        answer += FourShell(dp, q);
340                        }
341        }
342        // Already normalized
343        // add in the background
344        answer += pars->background;
345        return answer;
346}
347
348/**
349 * Function to evaluate 2D scattering function
350 * @param pars: parameters
351 * @param q: q-value
352 * @return: function value
353 */
354static double corefourshell_analytical_2DXY(CoreFourShellParameters *pars, double qx, double qy) {
355  double q;
356  q = sqrt(qx*qx+qy*qy);
357  return corefourshell_analytical_2D_scaled(pars, q, qx/q, qy/q);
358}
359
360
361double CoreFourShellModel :: operator()(double qx, double qy) {
362        CoreFourShellParameters dp;
363        dp.scale = scale();
364        dp.rad_core0 = rad_core0();
365        dp.sld_core0 = sld_core0();
366        dp.thick_shell1 = thick_shell1();
367        dp.sld_shell1 = sld_shell1();
368        dp.thick_shell2 = thick_shell2();
369        dp.sld_shell2 = sld_shell2();
370        dp.thick_shell3 = thick_shell3();
371        dp.sld_shell3 = sld_shell3();
372        dp.thick_shell4 = thick_shell4();
373        dp.sld_shell4 = sld_shell4();
374        dp.sld_solv = sld_solv();
375        dp.background = 0.0;
376        dp.M0_sld_shell1 = M0_sld_shell1();
377        dp.M_theta_shell1 = M_theta_shell1();
378        dp.M_phi_shell1 = M_phi_shell1();
379        dp.M0_sld_shell2 = M0_sld_shell2();
380        dp.M_theta_shell2 = M_theta_shell2();
381        dp.M_phi_shell2 = M_phi_shell2();
382        dp.M0_sld_shell3 = M0_sld_shell3();
383        dp.M_theta_shell3 = M_theta_shell3();
384        dp.M_phi_shell3 = M_phi_shell3();
385        dp.M0_sld_shell4 = M0_sld_shell4();
386        dp.M_theta_shell4 = M_theta_shell4();
387        dp.M_phi_shell4 = M_phi_shell4();
388        dp.M0_sld_core0 = M0_sld_core0();
389        dp.M_theta_core0 = M_theta_core0();
390        dp.M_phi_core0 = M_phi_core0();
391        dp.M0_sld_solv = M0_sld_solv();
392        dp.M_theta_solv = M_theta_solv();
393        dp.M_phi_solv = M_phi_solv();
394        dp.Up_frac_i = Up_frac_i();
395        dp.Up_frac_f = Up_frac_f();
396        dp.Up_theta = Up_theta();
397   
398        // Get the dispersion points for the radius
399        vector<WeightPoint> weights_rad;
400        rad_core0.get_weights(weights_rad);
401
402        // Get the dispersion points for the thick 1
403        vector<WeightPoint> weights_s1;
404        thick_shell1.get_weights(weights_s1);
405
406        // Get the dispersion points for the thick 2
407        vector<WeightPoint> weights_s2;
408        thick_shell2.get_weights(weights_s2);
409
410        // Get the dispersion points for the thick 3
411        vector<WeightPoint> weights_s3;
412        thick_shell3.get_weights(weights_s3);
413
414        // Get the dispersion points for the thick 4
415        vector<WeightPoint> weights_s4;
416        thick_shell4.get_weights(weights_s4);
417
418        // Perform the computation, with all weight points
419        double sum = 0.0;
420        double norm = 0.0;
421        double vol = 0.0;
422
423        // Loop over radius weight points
424        for(size_t i=0; i<weights_rad.size(); i++) {
425                dp.rad_core0 = weights_rad[i].value;
426                // Loop over radius weight points
427                for(size_t j=0; j<weights_s1.size(); j++) {
428                        dp.thick_shell1 = weights_s1[j].value;
429                        // Loop over radius weight points
430                        for(size_t k=0; k<weights_s2.size(); k++) {
431                                dp.thick_shell2 = weights_s2[k].value;
432                                // Loop over radius weight points
433                                for(size_t l=0; l<weights_s3.size(); l++) {
434                                        dp.thick_shell3 = weights_s3[l].value;
435                                        // Loop over radius weight points
436                                        for(size_t m=0; m<weights_s4.size(); m++) {
437                                                dp.thick_shell4 = weights_s4[m].value;
438                                                //Un-normalize FourShell by volume
439                                                sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
440                                                        * corefourshell_analytical_2DXY(&dp, qx, qy) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3);
441                                                //Find average volume
442                                                vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
443                                                        * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3);
444
445                                                norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight;
446                                        }
447                                }
448                        }
449                }
450        }
451
452        if (vol != 0.0 && norm != 0.0) {
453                //Re-normalize by avg volume
454                sum = sum/(vol/norm);}
455        return sum/norm + background();
456}
457
458/**
459 * Function to evaluate 2D scattering function
460 * @param pars: parameters of the sphere
461 * @param q: q-value
462 * @param phi: angle phi
463 * @return: function value
464 */
465double CoreFourShellModel :: evaluate_rphi(double q, double phi) {
466        double qx = q*cos(phi);
467        double qy = q*sin(phi);
468        return (*this).operator()(qx, qy);
469}
470
471/**
472 * Function to calculate effective radius
473 * @return: effective radius value
474 */
475double CoreFourShellModel :: calculate_ER() {
476        CoreFourShellParameters dp;
477
478        dp.scale = scale();
479        dp.rad_core0 = rad_core0();
480        dp.sld_core0 = sld_core0();
481        dp.thick_shell1 = thick_shell1();
482        dp.sld_shell1 = sld_shell1();
483        dp.thick_shell2 = thick_shell2();
484        dp.sld_shell2 = sld_shell2();
485        dp.thick_shell3 = thick_shell3();
486        dp.sld_shell3 = sld_shell3();
487        dp.thick_shell4 = thick_shell4();
488        dp.sld_shell4 = sld_shell4();
489        dp.sld_solv = sld_solv();
490        dp.background = 0.0;
491
492        // Get the dispersion points for the radius
493        vector<WeightPoint> weights_rad;
494        rad_core0.get_weights(weights_rad);
495
496        // Get the dispersion points for the thick 1
497        vector<WeightPoint> weights_s1;
498        thick_shell1.get_weights(weights_s1);
499
500        // Get the dispersion points for the thick 2
501        vector<WeightPoint> weights_s2;
502        thick_shell2.get_weights(weights_s2);
503
504        // Get the dispersion points for the thick 3
505        vector<WeightPoint> weights_s3;
506        thick_shell3.get_weights(weights_s3);
507
508        // Get the dispersion points for the thick 4
509        vector<WeightPoint> weights_s4;
510        thick_shell4.get_weights(weights_s4);
511
512        double rad_out = 0.0;
513        // Perform the computation, with all weight points
514        double sum = 0.0;
515        double norm = 0.0;
516
517        // Loop over radius weight points
518        for(size_t i=0; i<weights_rad.size(); i++) {
519                dp.rad_core0 = weights_rad[i].value;
520                // Loop over radius weight points
521                for(size_t j=0; j<weights_s1.size(); j++) {
522                        dp.thick_shell1 = weights_s1[j].value;
523                        // Loop over radius weight points
524                        for(size_t k=0; k<weights_s2.size(); k++) {
525                                dp.thick_shell2 = weights_s2[k].value;
526                                // Loop over radius weight points
527                                for(size_t l=0; l<weights_s3.size(); l++) {
528                                        dp.thick_shell3 = weights_s3[l].value;
529                                        // Loop over radius weight points
530                                        for(size_t m=0; m<weights_s4.size(); m++) {
531                                                dp.thick_shell4 = weights_s4[m].value;
532                                                //Un-normalize FourShell by volume
533                                                sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
534                                                        * (dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4);
535                                                norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight;
536                                        }
537                                }
538                        }
539                }
540        }
541        if (norm != 0){
542                //return the averaged value
543                rad_out =  sum/norm;}
544        else{
545                //return normal value
546                rad_out = dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4;}
547
548        return rad_out;
549}
550double CoreFourShellModel :: calculate_VR() {
551  return 1.0;
552}
Note: See TracBrowser for help on using the repository browser.