source: sasview/sansmodels/src/c_models/stackeddisks.cpp @ 6e10cff

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 6e10cff was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Reorganizing models in preparation of cpp cleanup

  • Property mode set to 100644
File size: 9.3 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 *      TODO: refactor so that we pull in the old sansmodels.c_extensions
21 *      TODO: add 2d
22 */
23
24#include <math.h>
25#include "models.hh"
26#include "parameters.hh"
27#include <stdio.h>
28using namespace std;
29
30extern "C" {
31        #include "libCylinder.h"
32        #include "libStructureFactor.h"
33        #include "stacked_disks.h"
34}
35
36StackedDisksModel :: StackedDisksModel() {
37        scale      = Parameter(1.0);
38        radius     = Parameter(3000.0, true);
39        radius.set_min(0.0);
40        core_thick  = Parameter(10.0, true);
41        core_thick.set_min(0.0);
42        layer_thick     = Parameter(15.0);
43        layer_thick.set_min(0.0);
44        core_sld = Parameter(4.0e-6);
45        layer_sld  = Parameter(-4.0e-7);
46        solvent_sld  = Parameter(5.0e-6);
47        n_stacking   = Parameter(1);
48        sigma_d   = Parameter(0);
49        background = Parameter(0.001);
50        axis_theta  = Parameter(0.0, true);
51        axis_phi    = Parameter(0.0, true);
52}
53
54/**
55 * Function to evaluate 1D scattering function
56 * The NIST IGOR library is used for the actual calculation.
57 * @param q: q-value
58 * @return: function value
59 */
60double StackedDisksModel :: operator()(double q) {
61        double dp[10];
62
63        // Fill parameter array for IGOR library
64        // Add the background after averaging
65        dp[0] = scale();
66        dp[1] = radius();
67        dp[2] = core_thick();
68        dp[3] = layer_thick();
69        dp[4] = core_sld();
70        dp[5] = layer_sld();
71        dp[6] = solvent_sld();
72        dp[7] = n_stacking();
73        dp[8] = sigma_d();
74        dp[9] = 0.0;
75
76        // Get the dispersion points for the radius
77        vector<WeightPoint> weights_radius;
78        radius.get_weights(weights_radius);
79
80        // Get the dispersion points for the core_thick
81        vector<WeightPoint> weights_core_thick;
82        core_thick.get_weights(weights_core_thick);
83
84        // Get the dispersion points for the layer_thick
85        vector<WeightPoint> weights_layer_thick;
86        layer_thick.get_weights(weights_layer_thick);
87
88        // Perform the computation, with all weight points
89        double sum = 0.0;
90        double norm = 0.0;
91        double vol = 0.0;
92
93        // Loop over length weight points
94        for(int i=0; i< (int)weights_radius.size(); i++) {
95                dp[1] = weights_radius[i].value;
96
97                // Loop over radius weight points
98                for(int j=0; j< (int)weights_core_thick.size(); j++) {
99                        dp[2] = weights_core_thick[j].value;
100
101                        // Loop over thickness weight points
102                        for(int k=0; k< (int)weights_layer_thick.size(); k++) {
103                                dp[3] = weights_layer_thick[k].value;
104                                //Un-normalize by volume
105                                sum += weights_radius[i].weight
106                                        * weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q)
107                                        *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
108                                //Find average volume
109                                vol += weights_radius[i].weight
110                                        * weights_core_thick[j].weight * weights_layer_thick[k].weight
111                                        *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
112                                norm += weights_radius[i].weight
113                                        * weights_core_thick[j].weight* weights_layer_thick[k].weight;
114                        }
115                }
116        }
117        if (vol != 0.0 && norm != 0.0) {
118                //Re-normalize by avg volume
119                sum = sum/(vol/norm);}
120
121        return sum/norm + background();
122}
123
124/**
125 * Function to evaluate 2D scattering function
126 * @param q_x: value of Q along x
127 * @param q_y: value of Q along y
128 * @return: function value
129 */
130double StackedDisksModel :: operator()(double qx, double qy) {
131        StackedDisksParameters dp;
132        // Fill parameter array
133        dp.scale      = scale();
134        dp.core_thick    = core_thick();
135        dp.radius         = radius();
136        dp.layer_thick  = layer_thick();
137        dp.core_sld   = core_sld();
138        dp.layer_sld  = layer_sld();
139        dp.solvent_sld= solvent_sld();
140        dp.n_stacking     = n_stacking();
141        dp.sigma_d   = sigma_d();
142        dp.background = 0.0;
143        dp.axis_theta = axis_theta();
144        dp.axis_phi   = axis_phi();
145
146        // Get the dispersion points for the length
147        vector<WeightPoint> weights_core_thick;
148        core_thick.get_weights(weights_core_thick);
149
150        // Get the dispersion points for the radius
151        vector<WeightPoint> weights_radius;
152        radius.get_weights(weights_radius);
153
154        // Get the dispersion points for the thickness
155        vector<WeightPoint> weights_layer_thick;
156        layer_thick.get_weights(weights_layer_thick);
157
158        // Get angular averaging for theta
159        vector<WeightPoint> weights_theta;
160        axis_theta.get_weights(weights_theta);
161
162        // Get angular averaging for phi
163        vector<WeightPoint> weights_phi;
164        axis_phi.get_weights(weights_phi);
165
166        // Perform the computation, with all weight points
167        double sum = 0.0;
168        double norm = 0.0;
169        double norm_vol = 0.0;
170        double vol = 0.0;
171
172        // Loop over length weight points
173        for(int i=0; i< (int)weights_core_thick.size(); i++) {
174                dp.core_thick = weights_core_thick[i].value;
175
176                // Loop over radius weight points
177                for(int j=0; j< (int)weights_radius.size(); j++) {
178                        dp.radius = weights_radius[j].value;
179
180                                // Loop over thickness weight points
181                                for(int k=0; k< (int)weights_layer_thick.size(); k++) {
182                                dp.layer_thick = weights_layer_thick[k].value;
183
184                                        for(int l=0; l< (int)weights_theta.size(); l++) {
185                                        dp.axis_theta = weights_theta[l].value;
186
187                                                // Average over phi distribution
188                                                for(int m=0; m <(int)weights_phi.size(); m++) {
189                                                        dp.axis_phi = weights_phi[m].value;
190
191                                                        //Un-normalize by volume
192                                                        double _ptvalue = weights_core_thick[i].weight
193                                                                * weights_radius[j].weight
194                                                                * weights_layer_thick[k].weight
195                                                                * weights_theta[l].weight
196                                                                * weights_phi[m].weight
197                                                                * stacked_disks_analytical_2DXY(&dp, qx, qy)
198                                                                *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
199                                                        if (weights_theta.size()>1) {
200                                                                _ptvalue *= fabs(sin(weights_theta[l].value));
201                                                        }
202                                                        sum += _ptvalue;
203                                                        //Find average volume
204                                                        vol += weights_radius[j].weight
205                                                                * weights_core_thick[i].weight * weights_layer_thick[k].weight
206                                                                *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
207                                                        //Find norm for volume
208                                                        norm_vol += weights_radius[j].weight
209                                                                * weights_core_thick[i].weight * weights_layer_thick[k].weight;
210
211                                                        norm += weights_core_thick[i].weight
212                                                                * weights_radius[j].weight
213                                                                * weights_layer_thick[k].weight
214                                                                * weights_theta[l].weight
215                                                                * weights_phi[m].weight;
216                                                }
217                                }
218                        }
219                }
220        }
221        // Averaging in theta needs an extra normalization
222        // factor to account for the sin(theta) term in the
223        // integration (see documentation).
224        if (weights_theta.size()>1) norm = norm / asin(1.0);
225        if (vol != 0.0 && norm_vol != 0.0) {
226                //Re-normalize by avg volume
227                sum = sum/(vol/norm_vol);}
228        return sum/norm + background();
229}
230
231/**
232 * Function to evaluate 2D scattering function
233 * @param pars: parameters of the triaxial ellipsoid
234 * @param q: q-value
235 * @param phi: angle phi
236 * @return: function value
237 */
238double StackedDisksModel :: evaluate_rphi(double q, double phi) {
239        double qx = q*cos(phi);
240        double qy = q*sin(phi);
241        return (*this).operator()(qx, qy);
242}
243/**
244 * Function to calculate effective radius
245 * @return: effective radius value
246 */
247double StackedDisksModel :: calculate_ER() {
248        StackedDisksParameters dp;
249
250        dp.core_thick    = core_thick();
251        dp.radius         = radius();
252        dp.layer_thick  = layer_thick();
253        dp.n_stacking     = n_stacking();
254
255        double rad_out = 0.0;
256        if (dp.n_stacking <= 0.0){
257                return rad_out;
258        }
259
260        // Perform the computation, with all weight points
261        double sum = 0.0;
262        double norm = 0.0;
263
264        // Get the dispersion points for the length
265        vector<WeightPoint> weights_core_thick;
266        core_thick.get_weights(weights_core_thick);
267
268        // Get the dispersion points for the radius
269        vector<WeightPoint> weights_radius;
270        radius.get_weights(weights_radius);
271
272        // Get the dispersion points for the thickness
273        vector<WeightPoint> weights_layer_thick;
274        layer_thick.get_weights(weights_layer_thick);
275
276        // Loop over major shell weight points
277        for(int i=0; i< (int)weights_core_thick.size(); i++) {
278                dp.core_thick = weights_core_thick[i].value;
279                for(int j=0; j< (int)weights_layer_thick.size(); j++) {
280                        dp.layer_thick = weights_layer_thick[j].value;
281                        for(int k=0; k< (int)weights_radius.size(); k++) {
282                                dp.radius = weights_radius[k].value;
283                                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
284                                sum +=weights_core_thick[i].weight*weights_layer_thick[j].weight
285                                        * weights_radius[k].weight*DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;
286                                norm += weights_core_thick[i].weight*weights_layer_thick[j].weight* weights_radius[k].weight;
287                        }
288                }
289        }
290        if (norm != 0){
291                //return the averaged value
292                rad_out =  sum/norm;}
293        else{
294                //return normal value
295                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
296                rad_out = DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;}
297
298        return rad_out;
299}
Note: See TracBrowser for help on using the repository browser.