source: sasview/sansmodels/src/sans/models/c_models/hollowcylinder.cpp @ 97603c0

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 97603c0 was c451be9, checked in by Jae Cho <jhjcho@…>, 15 years ago

corrections on the definition of polydispersity as suggested by steve K: should be normalized by average volume

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