source: sasview/src/sas/models/c_extension/c_models/sc.cpp @ fa6a8d1

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 fa6a8d1 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: 9.4 KB
RevLine 
[230f479]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 "sc.h"
27
28extern "C" {
29#include "libSphere.h"
30}
31
32// Convenience structure
33typedef struct {
34  double scale;
35  double dnn;
36  double d_factor;
37  double radius;
38  double sldSph;
39  double sldSolv;
40  double background;
41  double theta;
42  double phi;
43  double psi;
44} SCParameters;
45
46/**
47 * Function to evaluate 2D scattering function
48 * @param pars: parameters of the SCCrystalModel
49 * @param q: q-value
50 * @param q_x: q_x / q
51 * @param q_y: q_y / q
52 * @return: function value
53 */
54static double sc_analytical_2D_scaled(SCParameters *pars, double q, double q_x, double q_y) {
55  double a3_x, a3_y, a2_x, a2_y, a1_x, a1_y; //, a3_z
56  double q_z;
57  double cos_val_a3, cos_val_a2, cos_val_a1;
58  double a1_dot_q, a2_dot_q,a3_dot_q;
59  double answer;
60  double Pi = 4.0*atan(1.0);
61  double aa, Da, qDa_2, latticeScale, Zq;
62
63  double dp[5];
64  //convert angle degree to radian
65  double theta = pars->theta * Pi/180.0;
66  double phi = pars->phi * Pi/180.0;
67  double psi = pars->psi * Pi/180.0;
68  dp[0] = 1.0;
69  dp[1] = pars->radius;
70  dp[2] = pars->sldSph;
71  dp[3] = pars->sldSolv;
72  dp[4] = 0.0;
73
74
75  aa = pars->dnn;
76  Da = pars->d_factor*aa;
77  qDa_2 = pow(q*Da,2.0);
78
79  latticeScale = (4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa,3.0);
80  /// Angles here are respect to detector coordinate instead of against q coordinate(PRB 36, 3, 1754)
81  // a3 axis orientation
82  a3_x = cos(theta) * cos(phi);
83  a3_y = sin(theta);
84  //a3_z = -cos(theta) * sin(phi);
85  // q vector
86  q_z = 0.0;
87
88  // Compute the angle btw vector q and the a3 axis
89  cos_val_a3 = a3_x*q_x + a3_y*q_y;// + a3_z*q_z;
90
91  // a1 axis orientation
92  a1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
93  a1_y = sin(psi)*cos(theta);
94
95  cos_val_a1 = a1_x*q_x + a1_y*q_y;
96
97  // a2 axis orientation
98  a2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
99  a2_y = cos(theta)*cos(psi);
100  // a2 axis
101  cos_val_a2 =  a2_x*q_x + a2_y*q_y;
102
103  // The following test should always pass
104  if (fabs(cos_val_a3)>1.0) {
105    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
106    cos_val_a3 = 1.0;
107  }
108   if (fabs(cos_val_a1)>1.0) {
109    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
110    cos_val_a1 = 1.0;
111  }
112   if (fabs(cos_val_a2)>1.0) {
113    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
114    cos_val_a3 = 1.0;
115  }
116  a3_dot_q = aa*q*cos_val_a3;
117  a1_dot_q = aa*q*cos_val_a1;//*sin(alpha);
118  a2_dot_q = aa*q*cos_val_a2;
119 
120  // Call Zq=Z1*Z2*Z3
121  Zq = (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a1_dot_q)+exp(-qDa_2));
122  Zq *= (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a2_dot_q)+exp(-qDa_2));
123  Zq *= (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a3_dot_q)+exp(-qDa_2));
124
125  // Use SphereForm directly from libigor
126  answer = SphereForm(dp,q)*Zq;
127
128  //consider scales
129  answer *= latticeScale * pars->scale;
130
131  // This FIXES a singualrity the kernel in libigor.
132  if ( answer == INFINITY || answer == NAN){
133    answer = 0.0;
134  }
135
136  // add background
137  answer += pars->background;
138
139  return answer;
140}
141
142/**
143 * Function to evaluate 2D scattering function
144 * @param pars: parameters of the SC_ParaCrystal
145 * @param q: q-value
146 * @return: function value
147 */
148static double sc_analytical_2DXY(SCParameters *pars, double qx, double qy){
149  double q;
150  q = sqrt(qx*qx+qy*qy);
151  return sc_analytical_2D_scaled(pars, q, qx/q, qy/q);
152}
153
154SCCrystalModel :: SCCrystalModel() {
155  scale      = Parameter(1.0);
156  dnn           = Parameter(220.0);
157  d_factor = Parameter(0.06);
158  radius     = Parameter(40.0, true);
159  radius.set_min(0.0);
160  sldSph   = Parameter(3.0e-6);
161  sldSolv   = Parameter(6.3e-6);
162  background = Parameter(0.0);
163  theta  = Parameter(0.0, true);
164  phi    = Parameter(0.0, true);
165  psi    = Parameter(0.0, true);
166}
167
168/**
169 * Function to evaluate 1D scattering function
170 * The NIST IGOR library is used for the actual calculation.
171 * @param q: q-value
172 * @return: function value
173 */
174double SCCrystalModel :: operator()(double q) {
175  double dp[7];
176
177  // Fill parameter array for IGOR library
178  // Add the background after averaging
179  dp[0] = scale();
180  dp[1] = dnn();
181  dp[2] = d_factor();
182  dp[3] = radius();
183  dp[4] = sldSph();
184  dp[5] = sldSolv();
185  dp[6] = 0.0;
186
187  // Get the dispersion points for the radius
188  vector<WeightPoint> weights_rad;
189  radius.get_weights(weights_rad);
190
191  // Perform the computation, with all weight points
192  double sum = 0.0;
193  double norm = 0.0;
194  double vol = 0.0;
195  double result;
196
197  // Loop over radius weight points
198  for(size_t i=0; i<weights_rad.size(); i++) {
199    dp[3] = weights_rad[i].value;
200
201    //Un-normalize SphereForm by volume
202    result = SC_ParaCrystal(dp, q) * pow(weights_rad[i].value,3);
203    // This FIXES a singualrity the kernel in libigor.
204    if ( result == INFINITY || result == NAN){
205      result = 0.0;
206    }
207    sum += weights_rad[i].weight
208        * result;
209    //Find average volume
210    vol += weights_rad[i].weight
211        * pow(weights_rad[i].value,3);
212
213    norm += weights_rad[i].weight;
214  }
215
216  if (vol != 0.0 && norm != 0.0) {
217    //Re-normalize by avg volume
218    sum = sum/(vol/norm);}
219  return sum/norm + background();
220}
221
222/**
223 * Function to evaluate 2D scattering function
224 * @param q_x: value of Q along x
225 * @param q_y: value of Q along y
226 * @return: function value
227 */
228double SCCrystalModel :: operator()(double qx, double qy) {
229  SCParameters dp;
230  dp.scale      = scale();
231  dp.dnn   = dnn();
232  dp.d_factor   = d_factor();
233  dp.radius  = radius();
234  dp.sldSph   = sldSph();
235  dp.sldSolv   = sldSolv();
236  dp.background = 0.0;
237  dp.theta  = theta();
238  dp.phi    = phi();
239  dp.psi    = psi();
240
241  // Get the dispersion points for the radius
242  vector<WeightPoint> weights_rad;
243  radius.get_weights(weights_rad);
244
245  // Get angular averaging for theta
246  vector<WeightPoint> weights_theta;
247  theta.get_weights(weights_theta);
248
249  // Get angular averaging for phi
250  vector<WeightPoint> weights_phi;
251  phi.get_weights(weights_phi);
252
253  // Get angular averaging for psi
254  vector<WeightPoint> weights_psi;
255  psi.get_weights(weights_psi);
256
257  // Perform the computation, with all weight points
258  double sum = 0.0;
259  double norm = 0.0;
260  double norm_vol = 0.0;
261  double vol = 0.0;
262  double pi = 4.0*atan(1.0);
263  // Loop over radius weight points
264  for(size_t i=0; i<weights_rad.size(); i++) {
265    dp.radius = weights_rad[i].value;
266    // Average over theta distribution
267    for(size_t j=0; j< weights_theta.size(); j++) {
268      dp.theta = weights_theta[j].value;
269      // Average over phi distribution
270      for(size_t k=0; k< weights_phi.size(); k++) {
271        dp.phi = weights_phi[k].value;
272        // Average over phi distribution
273        for(size_t l=0; l< weights_psi.size(); l++) {
274          dp.psi = weights_psi[l].value;
275          //Un-normalize SphereForm by volume
276          double _ptvalue = weights_rad[i].weight
277              * weights_theta[j].weight
278              * weights_phi[k].weight
279              * weights_psi[l].weight
280              * sc_analytical_2DXY(&dp, qx, qy);
281          //* pow(weights_rad[i].value,3.0);
282          // Consider when there is infinte or nan.
283          if ( _ptvalue == INFINITY || _ptvalue == NAN){
284            _ptvalue = 0.0;
285          }
286          if (weights_theta.size()>1) {
287            _ptvalue *= fabs(cos(weights_theta[j].value*pi/180.0));
288          }
289          sum += _ptvalue;
290          // This model dose not need the volume of spheres correction!!!
291          //Find average volume
292          //vol += weights_rad[i].weight
293          //    * pow(weights_rad[i].value,3);
294          //Find norm for volume
295          //norm_vol += weights_rad[i].weight;
296
297          norm += weights_rad[i].weight
298              * weights_theta[j].weight
299              * weights_phi[k].weight
300              * weights_psi[l].weight;
301        }
302      }
303    }
304  }
305  // Averaging in theta needs an extra normalization
306  // factor to account for the sin(theta) term in the
307  // integration (see documentation).
308  if (weights_theta.size()>1) norm = norm / asin(1.0);
309
310  if (vol != 0.0 && norm_vol != 0.0) {
311    //Re-normalize by avg volume
312    sum = sum/(vol/norm_vol);}
313
314  return sum/norm + background();
315}
316
317/**
318 * Function to evaluate 2D scattering function
319 * @param pars: parameters of the SCCrystal
320 * @param q: q-value
321 * @param phi: angle phi
322 * @return: function value
323 */
324double SCCrystalModel :: evaluate_rphi(double q, double phi) {
325  return (*this).operator()(q);
326}
327
328/**
329 * Function to calculate effective radius
330 * @return: effective radius value
331 */
332double SCCrystalModel :: calculate_ER() {
333  //NOT implemented yet!!!
334  return 0.0;
335}
336double SCCrystalModel :: calculate_VR() {
337  return 1.0;
338}
Note: See TracBrowser for help on using the repository browser.