source: sasview/src/sas/models/c_extension/c_models/fcc.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: 9.4 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>
25using namespace std;
26#include "fcc.h"
27
28extern "C" {
29#include "libSphere.h"
30}
31
32// Convenience parameter 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} FCParameters;
45
46/**
47 * Function to evaluate 2D scattering function
48 * @param pars: parameters of the FCCCrystalModel
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 fc_analytical_2D_scaled(FCParameters *pars, double q, double q_x, double q_y) {
55  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //, b3_z
56  double q_z;
57  double cos_val_b3, cos_val_b2, cos_val_b1;
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, Fkq, Fkq_2;
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  aa = pars->dnn;
75  Da = pars->d_factor*aa;
76  qDa_2 = pow(q*Da,2.0);
77  //contrast = pars->sldSph - pars->sldSolv;
78
79  latticeScale = 4.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa*sqrt(2.0),3.0);
80  // q vector
81  q_z = 0.0; // for SANS; assuming qz is negligible
82  /// Angles here are respect to detector coordinate
83  ///  instead of against q coordinate in PRB 36(46), 3(6), 1754(3854)
84  // b3 axis orientation
85  b3_x = cos(theta) * cos(phi);
86  b3_y = sin(theta);
87  //b3_z = -cos(theta) * sin(phi);
88  cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
89 
90  // b1 axis orientation
91  b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
92  b1_y = sin(psi)*cos(theta);
93  cos_val_b1 = b1_x*q_x + b1_y*q_y;
94 
95  // b2 axis orientation
96  b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
97  b2_y = cos(theta)*cos(psi);
98  cos_val_b2 =  b2_x*q_x + b2_y*q_y;
99 
100  // The following test should always pass
101  if (fabs(cos_val_b3)>1.0) {
102    //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
103    cos_val_b3 = 1.0;
104  }
105  if (fabs(cos_val_b2)>1.0) {
106    //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
107    cos_val_b2 = 1.0;
108  }
109  if (fabs(cos_val_b1)>1.0) {
110    //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
111    cos_val_b1 = 1.0;
112  }
113  // Compute the angle btw vector q and the a3 axis
114  a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1);
115
116  // a1 axis
117  a1_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b3);
118
119  // a2 axis
120  a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1);
121
122
123  // Get Fkq and Fkq_2
124  Fkq = exp(-0.5*pow(Da/aa,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
125  Fkq_2 = Fkq*Fkq;
126  // Call Zq=Z1*Z2*Z3
127  Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
128  Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
129  Zq = Zq * (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
130
131  // Use SphereForm directly from libigor
132  answer = SphereForm(dp,q)*Zq;
133
134  //consider scales
135  answer *= latticeScale * pars->scale;
136
137  // This FIXES a singualrity the kernel in libigor.
138  if ( answer == INFINITY || answer == NAN){
139    answer = 0.0;
140  }
141
142  // add background
143  answer += pars->background;
144
145  return answer;
146}
147
148/**
149 * Function to evaluate 2D scattering function
150 * @param pars: parameters of the FCC_ParaCrystal
151 * @param q: q-value
152 * @return: function value
153 */
154static double fc_analytical_2DXY(FCParameters *pars, double qx, double qy){
155  double q;
156  q = sqrt(qx*qx+qy*qy);
157  return fc_analytical_2D_scaled(pars, q, qx/q, qy/q);
158}
159
160FCCrystalModel :: FCCrystalModel() {
161  scale      = Parameter(1.0);
162  dnn           = Parameter(220.0);
163  d_factor = Parameter(0.06);
164  radius     = Parameter(40.0, true);
165  radius.set_min(0.0);
166  sldSph   = Parameter(3.0e-6);
167  sldSolv   = Parameter(6.3e-6);
168  background = Parameter(0.0);
169  theta  = Parameter(0.0, true);
170  phi    = Parameter(0.0, true);
171  psi    = Parameter(0.0, true);
172}
173
174/**
175 * Function to evaluate 1D scattering function
176 * The NIST IGOR library is used for the actual calculation.
177 * @param q: q-value
178 * @return: function value
179 */
180double FCCrystalModel :: operator()(double q) {
181  double dp[7];
182
183  // Fill parameter array for IGOR library
184  // Add the background after averaging
185  dp[0] = scale();
186  dp[1] = dnn();
187  dp[2] = d_factor();
188  dp[3] = radius();
189  dp[4] = sldSph();
190  dp[5] = sldSolv();
191  dp[6] = 0.0;
192
193  // Get the dispersion points for the radius
194  vector<WeightPoint> weights_rad;
195  radius.get_weights(weights_rad);
196
197  // Perform the computation, with all weight points
198  double sum = 0.0;
199  double norm = 0.0;
200  double vol = 0.0;
201  double result;
202
203  // Loop over radius weight points
204  for(size_t i=0; i<weights_rad.size(); i++) {
205    dp[3] = weights_rad[i].value;
206
207    //Un-normalize SphereForm by volume
208    result = FCC_ParaCrystal(dp, q) * pow(weights_rad[i].value,3);
209    // This FIXES a singualrity the kernel in libigor.
210    if ( result == INFINITY || result == NAN){
211      result = 0.0;
212    }
213    sum += weights_rad[i].weight
214        * result;
215    //Find average volume
216    vol += weights_rad[i].weight
217        * pow(weights_rad[i].value,3);
218
219    norm += weights_rad[i].weight;
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 */
234double FCCrystalModel :: operator()(double qx, double qy) {
235  FCParameters dp;
236  dp.scale      = scale();
237  dp.dnn   = dnn();
238  dp.d_factor   = d_factor();
239  dp.radius  = radius();
240  dp.sldSph   = sldSph();
241  dp.sldSolv   = sldSolv();
242  dp.background = 0.0;
243  dp.theta  = theta();
244  dp.phi    = phi();
245  dp.psi    = psi();
246
247  // Get the dispersion points for the radius
248  vector<WeightPoint> weights_rad;
249  radius.get_weights(weights_rad);
250
251  // Get angular averaging for theta
252  vector<WeightPoint> weights_theta;
253  theta.get_weights(weights_theta);
254
255  // Get angular averaging for phi
256  vector<WeightPoint> weights_phi;
257  phi.get_weights(weights_phi);
258
259  // Get angular averaging for psi
260  vector<WeightPoint> weights_psi;
261  psi.get_weights(weights_psi);
262
263  // Perform the computation, with all weight points
264  double sum = 0.0;
265  double norm = 0.0;
266  double norm_vol = 0.0;
267  double vol = 0.0;
268  double pi = 4.0*atan(1.0);
269  // Loop over radius weight points
270  for(size_t i=0; i<weights_rad.size(); i++) {
271    dp.radius = weights_rad[i].value;
272    // Average over theta distribution
273    for(size_t j=0; j< weights_theta.size(); j++) {
274      dp.theta = weights_theta[j].value;
275      // Average over phi distribution
276      for(size_t k=0; k< weights_phi.size(); k++) {
277        dp.phi = weights_phi[k].value;
278        // Average over phi distribution
279        for(size_t l=0; l< weights_psi.size(); l++) {
280          dp.psi = weights_psi[l].value;
281          //Un-normalize SphereForm by volume
282          double _ptvalue = weights_rad[i].weight
283              * weights_theta[j].weight
284              * weights_phi[k].weight
285              * weights_psi[l].weight
286              * fc_analytical_2DXY(&dp, qx, qy);
287          //* pow(weights_rad[i].value,3.0);
288          // Consider when there is infinte or nan.
289          if ( _ptvalue == INFINITY || _ptvalue == NAN){
290            _ptvalue = 0.0;
291          }
292          if (weights_theta.size()>1) {
293            _ptvalue *= fabs(cos(weights_theta[j].value*pi/180.0));
294          }
295          sum += _ptvalue;
296          // This model dose not need the volume of spheres correction!!!
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 FCCCrystal
320 * @param q: q-value
321 * @param phi: angle phi
322 * @return: function value
323 */
324double FCCrystalModel :: 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 FCCrystalModel :: calculate_ER() {
333  //NOT implemented yet!!!
334  return 0.0;
335}
336double FCCrystalModel :: calculate_VR() {
337  return 1.0;
338}
Note: See TracBrowser for help on using the repository browser.