source: sasview/sansmodels/src/c_models/stackeddisks.cpp @ f368fd9

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 f368fd9 was cf7653d3, checked in by Jae Cho <jhjcho@…>, 12 years ago

set missing another pi

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