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