source: sasview/src/sas/models/c_extension/c_models/bcc.cpp @ 79492222

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