source: sasview/sansmodels/src/sans/models/c_models/cylinder.cpp @ e2f7b92

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 e2f7b92 was f9bf661, checked in by Jae Cho <jhjcho@…>, 15 years ago

updated documents

  • Property mode set to 100644
File size: 7.7 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 "libStructureFactor.h"
32        #include "cylinder.h"
33}
34
35CylinderModel :: CylinderModel() {
36        scale      = Parameter(1.0);
37        radius     = Parameter(20.0, true);
38        radius.set_min(0.0);
39        length     = Parameter(400.0, true);
40        length.set_min(0.0);
41        contrast   = Parameter(3.e-6);
42        background = Parameter(0.0);
43        cyl_theta  = Parameter(0.0, true);
44        cyl_phi    = Parameter(0.0, true);
45}
46
47/**
48 * Function to evaluate 1D scattering function
49 * The NIST IGOR library is used for the actual calculation.
50 * @param q: q-value
51 * @return: function value
52 */
53double CylinderModel :: operator()(double q) {
54        double dp[5];
55
56        // Fill parameter array for IGOR library
57        // Add the background after averaging
58        dp[0] = scale();
59        dp[1] = radius();
60        dp[2] = length();
61        dp[3] = contrast();
62        dp[4] = 0.0;
63
64        // Get the dispersion points for the radius
65        vector<WeightPoint> weights_rad;
66        radius.get_weights(weights_rad);
67
68        // Get the dispersion points for the length
69        vector<WeightPoint> weights_len;
70        length.get_weights(weights_len);
71
72        // Perform the computation, with all weight points
73        double sum = 0.0;
74        double norm = 0.0;
75
76        // Loop over radius weight points
77        for(int i=0; i<weights_rad.size(); i++) {
78                dp[1] = weights_rad[i].value;
79
80                // Loop over length weight points
81                for(int j=0; j<weights_len.size(); j++) {
82                        dp[2] = weights_len[j].value;
83
84                        sum += weights_rad[i].weight
85                                * weights_len[j].weight * CylinderForm(dp, q);
86                        norm += weights_rad[i].weight
87                                * weights_len[j].weight;
88                }
89        }
90        return sum/norm + background();
91}
92
93/**
94 * Function to evaluate 2D scattering function
95 * @param q_x: value of Q along x
96 * @param q_y: value of Q along y
97 * @return: function value
98 */
99double CylinderModel :: operator()(double qx, double qy) {
100        CylinderParameters dp;
101        // Fill parameter array
102        dp.scale      = scale();
103        dp.radius     = radius();
104        dp.length     = length();
105        dp.contrast   = contrast();
106        dp.background = 0.0;
107        dp.cyl_theta  = cyl_theta();
108        dp.cyl_phi    = cyl_phi();
109
110        // Get the dispersion points for the radius
111        vector<WeightPoint> weights_rad;
112        radius.get_weights(weights_rad);
113
114        // Get the dispersion points for the length
115        vector<WeightPoint> weights_len;
116        length.get_weights(weights_len);
117
118        // Get angular averaging for theta
119        vector<WeightPoint> weights_theta;
120        cyl_theta.get_weights(weights_theta);
121
122        // Get angular averaging for phi
123        vector<WeightPoint> weights_phi;
124        cyl_phi.get_weights(weights_phi);
125
126        // Perform the computation, with all weight points
127        double sum = 0.0;
128        double norm = 0.0;
129
130        // Loop over radius weight points
131        for(int i=0; i<weights_rad.size(); i++) {
132                dp.radius = weights_rad[i].value;
133
134
135                // Loop over length weight points
136                for(int j=0; j<weights_len.size(); j++) {
137                        dp.length = weights_len[j].value;
138
139                        // Average over theta distribution
140                        for(int k=0; k<weights_theta.size(); k++) {
141                                dp.cyl_theta = weights_theta[k].value;
142
143                                // Average over phi distribution
144                                for(int l=0; l<weights_phi.size(); l++) {
145                                        dp.cyl_phi = weights_phi[l].value;
146
147                                        double _ptvalue = weights_rad[i].weight
148                                                * weights_len[j].weight
149                                                * weights_theta[k].weight
150                                                * weights_phi[l].weight
151                                                * cylinder_analytical_2DXY(&dp, qx, qy);
152                                        if (weights_theta.size()>1) {
153                                                _ptvalue *= sin(weights_theta[k].value);
154                                        }
155                                        sum += _ptvalue;
156
157                                        norm += weights_rad[i].weight
158                                                * weights_len[j].weight
159                                                * weights_theta[k].weight
160                                                * weights_phi[l].weight;
161
162                                }
163                        }
164                }
165        }
166        // Averaging in theta needs an extra normalization
167        // factor to account for the sin(theta) term in the
168        // integration (see documentation).
169        if (weights_theta.size()>1) norm = norm / asin(1.0);
170        return sum/norm + background();
171}
172
173/**
174 * Function to evaluate 2D scattering function
175 * @param pars: parameters of the cylinder
176 * @param q: q-value
177 * @param phi: angle phi
178 * @return: function value
179 */
180double CylinderModel :: evaluate_rphi(double q, double phi) {
181        double qx = q*cos(phi);
182        double qy = q*sin(phi);
183        return (*this).operator()(qx, qy);
184}
185/**
186 * Function to calculate effective radius
187 * @return: effective radius value
188 */
189double CylinderModel :: calculate_ER() {
190        CylinderParameters dp;
191
192        dp.radius     = radius();
193        dp.length     = length();
194        double rad_out = 0.0;
195
196        // Perform the computation, with all weight points
197        double sum = 0.0;
198        double norm = 0.0;
199
200        // Get the dispersion points for the major shell
201        vector<WeightPoint> weights_length;
202        length.get_weights(weights_length);
203
204        // Get the dispersion points for the minor shell
205        vector<WeightPoint> weights_radius ;
206        radius.get_weights(weights_radius);
207
208        // Loop over major shell weight points
209        for(int i=0; i< (int)weights_length.size(); i++) {
210                dp.length = weights_length[i].value;
211                for(int k=0; k< (int)weights_radius.size(); k++) {
212                        dp.radius = weights_radius[k].value;
213                        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
214                        sum +=weights_length[i].weight
215                                * weights_radius[k].weight*DiamCyl(dp.length,dp.radius)/2.0;
216                        norm += weights_length[i].weight* weights_radius[k].weight;
217                }
218        }
219        if (norm != 0){
220                //return the averaged value
221                rad_out =  sum/norm;}
222        else{
223                //return normal value
224                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
225                rad_out = DiamCyl(dp.length,dp.radius)/2.0;}
226
227        return rad_out;
228}
229// Testing code
230int main(void)
231{
232        CylinderModel c = CylinderModel();
233
234        printf("Length = %g\n", c.length());
235        printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
236        printf("I(Q=%g) = %g\n", 0.001, c(0.001));
237        c.radius.dispersion = new GaussianDispersion();
238        c.radius.dispersion->npts = 100;
239        c.radius.dispersion->width = 5;
240
241        //c.length.dispersion = GaussianDispersion();
242        //c.length.dispersion.npts = 20;
243        //c.length.dispersion.width = 65;
244
245        printf("I(Q=%g) = %g\n", 0.001, c(0.001));
246        c.scale = 10.0;
247        printf("I(Q=%g) = %g\n", 0.001, c(0.001));
248        printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
249        printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854));
250
251        // Average over phi at theta=90 deg
252        c.cyl_theta = 1.57;
253        double values_th[100];
254        double values[100];
255        double weights[100];
256        double pi = acos(-1.0);
257        printf("pi=%g\n", pi);
258        for(int i=0; i<100; i++){
259                values[i] = (float)i*2.0*pi/99.0;
260                values_th[i] = (float)i*pi/99.0;
261                weights[i] = 1.0;
262        }
263        //c.radius.dispersion->width = 0;
264        c.cyl_phi.dispersion = new ArrayDispersion();
265        c.cyl_theta.dispersion = new ArrayDispersion();
266        (*c.cyl_phi.dispersion).set_weights(100, values, weights);
267        (*c.cyl_theta.dispersion).set_weights(100, values_th, weights);
268
269        double i_avg = c(0.01, 0.01);
270        double i_1d = c(sqrt(0.0002));
271
272        printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg);
273        printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d);
274        printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg);
275
276
277        return 0;
278}
279
Note: See TracBrowser for help on using the repository browser.