source: sasview/src/sas/models/c_extension/c_models/coreshellsphere.cpp @ 79492222

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 79492222 was 79492222, checked in by krzywon, 8 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 10.5 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 "core_shell.h"
27
28extern "C" {
29#include "libSphere.h"
30#include "libmultifunc/libfunc.h"
31}
32
33typedef struct {
34  double scale;
35  double radius;
36  double thickness;
37  double core_sld;
38  double shell_sld;
39  double solvent_sld;
40  double background;
41  double M0_sld_shell;
42  double M_theta_shell;
43  double M_phi_shell;
44  double M0_sld_core;
45  double M_theta_core;
46  double M_phi_core;
47  double M0_sld_solv;
48  double M_theta_solv;
49  double M_phi_solv;
50  double Up_frac_i;
51  double Up_frac_f;
52  double Up_theta;
53} CoreShellParameters;
54
55CoreShellModel :: CoreShellModel() {
56  scale      = Parameter(1.0);
57  radius     = Parameter(60.0, true);
58  radius.set_min(0.0);
59  thickness  = Parameter(10.0, true);
60  thickness.set_min(0.0);
61  core_sld   = Parameter(1.e-6);
62  shell_sld  = Parameter(2.e-6);
63  solvent_sld = Parameter(3.e-6);
64  background = Parameter(0.0);
65  M0_sld_shell = Parameter(0.0e-6);
66  M_theta_shell = Parameter(0.0);
67  M_phi_shell = Parameter(0.0);
68  M0_sld_core = Parameter(0.0e-6);
69  M_theta_core = Parameter(0.0);
70  M_phi_core = Parameter(0.0);
71  M0_sld_solv = Parameter(0.0e-6);
72  M_theta_solv = Parameter(0.0);
73  M_phi_solv = Parameter(0.0);
74  Up_frac_i = Parameter(0.5);
75  Up_frac_f = Parameter(0.5);
76  Up_theta = Parameter(0.0);
77}
78
79/**
80 * Function to evaluate 1D scattering function
81 * The NIST IGOR library is used for the actual calculation.
82 * @param q: q-value
83 * @return: function value
84 */
85double CoreShellModel :: operator()(double q) {
86  double dp[7];
87
88  // Fill parameter array for IGOR library
89  // Add the background after averaging
90
91  dp[0] = scale();
92  dp[1] = radius();
93  dp[2] = thickness();
94  dp[3] = core_sld();
95  dp[4] = shell_sld();
96  dp[5] = solvent_sld();
97  dp[6] = 0.0;
98
99  //im
100  ///dp[7] = 0.0;
101  ///dp[8] = 0.0;
102  ///dp[9] = 0.0;
103
104  // Get the dispersion points for the radius
105  vector<WeightPoint> weights_rad;
106  radius.get_weights(weights_rad);
107
108  // Get the dispersion points for the thickness
109  vector<WeightPoint> weights_thick;
110  thickness.get_weights(weights_thick);
111
112  // Perform the computation, with all weight points
113  double sum = 0.0;
114  double norm = 0.0;
115  double vol = 0.0;
116
117  // Loop over radius weight points
118  for(size_t i=0; i<weights_rad.size(); i++) {
119    dp[1] = weights_rad[i].value;
120
121    // Loop over thickness weight points
122    for(size_t j=0; j<weights_thick.size(); j++) {
123      dp[2] = weights_thick[j].value;
124      //Un-normalize SphereForm by volume
125      sum += weights_rad[i].weight
126          * weights_thick[j].weight * CoreShellForm(dp, q)* pow(weights_rad[i].value+weights_thick[j].value,3);
127
128      //Find average volume
129      vol += weights_rad[i].weight * weights_thick[j].weight
130          * pow(weights_rad[i].value+weights_thick[j].value,3);
131      norm += weights_rad[i].weight
132          * weights_thick[j].weight;
133    }
134  }
135
136  if (vol != 0.0 && norm != 0.0) {
137    //Re-normalize by avg volume
138    sum = sum/(vol/norm);}
139
140  return sum/norm + background();
141}
142
143
144
145/**
146 * Function to evaluate 2D scattering function
147 * @param pars: parameters
148 * @param q: q-value
149 * @param q_x: q_x / q
150 * @param q_y: q_y / q
151 * @return: function value
152 */
153
154static double coreshell_analytical_2D_scaled(CoreShellParameters *pars, double q, double q_x, double q_y) {
155        double dp[7];
156        //convert angle degree to radian
157        dp[0] = pars->scale;
158        dp[1] = pars->radius;
159        dp[2] = pars->thickness;
160        dp[3] = 0.0;
161        dp[4] = 0.0;
162        dp[5] = 0.0;
163        dp[6] = 0.0;
164        //im
165        ///dp[7] = 0.0;
166        ///dp[8] = 0.0;
167        ///dp[9] = 0.0;
168
169    double sld_core = pars->core_sld;
170    double sld_shell = pars->shell_sld;
171        double sld_solv = pars->solvent_sld;
172        double answer = 0.0;
173        double m_max = pars->M0_sld_core;
174        double m_max_shell = pars->M0_sld_shell;
175        double m_max_solv = pars->M0_sld_solv;
176
177        if (m_max < 1.0e-32 && m_max_solv < 1.0e-32 && m_max_shell < 1.0e-32){
178                dp[3] = sld_core;
179                dp[4] = sld_shell;
180                dp[5] = sld_solv;
181                answer = CoreShellForm(dp, q);
182        }
183        else{
184                double qx = q_x;
185                double qy = q_y;
186                double s_theta = pars->Up_theta;
187                double m_phi = pars->M_phi_core;
188                double m_theta = pars->M_theta_core;
189                double m_phi_shell = pars->M_phi_shell;
190                double m_theta_shell = pars->M_theta_shell;
191                double m_phi_solv = pars->M_phi_solv;
192                double m_theta_solv = pars->M_theta_solv;
193                double in_spin = pars->Up_frac_i;
194                double out_spin = pars->Up_frac_f;
195                polar_sld p_sld_core;
196                polar_sld p_sld_shell;
197                polar_sld p_sld_solv;
198                //Find (b+m) slds
199                p_sld_core = cal_msld(1, qx, qy, sld_core, m_max, m_theta, m_phi,
200                                                                in_spin, out_spin, s_theta);
201                p_sld_shell = cal_msld(1, qx, qy, sld_shell, m_max_shell, m_theta_shell, m_phi_shell,
202                                                                in_spin, out_spin, s_theta);
203                p_sld_solv = cal_msld(1, qx, qy, sld_solv, m_max_solv, m_theta_solv, m_phi_solv,
204                                                                in_spin, out_spin, s_theta);
205                //up_up
206                if (in_spin > 0.0 && out_spin > 0.0){
207                        dp[3] = p_sld_core.uu;
208                        dp[4] = p_sld_shell.uu;
209                        dp[5] = p_sld_solv.uu;
210                        answer += CoreShellForm(dp, q);
211                        }
212                //down_down
213                if (in_spin < 1.0 && out_spin < 1.0){
214                        dp[3] = p_sld_core.dd;
215                        dp[4] = p_sld_shell.dd;
216                        dp[5] = p_sld_solv.dd;
217                        answer += CoreShellForm(dp, q);
218                        }
219                //up_down
220                if (in_spin > 0.0 && out_spin < 1.0){
221                        dp[3] = p_sld_core.re_ud;
222                        dp[4] = p_sld_shell.re_ud;
223                        dp[5] = p_sld_solv.re_ud;
224                        answer += CoreShellForm(dp, q);
225                        dp[3] = p_sld_core.im_ud;
226                        dp[4] = p_sld_shell.im_ud;
227                        dp[5] = p_sld_solv.im_ud;
228                        answer += CoreShellForm(dp, q);
229                        }
230                //down_up
231                if (in_spin < 1.0 && out_spin > 0.0){
232                        dp[3] = p_sld_core.re_du;
233                        dp[4] = p_sld_shell.re_du;
234                        dp[5] = p_sld_solv.re_du;
235                        answer += CoreShellForm(dp, q);
236                        dp[3] = p_sld_core.im_du;
237                        dp[4] = p_sld_shell.im_du;
238                        dp[5] = p_sld_solv.im_du;
239                        answer += CoreShellForm(dp, q);
240                        }
241        }
242        // Already normalized
243        // add in the background
244        answer += pars->background;
245        return answer;
246}
247
248/**
249 * Function to evaluate 2D scattering function
250 * @param pars: parameters
251 * @param q: q-value
252 * @return: function value
253 */
254static double coreshell_analytical_2DXY(CoreShellParameters *pars, double qx, double qy) {
255  double q;
256  q = sqrt(qx*qx+qy*qy);
257  return coreshell_analytical_2D_scaled(pars, q, qx/q, qy/q);
258}
259
260
261/**
262 * Function to evaluate 2D scattering function
263 * @param q_x: value of Q along x
264 * @param q_y: value of Q along y
265 * @return: function value
266 */
267double CoreShellModel :: operator()(double qx, double qy) {
268        CoreShellParameters dp;
269    dp.scale = scale();
270    dp.radius = radius();
271    dp.thickness = thickness();
272    dp.core_sld = core_sld();
273    dp.shell_sld = shell_sld();
274    dp.solvent_sld = solvent_sld();
275    dp.background = 0.0;
276    dp.M0_sld_shell = M0_sld_shell();
277    dp.M_theta_shell = M_theta_shell();
278    dp.M_phi_shell = M_phi_shell();
279    dp.M0_sld_core = M0_sld_core();
280    dp.M_theta_core = M_theta_core();
281    dp.M_phi_core = M_phi_core();
282    dp.M0_sld_solv = M0_sld_solv();
283    dp.M_theta_solv = M_theta_solv();
284    dp.M_phi_solv = M_phi_solv();
285    dp.Up_frac_i = Up_frac_i();
286    dp.Up_frac_f = Up_frac_f();
287    dp.Up_theta = Up_theta();
288
289        // Get the dispersion points for the radius
290        vector<WeightPoint> weights_rad;
291        radius.get_weights(weights_rad);
292
293        // Get the dispersion points for the thickness
294        vector<WeightPoint> weights_thick;
295        thickness.get_weights(weights_thick);
296
297        // Perform the computation, with all weight points
298        double sum = 0.0;
299        double norm = 0.0;
300        double vol = 0.0;
301
302        // Loop over radius weight points
303        for(size_t i=0; i<weights_rad.size(); i++) {
304        dp.radius = weights_rad[i].value;
305
306        // Loop over thickness weight points
307        for(size_t j=0; j<weights_thick.size(); j++) {
308          dp.thickness = weights_thick[j].value;
309          //Un-normalize SphereForm by volume
310          sum += weights_rad[i].weight
311                  * weights_thick[j].weight * coreshell_analytical_2DXY(&dp, qx, qy) *
312                        pow(weights_rad[i].value+weights_thick[j].value,3);
313
314          //Find average volume
315          vol += weights_rad[i].weight * weights_thick[j].weight
316                  * pow(weights_rad[i].value+weights_thick[j].value,3);
317          norm += weights_rad[i].weight
318                  * weights_thick[j].weight;
319        }
320        }
321
322        if (vol != 0.0 && norm != 0.0) {
323                //Re-normalize by avg volume
324                sum = sum/(vol/norm);}
325        return sum/norm + background();
326}
327
328/**
329 * Function to evaluate 2D scattering function
330 * @param pars: parameters of the sphere
331 * @param q: q-value
332 * @param phi: angle phi
333 * @return: function value
334 */
335double CoreShellModel :: evaluate_rphi(double q, double phi) {
336        double qx = q*cos(phi);
337        double qy = q*sin(phi);
338        return (*this).operator()(qx, qy);
339}
340/**
341 * Function to calculate effective radius
342 * @return: effective radius value
343 */
344double CoreShellModel :: calculate_ER() {
345  CoreShellParameters dp;
346
347  dp.radius     = radius();
348  dp.thickness  = thickness();
349
350  double rad_out = 0.0;
351
352  // Perform the computation, with all weight points
353  double sum = 0.0;
354  double norm = 0.0;
355
356
357  // Get the dispersion points for the major shell
358  vector<WeightPoint> weights_thickness;
359  thickness.get_weights(weights_thickness);
360
361  // Get the dispersion points for the minor shell
362  vector<WeightPoint> weights_radius ;
363  radius.get_weights(weights_radius);
364
365  // Loop over major shell weight points
366  for(int j=0; j< (int)weights_thickness.size(); j++) {
367    dp.thickness = weights_thickness[j].value;
368    for(int k=0; k< (int)weights_radius.size(); k++) {
369      dp.radius = weights_radius[k].value;
370      sum += weights_thickness[j].weight
371          * weights_radius[k].weight*(dp.radius+dp.thickness);
372      norm += weights_thickness[j].weight* weights_radius[k].weight;
373    }
374  }
375  if (norm != 0){
376    //return the averaged value
377    rad_out =  sum/norm;}
378  else{
379    //return normal value
380    rad_out = (dp.radius+dp.thickness);}
381
382  return rad_out;
383}
384double CoreShellModel :: calculate_VR() {
385  return 1.0;
386}
Note: See TracBrowser for help on using the repository browser.