source: sasview/src/sas/models/c_extension/c_models/spheroidXT.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: 14.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// RKH 03Apr2014 a re-parametrised CoreShellEllipsoid with core axial ratio X and shell thickness T
22#include <math.h>
23#include "parameters.hh"
24#include <stdio.h>
25#include <stdlib.h>
26using namespace std;
27#include "spheroidXT.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 X_core;
38  double T_shell;
39  double XpolarShell;
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} SpheroidXTParameters;
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 spheroidXT_analytical_2D_scaled(SpheroidXTParameters *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 // was    answer = gfn4(cos_val,pars->equat_core,pars->polar_core,pars->equat_shell,pars->polar_shell,sldcs,sldss,q);
98  answer = gfn4(cos_val,pars->equat_core,  pars->equat_core*pars->X_core,   pars->equat_core + pars->T_shell,   pars->equat_core*pars->X_core + pars->T_shell*pars->XpolarShell ,sldcs,sldss,q);
99  //It seems that it should be normalized somehow. How???
100
101  //normalize by cylinder volume
102  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
103  // was vol = 4.0*Pi/3.0*pars->equat_shell*pars->equat_shell*pars->polar_shell;
104  vol = 4.0*Pi/3.0*(pars->equat_core + pars->T_shell)* (pars->equat_core + pars->T_shell) * ( pars->equat_core*pars->X_core + pars->T_shell*pars->XpolarShell);
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
119CoreShellEllipsoidXTModel :: CoreShellEllipsoidXTModel() {
120  scale      = Parameter(1.0);
121  equat_core     = Parameter(20.0, true);
122  equat_core.set_min(0.0);
123  X_core     = Parameter(3.0, true);
124  X_core.set_min(0.001);
125  T_shell   = Parameter(30.0, true);
126  T_shell.set_min(0.0);
127  XpolarShell    = Parameter(1.0, true);
128  XpolarShell.set_min(0.0);
129  sld_core   = Parameter(2e-6);
130  sld_shell  = Parameter(1e-6);
131  sld_solvent = Parameter(6.3e-6);
132  background = Parameter(0.0);
133  axis_theta  = Parameter(0.0, true);
134  axis_phi    = Parameter(0.0, true);
135
136}
137
138
139/**
140 * Function to evaluate 2D scattering function
141 * @param pars: parameters of the prolate
142 * @param q: q-value
143 * @return: function value
144 */
145
146static double spheroidXT_analytical_2DXY(SpheroidXTParameters *pars, double qx, double qy) {
147  double q;
148  q = sqrt(qx*qx+qy*qy);
149  return spheroidXT_analytical_2D_scaled(pars, q, qx/q, qy/q);
150}
151
152/**
153 * Function to evaluate 1D scattering function
154 * The NIST IGOR library is used for the actual calculation.
155 * @param q: q-value
156 * @return: function value
157 */
158double CoreShellEllipsoidXTModel :: operator()(double q) {
159  double dp[9];
160
161  // Fill parameter array for IGOR library
162  // Add the background after averaging
163  dp[0] = scale();
164  dp[1] = equat_core();
165  //dp[2] = polar_core();
166  //dp[3] = equat_shell();
167  //dp[4] = polar_shell();
168  dp[2] = equat_core()*X_core();
169  dp[3] = equat_core() + T_shell();
170  dp[4] = equat_core()*X_core() + T_shell()*XpolarShell();
171  dp[5] = sld_core();
172  dp[6] = sld_shell();
173  dp[7] = sld_solvent();
174  dp[8] = 0.0;
175
176  // Get the dispersion points for the major core
177  vector<WeightPoint> weights_equat_core;
178  equat_core.get_weights(weights_equat_core);
179
180  // Get the dispersion points
181  vector<WeightPoint> weights_X_core;
182  X_core.get_weights(weights_X_core);
183
184  // Get the dispersion points
185  vector<WeightPoint> weights_T_shell;
186  T_shell.get_weights(weights_T_shell);
187
188  // Get the dispersion points
189  vector<WeightPoint> weights_XpolarShell;
190  XpolarShell.get_weights(weights_XpolarShell);
191
192
193  // Perform the computation, with all weight points
194  double sum = 0.0;
195  double norm = 0.0;
196  double vol = 0.0;
197  double wproduct = 0.0;
198
199  // Loop over equat core weight points
200  for(int i=0; i<(int)weights_equat_core.size(); i++) {
201    dp[1] = weights_equat_core[i].value;
202
203    // Loop over polar core weight points
204    for(int j=0; j<(int)weights_X_core.size(); j++) {
205      dp[2] = dp[1]*weights_X_core[j].value;
206
207      // Loop over equat outer weight points
208      for(int k=0; k<(int)weights_T_shell.size(); k++) {
209        dp[3] = dp[1] + weights_T_shell[k].value;
210
211        // Loop over polar outer weight points
212        for(int l=0; l<(int)weights_XpolarShell.size(); l++) {
213          dp[4] = dp[2] + weights_XpolarShell[l].value*weights_T_shell[k].value;
214
215          //Un-normalize  by volume
216          wproduct =weights_equat_core[i].weight* weights_X_core[j].weight * weights_T_shell[k].weight
217                  * weights_XpolarShell[l].weight;
218          sum +=  wproduct * OblateForm(dp, q) * dp[3]*dp[3]*dp[4];
219          //Find average volume
220          vol += wproduct * dp[3]*dp[3]*dp[4];
221          // was pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
222          norm += wproduct;
223        }
224      }
225    }
226  }
227  if (vol != 0.0 && norm != 0.0) {
228    //Re-normalize by avg volume
229    sum = sum/(vol/norm);}
230  return sum/norm + background();
231}
232
233/**
234 * Function to evaluate 2D scattering function
235 * @param q_x: value of Q along x
236 * @param q_y: value of Q along y
237 * @return: function value
238 */
239/*
240double OblateModel :: operator()(double qx, double qy) {
241        double q = sqrt(qx*qx + qy*qy);
242
243        return (*this).operator()(q);
244}
245 */
246
247/**
248 * Function to evaluate 2D scattering function
249 * @param pars: parameters of the oblate
250 * @param q: q-value
251 * @param phi: angle phi
252 * @return: function value
253 */
254double CoreShellEllipsoidXTModel :: evaluate_rphi(double q, double phi) {
255  double qx = q*cos(phi);
256  double qy = q*sin(phi);
257  return (*this).operator()(qx, qy);
258}
259
260/**
261 * Function to evaluate 2D scattering function
262 * @param q_x: value of Q along x
263 * @param q_y: value of Q along y
264 * @return: function value
265 */
266double CoreShellEllipsoidXTModel :: operator()(double qx, double qy) {
267  SpheroidXTParameters dp;
268  // Fill parameter array
269  dp.scale      = scale();
270  dp.equat_core = equat_core();
271  dp.X_core = X_core();
272  dp.T_shell = T_shell();
273  dp.XpolarShell = XpolarShell();
274  dp.sld_core = sld_core();
275  dp.sld_shell = sld_shell();
276  dp.sld_solvent = sld_solvent();
277  dp.background = 0.0;
278  dp.axis_theta = axis_theta();
279  dp.axis_phi = axis_phi();
280
281  // Get the dispersion points for the major core
282  vector<WeightPoint> weights_equat_core;
283  equat_core.get_weights(weights_equat_core);
284
285  // Get the dispersion points
286  vector<WeightPoint> weights_X_core;
287  X_core.get_weights(weights_X_core);
288
289  // Get the dispersion points
290  vector<WeightPoint> weights_T_shell;
291  T_shell.get_weights(weights_T_shell);
292
293  // Get the dispersion points
294  vector<WeightPoint> weights_XpolarShell;
295  XpolarShell.get_weights(weights_XpolarShell);
296
297
298  // Get angular averaging for theta
299  vector<WeightPoint> weights_theta;
300  axis_theta.get_weights(weights_theta);
301
302  // Get angular averaging for phi
303  vector<WeightPoint> weights_phi;
304  axis_phi.get_weights(weights_phi);
305
306  // Perform the computation, with all weight points
307  double sum = 0.0;
308  double norm = 0.0;
309  double norm_vol = 0.0;
310  double vol = 0.0;
311  double pi = 4.0*atan(1.0);
312  double equat_outer = 0.0;
313  double polar_outer = 0.0;
314
315  // Loop over major core weight points
316  for(int i=0; i< (int)weights_equat_core.size(); i++) {
317    dp.equat_core = weights_equat_core[i].value;
318
319    // Loop over minor core weight points
320    for(int j=0; j< (int)weights_X_core.size(); j++) {
321      dp.X_core = weights_X_core[j].value;
322
323      // Loop over equat outer weight points
324      for(int k=0; k< (int)weights_T_shell.size(); k++) {
325        dp.T_shell = weights_T_shell[k].value;
326        equat_outer = weights_equat_core[i].value + weights_T_shell[k].value;
327
328        // Loop over polar outer weight points
329        for(int l=0; l< (int)weights_XpolarShell.size(); l++) {
330          dp.XpolarShell = weights_XpolarShell[l].value;
331          polar_outer = weights_equat_core[i].value*weights_X_core[j].value + weights_T_shell[k].value*weights_XpolarShell[l].value;
332
333          // Average over theta distribution
334          for(int m=0; m< (int)weights_theta.size(); m++) {
335            dp.axis_theta = weights_theta[m].value;
336
337            // Average over phi distribution
338            for(int n=0; n< (int)weights_phi.size(); n++) {
339              dp.axis_phi = weights_phi[n].value;
340              //Un-normalize by volume
341               double _ptvalue = weights_equat_core[i].weight *weights_X_core[j].weight
342                  * weights_T_shell[k].weight * weights_XpolarShell[l].weight
343                  * weights_theta[m].weight
344                  * weights_phi[n].weight
345                  // rkh NOTE this passes the NEW parameters
346                  * spheroidXT_analytical_2DXY(&dp, qx, qy) * pow(equat_outer,2)*polar_outer;
347              if (weights_theta.size()>1) {
348                _ptvalue *= fabs(cos(weights_theta[m].value*pi/180.0));
349              }
350              sum += _ptvalue;
351              //Find average volume
352              // rkh had to change this, original weighted by outer shell volume weights only, see spheroid.cpp, which looks odd,
353              // (as has to assume that weights of other loops sume to unity) and here we need all four loops to get the outer size.
354              vol += weights_equat_core[i].weight *weights_X_core[j].weight
355                      * weights_T_shell[k].weight * weights_XpolarShell[l].weight
356                      * pow(equat_outer,2)*polar_outer;
357              //Find norm for volume
358              norm_vol += weights_equat_core[i].weight *weights_X_core[j].weight
359                      * weights_T_shell[k].weight * weights_XpolarShell[l].weight;
360
361              norm += weights_equat_core[i].weight *weights_X_core[j].weight
362                  * weights_T_shell[k].weight * weights_XpolarShell[l].weight
363                  * weights_theta[m].weight * weights_phi[n].weight;
364            }
365          }
366        }
367      }
368    }
369  }
370  // Averaging in theta needs an extra normalization
371  // factor to account for the sin(theta) term in the
372  // integration (see documentation).
373  if (weights_theta.size()>1) norm = norm / asin(1.0);
374
375  if (vol != 0.0 && norm_vol != 0.0) {
376    //Re-normalize by avg volume
377    sum = sum/(vol/norm_vol);}
378
379  return sum/norm + background();
380}
381
382/**
383 * Function to calculate effective radius
384 * @return: effective radius value
385 * rkh This now needs to integrate over all four variables as above not just two
386 */
387double CoreShellEllipsoidXTModel :: calculate_ER() {
388  SpheroidXTParameters dp;
389
390  dp.equat_core = equat_core();
391  dp.X_core = X_core();
392  dp.T_shell = T_shell();
393  dp.XpolarShell = XpolarShell();
394
395  double rad_out = 0.0;
396  double equat_outer = 0.0;
397  double polar_outer = 0.0;
398
399  // Perform the computation, with all weight points
400  double sum = 0.0;
401  double norm = 0.0;
402
403  // Get the dispersion points for the core
404  vector<WeightPoint> weights_equat_core;
405  equat_core.get_weights(weights_equat_core);
406
407  // Get the dispersion points
408  vector<WeightPoint> weights_X_core;
409  X_core.get_weights(weights_X_core);
410
411  // Get the dispersion points
412  vector<WeightPoint> weights_T_shell;
413  T_shell.get_weights(weights_T_shell);
414
415  // Get the dispersion points
416  vector<WeightPoint> weights_XpolarShell;
417  XpolarShell.get_weights(weights_XpolarShell);
418
419
420  // Loop over core weight points
421  for(int i=0; i< (int)weights_equat_core.size(); i++) {
422    dp.equat_core = weights_equat_core[i].value;
423    // Loop over weight points
424    for(int j=0; j< (int)weights_X_core.size(); j++) {
425      // Loop over weight points
426      for(int k=0; k< (int)weights_T_shell.size(); k++) {
427        equat_outer = weights_equat_core[i].value + weights_T_shell[k].value;
428        // Loop over polar outer weight points
429        for(int l=0; l< (int)weights_XpolarShell.size(); l++) {
430          polar_outer = weights_equat_core[i].value*weights_X_core[j].value + weights_T_shell[k].value*weights_XpolarShell[l].value;
431        //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
432          sum +=weights_equat_core[i].weight *weights_X_core[j].weight
433                  * weights_T_shell[k].weight * weights_XpolarShell[l].weight*DiamEllip(polar_outer,equat_outer)/2.0;
434          norm += weights_equat_core[i].weight *weights_X_core[j].weight
435                  * weights_T_shell[k].weight * weights_XpolarShell[l].weight;
436        }
437      }
438    }
439  }
440  if (norm != 0){
441    //return the averaged value
442    rad_out =  sum/norm;}
443  else{
444    //return normal value
445    //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
446    rad_out = DiamEllip(dp.equat_core + dp.T_shell, dp.equat_core * dp.X_core + dp.T_shell*dp.XpolarShell)/2.0;}
447
448  return rad_out;
449}
450double CoreShellEllipsoidXTModel :: calculate_VR() {
451  return 1.0;
452}
Note: See TracBrowser for help on using the repository browser.