source: sasview/sansmodels/src/sans/models/c_models/csparallelepiped.cpp @ b1a65b6

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 b1a65b6 was 4628e31, checked in by Jae Cho <jhjcho@…>, 14 years ago

changed the unit of angles into degrees
  • Property mode set to 100644
File size: 10.2 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 2010, 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
21#include <math.h>
22#include "models.hh"
23#include "parameters.hh"
24#include <stdio.h>
25using namespace std;
26
27extern "C" {
28        #include "libCylinder.h"
29        #include "libStructureFactor.h"
30        #include "csparallelepiped.h"
31}
32
33CSParallelepipedModel :: CSParallelepipedModel() {
34        scale      = Parameter(1.0);
35        shortA     = Parameter(35.0, true);
36        shortA.set_min(1.0);
37        midB     = Parameter(75.0, true);
38        midB.set_min(1.0);
39        longC    = Parameter(400.0, true);
40        longC.set_min(1.0);
41        rimA     = Parameter(10.0, true);
42        rimB     = Parameter(10.0, true);
43        rimC     = Parameter(10.0, true);
44        sld_rimA     = Parameter(2.0e-6, true);
45        sld_rimB     = Parameter(4.0e-6, true);
46        sld_rimC    = Parameter(2.0e-6, true);
47        sld_pcore   = Parameter(1.0e-6);
48        sld_solv   = Parameter(6.0e-6);
49        background = Parameter(0.06);
50        parallel_theta  = Parameter(0.0, true);
51        parallel_phi    = Parameter(0.0, true);
52        parallel_psi    = Parameter(0.0, true);
53}
54
55/**
56 * Function to evaluate 1D scattering function
57 * The NIST IGOR library is used for the actual calculation.
58 * @param q: q-value
59 * @return: function value
60 */
61double CSParallelepipedModel :: operator()(double q) {
62        double dp[13];
63
64        // Fill parameter array for IGOR library
65        // Add the background after averaging
66        dp[0] = scale();
67        dp[1] = shortA();
68        dp[2] = midB();
69        dp[3] = longC();
70        dp[4] = rimA();
71        dp[5] = rimB();
72        dp[6] = rimC();
73        dp[7] = sld_rimA();
74        dp[8] = sld_rimB();
75        dp[9] = sld_rimC();
76        dp[10] = sld_pcore();
77        dp[11] = sld_solv();
78        dp[12] = 0.0;
79
80        // Get the dispersion points for the short_edgeA
81        vector<WeightPoint> weights_shortA;
82        shortA.get_weights(weights_shortA);
83
84        // Get the dispersion points for the longer_edgeB
85        vector<WeightPoint> weights_midB;
86        midB.get_weights(weights_midB);
87
88        // Get the dispersion points for the longuest_edgeC
89        vector<WeightPoint> weights_longC;
90        longC.get_weights(weights_longC);
91
92
93
94        // Perform the computation, with all weight points
95        double sum = 0.0;
96        double norm = 0.0;
97        double vol = 0.0;
98
99        // Loop over short_edgeA weight points
100        for(int i=0; i< (int)weights_shortA.size(); i++) {
101                dp[1] = weights_shortA[i].value;
102
103                // Loop over longer_edgeB weight points
104                for(int j=0; j< (int)weights_midB.size(); j++) {
105                        dp[2] = weights_midB[j].value;
106
107                        // Loop over longuest_edgeC weight points
108                        for(int k=0; k< (int)weights_longC.size(); k++) {
109                                dp[3] = weights_longC[k].value;
110                                //Un-normalize  by volume
111                                sum += weights_shortA[i].weight * weights_midB[j].weight
112                                        * weights_longC[k].weight * CSParallelepiped(dp, q)
113                                        * weights_shortA[i].value*weights_midB[j].value
114                                        * weights_longC[k].value;
115                                //Find average volume
116                                vol += weights_shortA[i].weight * weights_midB[j].weight
117                                        * weights_longC[k].weight
118                                        * weights_shortA[i].value * weights_midB[j].value
119                                        * weights_longC[k].value;
120
121                                norm += weights_shortA[i].weight
122                                         * weights_midB[j].weight * weights_longC[k].weight;
123                        }
124                }
125        }
126        if (vol != 0.0 && norm != 0.0) {
127                //Re-normalize by avg volume
128                sum = sum/(vol/norm);}
129
130        return sum/norm + background();
131}
132/**
133 * Function to evaluate 2D scattering function
134 * @param q_x: value of Q along x
135 * @param q_y: value of Q along y
136 * @return: function value
137 */
138double CSParallelepipedModel :: operator()(double qx, double qy) {
139        CSParallelepipedParameters dp;
140        // Fill parameter array
141        dp.scale      = scale();
142        dp.shortA   = shortA();
143        dp.midB   = midB();
144        dp.longC  = longC();
145        dp.rimA   = rimA();
146        dp.rimB   = rimB();
147        dp.rimC  = rimC();
148        dp.sld_rimA   = sld_rimA();
149        dp.sld_rimB   = sld_rimB();
150        dp.sld_rimC  = sld_rimC();
151        dp.sld_pcore   = sld_pcore();
152        dp.sld_solv   = sld_solv();
153        dp.background = 0.0;
154        //dp.background = background();
155        dp.parallel_theta  = parallel_theta();
156        dp.parallel_phi    = parallel_phi();
157        dp.parallel_psi    = parallel_psi();
158
159
160
161        // Get the dispersion points for the short_edgeA
162        vector<WeightPoint> weights_shortA;
163        shortA.get_weights(weights_shortA);
164
165        // Get the dispersion points for the longer_edgeB
166        vector<WeightPoint> weights_midB;
167        midB.get_weights(weights_midB);
168
169        // Get the dispersion points for the longuest_edgeC
170        vector<WeightPoint> weights_longC;
171        longC.get_weights(weights_longC);
172
173        // Get angular averaging for theta
174        vector<WeightPoint> weights_parallel_theta;
175        parallel_theta.get_weights(weights_parallel_theta);
176
177        // Get angular averaging for phi
178        vector<WeightPoint> weights_parallel_phi;
179        parallel_phi.get_weights(weights_parallel_phi);
180
181        // Get angular averaging for psi
182        vector<WeightPoint> weights_parallel_psi;
183        parallel_psi.get_weights(weights_parallel_psi);
184
185        // Perform the computation, with all weight points
186        double sum = 0.0;
187        double norm = 0.0;
188        double norm_vol = 0.0;
189        double vol = 0.0;
190        double pi = 4.0*atan(1.0);
191
192        // Loop over radius weight points
193        for(int i=0; i< (int)weights_shortA.size(); i++) {
194                dp.shortA = weights_shortA[i].value;
195
196                // Loop over longer_edgeB weight points
197                for(int j=0; j< (int)weights_midB.size(); j++) {
198                        dp.midB = weights_midB[j].value;
199
200                        // Average over longuest_edgeC distribution
201                        for(int k=0; k< (int)weights_longC.size(); k++) {
202                                dp.longC = weights_longC[k].value;
203
204                                // Average over theta distribution
205                                for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
206                                dp.parallel_theta = weights_parallel_theta[l].value;
207
208                                        // Average over phi distribution
209                                        for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
210                                                dp.parallel_phi = weights_parallel_phi[m].value;
211
212                                                // Average over phi distribution
213                                                for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
214                                                        dp.parallel_psi = weights_parallel_psi[n].value;
215                                                        //Un-normalize by volume
216                                                        double _ptvalue = weights_shortA[i].weight
217                                                                * weights_midB[j].weight
218                                                                * weights_longC[k].weight
219                                                                * weights_parallel_theta[l].weight
220                                                                * weights_parallel_phi[m].weight
221                                                                * weights_parallel_psi[n].weight
222                                                                * csparallelepiped_analytical_2DXY(&dp, qx, qy)
223                                                                * weights_shortA[i].value*weights_midB[j].value
224                                                                * weights_longC[k].value;
225
226                                                        if (weights_parallel_theta.size()>1) {
227                                                                _ptvalue *= fabs(sin(weights_parallel_theta[l].value*pi/180.0));
228                                                        }
229                                                        sum += _ptvalue;
230                                                        //Find average volume
231                                                        vol += weights_shortA[i].weight
232                                                                * weights_midB[j].weight
233                                                                * weights_longC[k].weight
234                                                                * weights_shortA[i].value*weights_midB[j].value
235                                                                * weights_longC[k].value;
236                                                        //Find norm for volume
237                                                        norm_vol += weights_shortA[i].weight
238                                                                * weights_midB[j].weight
239                                                                * weights_longC[k].weight;
240
241                                                        norm += weights_shortA[i].weight
242                                                                * weights_midB[j].weight
243                                                                * weights_longC[k].weight
244                                                                * weights_parallel_theta[l].weight
245                                                                * weights_parallel_phi[m].weight
246                                                                * weights_parallel_psi[n].weight;
247                                                }
248                                        }
249
250                                }
251                        }
252                }
253        }
254        // Averaging in theta needs an extra normalization
255        // factor to account for the sin(theta) term in the
256        // integration (see documentation).
257        if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
258
259        if (vol != 0.0 && norm_vol != 0.0) {
260                //Re-normalize by avg volume
261                sum = sum/(vol/norm_vol);}
262
263        return sum/norm + background();
264}
265
266
267/**
268 * Function to evaluate 2D scattering function
269 * @param pars: parameters of the cylinder
270 * @param q: q-value
271 * @param phi: angle phi
272 * @return: function value
273 */
274double CSParallelepipedModel :: evaluate_rphi(double q, double phi) {
275        double qx = q*cos(phi);
276        double qy = q*sin(phi);
277        return (*this).operator()(qx, qy);
278}
279/**
280 * Function to calculate effective radius
281 * @return: effective radius value
282 */
283double CSParallelepipedModel :: calculate_ER() {
284        CSParallelepipedParameters dp;
285        dp.shortA   = shortA();
286        dp.midB   = midB();
287        dp.longC  = longC();
288        dp.rimA   = rimA();
289        dp.rimB   = rimB();
290        dp.rimC  = rimC();
291
292        double rad_out = 0.0;
293        double pi = 4.0*atan(1.0);
294        double suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
295        double height =(dp.longC + 2.0*dp.rimC);
296        // Perform the computation, with all weight points
297        double sum = 0.0;
298        double norm = 0.0;
299
300        // Get the dispersion points for the short_edgeA
301        vector<WeightPoint> weights_shortA;
302        shortA.get_weights(weights_shortA);
303
304        // Get the dispersion points for the longer_edgeB
305        vector<WeightPoint> weights_midB;
306        midB.get_weights(weights_midB);
307
308        // Get angular averaging for the longuest_edgeC
309        vector<WeightPoint> weights_longC;
310        longC.get_weights(weights_longC);
311
312        // Loop over radius weight points
313        for(int i=0; i< (int)weights_shortA.size(); i++) {
314                dp.shortA = weights_shortA[i].value;
315
316                // Loop over longer_edgeB weight points
317                for(int j=0; j< (int)weights_midB.size(); j++) {
318                        dp.midB = weights_midB[j].value;
319
320                        // Average over longuest_edgeC distribution
321                        for(int k=0; k< (int)weights_longC.size(); k++) {
322                                dp.longC = weights_longC[k].value;
323                                //Calculate surface averaged radius
324                                //This is rough approximation.
325                                suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
326                                height =(dp.longC + 2.0*dp.rimC);
327                                //Note: output of "DiamCyl(dp.length,dp.radius)" is a DIAMETER.
328                                sum +=weights_shortA[i].weight* weights_midB[j].weight
329                                        * weights_longC[k].weight*DiamCyl(height, suf_rad)/2.0;
330                                norm += weights_shortA[i].weight* weights_midB[j].weight*weights_longC[k].weight;
331                        }
332                }
333        }
334
335        if (norm != 0){
336                //return the averaged value
337                rad_out =  sum/norm;}
338        else{
339                //return normal value
340                //Note: output of "DiamCyl(length,radius)" is DIAMETER.
341                rad_out = DiamCyl(dp.longC, suf_rad)/2.0;}
342        return rad_out;
343
344}
Note: See TracBrowser for help on using the repository browser.