source: sasview/src/sas/models/c_extension/c_models/ellipsoid.cpp @ 79492222

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 79492222 was 79492222, checked in by krzywon, 8 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 9.6 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 */
21
22#include <math.h>
23#include "parameters.hh"
24#include <stdio.h>
25using namespace std;
26#include "ellipsoid.h"
27
28extern "C" {
29#include "libCylinder.h"
30#include "libStructureFactor.h"
31}
32
33typedef struct {
34  double scale;
35  double radius_a;
36  double radius_b;
37  double sldEll;
38  double sldSolv;
39  double background;
40  double axis_theta;
41  double axis_phi;
42} EllipsoidParameters;
43
44/**
45 * Function to evaluate 2D scattering function
46 * @param pars: parameters of the ellipsoid
47 * @param q: q-value
48 * @param q_x: q_x / q
49 * @param q_y: q_y / q
50 * @return: function value
51 */
52static double ellipsoid_analytical_2D_scaled(EllipsoidParameters *pars, double q, double q_x, double q_y) {
53  double cyl_x, cyl_y;//, cyl_z;
54  //double q_z;
55  double alpha, vol, cos_val;
56  double answer;
57  //convert angle degree to radian
58  double pi = 4.0*atan(1.0);
59  double theta = pars->axis_theta * pi/180.0;
60  double phi = pars->axis_phi * pi/180.0;
61
62  // Ellipsoid orientation
63  cyl_x = cos(theta) * cos(phi);
64  cyl_y = sin(theta);
65  //cyl_z = -cos(theta) * sin(phi);
66
67  // q vector
68  //q_z = 0.0;
69
70  // Compute the angle btw vector q and the
71  // axis of the cylinder
72  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
73
74  // The following test should always pass
75  if (fabs(cos_val)>1.0) {
76    printf("ellipsoid_ana_2D: Unexpected error: cos(alpha)>1\n");
77    return 0;
78  }
79
80  // Angle between rotation axis and q vector
81  alpha = acos( cos_val );
82
83  // Call the IGOR library function to get the kernel
84  answer = EllipsoidKernel(q, pars->radius_b, pars->radius_a, cos_val);
85
86  // Multiply by contrast^2
87  answer *= (pars->sldEll - pars->sldSolv) * (pars->sldEll - pars->sldSolv);
88
89  //normalize by cylinder volume
90  vol = 4.0/3.0 * acos(-1.0) * pars->radius_b * pars->radius_b * pars->radius_a;
91  answer *= vol;
92
93  //convert to [cm-1]
94  answer *= 1.0e8;
95
96  //Scale
97  answer *= pars->scale;
98
99  // add in the background
100  answer += pars->background;
101
102  return answer;
103}
104
105/**
106 * Function to evaluate 2D scattering function
107 * @param pars: parameters of the ellipsoid
108 * @param q: q-value
109 * @return: function value
110 */
111static double ellipsoid_analytical_2DXY(EllipsoidParameters *pars, double qx, double qy) {
112  double q;
113  q = sqrt(qx*qx+qy*qy);
114  return ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q);
115}
116
117EllipsoidModel :: EllipsoidModel() {
118  scale      = Parameter(1.0);
119  radius_a   = Parameter(20.0, true);
120  radius_a.set_min(0.0);
121  radius_b   = Parameter(400.0, true);
122  radius_b.set_min(0.0);
123  sldEll   = Parameter(4.e-6);
124  sldSolv   = Parameter(1.e-6);
125  background = Parameter(0.0);
126  axis_theta  = Parameter(57.325, true);
127  axis_phi    = Parameter(0.0, true);
128}
129
130/**
131 * Function to evaluate 1D scattering function
132 * The NIST IGOR library is used for the actual calculation.
133 * @param q: q-value
134 * @return: function value
135 */
136double EllipsoidModel :: operator()(double q) {
137  double dp[6];
138
139  // Fill parameter array for IGOR library
140  // Add the background after averaging
141  dp[0] = scale();
142  dp[1] = radius_a();
143  dp[2] = radius_b();
144  dp[3] = sldEll();
145  dp[4] = sldSolv();
146  dp[5] = 0.0;
147
148  // Get the dispersion points for the radius_a
149  vector<WeightPoint> weights_rad_a;
150  radius_a.get_weights(weights_rad_a);
151
152  // Get the dispersion points for the radius_b
153  vector<WeightPoint> weights_rad_b;
154  radius_b.get_weights(weights_rad_b);
155
156  // Perform the computation, with all weight points
157  double sum = 0.0;
158  double norm = 0.0;
159  double vol = 0.0;
160
161  // Loop over radius_a weight points
162  for(size_t i=0; i<weights_rad_a.size(); i++) {
163    dp[1] = weights_rad_a[i].value;
164
165    // Loop over radius_b weight points
166    for(size_t j=0; j<weights_rad_b.size(); j++) {
167      dp[2] = weights_rad_b[j].value;
168      //Un-normalize  by volume
169      sum += weights_rad_a[i].weight
170          * weights_rad_b[j].weight * EllipsoidForm(dp, q)
171      * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
172
173      //Find average volume
174      vol += weights_rad_a[i].weight
175          * weights_rad_b[j].weight
176          * pow(weights_rad_b[j].value,2)
177      * weights_rad_a[i].value;
178      norm += weights_rad_a[i].weight
179          * weights_rad_b[j].weight;
180    }
181  }
182
183  if (vol != 0.0 && norm != 0.0) {
184    //Re-normalize by avg volume
185    sum = sum/(vol/norm);}
186
187  return sum/norm + background();
188}
189
190/**
191 * Function to evaluate 2D scattering function
192 * @param q_x: value of Q along x
193 * @param q_y: value of Q along y
194 * @return: function value
195 */
196double EllipsoidModel :: operator()(double qx, double qy) {
197  EllipsoidParameters dp;
198  // Fill parameter array
199  dp.scale      = scale();
200  dp.radius_a   = radius_a();
201  dp.radius_b   = radius_b();
202  dp.sldEll   = sldEll();
203  dp.sldSolv   = sldSolv();
204  dp.background = 0.0;
205  dp.axis_theta = axis_theta();
206  dp.axis_phi   = axis_phi();
207
208  // Get the dispersion points for the radius_a
209  vector<WeightPoint> weights_rad_a;
210  radius_a.get_weights(weights_rad_a);
211
212  // Get the dispersion points for the radius_b
213  vector<WeightPoint> weights_rad_b;
214  radius_b.get_weights(weights_rad_b);
215
216  // Get angular averaging for theta
217  vector<WeightPoint> weights_theta;
218  axis_theta.get_weights(weights_theta);
219
220  // Get angular averaging for phi
221  vector<WeightPoint> weights_phi;
222  axis_phi.get_weights(weights_phi);
223
224  // Perform the computation, with all weight points
225  double sum = 0.0;
226  double norm = 0.0;
227  double norm_vol = 0.0;
228  double vol = 0.0;
229  double pi = 4.0*atan(1.0);
230  // Loop over radius weight points
231  for(size_t i=0; i<weights_rad_a.size(); i++) {
232    dp.radius_a = weights_rad_a[i].value;
233
234
235    // Loop over length weight points
236    for(size_t j=0; j<weights_rad_b.size(); j++) {
237      dp.radius_b = weights_rad_b[j].value;
238
239      // Average over theta distribution
240      for(size_t k=0; k<weights_theta.size(); k++) {
241        dp.axis_theta = weights_theta[k].value;
242
243        // Average over phi distribution
244        for(size_t l=0; l<weights_phi.size(); l++) {
245          dp.axis_phi = weights_phi[l].value;
246          //Un-normalize by volume
247          double _ptvalue = weights_rad_a[i].weight
248              * weights_rad_b[j].weight
249              * weights_theta[k].weight
250              * weights_phi[l].weight
251              * ellipsoid_analytical_2DXY(&dp, qx, qy)
252          * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
253          if (weights_theta.size()>1) {
254            _ptvalue *= fabs(cos(weights_theta[k].value*pi/180.0));
255          }
256          sum += _ptvalue;
257          //Find average volume
258          vol += weights_rad_a[i].weight
259              * weights_rad_b[j].weight
260              * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
261          //Find norm for volume
262          norm_vol += weights_rad_a[i].weight
263              * weights_rad_b[j].weight;
264
265          norm += weights_rad_a[i].weight
266              * weights_rad_b[j].weight
267              * weights_theta[k].weight
268              * weights_phi[l].weight;
269
270        }
271      }
272    }
273  }
274  // Averaging in theta needs an extra normalization
275  // factor to account for the sin(theta) term in the
276  // integration (see documentation).
277  if (weights_theta.size()>1) norm = norm / asin(1.0);
278
279  if (vol != 0.0 && norm_vol != 0.0) {
280    //Re-normalize by avg volume
281    sum = sum/(vol/norm_vol);}
282
283  return sum/norm + background();
284}
285
286/**
287 * Function to evaluate 2D scattering function
288 * @param pars: parameters of the cylinder
289 * @param q: q-value
290 * @param phi: angle phi
291 * @return: function value
292 */
293double EllipsoidModel :: evaluate_rphi(double q, double phi) {
294  double qx = q*cos(phi);
295  double qy = q*sin(phi);
296  return (*this).operator()(qx, qy);
297}
298
299/**
300 * Function to calculate effective radius
301 * @return: effective radius value
302 */
303double EllipsoidModel :: calculate_ER() {
304  EllipsoidParameters dp;
305
306  dp.radius_a = radius_a();
307  dp.radius_b = radius_b();
308
309  double rad_out = 0.0;
310
311  // Perform the computation, with all weight points
312  double sum = 0.0;
313  double norm = 0.0;
314
315  // Get the dispersion points for the major shell
316  vector<WeightPoint> weights_radius_a;
317  radius_a.get_weights(weights_radius_a);
318
319  // Get the dispersion points for the minor shell
320  vector<WeightPoint> weights_radius_b;
321  radius_b.get_weights(weights_radius_b);
322
323  // Loop over major shell weight points
324  for(int i=0; i< (int)weights_radius_b.size(); i++) {
325    dp.radius_b = weights_radius_b[i].value;
326    for(int k=0; k< (int)weights_radius_a.size(); k++) {
327      dp.radius_a = weights_radius_a[k].value;
328      sum +=weights_radius_b[i].weight
329          * weights_radius_a[k].weight*DiamEllip(dp.radius_a,dp.radius_b)/2.0;
330      norm += weights_radius_b[i].weight* weights_radius_a[k].weight;
331    }
332  }
333  if (norm != 0){
334    //return the averaged value
335    rad_out =  sum/norm;}
336  else{
337    //return normal value
338    rad_out = DiamEllip(dp.radius_a,dp.radius_b)/2.0;}
339
340  return rad_out;
341}
342double EllipsoidModel :: calculate_VR() {
343  return 1.0;
344}
Note: See TracBrowser for help on using the repository browser.