source: sasview/src/sas/models/c_extension/c_models/cylinder.cpp @ 7bcbeb7

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 7bcbeb7 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.1 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
27extern "C" {
28        #include "libCylinder.h"
29        #include "libStructureFactor.h"
30        #include "libmultifunc/libfunc.h"
31}
32#include "cylinder.h"
33
34// Convenience parameter structure
35typedef struct {
36    double scale;
37    double radius;
38    double length;
39    double sldCyl;
40    double sldSolv;
41    double background;
42    double cyl_theta;
43    double cyl_phi;
44    double M0_sld_cyl;
45    double M_theta_cyl;
46    double M_phi_cyl;
47    double M0_sld_solv;
48    double M_theta_solv;
49    double M_phi_solv;
50    double Up_frac_i;
51        double Up_frac_f;
52        double Up_theta;
53} CylinderParameters;
54
55CylinderModel :: CylinderModel() {
56        scale      = Parameter(1.0);
57        radius     = Parameter(20.0, true);
58        radius.set_min(0.0);
59        length     = Parameter(400.0, true);
60        length.set_min(0.0);
61        sldCyl   = Parameter(4.e-6);
62        sldSolv   = Parameter(1.e-6);
63        background = Parameter(0.0);
64        cyl_theta  = Parameter(0.0, true);
65        cyl_phi    = Parameter(0.0, true);
66        M0_sld_cyl = Parameter(0.0e-6);
67        M_theta_cyl = Parameter(0.0);
68        M_phi_cyl = Parameter(0.0); 
69        M0_sld_solv = Parameter(0.0e-6);
70        M_theta_solv = Parameter(0.0);
71        M_phi_solv = Parameter(0.0); 
72        Up_frac_i = Parameter(0.5); 
73        Up_frac_f = Parameter(0.5);
74        Up_theta = Parameter(0.0);
75}
76
77/**
78 * Function to evaluate 1D scattering function
79 * The NIST IGOR library is used for the actual calculation.
80 * @param q: q-value
81 * @return: function value
82 */
83double CylinderModel :: operator()(double q) {
84        double dp[6];
85
86        // Fill parameter array for IGOR library
87        // Add the background after averaging
88        dp[0] = scale();
89        dp[1] = radius();
90        dp[2] = length();
91        dp[3] = sldCyl();
92        dp[4] = sldSolv();
93        dp[5] = 0.0;
94
95        // Get the dispersion points for the radius
96        vector<WeightPoint> weights_rad;
97        radius.get_weights(weights_rad);
98
99        // Get the dispersion points for the length
100        vector<WeightPoint> weights_len;
101        length.get_weights(weights_len);
102
103        // Perform the computation, with all weight points
104        double sum = 0.0;
105        double norm = 0.0;
106        double vol = 0.0;
107
108        // Loop over radius weight points
109        for(size_t i=0; i<weights_rad.size(); i++) {
110                dp[1] = weights_rad[i].value;
111
112                // Loop over length weight points
113                for(size_t j=0; j<weights_len.size(); j++) {
114                        dp[2] = weights_len[j].value;
115                        //Un-normalize by volume
116                        sum += weights_rad[i].weight
117                                * weights_len[j].weight * CylinderForm(dp, q)
118                                *pow(weights_rad[i].value,2)*weights_len[j].value;
119
120                        //Find average volume
121                        vol += weights_rad[i].weight
122                                * weights_len[j].weight *pow(weights_rad[i].value,2)*weights_len[j].value;
123                        norm += weights_rad[i].weight
124                                * weights_len[j].weight;
125                }
126        }
127        if (vol != 0.0 && norm != 0.0) {
128                //Re-normalize by avg volume
129                sum = sum/(vol/norm);}
130
131        return sum/norm + background();
132}
133
134/**
135 * Function to evaluate 2D scattering function
136 * @param pars: parameters of the cylinder
137 * @param q: q-value
138 * @param q_x: q_x / q
139 * @param q_y: q_y / q
140 * @return: function value
141 */
142static double cylinder_analytical_2D_scaled(CylinderParameters *pars, double q, double q_x, double q_y) {
143    double cyl_x, cyl_y;//, cyl_z;
144    //double q_z;
145    double alpha, vol, cos_val;
146    double answer = 0.0;
147    double form = 0.0;
148    //convert angle degree to radian
149    double pi = 4.0*atan(1.0);
150    double theta = pars->cyl_theta * pi/180.0;
151    double phi = pars->cyl_phi * pi/180.0;
152    double sld_solv = pars->sldSolv;
153    double sld_cyl = pars->sldCyl;
154    double m_max = pars->M0_sld_cyl;
155    double m_max_solv = pars->M0_sld_solv;
156    double contrast = 0.0;
157
158    // Cylinder orientation
159        cyl_x = cos(theta) * cos(phi);
160    cyl_y = sin(theta);
161    //cyl_z = -cos(theta) * sin(phi);
162    // q vector
163    //q_z = 0.0;
164
165    // Compute the angle btw vector q and the
166    // axis of the cylinder
167    cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
168
169    // The following test should always pass
170    if (fabs(cos_val)>1.0) {
171      printf("cyl_ana_2D: Unexpected error: |cos(alpha)=%g|>1\n", cos_val);
172      printf("cyl_ana_2D: at theta=%g and phi=%g.", theta, phi);
173      return 1.0;
174    }
175
176    // Note: cos(alpha) = 0 and 1 will get an
177    // undefined value from CylKernel
178  alpha = acos( cos_val );
179  if (alpha == 0.0){
180        alpha = 1.0e-26;
181        }
182  // Call the IGOR library function to get the kernel
183  //answer = CylKernel(q, pars->radius, pars->length/2.0, alpha) / sin(alpha);
184
185    // Call the IGOR library function to get the kernel
186    form = CylKernel(q, pars->radius, pars->length/2.0, alpha) / sin(alpha);
187
188        if (m_max < 1.0e-32 && m_max_solv < 1.0e-32){
189                // Multiply by contrast^2
190                contrast = (pars->sldCyl - pars->sldSolv);
191                answer = contrast * contrast * form;
192        }
193        else{
194                double qx = q_x;
195                double qy = q_y;
196                double s_theta = pars->Up_theta;
197                double m_phi = pars->M_phi_cyl;
198                double m_theta = pars->M_theta_cyl;
199                double m_phi_solv = pars->M_phi_solv;
200                double m_theta_solv = pars->M_theta_solv;
201                double in_spin = pars->Up_frac_i;
202                double out_spin = pars->Up_frac_f;
203                polar_sld p_sld;
204                polar_sld p_sld_solv;
205                p_sld = cal_msld(1, qx, qy, sld_cyl, m_max, m_theta, m_phi, 
206                                                in_spin, out_spin, s_theta);
207                p_sld_solv = cal_msld(1, qx, qy, sld_solv, m_max_solv, m_theta_solv, m_phi_solv, 
208                                                in_spin, out_spin, s_theta);
209                //up_up
210                if (in_spin > 0.0 && out_spin > 0.0){                   
211                        answer += ((p_sld.uu- p_sld_solv.uu) * (p_sld.uu- p_sld_solv.uu) * form);
212                        }
213                //down_down
214                if (in_spin < 1.0 && out_spin < 1.0){
215                        answer += ((p_sld.dd - p_sld_solv.dd) * (p_sld.dd - p_sld_solv.dd) * form);
216                        }
217                //up_down
218                if (in_spin > 0.0 && out_spin < 1.0){
219                        answer += ((p_sld.re_ud - p_sld_solv.re_ud) * (p_sld.re_ud - p_sld_solv.re_ud) * form);
220                        answer += ((p_sld.im_ud - p_sld_solv.im_ud) * (p_sld.im_ud - p_sld_solv.im_ud) * form);
221                        }
222                //down_up       
223                if (in_spin < 1.0 && out_spin > 0.0){
224                        answer += ((p_sld.re_du - p_sld_solv.re_du) * (p_sld.re_du - p_sld_solv.re_du) * form);
225                        answer += ((p_sld.im_du - p_sld_solv.im_du) * (p_sld.im_du - p_sld_solv.im_du) * form);
226                        }
227        }
228
229    //normalize by cylinder volume
230    //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
231    vol = acos(-1.0) * pars->radius * pars->radius * pars->length;
232    answer *= vol;
233
234    //convert to [cm-1]
235    answer *= 1.0e8;
236
237    //Scale
238    answer *= pars->scale;
239
240    // add in the background
241    answer += pars->background;
242
243    return answer;
244}
245
246/**
247 * Function to evaluate 2D scattering function
248 * @param pars: parameters of the cylinder
249 * @param q: q-value
250 * @return: function value
251 */
252static double cylinder_analytical_2DXY(CylinderParameters *pars, double qx, double qy) {
253  double q;
254  q = sqrt(qx*qx+qy*qy);
255  return cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
256}
257
258/**
259 * Function to evaluate 2D scattering function
260 * @param q_x: value of Q along x
261 * @param q_y: value of Q along y
262 * @return: function value
263 */
264double CylinderModel :: operator()(double qx, double qy) {
265        CylinderParameters dp;
266        // Fill parameter array
267        dp.scale      = scale();
268        dp.radius     = radius();
269        dp.length     = length();
270        dp.sldCyl   = sldCyl();
271        dp.sldSolv   = sldSolv();
272        dp.background = 0.0;
273        dp.cyl_theta  = cyl_theta();
274        dp.cyl_phi    = cyl_phi();
275        dp.Up_theta =  Up_theta();
276        dp.M_phi_cyl =  M_phi_cyl();
277        dp.M_theta_cyl =  M_theta_cyl();
278        dp.M0_sld_cyl =  M0_sld_cyl();
279        dp.M_phi_solv =  M_phi_solv();
280        dp.M_theta_solv =  M_theta_solv();
281        dp.M0_sld_solv =  M0_sld_solv();
282        dp.Up_frac_i =  Up_frac_i();
283        dp.Up_frac_f =  Up_frac_f();
284       
285        // Get the dispersion points for the radius
286        vector<WeightPoint> weights_rad;
287        radius.get_weights(weights_rad);
288
289        // Get the dispersion points for the length
290        vector<WeightPoint> weights_len;
291        length.get_weights(weights_len);
292
293        // Get angular averaging for theta
294        vector<WeightPoint> weights_theta;
295        cyl_theta.get_weights(weights_theta);
296
297        // Get angular averaging for phi
298        vector<WeightPoint> weights_phi;
299        cyl_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 radius weight points
308        for(size_t i=0; i<weights_rad.size(); i++) {
309                dp.radius = weights_rad[i].value;
310
311
312                // Loop over length weight points
313                for(size_t j=0; j<weights_len.size(); j++) {
314                        dp.length = weights_len[j].value;
315
316                        // Average over theta distribution
317                        for(size_t k=0; k<weights_theta.size(); k++) {
318                                dp.cyl_theta = weights_theta[k].value;
319
320                                // Average over phi distribution
321                                for(size_t l=0; l<weights_phi.size(); l++) {
322                                        dp.cyl_phi = weights_phi[l].value;
323                                        //Un-normalize by volume
324                                        double _ptvalue = weights_rad[i].weight
325                                                * weights_len[j].weight
326                                                * weights_theta[k].weight
327                                                * weights_phi[l].weight
328                                                * cylinder_analytical_2DXY(&dp, qx, qy)
329                                                *pow(weights_rad[i].value,2)*weights_len[j].value;
330                                        if (weights_theta.size()>1) {
331                                                _ptvalue *= fabs(cos(weights_theta[k].value*pi/180.0));
332                                        }
333                                        sum += _ptvalue;
334                                        //Find average volume
335                                        vol += weights_rad[i].weight
336                                                        * weights_len[j].weight
337                                                        * pow(weights_rad[i].value,2)*weights_len[j].value;
338                                        //Find norm for volume
339                                        norm_vol += weights_rad[i].weight
340                                                        * weights_len[j].weight;
341
342                                        norm += weights_rad[i].weight
343                                                * weights_len[j].weight
344                                                * weights_theta[k].weight
345                                                * weights_phi[l].weight;
346
347                                }
348                        }
349                }
350        }
351        // Averaging in theta needs an extra normalization
352        // factor to account for the sin(theta) term in the
353        // integration (see documentation).
354        if (weights_theta.size()>1) norm = norm / asin(1.0);
355        if (vol != 0.0 && norm_vol != 0.0) {
356                //Re-normalize by avg volume
357                sum = sum/(vol/norm_vol);}
358
359        return sum/norm + background();
360}
361
362/**
363 * Function to evaluate 2D scattering function
364 * @param pars: parameters of the cylinder
365 * @param q: q-value
366 * @param phi: angle phi
367 * @return: function value
368 */
369double CylinderModel :: evaluate_rphi(double q, double phi) {
370        double qx = q*cos(phi);
371        double qy = q*sin(phi);
372        return (*this).operator()(qx, qy);
373}
374/**
375 * Function to calculate effective radius
376 * @return: effective radius value
377 */
378double CylinderModel :: calculate_ER() {
379        CylinderParameters dp;
380
381        dp.radius     = radius();
382        dp.length     = length();
383        double rad_out = 0.0;
384
385        // Perform the computation, with all weight points
386        double sum = 0.0;
387        double norm = 0.0;
388
389        // Get the dispersion points for the major shell
390        vector<WeightPoint> weights_length;
391        length.get_weights(weights_length);
392
393        // Get the dispersion points for the minor shell
394        vector<WeightPoint> weights_radius ;
395        radius.get_weights(weights_radius);
396
397        // Loop over major shell weight points
398        for(int i=0; i< (int)weights_length.size(); i++) {
399                dp.length = weights_length[i].value;
400                for(int k=0; k< (int)weights_radius.size(); k++) {
401                        dp.radius = weights_radius[k].value;
402                        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
403                        sum +=weights_length[i].weight
404                                * weights_radius[k].weight*DiamCyl(dp.length,dp.radius)/2.0;
405                        norm += weights_length[i].weight* weights_radius[k].weight;
406                }
407        }
408        if (norm != 0){
409                //return the averaged value
410                rad_out =  sum/norm;}
411        else{
412                //return normal value
413                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
414                rad_out = DiamCyl(dp.length,dp.radius)/2.0;}
415
416        return rad_out;
417}
418double CylinderModel :: calculate_VR() {
419  return 1.0;
420}
Note: See TracBrowser for help on using the repository browser.