source: sasview/sansmodels/src/sans/models/c_models/pearlnecklace.cpp @ 8bac371

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

Re #4 Fixed warnings

  • Property mode set to 100644
File size: 5.4 KB
Line 
1
2#include <math.h>
3#include "models.hh"
4#include "parameters.hh"
5#include <stdio.h>
6using namespace std;
7
8extern "C" {
9        #include "pearlnecklace.h"
10}
11
12PearlNecklaceModel :: PearlNecklaceModel() {
13        scale = Parameter(1.0);
14        radius = Parameter(80.0, true);
15        radius.set_min(0.0);
16        edge_separation = Parameter(350.0, true);
17        edge_separation.set_min(0.0);
18        thick_string = Parameter(2.5, true);
19        thick_string.set_min(0.0);
20        num_pearls = Parameter(3);
21        num_pearls.set_min(0.0);
22        sld_pearl = Parameter(1.0e-06);
23        sld_string = Parameter(1.0e-06);
24        sld_solv = Parameter(6.3e-06);
25    background = Parameter(0.0);
26
27}
28
29/**
30 * Function to evaluate 1D PearlNecklaceModel function
31 * @param q: q-value
32 * @return: function value
33 */
34double PearlNecklaceModel :: operator()(double q) {
35        double dp[9];
36        // Fill parameter array for IGOR library
37        // Add the background after averaging
38        dp[0] = scale();
39        dp[1] = radius();
40        dp[2] = edge_separation();
41        dp[3] = thick_string();
42        dp[4] = num_pearls();
43        dp[5] = sld_pearl();
44        dp[6] = sld_string();
45        dp[7] = sld_solv();
46        dp[8] = 0.0;
47        double pi = 4.0*atan(1.0);
48        // No polydispersion supported in this model.
49        // Get the dispersion points for the radius
50        vector<WeightPoint> weights_radius;
51        radius.get_weights(weights_radius);
52        vector<WeightPoint> weights_edge_separation;
53        edge_separation.get_weights(weights_edge_separation);
54        // Perform the computation, with all weight points
55        double sum = 0.0;
56        double norm = 0.0;
57        double vol = 0.0;
58        double string_vol = 0.0;
59        double pearl_vol = 0.0;
60        double tot_vol = 0.0;
61        // Loop over core weight points
62        for(size_t i=0; i<weights_radius.size(); i++) {
63                dp[1] = weights_radius[i].value;
64                // Loop over thick_inter0 weight points
65                for(size_t j=0; j<weights_edge_separation.size(); j++) {
66                        dp[2] = weights_edge_separation[j].value;
67                        pearl_vol = 4.0 /3.0 * pi * pow(dp[1], 3);
68                        string_vol =dp[2] * pi * pow((dp[3] / 2.0), 2);
69                        tot_vol = (dp[4] - 1.0) * string_vol;
70                        tot_vol += dp[4] * pearl_vol;
71                        //Un-normalize Sphere by volume
72                        sum += weights_radius[i].weight * weights_edge_separation[j].weight
73                                * pearl_necklace_kernel(dp,q) * tot_vol;
74                        //Find average volume
75                        vol += weights_radius[i].weight * weights_edge_separation[j].weight
76                                * tot_vol;
77                        norm += weights_radius[i].weight * weights_edge_separation[j].weight;
78                }
79        }
80
81        if (vol != 0.0 && norm != 0.0) {
82                //Re-normalize by avg volume
83                sum = sum/(vol/norm);}
84
85        return sum/norm + background();
86}
87
88/**
89 * Function to evaluate 2D PearlNecklaceModel function
90 * @param q_x: value of Q along x
91 * @param q_y: value of Q along y
92 * @return: function value
93 */
94double PearlNecklaceModel :: operator()(double qx, double qy) {
95        double q = sqrt(qx*qx + qy*qy);
96        return (*this).operator()(q);
97}
98
99/**
100 * Function to evaluate PearlNecklaceModel function
101 * @param pars: parameters of the PearlNecklaceModel
102 * @param q: q-value
103 * @param phi: angle phi
104 * @return: function value
105 */
106double PearlNecklaceModel :: evaluate_rphi(double q, double phi) {
107        return (*this).operator()(q);
108}
109
110/**
111 * Function to calculate TOTAL radius
112 * Todo: decide whether or not we keep this calculation
113 * @return: effective radius value
114 */
115// No polydispersion supported in this model.
116// Calculate max radius assumming max_radius = effective radius
117// Note that this max radius is not affected by sld of layer, sld of interface, or
118// sld of solvent.
119double PearlNecklaceModel :: calculate_ER() {
120        PeralNecklaceParameters dp;
121
122        dp.scale = scale();
123        dp.radius = radius();
124        dp.edge_separation = edge_separation();
125        dp.thick_string = thick_string();
126        dp.num_pearls = num_pearls();
127
128        double rad_out = 0.0;
129        // Perform the computation, with all weight points
130        double sum = 0.0;
131        double norm = 0.0;
132        double pi = 4.0*atan(1.0);
133        // No polydispersion supported in this model.
134        // Get the dispersion points for the radius
135        vector<WeightPoint> weights_radius;
136        radius.get_weights(weights_radius);
137        vector<WeightPoint> weights_edge_separation;
138        edge_separation.get_weights(weights_edge_separation);
139        // Perform the computation, with all weight points
140        double string_vol = 0.0;
141        double pearl_vol = 0.0;
142        double tot_vol = 0.0;
143        // Loop over core weight points
144        for(size_t i=0; i<weights_radius.size(); i++) {
145                dp.radius = weights_radius[i].value;
146                // Loop over thick_inter0 weight points
147                for(size_t j=0; j<weights_edge_separation.size(); j++) {
148                        dp.edge_separation = weights_edge_separation[j].value;
149                        pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
150                        string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
151                        tot_vol = (dp.num_pearls - 1.0) * string_vol;
152                        tot_vol += dp.num_pearls * pearl_vol;
153                        //Find  volume
154                        // This may be a too much approximation
155                        //Todo: decided whether or not we keep this calculation
156                        sum += weights_radius[i].weight * weights_edge_separation[j].weight
157                                * pow(3.0*tot_vol/4.0/pi,0.333333);
158                        norm += weights_radius[i].weight * weights_edge_separation[j].weight;
159                }
160        }
161
162        if (norm != 0){
163                //return the averaged value
164                rad_out =  sum/norm;}
165        else{
166                //return normal value
167                pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
168                string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
169                tot_vol = (dp.num_pearls - 1.0) * string_vol;
170                tot_vol += dp.num_pearls * pearl_vol;
171
172                rad_out =  pow((3.0*tot_vol/4.0/pi), 0.33333);
173                }
174
175        return rad_out;
176
177}
Note: See TracBrowser for help on using the repository browser.