source: sasview/sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp @ 0a9871c

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