source: sasview/src/sas/models/c_extension/c_models/hollowcylinder.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 79492222, checked in by krzywon, 10 years ago

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

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