source: sasview/sansmodels/src/sans/models/c_models/capcyl.cpp @ 0cbfaeb

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 0cbfaeb was 4628e31, checked in by Jae Cho <jhjcho@…>, 14 years ago

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