source: sasview/sansmodels/src/sans/models/c_models/ellipticalcylinder.cpp @ 27953d1

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 27953d1 was 0f5bc9f, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

Update of all C models to the new style of C++ models

  • Property mode set to 100644
File size: 5.9 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 *      TODO: refactor so that we pull in the old sansmodels.c_extensions
21 */
22
23#include <math.h>
24#include "models.hh"
25#include "parameters.hh"
26#include <stdio.h>
27using namespace std;
28
29extern "C" {
30        #include "libCylinder.h"
31        #include "elliptical_cylinder.h"
32}
33
34EllipticalCylinderModel :: EllipticalCylinderModel() {
35        scale      = Parameter(1.0);
36        r_minor    = Parameter(20.0, true);
37        r_minor.set_min(0.0);
38        r_ratio    = Parameter(1.5, true);
39        r_ratio.set_min(0.0);
40        length     = Parameter(400.0, true);
41        length.set_min(0.0);
42        contrast   = Parameter(3.e-6);
43        background = Parameter(0.0);
44        cyl_theta  = Parameter(1.57, true);
45        cyl_phi    = Parameter(0.0, true);
46        cyl_psi    = Parameter(0.0, true);
47}
48
49/**
50 * Function to evaluate 1D scattering function
51 * The NIST IGOR library is used for the actual calculation.
52 * @param q: q-value
53 * @return: function value
54 */
55double EllipticalCylinderModel :: operator()(double q) {
56        double dp[6];
57
58        dp[0] = scale();
59        dp[1] = r_minor();
60        dp[2] = r_ratio();
61        dp[3] = length();
62        dp[4] = contrast();
63        dp[5] = 0.0;
64
65        // Get the dispersion points for the r_minor
66        vector<WeightPoint> weights_rad;
67        r_minor.get_weights(weights_rad);
68
69        // Get the dispersion points for the r_ratio
70        vector<WeightPoint> weights_rat;
71        r_ratio.get_weights(weights_rat);
72
73        // Get the dispersion points for the length
74        vector<WeightPoint> weights_len;
75        length.get_weights(weights_len);
76
77        // Perform the computation, with all weight points
78        double sum = 0.0;
79        double norm = 0.0;
80
81        // Loop over r_minor weight points
82        for(int i=0; i<weights_rad.size(); i++) {
83                dp[1] = weights_rad[i].value;
84
85                // Loop over r_ratio weight points
86                for(int j=0; j<weights_rat.size(); j++) {
87                        dp[2] = weights_rat[j].value;
88
89                        // Loop over length weight points
90                        for(int k=0; k<weights_len.size(); k++) {
91                                dp[3] = weights_len[k].value;
92
93                                sum += weights_rad[i].weight
94                                        * weights_len[k].weight
95                                        * weights_rat[j].weight
96                                        * EllipCyl20(dp, q);
97                                norm += weights_rad[i].weight
98                                * weights_len[k].weight
99                                * weights_rat[j].weight;
100                        }
101                }
102        }
103        return sum/norm + background();
104}
105
106/**
107 * Function to evaluate 2D scattering function
108 * @param q_x: value of Q along x
109 * @param q_y: value of Q along y
110 * @return: function value
111 */
112double EllipticalCylinderModel :: operator()(double qx, double qy) {
113        EllipticalCylinderParameters dp;
114        // Fill parameter array
115        dp.scale      = scale();
116        dp.r_minor    = r_minor();
117        dp.r_ratio    = r_ratio();
118        dp.length     = length();
119        dp.contrast   = contrast();
120        dp.background = 0.0;
121        dp.cyl_theta  = cyl_theta();
122        dp.cyl_phi    = cyl_phi();
123        dp.cyl_psi    = cyl_psi();
124
125        // Get the dispersion points for the r_minor
126        vector<WeightPoint> weights_rad;
127        r_minor.get_weights(weights_rad);
128
129        // Get the dispersion points for the r_ratio
130        vector<WeightPoint> weights_rat;
131        r_ratio.get_weights(weights_rat);
132
133        // Get the dispersion points for the length
134        vector<WeightPoint> weights_len;
135        length.get_weights(weights_len);
136
137        // Get angular averaging for theta
138        vector<WeightPoint> weights_theta;
139        cyl_theta.get_weights(weights_theta);
140
141        // Get angular averaging for phi
142        vector<WeightPoint> weights_phi;
143        cyl_phi.get_weights(weights_phi);
144
145        // Get angular averaging for psi
146        vector<WeightPoint> weights_psi;
147        cyl_psi.get_weights(weights_psi);
148
149        // Perform the computation, with all weight points
150        double sum = 0.0;
151        double norm = 0.0;
152
153        // Loop over minor radius weight points
154        for(int i=0; i<weights_rad.size(); i++) {
155                dp.r_minor = weights_rad[i].value;
156
157
158                // Loop over length weight points
159                for(int j=0; j<weights_len.size(); j++) {
160                        dp.length = weights_len[j].value;
161
162                        // Loop over r_ration weight points
163                        for(int m=0; m<weights_rat.size(); m++) {
164                                dp.r_ratio = weights_rat[m].value;
165
166                        // Average over theta distribution
167                        for(int k=0; k<weights_theta.size(); k++) {
168                                dp.cyl_theta = weights_theta[k].value;
169
170                                // Average over phi distribution
171                                for(int l=0; l<weights_phi.size(); l++) {
172                                        dp.cyl_phi = weights_phi[l].value;
173
174                                // Average over phi distribution
175                                for(int o=0; o<weights_psi.size(); o++) {
176                                        dp.cyl_psi = weights_psi[o].value;
177
178                                        double _ptvalue = weights_rad[i].weight
179                                                * weights_len[j].weight
180                                                * weights_rat[m].weight
181                                                * weights_theta[k].weight
182                                                * weights_phi[l].weight
183                                                * weights_psi[o].weight
184                                                * elliptical_cylinder_analytical_2DXY(&dp, qx, qy);
185                                        if (weights_theta.size()>1) {
186                                                _ptvalue *= sin(weights_theta[k].value);
187                                        }
188                                        sum += _ptvalue;
189
190                                        norm += weights_rad[i].weight
191                                                * weights_len[j].weight
192                                                * weights_rat[m].weight
193                                                * weights_theta[k].weight
194                                                * weights_phi[l].weight
195                                                * weights_psi[o].weight;
196
197                                }
198                                }
199                        }
200                        }
201                }
202        }
203        // Averaging in theta needs an extra normalization
204        // factor to account for the sin(theta) term in the
205        // integration (see documentation).
206        if (weights_theta.size()>1) norm = norm / asin(1.0);
207        return sum/norm + background();
208}
209
210/**
211 * Function to evaluate 2D scattering function
212 * @param pars: parameters of the cylinder
213 * @param q: q-value
214 * @param phi: angle phi
215 * @return: function value
216 */
217double EllipticalCylinderModel :: evaluate_rphi(double q, double phi) {
218        double qx = q*cos(phi);
219        double qy = q*sin(phi);
220        return (*this).operator()(qx, qy);
221}
Note: See TracBrowser for help on using the repository browser.