source: sasview/sansmodels/src/c_models/parallelepiped.cpp @ acd0fd10

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 acd0fd10 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.7 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 function
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 "parallelepiped.h"
34}
35
36ParallelepipedModel :: ParallelepipedModel() {
37        scale      = Parameter(1.0);
38        short_a     = Parameter(35.0, true);
39        short_a.set_min(1.0);
40        short_b     = Parameter(75.0, true);
41        short_b.set_min(1.0);
42        long_c     = Parameter(400.0, true);
43        long_c.set_min(1.0);
44        sldPipe   = Parameter(6.3e-6);
45        sldSolv   = Parameter(1.0e-6);
46        background = Parameter(0.0);
47        parallel_theta  = Parameter(0.0, true);
48        parallel_phi    = Parameter(0.0, true);
49        parallel_psi    = Parameter(0.0, true);
50}
51
52/**
53 * Function to evaluate 1D scattering function
54 * The NIST IGOR library is used for the actual calculation.
55 * @param q: q-value
56 * @return: function value
57 */
58double ParallelepipedModel :: operator()(double q) {
59        double dp[7];
60
61        // Fill parameter array for IGOR library
62        // Add the background after averaging
63        dp[0] = scale();
64        dp[1] = short_a();
65        dp[2] = short_b();
66        dp[3] = long_c();
67        dp[4] = sldPipe();
68        dp[5] = sldSolv();
69        dp[6] = 0.0;
70
71        // Get the dispersion points for the short_edgeA
72        vector<WeightPoint> weights_short_a;
73        short_a.get_weights(weights_short_a);
74
75        // Get the dispersion points for the longer_edgeB
76        vector<WeightPoint> weights_short_b;
77        short_b.get_weights(weights_short_b);
78
79        // Get the dispersion points for the longuest_edgeC
80        vector<WeightPoint> weights_long_c;
81        long_c.get_weights(weights_long_c);
82
83
84
85        // Perform the computation, with all weight points
86        double sum = 0.0;
87        double norm = 0.0;
88        double vol = 0.0;
89
90        // Loop over short_edgeA weight points
91        for(int i=0; i< (int)weights_short_a.size(); i++) {
92                dp[1] = weights_short_a[i].value;
93
94                // Loop over longer_edgeB weight points
95                for(int j=0; j< (int)weights_short_b.size(); j++) {
96                        dp[2] = weights_short_b[j].value;
97
98                        // Loop over longuest_edgeC weight points
99                        for(int k=0; k< (int)weights_long_c.size(); k++) {
100                                dp[3] = weights_long_c[k].value;
101                                //Un-normalize  by volume
102                                sum += weights_short_a[i].weight * weights_short_b[j].weight
103                                        * weights_long_c[k].weight * Parallelepiped(dp, q)
104                                        * weights_short_a[i].value*weights_short_b[j].value
105                                        * weights_long_c[k].value;
106                                //Find average volume
107                                vol += weights_short_a[i].weight * weights_short_b[j].weight
108                                        * weights_long_c[k].weight
109                                        * weights_short_a[i].value * weights_short_b[j].value
110                                        * weights_long_c[k].value;
111
112                                norm += weights_short_a[i].weight
113                                         * weights_short_b[j].weight * weights_long_c[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 * 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 ParallelepipedModel :: operator()(double qx, double qy) {
130        ParallelepipedParameters dp;
131        // Fill parameter array
132        dp.scale      = scale();
133        dp.short_a   = short_a();
134        dp.short_b   = short_b();
135        dp.long_c  = long_c();
136        dp.sldPipe   = sldPipe();
137        dp.sldSolv   = sldSolv();
138        dp.background = 0.0;
139        //dp.background = background();
140        dp.parallel_theta  = parallel_theta();
141        dp.parallel_phi    = parallel_phi();
142        dp.parallel_psi    = parallel_psi();
143
144
145        // Get the dispersion points for the short_edgeA
146        vector<WeightPoint> weights_short_a;
147        short_a.get_weights(weights_short_a);
148
149        // Get the dispersion points for the longer_edgeB
150        vector<WeightPoint> weights_short_b;
151        short_b.get_weights(weights_short_b);
152
153        // Get angular averaging for the longuest_edgeC
154        vector<WeightPoint> weights_long_c;
155        long_c.get_weights(weights_long_c);
156
157        // Get angular averaging for theta
158        vector<WeightPoint> weights_parallel_theta;
159        parallel_theta.get_weights(weights_parallel_theta);
160
161        // Get angular averaging for phi
162        vector<WeightPoint> weights_parallel_phi;
163        parallel_phi.get_weights(weights_parallel_phi);
164
165        // Get angular averaging for psi
166        vector<WeightPoint> weights_parallel_psi;
167        parallel_psi.get_weights(weights_parallel_psi);
168
169        // Perform the computation, with all weight points
170        double sum = 0.0;
171        double norm = 0.0;
172        double norm_vol = 0.0;
173        double vol = 0.0;
174        double pi = 4.0*atan(1.0);
175        // Loop over radius weight points
176        for(int i=0; i< (int)weights_short_a.size(); i++) {
177                dp.short_a = weights_short_a[i].value;
178
179                // Loop over longer_edgeB weight points
180                for(int j=0; j< (int)weights_short_b.size(); j++) {
181                        dp.short_b = weights_short_b[j].value;
182
183                        // Average over longuest_edgeC distribution
184                        for(int k=0; k< (int)weights_long_c.size(); k++) {
185                                dp.long_c = weights_long_c[k].value;
186
187                                // Average over theta distribution
188                                for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
189                                dp.parallel_theta = weights_parallel_theta[l].value;
190
191                                        // Average over phi distribution
192                                        for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
193                                                dp.parallel_phi = weights_parallel_phi[m].value;
194
195                                                // Average over phi distribution
196                                                for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
197                                                        dp.parallel_psi = weights_parallel_psi[n].value;
198                                                        //Un-normalize by volume
199                                                        double _ptvalue = weights_short_a[i].weight
200                                                                * weights_short_b[j].weight
201                                                                * weights_long_c[k].weight
202                                                                * weights_parallel_theta[l].weight
203                                                                * weights_parallel_phi[m].weight
204                                                                * weights_parallel_psi[n].weight
205                                                                * parallelepiped_analytical_2DXY(&dp, qx, qy)
206                                                                * weights_short_a[i].value*weights_short_b[j].value
207                                                                * weights_long_c[k].value;
208
209                                                        if (weights_parallel_theta.size()>1) {
210                                                                _ptvalue *= fabs(sin(weights_parallel_theta[l].value*pi/180.0));
211                                                        }
212                                                        sum += _ptvalue;
213                                                        //Find average volume
214                                                        vol += weights_short_a[i].weight
215                                                                * weights_short_b[j].weight
216                                                                * weights_long_c[k].weight
217                                                                * weights_short_a[i].value*weights_short_b[j].value
218                                                                * weights_long_c[k].value;
219                                                        //Find norm for volume
220                                                        norm_vol += weights_short_a[i].weight
221                                                                * weights_short_b[j].weight
222                                                                * weights_long_c[k].weight;
223
224                                                        norm += weights_short_a[i].weight
225                                                                * weights_short_b[j].weight
226                                                                * weights_long_c[k].weight
227                                                                * weights_parallel_theta[l].weight
228                                                                * weights_parallel_phi[m].weight
229                                                                * weights_parallel_psi[n].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_parallel_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 ParallelepipedModel :: 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 ParallelepipedModel :: calculate_ER() {
267        ParallelepipedParameters dp;
268        dp.short_a   = short_a();
269        dp.short_b   = short_b();
270        dp.long_c  = long_c();
271        double rad_out = 0.0;
272        double pi = 4.0*atan(1.0);
273        double suf_rad = sqrt(dp.short_a*dp.short_b/pi);
274
275        // Perform the computation, with all weight points
276        double sum = 0.0;
277        double norm = 0.0;
278
279        // Get the dispersion points for the short_edgeA
280        vector<WeightPoint> weights_short_a;
281        short_a.get_weights(weights_short_a);
282
283        // Get the dispersion points for the longer_edgeB
284        vector<WeightPoint> weights_short_b;
285        short_b.get_weights(weights_short_b);
286
287        // Get angular averaging for the longuest_edgeC
288        vector<WeightPoint> weights_long_c;
289        long_c.get_weights(weights_long_c);
290
291        // Loop over radius weight points
292        for(int i=0; i< (int)weights_short_a.size(); i++) {
293                dp.short_a = weights_short_a[i].value;
294
295                // Loop over longer_edgeB weight points
296                for(int j=0; j< (int)weights_short_b.size(); j++) {
297                        dp.short_b = weights_short_b[j].value;
298
299                        // Average over longuest_edgeC distribution
300                        for(int k=0; k< (int)weights_long_c.size(); k++) {
301                                dp.long_c = weights_long_c[k].value;
302                                //Calculate surface averaged radius
303                                //This is rough approximation.
304                                suf_rad = sqrt(dp.short_a*dp.short_b/pi);
305
306                                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
307                                sum +=weights_short_a[i].weight* weights_short_b[j].weight
308                                        * weights_long_c[k].weight*DiamCyl(dp.long_c, suf_rad)/2.0;
309                                norm += weights_short_a[i].weight* weights_short_b[j].weight*weights_long_c[k].weight;
310                        }
311                }
312        }
313
314        if (norm != 0){
315                //return the averaged value
316                rad_out =  sum/norm;}
317        else{
318                //return normal value
319                //Note: output of "DiamCyl(length,radius)" is DIAMETER.
320                rad_out = DiamCyl(dp.long_c, suf_rad)/2.0;}
321        return rad_out;
322
323}
Note: See TracBrowser for help on using the repository browser.