source: sasview/src/sas/models/c_extension/c_models/barbell.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, 9 years ago

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

  • Property mode set to 100644
File size: 10.0 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
27extern "C" {
28        #include "libCylinder.h"
29  #include "GaussWeights.h"
30        #include "barbell.h"
31}
32
33BarBellModel :: BarBellModel() {
34        scale      = Parameter(1.0);
35        rad_bar         = Parameter(20.0);
36        rad_bar.set_min(0.0);
37        len_bar     = Parameter(400.0, true);
38        len_bar.set_min(0.0);
39        rad_bell = Parameter(40.0);
40        rad_bell.set_min(0.0);
41        sld_barbell   = Parameter(1.0e-6);
42        sld_solv   = Parameter(6.3e-6);
43        background = Parameter(0.0);
44        theta  = Parameter(0.0, true);
45        phi    = Parameter(0.0, true);
46}
47
48double bar2d_kernel(double dp[], double q, double alpha) {
49  int j;
50  double Pi;
51  double scale,contr,bkg,sldc,slds;
52  double len,rad,hDist,endRad;
53  int nordj=76;
54  double zi=alpha,yyy,answer;     //running tally of integration
55  double summj,vaj,vbj,zij;     //for the inner integration
56  double arg1,arg2,inner,be;
57
58
59  scale = dp[0];
60  rad = dp[1];
61  len = dp[2];
62  endRad = dp[3];
63  sldc = dp[4];
64  slds = dp[5];
65  bkg = dp[6];
66
67  hDist = sqrt(fabs(endRad*endRad-rad*rad));    //by definition for this model
68
69  contr = sldc-slds;
70
71  Pi = 4.0*atan(1.0);
72  vaj = -1.0*hDist/endRad;
73  vbj = 1.0;    //endpoints of inner integral
74
75  summj=0.0;
76
77  for(j=0;j<nordj;j++) {
78    //20 gauss points for the inner integral
79    zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;    //the "t" dummy
80    yyy = Gauss76Wt[j] * Dumb_kernel(dp,q,zij,zi);    //uses the same Kernel as the Dumbbell, here L>0
81    summj += yyy;
82  }
83  //now calculate the value of the inner integral
84  inner = (vbj-vaj)/2.0*summj;
85  inner *= 4.0*Pi*endRad*endRad*endRad;
86
87  //now calculate outer integrand
88  arg1 = q*len/2.0*cos(zi);
89  arg2 = q*rad*sin(zi);
90  yyy = inner;
91
92  if(arg2 == 0) {
93    be = 0.5;
94  } else {
95    be = NR_BessJ1(arg2)/arg2;
96  }
97
98  if(arg1 == 0.0) {   //limiting value of sinc(0) is 1; sinc is not defined in math.h
99    yyy += Pi*rad*rad*len*2.0*be;
100  } else {
101    yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
102  }
103  yyy *= yyy;   //sin(zi);
104  answer = yyy;
105
106
107  answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);   //divide by volume
108  answer *= 1.0e8;    //convert to cm^-1
109  answer *= contr*contr;
110  answer *= scale;
111  answer += bkg;
112
113  return answer;
114}
115/**
116 * Function to evaluate 1D scattering function
117 * The NIST IGOR library is used for the actual calculation.
118 * @param q: q-value
119 * @return: function value
120 */
121double BarBellModel :: operator()(double q) {
122        double dp[7];
123
124        // Fill parameter array for IGOR library
125        // Add the background after averaging
126        dp[0] = scale();
127        dp[1] = rad_bar();
128        dp[2] = len_bar();
129        dp[3] = rad_bell();
130        dp[4] = sld_barbell();
131        dp[5] = sld_solv();
132        dp[6] = 0.0;
133
134        // Get the dispersion points for the rad_bar
135        vector<WeightPoint> weights_rad_bar;
136        rad_bar.get_weights(weights_rad_bar);
137        // Get the dispersion points for the len_bar
138        vector<WeightPoint> weights_len_bar;
139        len_bar.get_weights(weights_len_bar);
140        // Get the dispersion points for the rad_bell
141        vector<WeightPoint> weights_rad_bell;
142        rad_bell.get_weights(weights_rad_bell);
143
144        // Perform the computation, with all weight points
145        double sum = 0.0;
146        double norm = 0.0;
147        double vol = 0.0;
148        double pi,hDist,result;
149        double vol_i = 0.0;
150        pi = 4.0*atan(1.0);
151        // Loop over radius weight points
152        for(size_t i=0; i<weights_rad_bar.size(); i++) {
153                dp[1] = weights_rad_bar[i].value;
154                for(size_t j=0; j<weights_len_bar.size(); j++) {
155                        dp[2] = weights_len_bar[j].value;
156                        for(size_t k=0; k<weights_rad_bell.size(); k++) {
157                                dp[3] = weights_rad_bell[k].value;
158
159                                //Un-normalize SphereForm by volume
160                                hDist = sqrt(fabs(dp[3]*dp[3]-dp[1]*dp[1]));
161                                vol_i = pi*dp[1]*dp[1]*dp[2]+2.0*pi*(2.0*dp[3]*dp[3]*dp[3]/3.0
162                                                                +dp[3]*dp[3]*hDist-hDist*hDist*hDist/3.0);
163                                result =  Barbell(dp, q) * vol_i;
164                                // This FIXES a singualrity the kernel in libigor.
165                                if ( result == INFINITY || result == NAN){
166                                        result = 0.0;
167                                }
168                                sum += weights_rad_bar[i].weight*weights_len_bar[j].weight*weights_rad_bell[k].weight
169                                        * result;
170                                //Find average volume
171                                vol += weights_rad_bar[i].weight*weights_len_bar[j].weight*weights_rad_bell[k].weight
172                                        * vol_i;
173
174                                norm += weights_rad_bar[i].weight*weights_len_bar[j].weight*weights_rad_bell[k].weight;
175                        }
176                }
177        }
178
179        if (vol != 0.0 && norm != 0.0) {
180                //Re-normalize by avg volume
181                sum = sum/(vol/norm);}
182        return sum/norm + background();
183}
184
185/**
186 * Function to evaluate 2D scattering function
187 * @param q_x: value of Q along x
188 * @param q_y: value of Q along y
189 * @return: function value
190 */
191double BarBellModel :: operator()(double qx, double qy) {
192  double dp[7];
193
194  // Fill parameter array for IGOR library
195  // Add the background after averaging
196  dp[0] = scale();
197  dp[1] = rad_bar();
198  dp[2] = len_bar();
199  dp[3] = rad_bell();
200  dp[4] = sld_barbell();
201  dp[5] = sld_solv();
202  dp[6] = 0.0;
203
204  double _theta = theta();
205  double _phi = phi();
206
207        // Get the dispersion points for the rad_bar
208        vector<WeightPoint> weights_rad_bar;
209        rad_bar.get_weights(weights_rad_bar);
210
211        // Get the dispersion points for the len_bar
212        vector<WeightPoint> weights_len_bar;
213        len_bar.get_weights(weights_len_bar);
214
215        // Get the dispersion points for the rad_bell
216        vector<WeightPoint> weights_rad_bell;
217        rad_bell.get_weights(weights_rad_bell);
218
219        // Get angular averaging for theta
220        vector<WeightPoint> weights_theta;
221        theta.get_weights(weights_theta);
222
223        // Get angular averaging for phi
224        vector<WeightPoint> weights_phi;
225        phi.get_weights(weights_phi);
226
227
228        // Perform the computation, with all weight points
229        double sum = 0.0;
230        double norm = 0.0;
231        double norm_vol = 0.0;
232        double vol = 0.0;
233        double pi,hDist;
234        double vol_i = 0.0;
235        pi = 4.0*atan(1.0);
236
237        // Loop over radius weight points
238        for(size_t i=0; i<weights_rad_bar.size(); i++) {
239                dp[1] = weights_rad_bar[i].value;
240                for(size_t j=0; j<weights_len_bar.size(); j++) {
241                        dp[2] = weights_len_bar[j].value;
242                        for(size_t k=0; k<weights_rad_bell.size(); k++) {
243                                dp[3] = weights_rad_bell[k].value;
244                                // Average over theta distribution
245                                for(size_t l=0; l< weights_theta.size(); l++) {
246                                        _theta = weights_theta[l].value;
247                                        // Average over phi distribution
248                                        for(size_t m=0; m< weights_phi.size(); m++) {
249                                                _phi = weights_phi[m].value;
250                                                //Un-normalize Form by volume
251                                                hDist = sqrt(fabs(dp[3]*dp[3]-dp[1]*dp[1]));
252                                                vol_i = pi*dp[1]*dp[1]*dp[2]+2.0*pi*(2.0*dp[3]*dp[3]*dp[3]/3.0
253                                                                                +dp[3]*dp[3]*hDist-hDist*hDist*hDist/3.0);
254
255                                          const double q = sqrt(qx*qx+qy*qy);
256                                          //convert angle degree to radian
257                                          const double pi = 4.0*atan(1.0);
258
259                                          // Cylinder orientation
260                                    const double cyl_x = cos(_theta * pi/180.0) * cos(_phi * pi/180.0);
261                                    const double cyl_y = sin(_theta * pi/180.0);
262
263                                    // Compute the angle btw vector q and the
264                                    // axis of the cylinder (assume qz = 0)
265                                    const double cos_val = cyl_x*qx + cyl_y*qy;
266
267                                    // The following test should always pass
268                                    if (fabs(cos_val)>1.0) {
269                                      return 0;
270                                    }
271
272                                    // Note: cos(alpha) = 0 and 1 will get an
273                                    // undefined value from CylKernel
274                                    const double alpha = acos( cos_val );
275
276            // Call the IGOR library function to get the kernel
277            const double output = bar2d_kernel(dp, q, alpha)/sin(alpha);
278
279                                                double _ptvalue = weights_rad_bar[i].weight
280                                                                                        * weights_len_bar[j].weight
281                                                                                        * weights_rad_bell[k].weight
282                                                                                        * weights_theta[l].weight
283                                                                                        * weights_phi[m].weight
284                                                                                        * vol_i
285                                                                                        * output;
286                                                                                        //* pow(weights_rad[i].value,3.0);
287
288                                                // Consider when there is infinte or nan.
289                                                if ( _ptvalue == INFINITY || _ptvalue == NAN){
290                                                        _ptvalue = 0.0;
291                                                }
292                                                if (weights_theta.size()>1) {
293                                                        _ptvalue *= fabs(cos(weights_theta[l].value*pi/180.0));
294                                                }
295                                                sum += _ptvalue;
296                                                // This model dose not need the volume of spheres correction!!!
297                                                //Find average volume
298                                                vol += weights_rad_bar[i].weight
299                                                                * weights_len_bar[j].weight
300                                                                * weights_rad_bell[k].weight
301                                                                * vol_i;
302                                                //Find norm for volume
303                                                norm_vol += weights_rad_bar[i].weight
304                                                                * weights_len_bar[j].weight
305                                                                * weights_rad_bell[k].weight;
306
307                                                norm += weights_rad_bar[i].weight
308                                                                * weights_len_bar[j].weight
309                                                                * weights_rad_bell[k].weight
310                                                                * weights_theta[l].weight
311                                                                * weights_phi[m].weight;
312                                        }
313                                }
314                        }
315                }
316        }
317        // Averaging in theta needs an extra normalization
318        // factor to account for the sin(theta) term in the
319        // integration (see documentation).
320        if (weights_theta.size()>1) norm = norm / asin(1.0);
321
322        if (vol != 0.0 && norm_vol != 0.0) {
323                //Re-normalize by avg volume
324                sum = sum/(vol/norm_vol);}
325
326        return sum/norm + background();
327}
328
329/**
330 * Function to evaluate 2D scattering function
331 * @param pars: parameters of the SCCrystal
332 * @param q: q-value
333 * @param phi: angle phi
334 * @return: function value
335 */
336double BarBellModel :: evaluate_rphi(double q, double phi) {
337        return (*this).operator()(q);
338}
339
340/**
341 * Function to calculate effective radius
342 * @return: effective radius value
343 */
344double BarBellModel :: calculate_ER() {
345        //NOT implemented yet!!!
346        return 0.0;
347}
348/**
349 * Function to calculate volf_ratio for shell/tot
350 * @return: volf_ratio value
351 */
352double BarBellModel :: calculate_VR() {
353  return 1.0;
354}
Note: See TracBrowser for help on using the repository browser.