source: sasview/sansmodels/src/sans/models/c_models/capcyl.cpp @ 67e258c

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 67e258c was 85bb870, checked in by Jae Cho <jhjcho@…>, 14 years ago

minor correcions of orientational order,etc
  • 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(int i=0; i<weights_rad_cyl.size(); i++) {
86                dp[1] = weights_rad_cyl[i].value;
87                for(int j=0; j<weights_len_cyl.size(); j++) {
88                        dp[2] = weights_len_cyl[j].value;
89                        for(int 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        double q = sqrt(qx*qx + qy*qy);
128        dp.scale = scale();
129        dp.rad_cyl = rad_cyl();
130        dp.len_cyl = len_cyl();
131        dp.rad_cap = rad_cap();
132        dp.sld_capcyl = sld_capcyl();
133        dp.sld_solv = sld_solv();
134        dp.background = 0.0;
135        dp.theta  = theta();
136        dp.phi    = phi();
137
138
139        // Get the dispersion points for the rad_bar
140        vector<WeightPoint> weights_rad_cyl;
141        rad_cyl.get_weights(weights_rad_cyl);
142
143        // Get the dispersion points for the len_bar
144        vector<WeightPoint> weights_len_cyl;
145        len_cyl.get_weights(weights_len_cyl);
146
147        // Get the dispersion points for the rad_bell
148        vector<WeightPoint> weights_rad_cap;
149        rad_cap.get_weights(weights_rad_cap);
150
151        // Get angular averaging for theta
152        vector<WeightPoint> weights_theta;
153        theta.get_weights(weights_theta);
154
155        // Get angular averaging for phi
156        vector<WeightPoint> weights_phi;
157        phi.get_weights(weights_phi);
158
159
160        // Perform the computation, with all weight points
161        double sum = 0.0;
162        double norm = 0.0;
163        double norm_vol = 0.0;
164        double vol = 0.0;
165        double pi,hDist,result;
166        double vol_i = 0.0;
167        pi = 4.0*atan(1.0);
168
169        // Loop over radius weight points
170        for(int i=0; i<weights_rad_cyl.size(); i++) {
171                dp.rad_cyl = weights_rad_cyl[i].value;
172                for(int j=0; j<weights_len_cyl.size(); j++) {
173                        dp.len_cyl = weights_len_cyl[j].value;
174                        for(int k=0; k<weights_rad_cap.size(); k++) {
175                                dp.rad_cap = weights_rad_cap[k].value;
176                                // Average over theta distribution
177                                for(int l=0; l< weights_theta.size(); l++) {
178                                        dp.theta = weights_theta[l].value;
179                                        // Average over phi distribution
180                                        for(int m=0; m< weights_phi.size(); m++) {
181                                                dp.phi = weights_phi[m].value;
182                                                //Un-normalize Form by volume
183                                                hDist = -1.0*sqrt(fabs(dp.rad_cap*dp.rad_cap-dp.rad_cyl*dp.rad_cyl));
184                                                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)*
185                                                                                (2*dp.rad_cap+hDist));
186
187                                                double _ptvalue = weights_rad_cyl[i].weight
188                                                                                        * weights_len_cyl[j].weight
189                                                                                        * weights_rad_cap[k].weight
190                                                                                        * weights_theta[l].weight
191                                                                                        * weights_phi[m].weight
192                                                                                        * vol_i
193                                                                                        * capcyl_analytical_2DXY(&dp, qx, qy);
194                                                                                        //* pow(weights_rad[i].value,3.0);
195                                                // Consider when there is infinte or nan.
196                                                if ( _ptvalue == INFINITY || _ptvalue == NAN){
197                                                        _ptvalue = 0.0;
198                                                }
199                                                if (weights_theta.size()>1) {
200                                                        _ptvalue *= fabs(sin(weights_theta[l].value));
201                                                }
202                                                sum += _ptvalue;
203                                                // This model dose not need the volume of spheres correction!!!
204                                                //Find average volume
205                                                vol += weights_rad_cyl[i].weight
206                                                                * weights_len_cyl[j].weight
207                                                                * weights_rad_cap[k].weight
208                                                                * vol_i;
209                                                //Find norm for volume
210                                                norm_vol += weights_rad_cyl[i].weight
211                                                                * weights_len_cyl[j].weight
212                                                                * weights_rad_cap[k].weight;
213
214                                                norm += weights_rad_cyl[i].weight
215                                                                * weights_len_cyl[j].weight
216                                                                * weights_rad_cap[k].weight
217                                                                * weights_theta[l].weight
218                                                                * weights_phi[m].weight;
219                                        }
220                                }
221                        }
222                }
223        }
224        // Averaging in theta needs an extra normalization
225        // factor to account for the sin(theta) term in the
226        // integration (see documentation).
227        if (weights_theta.size()>1) norm = norm / asin(1.0);
228
229        if (vol != 0.0 && norm_vol != 0.0) {
230                //Re-normalize by avg volume
231                sum = sum/(vol/norm_vol);}
232
233        return sum/norm + background();
234}
235
236/**
237 * Function to evaluate 2D scattering function
238 * @param pars: parameters of the SCCrystal
239 * @param q: q-value
240 * @param phi: angle phi
241 * @return: function value
242 */
243double CappedCylinderModel :: evaluate_rphi(double q, double phi) {
244        return (*this).operator()(q);
245}
246
247/**
248 * Function to calculate effective radius
249 * @return: effective radius value
250 */
251double CappedCylinderModel :: calculate_ER() {
252        //NOT implemented yet!!!
253}
Note: See TracBrowser for help on using the repository browser.