source: sasview/sansmodels/src/sans/models/c_models/capcyl.cpp @ 916f5c0

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 916f5c0 was 34c2649, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Re #4 Fixed warnings
  • Property mode set to 100644
File size: 7.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 "models.hh"
24#include "parameters.hh"
25#include <stdio.h>
26using namespace std;
27
28extern "C" {
29        #include "libCylinder.h"
30        #include "capcyl.h"
31}
32
33CappedCylinderModel :: CappedCylinderModel() {
34        scale      = Parameter(1.0);
35        rad_cyl         = Parameter(20.0);
36        rad_cyl.set_min(0.0);
37        len_cyl     = Parameter(400.0, true);
38        len_cyl.set_min(0.0);
39        rad_cap = Parameter(40.0);
40        rad_cap.set_min(0.0);
41        sld_capcyl  = Parameter(1.0e-6);
42        sld_solv   = Parameter(6.3e-6);
43        background = Parameter(0.0);
44        theta  = Parameter(0.0, true);
45        phi    = Parameter(0.0, true);
46}
47
48/**
49 * Function to evaluate 1D scattering function
50 * The NIST IGOR library is used for the actual calculation.
51 * @param q: q-value
52 * @return: function value
53 */
54double CappedCylinderModel :: operator()(double q) {
55        double dp[7];
56
57        // Fill parameter array for IGOR library
58        // Add the background after averaging
59        dp[0] = scale();
60        dp[1] = rad_cyl();
61        dp[2] = len_cyl();
62        dp[3] = rad_cap();
63        dp[4] = sld_capcyl();
64        dp[5] = sld_solv();
65        dp[6] = 0.0;
66
67        // Get the dispersion points for the rad_cyl
68        vector<WeightPoint> weights_rad_cyl;
69        rad_cyl.get_weights(weights_rad_cyl);
70        // Get the dispersion points for the len_cyl
71        vector<WeightPoint> weights_len_cyl;
72        len_cyl.get_weights(weights_len_cyl);
73        // Get the dispersion points for the rad_cap
74        vector<WeightPoint> weights_rad_cap;
75        rad_cap.get_weights(weights_rad_cap);
76
77        // Perform the computation, with all weight points
78        double sum = 0.0;
79        double norm = 0.0;
80        double vol = 0.0;
81        double pi,hDist,result;
82        double vol_i = 0.0;
83        pi = 4.0*atan(1.0);
84        // Loop over radius weight points
85        for(size_t i=0; i<weights_rad_cyl.size(); i++) {
86                dp[1] = weights_rad_cyl[i].value;
87                for(size_t j=0; j<weights_len_cyl.size(); j++) {
88                        dp[2] = weights_len_cyl[j].value;
89                        for(size_t k=0; k<weights_rad_cap.size(); k++) {
90                                dp[3] = weights_rad_cap[k].value;
91
92                                //Un-normalize SphereForm by volume
93                                hDist = -1.0*sqrt(fabs(dp[3]*dp[3]-dp[1]*dp[1]));
94                                vol_i = pi*dp[1]*dp[1]*dp[2]+2.0*pi/3.0*((dp[3]-hDist)*(dp[3]-hDist)*
95                                                                (2.0*(dp[3]+hDist)));
96                                result =  CappedCylinder(dp, q) * vol_i;
97                                // This FIXES a singualrity the kernel in libigor.
98                                if ( result == INFINITY || result == NAN){
99                                        result = 0.0;
100                                }
101                                sum += weights_rad_cyl[i].weight*weights_len_cyl[j].weight*weights_rad_cap[k].weight
102                                        * result;
103                                //Find average volume
104                                vol += weights_rad_cyl[i].weight*weights_len_cyl[j].weight*weights_rad_cap[k].weight
105                                        * vol_i;
106
107                                norm += weights_rad_cyl[i].weight*weights_len_cyl[j].weight*weights_rad_cap[k].weight;
108                        }
109                }
110        }
111
112        if (vol != 0.0 && norm != 0.0) {
113                //Re-normalize by avg volume
114                sum = sum/(vol/norm);}
115        return sum/norm + background();
116}
117
118/**
119 * Function to evaluate 2D scattering function
120 * @param q_x: value of Q along x
121 * @param q_y: value of Q along y
122 * @return: function value
123 */
124double CappedCylinderModel :: operator()(double qx, double qy) {
125        CapCylParameters dp;
126
127        dp.scale = scale();
128        dp.rad_cyl = rad_cyl();
129        dp.len_cyl = len_cyl();
130        dp.rad_cap = rad_cap();
131        dp.sld_capcyl = sld_capcyl();
132        dp.sld_solv = sld_solv();
133        dp.background = 0.0;
134        dp.theta  = theta();
135        dp.phi    = phi();
136
137        // Get the dispersion points for the rad_bar
138        vector<WeightPoint> weights_rad_cyl;
139        rad_cyl.get_weights(weights_rad_cyl);
140
141        // Get the dispersion points for the len_bar
142        vector<WeightPoint> weights_len_cyl;
143        len_cyl.get_weights(weights_len_cyl);
144
145        // Get the dispersion points for the rad_bell
146        vector<WeightPoint> weights_rad_cap;
147        rad_cap.get_weights(weights_rad_cap);
148
149        // Get angular averaging for theta
150        vector<WeightPoint> weights_theta;
151        theta.get_weights(weights_theta);
152
153        // Get angular averaging for phi
154        vector<WeightPoint> weights_phi;
155        phi.get_weights(weights_phi);
156
157
158        // Perform the computation, with all weight points
159        double sum = 0.0;
160        double norm = 0.0;
161        double norm_vol = 0.0;
162        double vol = 0.0;
163        double pi,hDist;
164        double vol_i = 0.0;
165        pi = 4.0*atan(1.0);
166
167        // Loop over radius weight points
168        for(size_t i=0; i<weights_rad_cyl.size(); i++) {
169                dp.rad_cyl = weights_rad_cyl[i].value;
170                for(size_t j=0; j<weights_len_cyl.size(); j++) {
171                        dp.len_cyl = weights_len_cyl[j].value;
172                        for(size_t k=0; k<weights_rad_cap.size(); k++) {
173                                dp.rad_cap = weights_rad_cap[k].value;
174                                // Average over theta distribution
175                                for(size_t l=0; l< weights_theta.size(); l++) {
176                                        dp.theta = weights_theta[l].value;
177                                        // Average over phi distribution
178                                        for(size_t m=0; m< weights_phi.size(); m++) {
179                                                dp.phi = weights_phi[m].value;
180                                                //Un-normalize Form by volume
181                                                hDist = -1.0*sqrt(fabs(dp.rad_cap*dp.rad_cap-dp.rad_cyl*dp.rad_cyl));
182                                                vol_i = pi*dp.rad_cyl*dp.rad_cyl*dp.len_cyl+2.0*pi/3.0*((dp.rad_cap-hDist)*(dp.rad_cap-hDist)*
183                                                                                (2*dp.rad_cap+hDist));
184
185                                                double _ptvalue = weights_rad_cyl[i].weight
186                                                                                        * weights_len_cyl[j].weight
187                                                                                        * weights_rad_cap[k].weight
188                                                                                        * weights_theta[l].weight
189                                                                                        * weights_phi[m].weight
190                                                                                        * vol_i
191                                                                                        * capcyl_analytical_2DXY(&dp, qx, qy);
192                                                                                        //* pow(weights_rad[i].value,3.0);
193                                                // Consider when there is infinte or nan.
194                                                if ( _ptvalue == INFINITY || _ptvalue == NAN){
195                                                        _ptvalue = 0.0;
196                                                }
197                                                if (weights_theta.size()>1) {
198                                                        _ptvalue *= fabs(sin(weights_theta[l].value*pi/180.0));
199                                                }
200                                                sum += _ptvalue;
201                                                // This model dose not need the volume of spheres correction!!!
202                                                //Find average volume
203                                                vol += weights_rad_cyl[i].weight
204                                                                * weights_len_cyl[j].weight
205                                                                * weights_rad_cap[k].weight
206                                                                * vol_i;
207                                                //Find norm for volume
208                                                norm_vol += weights_rad_cyl[i].weight
209                                                                * weights_len_cyl[j].weight
210                                                                * weights_rad_cap[k].weight;
211
212                                                norm += weights_rad_cyl[i].weight
213                                                                * weights_len_cyl[j].weight
214                                                                * weights_rad_cap[k].weight
215                                                                * weights_theta[l].weight
216                                                                * weights_phi[m].weight;
217                                        }
218                                }
219                        }
220                }
221        }
222        // Averaging in theta needs an extra normalization
223        // factor to account for the sin(theta) term in the
224        // integration (see documentation).
225        if (weights_theta.size()>1) norm = norm / asin(1.0);
226
227        if (vol != 0.0 && norm_vol != 0.0) {
228                //Re-normalize by avg volume
229                sum = sum/(vol/norm_vol);}
230
231        return sum/norm + background();
232}
233
234/**
235 * Function to evaluate 2D scattering function
236 * @param pars: parameters of the SCCrystal
237 * @param q: q-value
238 * @param phi: angle phi
239 * @return: function value
240 */
241double CappedCylinderModel :: evaluate_rphi(double q, double phi) {
242        return (*this).operator()(q);
243}
244
245/**
246 * Function to calculate effective radius
247 * @return: effective radius value
248 */
249double CappedCylinderModel :: calculate_ER() {
250        //NOT implemented yet!!!
251        return 0.0;
252}
Note: See TracBrowser for help on using the repository browser.