source: sasview/src/sas/models/c_extension/c_models/pringles.cpp @ 79492222

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 79492222 was 79492222, checked in by krzywon, 8 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 8.2 KB
Line 
1/**
2 * PringlesModel
3 *
4 * (c) 2013 / Andrew J Jackson / andrew.jackson@esss.se
5 *
6 * Model for "pringles" particle from K. Edler @ Bath University
7 *
8 */
9#include <math.h>
10#include "parameters.hh"
11#include "pringles.h"
12#include "cephes.h"
13using namespace std;
14
15extern "C"
16{
17#include "GaussWeights.h"
18#include "libStructureFactor.h"
19}
20
21// Convenience parameter structure
22typedef struct {
23    double scale;
24    double radius;
25    double thickness;
26    double alpha;
27    double beta;
28    double sldCyl;
29    double sldSolv;
30    double background;
31    double cyl_theta;
32    double cyl_phi;
33} PringleParameters;
34
35
36PringlesModel::PringlesModel() {
37        scale = Parameter(1.0);
38        radius = Parameter(60.0,true);
39        radius.set_min(0.0);
40        thickness = Parameter(10.0,true);
41        thickness.set_min(0.0);
42        alpha = Parameter(0.001,true);
43        beta = Parameter(0.02,true);
44        sld_pringle = Parameter(1.0e-6);
45        sld_solvent = Parameter(6.35e-6);
46        background = Parameter(0.0);
47}
48
49/**
50 * Function to evaluate 1D scattering function
51 * @param q: q-value
52 * @return: function value
53 */
54double PringlesModel::operator()(double q) {
55        double dp[8];
56        // Fill parameter array
57        // Add the background after averaging
58        dp[0] = scale();
59        dp[1] = radius();
60        dp[2] = thickness();
61        dp[3] = alpha();
62        dp[4] = beta();
63        dp[5] = sld_pringle();
64        dp[6] = sld_solvent();
65        dp[7] = 0.0;
66
67        // Get the dispersion points for the radius
68        vector<WeightPoint> weights_rad;
69        radius.get_weights(weights_rad);
70
71        // Get the dispersion points for the thickness
72        vector<WeightPoint> weights_thick;
73        thickness.get_weights(weights_thick);
74
75        // Get the dispersion points for alpha
76        vector<WeightPoint> weights_alpha;
77        alpha.get_weights(weights_alpha);
78
79        // Get the dispersion points for beta
80        vector<WeightPoint> weights_beta;
81        beta.get_weights(weights_beta);
82
83        // Perform the computation, with all weight points
84        double sum = 0.0;
85        double norm = 0.0;
86        double volnorm = 0.0;
87        double vol = 0.0;
88        double Pi;
89
90        Pi = 4.0 * atan(1.0);
91
92        // Loop over alpha weight points
93        for (size_t i = 0; i < weights_alpha.size(); i++) {
94                dp[3] = weights_alpha[i].value;
95
96                //Loop over beta weight points
97                for (size_t j = 0; j < weights_beta.size(); j++) {
98                        dp[4] = weights_beta[j].value;
99
100                        // Loop over thickness weight points
101                        for (size_t k = 0; k < weights_thick.size(); k++) {
102                                dp[2] = weights_thick[k].value;
103
104                                // Loop over radius weight points
105                                for (size_t l = 0; l < weights_rad.size(); l++) {
106                                        dp[1] = weights_rad[l].value;
107                                        sum += weights_rad[l].weight * weights_thick[k].weight
108                                                        * weights_alpha[i].weight * weights_beta[j].weight
109                                                        * pringle_form(dp, q);
110                                        //Find average volume
111                                        vol += weights_rad[l].weight * weights_thick[k].weight * Pi * pow(weights_rad[l].value, 2) * weights_thick[k].value;
112                                        volnorm += weights_rad[l].weight * weights_thick[k].weight;
113                                        norm +=  weights_rad[l].weight * weights_thick[k].weight * weights_alpha[i].weight * weights_beta[j].weight;
114
115                                }
116                        }
117                }
118        }
119
120        if (vol > 0.0 && norm > 0.0) {
121                //normalize by avg volume
122                sum = sum * (vol/volnorm);
123                return sum/norm + background();
124        } else {
125                return 0.0;
126        }
127}
128
129/**
130 * Function to evaluate 2D scattering function
131 * @param q_x: value of Q along x
132 * @param q_y: value of Q along y
133 * @return: function value
134 */
135double PringlesModel::operator()(double qx, double qy) {
136        double q = sqrt(qx * qx + qy * qy);
137        return (*this).operator()(q);
138}
139/**
140 * Function to evaluate 2D scattering function
141 * @param pars: parameters of the model
142 * @param q: q-value
143 * @param phi: angle phi
144 * @return: function value
145 */
146double PringlesModel::evaluate_rphi(double q, double phi) {
147        return (*this).operator()(q);
148}
149
150/**
151 * Function to calculate effective radius
152 * @return: effective radius value
153 */
154double PringlesModel::calculate_ER() {
155        PringleParameters dp;
156
157        dp.radius     = radius();
158        dp.thickness  = thickness();
159        double rad_out = 0.0;
160
161        // Perform the computation, with all weight points
162        double sum = 0.0;
163        double norm = 0.0;
164
165        // Get the dispersion points for the major shell
166        vector<WeightPoint> weights_thick;
167        thickness.get_weights(weights_thick);
168
169        // Get the dispersion points for the minor shell
170        vector<WeightPoint> weights_radius ;
171        radius.get_weights(weights_radius);
172
173        // Loop over major shell weight points
174        for(int i=0; i< (int)weights_thick.size(); i++) {
175                dp.thickness = weights_thick[i].value;
176                for(int k=0; k< (int)weights_radius.size(); k++) {
177                        dp.radius = weights_radius[k].value;
178                        //Note: output of "DiamCyl(dp.thick,dp.radius)" is DIAMETER.
179                        sum +=weights_thick[i].weight
180                                * weights_radius[k].weight*DiamCyl(dp.thickness,dp.radius)/2.0;
181                        norm += weights_thick[i].weight* weights_radius[k].weight;
182                }
183        }
184        if (norm != 0){
185                //return the averaged value
186                rad_out =  sum/norm;}
187        else{
188                //return normal value
189                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
190                rad_out = DiamCyl(dp.thickness,dp.radius)/2.0;}
191
192        return rad_out;
193}
194/**
195 * Function to calculate particle volume/total volume for shape models:
196 *    Most case returns 1 but for example for the vesicle model it is
197 *    (total volume - core volume)/total volume
198 *    (< 1 depending on the thickness).
199 * @return: effective radius value
200 */
201double PringlesModel::calculate_VR() {
202        return 1.0;
203}
204
205/*
206 * Useful work functions start here!
207 *
208 */
209
210static double pringle_form(double dp[], double q) {
211
212        double Pi;
213        int nord = 76, i=0;                     //order of integration
214        double uplim, lolim;            //upper and lower integration limits
215        double summ, phi, yyy, answer, vcyl;                    //running tally of integration
216        double delrho;
217
218        Pi = 4.0 * atan(1.0);
219        lolim = 0.0;
220        uplim = Pi / 2.0;
221
222        summ = 0.0;                     //initialize integral
223
224        delrho = dp[5] - dp[6] ; //make contrast term
225
226        for (i = 0; i < nord; i++) {
227                phi = (Gauss76Z[i] * (uplim - lolim) + uplim + lolim) / 2.0;
228                yyy = Gauss76Wt[i] * pringle_kernel(dp, q, phi);
229                summ += yyy;
230        }
231
232        answer = (uplim - lolim) / 2.0 * summ;
233
234        answer *= delrho*delrho;
235
236    //normalize by cylinder volume
237        //vcyl=Pi*dp[1]*dp[1]*dp[2];
238        //answer *= vcyl;
239
240    //convert to [cm-1]
241        answer *= 1.0e8;
242
243    //Scale by volume fraction
244        answer *= dp[0];
245
246        return answer;
247}
248
249static double pringle_kernel(double dp[], double q, double phi) {
250
251        double sumterm, sincarg, sincterm, nn, retval;
252
253        sincarg = q * dp[2] * cos(phi) / 2.0; //dp[2] = thickness
254        sincterm = pow(sin(sincarg) / sincarg, 2.0);
255
256        //calculate sum term from n = -3 to 3
257        sumterm = 0;
258        for (nn = -3; nn <= 3; nn = nn + 1) {
259                sumterm = sumterm
260                                + (pow(pringleC(dp, q, phi, nn), 2.0)
261                                                + pow(pringleS(dp, q, phi, nn), 2.0));
262        }
263
264        retval = 4.0 * sin(phi) * sumterm * sincterm;
265
266        return retval;
267
268}
269
270static double pringleC(double dp[], double q, double phi, double n) {
271
272        double nord, va, vb, summ;
273        double bessargs, cosarg, bessargcb;
274        double r, retval, yyy;
275        int ii;
276        // set up the integration
277        // end points and weights
278        nord = 76;
279        va = 0;
280        vb = dp[1]; //radius
281
282        // evaluate at Gauss points
283        // remember to index from 0,size-1
284
285        summ = 0.0;             // initialize integral
286        ii = 0;
287        do {
288                // Using 76 Gauss points
289                r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0;
290
291                bessargs = q * r * sin(phi);
292                cosarg = q * r * r * dp[3] * cos(phi);
293                bessargcb = q * r * r * dp[4] * cos(phi);
294
295                yyy = Gauss76Wt[ii] * r * cos(cosarg) * jn(n, bessargcb)
296                                * jn(2 * n, bessargs);
297                summ += yyy;
298
299                ii += 1;
300        } while (ii < nord);                    // end of loop over quadrature points
301        //
302        // calculate value of integral to return
303
304        retval = (vb - va) / 2.0 * summ;
305
306        retval = retval / pow(r, 2.0);
307
308        return retval;
309}
310
311static double pringleS(double dp[], double q, double phi, double n) {
312
313        double nord, va, vb, summ;
314        double bessargs, sinarg, bessargcb;
315        double r, retval, yyy;
316        int ii;
317        // set up the integration
318        // end points and weights
319        nord = 76;
320        va = 0;
321        vb = dp[1]; //radius
322
323        // evaluate at Gauss points
324        // remember to index from 0,size-1
325
326        summ = 0.0;             // initialize integral
327        ii = 0;
328        do {
329                // Using 76 Gauss points
330                r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0;
331
332                bessargs = q * r * sin(phi);
333                sinarg = q * r * r * dp[3] * cos(phi);
334                bessargcb = q * r * r * dp[4] * cos(phi);
335
336                yyy = Gauss76Wt[ii] * r * sin(sinarg) * jn(n, bessargcb)
337                                * jn(2 * n, bessargs);
338
339                summ += yyy;
340
341                ii += 1;
342        } while (ii < nord);                    // end of loop over quadrature points
343        //
344        // calculate value of integral to return
345
346        retval = (vb - va) / 2.0 * summ;
347
348        retval = retval / pow(r, 2.0);
349
350        return retval;
351}
Note: See TracBrowser for help on using the repository browser.