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

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 d67fc8d was f9bf661, checked in by Jae Cho <jhjcho@…>, 15 years ago

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