source: sasview/sansmodels/src/sans/models/c_models/ellipticalcylinder.cpp @ dfa8832

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