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

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 8f5b34a was 34c2649, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Re #4 Fixed warnings

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