source: sasview/src/sas/models/c_extension/c_models/barbell.cpp @ fc18690

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 fc18690 was db46d13, checked in by Paul Kienzle <pkienzle@…>, 9 years ago

fix 2D barbell model and add it to the fit page

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