source: sasview/src/sas/models/c_extension/c_models/stackeddisks.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: 12.2 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 *      TODO: refactor so that we pull in the old sansmodels.c_extensions
21 *      TODO: add 2d
22 */
23
24#include <math.h>
25#include "parameters.hh"
26#include <stdio.h>
27using namespace std;
28#include "stacked_disks.h"
29
30extern "C" {
31#include "libCylinder.h"
32#include "libStructureFactor.h"
33}
34
35typedef struct {
36  double scale;
37  double radius;
38  double core_thick;
39  double layer_thick;
40  double core_sld;
41  double layer_sld;
42  double solvent_sld;
43  double n_stacking;
44  double sigma_d;
45  double background;
46  double axis_theta;
47  double axis_phi;
48} StackedDisksParameters;
49
50
51/**
52 * Function to evaluate 2D scattering function
53 * @param pars: parameters of the staked disks
54 * @param q: q-value
55 * @param q_x: q_x / q
56 * @param q_y: q_y / q
57 * @return: function value
58 */
59static double stacked_disks_analytical_2D_scaled(StackedDisksParameters *pars, double q, double q_x, double q_y) {
60  double cyl_x, cyl_y;//, cyl_z;
61  //double q_z;
62  double alpha, vol, cos_val;
63  double d, dum, halfheight;
64  double answer;
65  double pi = 4.0*atan(1.0);
66  double theta = pars->axis_theta * pi/180.0;
67  double phi = pars->axis_phi * pi/180.0;
68
69
70  // parallelepiped orientation
71  cyl_x = cos(theta) * cos(phi);
72  cyl_y = sin(theta);
73
74  // q vector
75  //q_z = 0;
76
77  // Compute the angle btw vector q and the
78  // axis of the parallelepiped
79  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
80
81  // The following test should always pass
82  if (fabs(cos_val)>1.0) {
83    printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
84    return 0;
85  }
86
87  // Note: cos(alpha) = 0 and 1 will get an
88  // undefined value from Stackdisc_kern
89  alpha = acos( cos_val );
90
91  // Call the IGOR library function to get the kernel
92  d = 2 * pars->layer_thick + pars->core_thick;
93  halfheight = pars->core_thick/2.0;
94  dum =alpha ;
95  answer = Stackdisc_kern(q, pars->radius, pars->core_sld,pars->layer_sld,
96      pars->solvent_sld, halfheight, pars->layer_thick, dum, pars->sigma_d, d, pars->n_stacking)/sin(alpha);
97
98  // Multiply by contrast^2
99  //answer *= pars->contrast*pars->contrast;
100
101  //normalize by staked disks volume
102  vol = acos(-1.0) * pars->radius * pars->radius * d * pars->n_stacking;
103  answer /= vol;
104
105  //convert to [cm-1]
106  answer *= 1.0e8;
107
108  //Scale
109  answer *= pars->scale;
110
111  // add in the background
112  answer += pars->background;
113
114  return answer;
115}
116
117/**
118 * Function to evaluate 2D scattering function
119 * @param pars: parameters of the staked disks
120 * @param q: q-value
121 * @return: function value
122 */
123static double stacked_disks_analytical_2DXY(StackedDisksParameters *pars, double qx, double qy) {
124  double q;
125  q = sqrt(qx*qx+qy*qy);
126  return stacked_disks_analytical_2D_scaled(pars, q, qx/q, qy/q);
127}
128
129StackedDisksModel :: StackedDisksModel() {
130  scale      = Parameter(1.0);
131  radius     = Parameter(3000.0, true);
132  radius.set_min(0.0);
133  core_thick  = Parameter(10.0, true);
134  core_thick.set_min(0.0);
135  layer_thick     = Parameter(15.0);
136  layer_thick.set_min(0.0);
137  core_sld = Parameter(4.0e-6);
138  layer_sld  = Parameter(-4.0e-7);
139  solvent_sld  = Parameter(5.0e-6);
140  n_stacking   = Parameter(1);
141  sigma_d   = Parameter(0);
142  background = Parameter(0.001);
143  axis_theta  = Parameter(0.0, true);
144  axis_phi    = Parameter(0.0, true);
145}
146
147/**
148 * Function to evaluate 1D scattering function
149 * The NIST IGOR library is used for the actual calculation.
150 * @param q: q-value
151 * @return: function value
152 */
153double StackedDisksModel :: operator()(double q) {
154  double dp[10];
155
156  // Fill parameter array for IGOR library
157  // Add the background after averaging
158  dp[0] = scale();
159  dp[1] = radius();
160  dp[2] = core_thick();
161  dp[3] = layer_thick();
162  dp[4] = core_sld();
163  dp[5] = layer_sld();
164  dp[6] = solvent_sld();
165  dp[7] = n_stacking();
166  dp[8] = sigma_d();
167  dp[9] = 0.0;
168
169  // Get the dispersion points for the radius
170  vector<WeightPoint> weights_radius;
171  radius.get_weights(weights_radius);
172
173  // Get the dispersion points for the core_thick
174  vector<WeightPoint> weights_core_thick;
175  core_thick.get_weights(weights_core_thick);
176
177  // Get the dispersion points for the layer_thick
178  vector<WeightPoint> weights_layer_thick;
179  layer_thick.get_weights(weights_layer_thick);
180
181  // Perform the computation, with all weight points
182  double sum = 0.0;
183  double norm = 0.0;
184  double vol = 0.0;
185
186  // Loop over length weight points
187  for(int i=0; i< (int)weights_radius.size(); i++) {
188    dp[1] = weights_radius[i].value;
189
190    // Loop over radius weight points
191    for(int j=0; j< (int)weights_core_thick.size(); j++) {
192      dp[2] = weights_core_thick[j].value;
193
194      // Loop over thickness weight points
195      for(int k=0; k< (int)weights_layer_thick.size(); k++) {
196        dp[3] = weights_layer_thick[k].value;
197        //Un-normalize by volume
198        sum += weights_radius[i].weight
199            * weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q)
200        *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
201        //Find average volume
202        vol += weights_radius[i].weight
203            * weights_core_thick[j].weight * weights_layer_thick[k].weight
204            *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
205        norm += weights_radius[i].weight
206            * weights_core_thick[j].weight* weights_layer_thick[k].weight;
207      }
208    }
209  }
210  if (vol != 0.0 && norm != 0.0) {
211    //Re-normalize by avg volume
212    sum = sum/(vol/norm);}
213
214  return sum/norm + background();
215}
216
217/**
218 * Function to evaluate 2D scattering function
219 * @param q_x: value of Q along x
220 * @param q_y: value of Q along y
221 * @return: function value
222 */
223double StackedDisksModel :: operator()(double qx, double qy) {
224  StackedDisksParameters dp;
225  // Fill parameter array
226  dp.scale      = scale();
227  dp.core_thick    = core_thick();
228  dp.radius       = radius();
229  dp.layer_thick  = layer_thick();
230  dp.core_sld   = core_sld();
231  dp.layer_sld  = layer_sld();
232  dp.solvent_sld= solvent_sld();
233  dp.n_stacking   = n_stacking();
234  dp.sigma_d   = sigma_d();
235  dp.background = 0.0;
236  dp.axis_theta = axis_theta();
237  dp.axis_phi   = axis_phi();
238
239  // Get the dispersion points for the length
240  vector<WeightPoint> weights_core_thick;
241  core_thick.get_weights(weights_core_thick);
242
243  // Get the dispersion points for the radius
244  vector<WeightPoint> weights_radius;
245  radius.get_weights(weights_radius);
246
247  // Get the dispersion points for the thickness
248  vector<WeightPoint> weights_layer_thick;
249  layer_thick.get_weights(weights_layer_thick);
250
251  // Get angular averaging for theta
252  vector<WeightPoint> weights_theta;
253  axis_theta.get_weights(weights_theta);
254
255  // Get angular averaging for phi
256  vector<WeightPoint> weights_phi;
257  axis_phi.get_weights(weights_phi);
258
259  // Perform the computation, with all weight points
260  double sum = 0.0;
261  double norm = 0.0;
262  double norm_vol = 0.0;
263  double vol = 0.0;
264  double pi = 4.0*atan(1.0);
265 
266  // Loop over length weight points
267  for(int i=0; i< (int)weights_core_thick.size(); i++) {
268    dp.core_thick = weights_core_thick[i].value;
269
270    // Loop over radius weight points
271    for(int j=0; j< (int)weights_radius.size(); j++) {
272      dp.radius = weights_radius[j].value;
273
274      // Loop over thickness weight points
275      for(int k=0; k< (int)weights_layer_thick.size(); k++) {
276        dp.layer_thick = weights_layer_thick[k].value;
277
278        for(int l=0; l< (int)weights_theta.size(); l++) {
279          dp.axis_theta = weights_theta[l].value;
280
281          // Average over phi distribution
282          for(int m=0; m <(int)weights_phi.size(); m++) {
283            dp.axis_phi = weights_phi[m].value;
284
285            //Un-normalize by volume
286            double _ptvalue = weights_core_thick[i].weight
287                * weights_radius[j].weight
288                * weights_layer_thick[k].weight
289                * weights_theta[l].weight
290                * weights_phi[m].weight
291                * stacked_disks_analytical_2DXY(&dp, qx, qy)
292            *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
293            if (weights_theta.size()>1) {
294              _ptvalue *= fabs(cos(weights_theta[l].value*pi/180.0));
295            }
296            sum += _ptvalue;
297            //Find average volume
298            vol += weights_radius[j].weight
299                * weights_core_thick[i].weight * weights_layer_thick[k].weight
300                *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
301            //Find norm for volume
302            norm_vol += weights_radius[j].weight
303                * weights_core_thick[i].weight * weights_layer_thick[k].weight;
304
305            norm += weights_core_thick[i].weight
306                * weights_radius[j].weight
307                * weights_layer_thick[k].weight
308                * weights_theta[l].weight
309                * weights_phi[m].weight;
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  if (vol != 0.0 && norm_vol != 0.0) {
320    //Re-normalize by avg volume
321    sum = sum/(vol/norm_vol);}
322  return sum/norm + background();
323}
324
325/**
326 * Function to evaluate 2D scattering function
327 * @param pars: parameters of the triaxial ellipsoid
328 * @param q: q-value
329 * @param phi: angle phi
330 * @return: function value
331 */
332double StackedDisksModel :: evaluate_rphi(double q, double phi) {
333  double qx = q*cos(phi);
334  double qy = q*sin(phi);
335  return (*this).operator()(qx, qy);
336}
337/**
338 * Function to calculate effective radius
339 * @return: effective radius value
340 */
341double StackedDisksModel :: calculate_ER() {
342  StackedDisksParameters dp;
343
344  dp.core_thick    = core_thick();
345  dp.radius       = radius();
346  dp.layer_thick  = layer_thick();
347  dp.n_stacking   = n_stacking();
348
349  double rad_out = 0.0;
350  if (dp.n_stacking <= 0.0){
351    return rad_out;
352  }
353
354  // Perform the computation, with all weight points
355  double sum = 0.0;
356  double norm = 0.0;
357
358  // Get the dispersion points for the length
359  vector<WeightPoint> weights_core_thick;
360  core_thick.get_weights(weights_core_thick);
361
362  // Get the dispersion points for the radius
363  vector<WeightPoint> weights_radius;
364  radius.get_weights(weights_radius);
365
366  // Get the dispersion points for the thickness
367  vector<WeightPoint> weights_layer_thick;
368  layer_thick.get_weights(weights_layer_thick);
369
370  // Loop over major shell weight points
371  for(int i=0; i< (int)weights_core_thick.size(); i++) {
372    dp.core_thick = weights_core_thick[i].value;
373    for(int j=0; j< (int)weights_layer_thick.size(); j++) {
374      dp.layer_thick = weights_layer_thick[j].value;
375      for(int k=0; k< (int)weights_radius.size(); k++) {
376        dp.radius = weights_radius[k].value;
377        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
378        sum +=weights_core_thick[i].weight*weights_layer_thick[j].weight
379            * weights_radius[k].weight*DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;
380        norm += weights_core_thick[i].weight*weights_layer_thick[j].weight* weights_radius[k].weight;
381      }
382    }
383  }
384  if (norm != 0){
385    //return the averaged value
386    rad_out =  sum/norm;}
387  else{
388    //return normal value
389    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
390    rad_out = DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;}
391
392  return rad_out;
393}
394double StackedDisksModel :: calculate_VR() {
395  return 1.0;
396}
Note: See TracBrowser for help on using the repository browser.