source: sasview/src/sas/models/c_extension/c_models/coreshellbicelle.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: 12.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 "core_shell_bicelle.h"
27
28extern "C" {
29#include "libCylinder.h"
30#include "libStructureFactor.h"
31}
32
33typedef struct {
34  double scale;
35  double radius;
36  double rim_thick;
37  double face_thick;
38  double length;
39  double core_sld;
40  double rim_sld;
41  double face_sld;
42  double solvent_sld;
43  double background;
44  double axis_theta;
45  double axis_phi;
46} CoreShellBicelleParameters;
47
48
49/**
50 * Function to evaluate 2D scattering function
51 * @param pars: parameters of the core-shell Bicelle
52 * @param q: q-value
53 * @param q_x: q_x / q
54 * @param q_y: q_y / q
55 * @return: function value
56 */
57static double core_shell_bicelle_analytical_2D_scaled(CoreShellBicelleParameters *pars, double q, double q_x, double q_y) {
58  double cyl_x, cyl_y;//, cyl_z;
59  //double q_z;
60  double alpha, vol, cos_val;
61  double answer;
62  //convert angle degree to radian
63  double pi = 4.0*atan(1.0);
64  double theta = pars->axis_theta * pi/180.0;
65  double phi = pars->axis_phi * pi/180.0;
66
67  // Cylinder orientation
68  cyl_x = cos(theta) * cos(phi);
69  cyl_y = sin(theta);
70  //cyl_z = -cos(theta) * sin(phi);
71
72  // q vector
73  //q_z = 0;
74
75  // Compute the angle btw vector q and the
76  // axis of the cylinder
77  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
78
79  // The following test should always pass
80  if (fabs(cos_val)>1.0) {
81    printf("core_shell_bicelle_analytical_2D: Unexpected error: cos(alpha)=%g\n", cos_val);
82    return 0;
83  }
84
85  alpha = acos( cos_val );
86
87  // Call the IGOR library function to get the kernel
88  answer = BicelleKernel(q, pars->radius, pars->rim_thick, pars->face_thick,
89      pars->core_sld,pars->face_sld,pars->rim_sld,
90      pars->solvent_sld, pars->length/2.0, alpha) / fabs(sin(alpha));
91
92  //normalize by cylinder volume
93  vol=pi*(pars->radius+pars->rim_thick)
94          *(pars->radius+pars->rim_thick)
95          *(pars->length+2.0*pars->face_thick);
96  answer /= vol;
97
98  //convert to [cm-1]
99  answer *= 1.0e8;
100
101  //Scale
102  answer *= pars->scale;
103
104  // add in the background
105  answer += pars->background;
106
107  return answer;
108}
109
110/**
111 * Function to evaluate 2D scattering function
112 * @param pars: parameters of the core-shell cylinder
113 * @param q: q-value
114 * @return: function value
115 */
116static double core_shell_bicelle_analytical_2DXY(CoreShellBicelleParameters *pars, double qx, double qy) {
117  double q;
118  q = sqrt(qx*qx+qy*qy);
119  return core_shell_bicelle_analytical_2D_scaled(pars, q, qx/q, qy/q);
120}
121
122
123CoreShellBicelleModel :: CoreShellBicelleModel() {
124  scale      = Parameter(1.0);
125  radius     = Parameter(20.0, true);
126  radius.set_min(0.0);
127  face_thick  = Parameter(10.0, true);
128  face_thick.set_min(0.0);
129  rim_thick  = Parameter(10.0, true);
130  rim_thick.set_min(0.0);
131  length     = Parameter(400.0, true);
132  length.set_min(0.0);
133  core_sld   = Parameter(1.e-6);
134  face_sld  = Parameter(4.e-6);
135  rim_sld  = Parameter(4.e-6);
136  solvent_sld= Parameter(1.e-6);
137  background = Parameter(0.0);
138  axis_theta = Parameter(90.0, true);
139  axis_phi   = Parameter(0.0, true);
140}
141
142/**
143 * Function to evaluate 1D scattering function
144 * The NIST IGOR library is used for the actual calculation.
145 * @param q: q-value
146 * @return: function value
147 */
148double CoreShellBicelleModel :: operator()(double q) {
149  double dp[10];
150
151  dp[0] = scale();
152  dp[1] = radius();
153  dp[2] = rim_thick();
154  dp[3] = face_thick();
155  dp[4] = length();
156  dp[5] = core_sld();
157  dp[6] = face_sld();
158  dp[7] = rim_sld();
159  dp[8] = solvent_sld();
160  dp[9] = 0.0;
161  double pi = 4.0*atan(1.0);
162  // Get the dispersion points for the radius
163  vector<WeightPoint> weights_rad;
164  radius.get_weights(weights_rad);
165
166  // Get the dispersion points for the thickness
167  vector<WeightPoint> weights_rthick;
168  rim_thick.get_weights(weights_rthick);
169
170  // Get the dispersion points for the thickness
171  vector<WeightPoint> weights_fthick;
172  face_thick.get_weights(weights_fthick);
173
174  // Get the dispersion points for the length
175  vector<WeightPoint> weights_len;
176  length.get_weights(weights_len);
177
178  // Perform the computation, with all weight points
179  double sum = 0.0;
180  double norm = 0.0;
181  double vol = 0.0;
182
183  // Loop over radius weight points
184  for(size_t i=0; i<weights_rad.size(); i++) {
185    dp[1] = weights_rad[i].value;
186
187    // Loop over length weight points
188    for(size_t j=0; j<weights_len.size(); j++) {
189      dp[4] = weights_len[j].value;
190
191      // Loop over thickness weight points
192      for(size_t k=0; k<weights_rthick.size(); k++) {
193        dp[2] = weights_rthick[k].value;
194        for(size_t l=0; l<weights_fthick.size(); l++) {
195            dp[3] = weights_fthick[l].value;
196                        //Un-normalize by volume
197                        sum += weights_rad[i].weight
198                                * weights_len[j].weight
199                                * weights_rthick[k].weight
200                                * weights_fthick[l].weight
201                                * BicelleIntegration(q, dp[1],dp[2], dp[3],
202                                                dp[5],dp[6],dp[7],dp[8], dp[4]);
203                        //* pow(weights_rad[i].value+weights_rthick[k].value,2)
204                        //*(weights_len[j].value+2.0*weights_fthick[k].value);
205                        //Find average volume
206                        vol += weights_rad[i].weight
207                                * weights_len[j].weight
208                                * weights_rthick[k].weight
209                                * weights_fthick[l].weight
210                                * pi * pow(weights_rad[i].value+weights_rthick[k].value,2)
211                        *(weights_len[j].value+2.0*weights_fthick[l].value);
212                        norm += weights_rad[i].weight
213                                * weights_len[j].weight
214                                * weights_rthick[k].weight
215                                * weights_fthick[l].weight;
216        }
217      }
218    }
219  }
220
221  if (vol != 0.0 && norm != 0.0) {
222    //Re-normalize by avg volume
223    sum = sum/(vol/norm);}
224  //convert to [cm-1]
225  sum *= 1.0e8;
226  sum *= dp[0];
227  return sum/norm + background();
228}
229
230/**
231 * Function to evaluate 2D scattering function
232 * @param q_x: value of Q along x
233 * @param q_y: value of Q along y
234 * @return: function value
235 */
236double CoreShellBicelleModel :: operator()(double qx, double qy) {
237  CoreShellBicelleParameters dp;
238  // Fill parameter array
239  dp.scale      = scale();
240  dp.radius     = radius();
241  dp.rim_thick  = rim_thick();
242  dp.face_thick  = face_thick();
243  dp.length     = length();
244  dp.core_sld   = core_sld();
245  dp.rim_sld  = rim_sld();
246  dp.face_sld  = face_sld();
247  dp.solvent_sld= solvent_sld();
248  dp.background = 0.0;
249  dp.axis_theta = axis_theta();
250  dp.axis_phi   = axis_phi();
251
252  // Get the dispersion points for the radius
253  vector<WeightPoint> weights_rad;
254  radius.get_weights(weights_rad);
255
256  // Get the dispersion points for the thickness
257  vector<WeightPoint> weights_rthick;
258  rim_thick.get_weights(weights_rthick);
259
260  // Get the dispersion points for the thickness
261  vector<WeightPoint> weights_fthick;
262  face_thick.get_weights(weights_fthick);
263
264  // Get the dispersion points for the length
265  vector<WeightPoint> weights_len;
266  length.get_weights(weights_len);
267
268  // Get angular averaging for theta
269  vector<WeightPoint> weights_theta;
270  axis_theta.get_weights(weights_theta);
271
272  // Get angular averaging for phi
273  vector<WeightPoint> weights_phi;
274  axis_phi.get_weights(weights_phi);
275
276  // Perform the computation, with all weight points
277  double sum = 0.0;
278  double norm = 0.0;
279  double norm_vol = 0.0;
280  double vol = 0.0;
281  double pi = 4.0*atan(1.0);
282  // Loop over radius weight points
283  for(size_t i=0; i<weights_rad.size(); i++) {
284    dp.radius = weights_rad[i].value;
285
286    // Loop over length weight points
287    for(size_t j=0; j<weights_len.size(); j++) {
288      dp.length = weights_len[j].value;
289
290      // Loop over thickness weight points
291      for(size_t m=0; m<weights_rthick.size(); m++) {
292        dp.rim_thick = weights_rthick[m].value;
293
294        // Loop over thickness weight points
295        for(size_t n=0; n<weights_fthick.size(); n++) {
296          dp.face_thick = weights_fthick[n].value;
297
298                        // Average over theta distribution
299                        for(size_t k=0; k<weights_theta.size(); k++) {
300                          dp.axis_theta = weights_theta[k].value;
301
302                          // Average over phi distribution
303                          for(size_t l=0; l<weights_phi.size(); l++) {
304                                dp.axis_phi = weights_phi[l].value;
305                                //Un-normalize by volume
306                                double _ptvalue = weights_rad[i].weight
307                                        * weights_len[j].weight
308                                        * weights_rthick[m].weight
309                                        * weights_fthick[n].weight
310                                        * weights_theta[k].weight
311                                        * weights_phi[l].weight
312                                        * core_shell_bicelle_analytical_2DXY(&dp, qx, qy)
313                                * pow(weights_rad[i].value+weights_rthick[m].value,2)
314                                *(weights_len[j].value+2.0*weights_fthick[n].value);
315
316                                if (weights_theta.size()>1) {
317                                  _ptvalue *= fabs(cos(weights_theta[k].value*pi/180.0));
318                                }
319                                sum += _ptvalue;
320
321                                //Find average volume
322                                vol += weights_rad[i].weight
323                                        * weights_len[j].weight
324                                        * weights_rthick[m].weight
325                                        * weights_fthick[n].weight
326                                        * pow(weights_rad[i].value+weights_rthick[m].value,2)
327                                *(weights_len[j].value+2.0*weights_fthick[n].value);
328                                //Find norm for volume
329                                norm_vol += weights_rad[i].weight
330                                        * weights_len[j].weight
331                                        * weights_rthick[m].weight
332                                        * weights_fthick[n].weight;
333
334                                norm += weights_rad[i].weight
335                                        * weights_len[j].weight
336                                        * weights_rthick[m].weight
337                                        * weights_fthick[n].weight
338                                        * weights_theta[k].weight
339                                        * weights_phi[l].weight;
340
341                          }
342                        }
343        }
344      }
345    }
346  }
347  // Averaging in theta needs an extra normalization
348  // factor to account for the sin(theta) term in the
349  // integration (see documentation).
350  if (weights_theta.size()>1) norm = norm / asin(1.0);
351
352  if (vol != 0.0 && norm_vol != 0.0) {
353    //Re-normalize by avg volume
354    sum = sum/(vol/norm_vol);}
355
356  return sum/norm + background();
357}
358
359/**
360 * Function to evaluate 2D scattering function
361 * @param pars: parameters of the cylinder
362 * @param q: q-value
363 * @param phi: angle phi
364 * @return: function value
365 */
366double CoreShellBicelleModel :: evaluate_rphi(double q, double phi) {
367  double qx = q*cos(phi);
368  double qy = q*sin(phi);
369  return (*this).operator()(qx, qy);
370}
371/**
372 * Function to calculate effective radius
373 * @return: effective radius value
374 */
375double CoreShellBicelleModel :: calculate_ER() {
376  CoreShellBicelleParameters dp;
377
378  dp.radius     = radius();
379  dp.rim_thick  = rim_thick();
380  dp.face_thick  = face_thick();
381  dp.length     = length();
382  double rad_out = 0.0;
383
384  // Perform the computation, with all weight points
385  double sum = 0.0;
386  double norm = 0.0;
387
388  // Get the dispersion points for the length
389  vector<WeightPoint> weights_length;
390  length.get_weights(weights_length);
391
392  // Get the dispersion points for the thickness
393  vector<WeightPoint> weights_rthick;
394  rim_thick.get_weights(weights_rthick);
395
396  // Get the dispersion points for the thickness
397  vector<WeightPoint> weights_fthick;
398  face_thick.get_weights(weights_fthick);
399
400  // Get the dispersion points for the radius
401  vector<WeightPoint> weights_radius ;
402  radius.get_weights(weights_radius);
403
404  // Loop over major shell weight points
405  for(int i=0; i< (int)weights_length.size(); i++) {
406    dp.length = weights_length[i].value;
407    for(int j=0; j< (int)weights_rthick.size(); j++) {
408      dp.rim_thick = weights_rthick[j].value;
409      for(int l=0; l< (int)weights_fthick.size(); l++) {
410          dp.face_thick = weights_fthick[l].value;
411                  for(int k=0; k< (int)weights_radius.size(); k++) {
412                        dp.radius = weights_radius[k].value;
413                        //Note: output of "DiamCyl( )" is DIAMETER.
414                        sum +=weights_length[i].weight * weights_rthick[j].weight * weights_fthick[l].weight
415                                * weights_radius[k].weight*DiamCyl(dp.length+2.0*dp.face_thick,dp.radius+dp.rim_thick)/2.0;
416                        norm += weights_length[i].weight* weights_rthick[j].weight * weights_fthick[l].weight* weights_radius[k].weight;
417                  }
418      }
419    }
420  }
421  if (norm != 0){
422    //return the averaged value
423    rad_out =  sum/norm;}
424  else{
425    //return normal value
426    //Note: output of "DiamCyl()" is DIAMETER.
427    rad_out = DiamCyl(dp.length+2.0*dp.face_thick,dp.radius+dp.rim_thick)/2.0;}
428
429  return rad_out;
430}
431double CoreShellBicelleModel :: calculate_VR() {
432  return 1.0;
433}
Note: See TracBrowser for help on using the repository browser.