source: sasview/src/sas/models/c_extension/c_models/spheroid.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.9 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>
25#include <stdlib.h>
26using namespace std;
27#include "spheroid.h"
28
29extern "C" {
30#include "libCylinder.h"
31#include "libStructureFactor.h"
32}
33
34typedef struct {
35  double scale;
36  double equat_core;
37  double polar_core;
38  double equat_shell;
39  double polar_shell;
40  double sld_core;
41  double sld_shell;
42  double sld_solvent;
43  double background;
44  double axis_theta;
45  double axis_phi;
46
47} SpheroidParameters;
48
49/**
50 * Function to evaluate 2D scattering function
51 * @param pars: parameters of the prolate
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 spheroid_analytical_2D_scaled(SpheroidParameters *pars, double q, double q_x, double q_y) {
58
59  double cyl_x, cyl_y;//, cyl_z;
60  //double q_z;
61  double alpha, vol, cos_val;
62  double answer;
63  double Pi = 4.0*atan(1.0);
64  double sldcs,sldss;
65
66  //convert angle degree to radian
67  double theta = pars->axis_theta * Pi/180.0;
68  double phi = pars->axis_phi * Pi/180.0;
69
70
71  // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
72  cyl_x = cos(theta) * cos(phi);
73  cyl_y = sin(theta);
74  //cyl_z = -cos(theta) * sin(phi);
75  //del sld
76  sldcs = pars->sld_core - pars->sld_shell;
77  sldss = pars->sld_shell- pars->sld_solvent;
78
79  // q vector
80  //q_z = 0;
81
82  // Compute the angle btw vector q and the
83  // axis of the cylinder
84  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
85
86  // The following test should always pass
87  if (fabs(cos_val)>1.0) {
88    printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
89    return 0;
90  }
91
92  // Note: cos(alpha) = 0 and 1 will get an
93  // undefined value from CylKernel
94  alpha = acos( cos_val );
95
96  // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params.
97  answer = gfn4(cos_val,pars->equat_core,pars->polar_core,pars->equat_shell,pars->polar_shell,sldcs,sldss,q);
98  //It seems that it should be normalized somehow. How???
99
100  //normalize by cylinder volume
101  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
102  vol = 4.0*Pi/3.0*pars->equat_shell*pars->equat_shell*pars->polar_shell;
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
117CoreShellEllipsoidModel :: CoreShellEllipsoidModel() {
118  scale      = Parameter(1.0);
119  equat_core     = Parameter(200.0, true);
120  equat_core.set_min(0.0);
121  polar_core     = Parameter(20.0, true);
122  polar_core.set_min(0.0);
123  equat_shell   = Parameter(250.0, true);
124  equat_shell.set_min(0.0);
125  polar_shell    = Parameter(30.0, true);
126  polar_shell.set_min(0.0);
127  sld_core   = Parameter(2e-6);
128  sld_shell  = Parameter(1e-6);
129  sld_solvent = Parameter(6.3e-6);
130  background = Parameter(0.0);
131  axis_theta  = Parameter(0.0, true);
132  axis_phi    = Parameter(0.0, true);
133
134}
135
136
137/**
138 * Function to evaluate 2D scattering function
139 * @param pars: parameters of the prolate
140 * @param q: q-value
141 * @return: function value
142 */
143static double spheroid_analytical_2DXY(SpheroidParameters *pars, double qx, double qy) {
144  double q;
145  q = sqrt(qx*qx+qy*qy);
146  return spheroid_analytical_2D_scaled(pars, q, qx/q, qy/q);
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 CoreShellEllipsoidModel :: operator()(double q) {
156  double dp[9];
157
158  // Fill parameter array for IGOR library
159  // Add the background after averaging
160  dp[0] = scale();
161  dp[1] = equat_core();
162  dp[2] = polar_core();
163  dp[3] = equat_shell();
164  dp[4] = polar_shell();
165  dp[5] = sld_core();
166  dp[6] = sld_shell();
167  dp[7] = sld_solvent();
168  dp[8] = 0.0;
169
170  // Get the dispersion points for the major core
171  vector<WeightPoint> weights_equat_core;
172  equat_core.get_weights(weights_equat_core);
173
174  // Get the dispersion points for the minor core
175  vector<WeightPoint> weights_polar_core;
176  polar_core.get_weights(weights_polar_core);
177
178  // Get the dispersion points for the major shell
179  vector<WeightPoint> weights_equat_shell;
180  equat_shell.get_weights(weights_equat_shell);
181
182  // Get the dispersion points for the minor_shell
183  vector<WeightPoint> weights_polar_shell;
184  polar_shell.get_weights(weights_polar_shell);
185
186
187  // Perform the computation, with all weight points
188  double sum = 0.0;
189  double norm = 0.0;
190  double vol = 0.0;
191
192  // Loop over major core weight points
193  for(int i=0; i<(int)weights_equat_core.size(); i++) {
194    dp[1] = weights_equat_core[i].value;
195
196    // Loop over minor core weight points
197    for(int j=0; j<(int)weights_polar_core.size(); j++) {
198      dp[2] = weights_polar_core[j].value;
199
200      // Loop over major shell weight points
201      for(int k=0; k<(int)weights_equat_shell.size(); k++) {
202        dp[3] = weights_equat_shell[k].value;
203
204        // Loop over minor shell weight points
205        for(int l=0; l<(int)weights_polar_shell.size(); l++) {
206          dp[4] = weights_polar_shell[l].value;
207          //Un-normalize  by volume
208          sum += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
209              * weights_polar_shell[l].weight * OblateForm(dp, q)
210          * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
211          //Find average volume
212          vol += weights_equat_core[i].weight* weights_polar_core[j].weight
213              * weights_equat_shell[k].weight
214              * weights_polar_shell[l].weight
215              * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
216          norm += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
217              * weights_polar_shell[l].weight;
218        }
219      }
220    }
221  }
222  if (vol != 0.0 && norm != 0.0) {
223    //Re-normalize by avg volume
224    sum = sum/(vol/norm);}
225  return sum/norm + background();
226}
227
228/**
229 * Function to evaluate 2D scattering function
230 * @param q_x: value of Q along x
231 * @param q_y: value of Q along y
232 * @return: function value
233 */
234/*
235double OblateModel :: operator()(double qx, double qy) {
236        double q = sqrt(qx*qx + qy*qy);
237
238        return (*this).operator()(q);
239}
240 */
241
242/**
243 * Function to evaluate 2D scattering function
244 * @param pars: parameters of the oblate
245 * @param q: q-value
246 * @param phi: angle phi
247 * @return: function value
248 */
249double CoreShellEllipsoidModel :: evaluate_rphi(double q, double phi) {
250  double qx = q*cos(phi);
251  double qy = q*sin(phi);
252  return (*this).operator()(qx, qy);
253}
254
255/**
256 * Function to evaluate 2D scattering function
257 * @param q_x: value of Q along x
258 * @param q_y: value of Q along y
259 * @return: function value
260 */
261double CoreShellEllipsoidModel :: operator()(double qx, double qy) {
262  SpheroidParameters dp;
263  // Fill parameter array
264  dp.scale      = scale();
265  dp.equat_core = equat_core();
266  dp.polar_core = polar_core();
267  dp.equat_shell = equat_shell();
268  dp.polar_shell = polar_shell();
269  dp.sld_core = sld_core();
270  dp.sld_shell = sld_shell();
271  dp.sld_solvent = sld_solvent();
272  dp.background = 0.0;
273  dp.axis_theta = axis_theta();
274  dp.axis_phi = axis_phi();
275
276  // Get the dispersion points for the major core
277  vector<WeightPoint> weights_equat_core;
278  equat_core.get_weights(weights_equat_core);
279
280  // Get the dispersion points for the minor core
281  vector<WeightPoint> weights_polar_core;
282  polar_core.get_weights(weights_polar_core);
283
284  // Get the dispersion points for the major shell
285  vector<WeightPoint> weights_equat_shell;
286  equat_shell.get_weights(weights_equat_shell);
287
288  // Get the dispersion points for the minor shell
289  vector<WeightPoint> weights_polar_shell;
290  polar_shell.get_weights(weights_polar_shell);
291
292
293  // Get angular averaging for theta
294  vector<WeightPoint> weights_theta;
295  axis_theta.get_weights(weights_theta);
296
297  // Get angular averaging for phi
298  vector<WeightPoint> weights_phi;
299  axis_phi.get_weights(weights_phi);
300
301  // Perform the computation, with all weight points
302  double sum = 0.0;
303  double norm = 0.0;
304  double norm_vol = 0.0;
305  double vol = 0.0;
306  double pi = 4.0*atan(1.0);
307  // Loop over major core weight points
308  for(int i=0; i< (int)weights_equat_core.size(); i++) {
309    dp.equat_core = weights_equat_core[i].value;
310
311    // Loop over minor core weight points
312    for(int j=0; j< (int)weights_polar_core.size(); j++) {
313      dp.polar_core = weights_polar_core[j].value;
314
315      // Loop over major shell weight points
316      for(int k=0; k< (int)weights_equat_shell.size(); k++) {
317        dp.equat_shell = weights_equat_shell[i].value;
318
319        // Loop over minor shell weight points
320        for(int l=0; l< (int)weights_polar_shell.size(); l++) {
321          dp.polar_shell = weights_polar_shell[l].value;
322
323          // Average over theta distribution
324          for(int m=0; m< (int)weights_theta.size(); m++) {
325            dp.axis_theta = weights_theta[m].value;
326
327            // Average over phi distribution
328            for(int n=0; n< (int)weights_phi.size(); n++) {
329              dp.axis_phi = weights_phi[n].value;
330              //Un-normalize by volume
331              double _ptvalue = weights_equat_core[i].weight *weights_polar_core[j].weight
332                  * weights_equat_shell[k].weight * weights_polar_shell[l].weight
333                  * weights_theta[m].weight
334                  * weights_phi[n].weight
335                  * spheroid_analytical_2DXY(&dp, qx, qy)
336              * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
337              if (weights_theta.size()>1) {
338                _ptvalue *= fabs(cos(weights_theta[m].value*pi/180.0));
339              }
340              sum += _ptvalue;
341              //Find average volume
342              vol += weights_equat_shell[k].weight
343                  * weights_polar_shell[l].weight
344                  * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
345              //Find norm for volume
346              norm_vol += weights_equat_shell[k].weight
347                  * weights_polar_shell[l].weight;
348
349              norm += weights_equat_core[i].weight *weights_polar_core[j].weight
350                  * weights_equat_shell[k].weight * weights_polar_shell[l].weight
351                  * weights_theta[m].weight * weights_phi[n].weight;
352            }
353          }
354        }
355      }
356    }
357  }
358  // Averaging in theta needs an extra normalization
359  // factor to account for the sin(theta) term in the
360  // integration (see documentation).
361  if (weights_theta.size()>1) norm = norm / asin(1.0);
362
363  if (vol != 0.0 && norm_vol != 0.0) {
364    //Re-normalize by avg volume
365    sum = sum/(vol/norm_vol);}
366
367  return sum/norm + background();
368}
369
370/**
371 * Function to calculate effective radius
372 * @return: effective radius value
373 */
374double CoreShellEllipsoidModel :: calculate_ER() {
375  SpheroidParameters dp;
376
377  dp.equat_shell = equat_shell();
378  dp.polar_shell = polar_shell();
379
380  double rad_out = 0.0;
381
382  // Perform the computation, with all weight points
383  double sum = 0.0;
384  double norm = 0.0;
385
386  // Get the dispersion points for the major shell
387  vector<WeightPoint> weights_equat_shell;
388  equat_shell.get_weights(weights_equat_shell);
389
390  // Get the dispersion points for the minor shell
391  vector<WeightPoint> weights_polar_shell;
392  polar_shell.get_weights(weights_polar_shell);
393
394  // Loop over major shell weight points
395  for(int i=0; i< (int)weights_equat_shell.size(); i++) {
396    dp.equat_shell = weights_equat_shell[i].value;
397    for(int k=0; k< (int)weights_polar_shell.size(); k++) {
398      dp.polar_shell = weights_polar_shell[k].value;
399      //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
400      sum +=weights_equat_shell[i].weight
401          * weights_polar_shell[k].weight*DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;
402      norm += weights_equat_shell[i].weight* weights_polar_shell[k].weight;
403    }
404  }
405  if (norm != 0){
406    //return the averaged value
407    rad_out =  sum/norm;}
408  else{
409    //return normal value
410    //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
411    rad_out = DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;}
412
413  return rad_out;
414}
415double CoreShellEllipsoidModel :: calculate_VR() {
416  return 1.0;
417}
Note: See TracBrowser for help on using the repository browser.