source: sasview/src/sas/models/c_extension/c_models/pearlnecklace.cpp @ fc18690

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 fc18690 was 422e7e8, checked in by Paul Kienzle <pkienzle@…>, 9 years ago

allow non-integer values for number of pearls in model by rounding to the nearest integer in the computation.

  • Property mode set to 100644
File size: 8.5 KB
Line 
1
2#include <math.h>
3#include "parameters.hh"
4#include <stdio.h>
5using namespace std;
6#include "pearlnecklace.h"
7
8extern "C" {
9#include "libmultifunc/libfunc.h"
10}
11
12static double pearl_necklace_kernel(double dp[], double q) {
13  // fit parameters
14  double scale = dp[0];
15  double radius = dp[1];
16  double edge_separation = dp[2];
17  double thick_string = dp[3];
18  double num_pearls = floor(dp[4]+0.5);  // Force integer number of pearls
19  double sld_pearl = dp[5];
20  double sld_string = dp[6];
21  double sld_solv = dp[7];
22  double background = dp[8];
23
24  //relative slds
25  double contrast_pearl = sld_pearl - sld_solv;
26  double contrast_string = sld_string - sld_solv;
27
28  // number of string segments
29  double num_strings = num_pearls - 1.0;
30
31  //Pi
32  double pi = 4.0*atan(1.0);
33
34  // each volumes
35  double string_vol = edge_separation * pi * pow((thick_string / 2.0), 2);
36  double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3);
37
38  //total volume
39  double tot_vol;
40  //each masses
41  double m_r= contrast_string * string_vol;
42  double m_s= contrast_pearl * pearl_vol;
43  double psi;
44  double gamma;
45  double beta;
46  //form factors
47  double sss; //pearls
48  double srr; //strings
49  double srs; //cross
50  double A_s;
51  double srr_1;
52  double srr_2;
53  double srr_3;
54  double form_factor;
55
56  tot_vol = num_strings * string_vol;
57  tot_vol += num_pearls * pearl_vol;
58
59  //sine functions of a pearl
60  psi = sin(q*radius);
61  psi -= q * radius * cos(q * radius);
62  psi /= pow((q * radius), 3);
63  psi *= 3.0;
64
65  // Note take only 20 terms in Si series: 10 terms may be enough though.
66  gamma = Si(q* edge_separation);
67  gamma /= (q* edge_separation);
68  beta = Si(q * (edge_separation + radius));
69  beta -= Si(q * radius);
70  beta /= (q* edge_separation);
71
72  // center to center distance between the neighboring pearls
73  A_s = edge_separation + 2.0 * radius;
74
75  // form factor for num_pearls
76  sss = 1.0 - pow(sinc(q*A_s), num_pearls );
77  sss /= pow((1.0-sinc(q*A_s)), 2);
78  sss *= -sinc(q*A_s);
79  sss -= num_pearls/2.0;
80  sss += num_pearls/(1.0-sinc(q*A_s));
81  sss *= 2.0 * pow((m_s*psi), 2);
82
83  // form factor for num_strings (like thin rods)
84  srr_1 = -pow(sinc(q*edge_separation/2.0), 2);
85
86  srr_1 += 2.0 * gamma;
87  srr_1 *= num_strings;
88  srr_2 = 2.0/(1.0-sinc(q*A_s));
89  srr_2 *= num_strings;
90  srr_2 *= pow(beta, 2);
91  srr_3 = 1.0 - pow(sinc(q*A_s), num_strings);
92  srr_3 /= pow((1.0-sinc(q*A_s)), 2);
93  srr_3 *= pow(beta, 2);
94  srr_3 *= -2.0;
95
96  // total srr
97  srr = srr_1 + srr_2 + srr_3;
98  srr *= pow(m_r, 2);
99
100  // form factor for correlations
101  srs = 1.0;
102  srs -= pow(sinc(q*A_s), num_strings);
103  srs /= pow((1.0-sinc(q*A_s)), 2);
104  srs *= -sinc(q*A_s);
105  srs += (num_strings/(1.0-sinc(q*A_s)));
106  srs *= 4.0;
107  srs *= (m_r * m_s * beta * psi);
108
109  form_factor = sss + srr + srs;
110  form_factor /= (tot_vol * 1.0e-8); // norm by volume and A^-1 to cm^-1
111
112  // scale and background
113  form_factor *= scale;
114  form_factor += background;
115  return (form_factor);
116}
117
118PearlNecklaceModel :: PearlNecklaceModel() {
119  scale = Parameter(1.0);
120  radius = Parameter(80.0, true);
121  radius.set_min(0.0);
122  edge_separation = Parameter(350.0, true);
123  edge_separation.set_min(0.0);
124  thick_string = Parameter(2.5, true);
125  thick_string.set_min(0.0);
126  num_pearls = Parameter(3);
127  num_pearls.set_min(0.0);
128  sld_pearl = Parameter(1.0e-06);
129  sld_string = Parameter(1.0e-06);
130  sld_solv = Parameter(6.3e-06);
131  background = Parameter(0.0);
132
133}
134
135/**
136 * Function to evaluate 1D PearlNecklaceModel function
137 * @param q: q-value
138 * @return: function value
139 */
140double PearlNecklaceModel :: operator()(double q) {
141  double dp[9];
142  // Fill parameter array for IGOR library
143  // Add the background after averaging
144  dp[0] = scale();
145  dp[1] = radius();
146  dp[2] = edge_separation();
147  dp[3] = thick_string();
148  dp[4] = num_pearls();
149  dp[5] = sld_pearl();
150  dp[6] = sld_string();
151  dp[7] = sld_solv();
152  dp[8] = 0.0;
153  double pi = 4.0*atan(1.0);
154  // No polydispersion supported in this model.
155  // Get the dispersion points for the radius
156  vector<WeightPoint> weights_radius;
157  radius.get_weights(weights_radius);
158  vector<WeightPoint> weights_edge_separation;
159  edge_separation.get_weights(weights_edge_separation);
160  // Perform the computation, with all weight points
161  double sum = 0.0;
162  double norm = 0.0;
163  double vol = 0.0;
164  double string_vol = 0.0;
165  double pearl_vol = 0.0;
166  double tot_vol = 0.0;
167  // Loop over core weight points
168  for(size_t i=0; i<weights_radius.size(); i++) {
169    dp[1] = weights_radius[i].value;
170    // Loop over thick_inter0 weight points
171    for(size_t j=0; j<weights_edge_separation.size(); j++) {
172      dp[2] = weights_edge_separation[j].value;
173      pearl_vol = 4.0 /3.0 * pi * pow(dp[1], 3);
174      string_vol =dp[2] * pi * pow((dp[3] / 2.0), 2);
175      tot_vol = (dp[4] - 1.0) * string_vol;
176      tot_vol += dp[4] * pearl_vol;
177      //Un-normalize Sphere by volume
178      sum += weights_radius[i].weight * weights_edge_separation[j].weight
179          * pearl_necklace_kernel(dp,q) * tot_vol;
180      //Find average volume
181      vol += weights_radius[i].weight * weights_edge_separation[j].weight
182          * tot_vol;
183      norm += weights_radius[i].weight * weights_edge_separation[j].weight;
184    }
185  }
186
187  if (vol != 0.0 && norm != 0.0) {
188    //Re-normalize by avg volume
189    sum = sum/(vol/norm);}
190
191  return sum/norm + background();
192}
193
194/**
195 * Function to evaluate 2D PearlNecklaceModel function
196 * @param q_x: value of Q along x
197 * @param q_y: value of Q along y
198 * @return: function value
199 */
200double PearlNecklaceModel :: operator()(double qx, double qy) {
201  double q = sqrt(qx*qx + qy*qy);
202  return (*this).operator()(q);
203}
204
205/**
206 * Function to evaluate PearlNecklaceModel function
207 * @param pars: parameters of the PearlNecklaceModel
208 * @param q: q-value
209 * @param phi: angle phi
210 * @return: function value
211 */
212double PearlNecklaceModel :: evaluate_rphi(double q, double phi) {
213  return (*this).operator()(q);
214}
215
216/**
217 * Function to calculate TOTAL radius
218 * Todo: decide whether or not we keep this calculation
219 * @return: effective radius value
220 */
221// No polydispersion supported in this model.
222// Calculate max radius assumming max_radius = effective radius
223// Note that this max radius is not affected by sld of layer, sld of interface, or
224// sld of solvent.
225double PearlNecklaceModel :: calculate_ER() {
226  PeralNecklaceParameters dp;
227
228  dp.scale = scale();
229  dp.radius = radius();
230  dp.edge_separation = edge_separation();
231  dp.thick_string = thick_string();
232  dp.num_pearls = num_pearls();
233
234  double rad_out = 0.0;
235  // Perform the computation, with all weight points
236  double sum = 0.0;
237  double norm = 0.0;
238  double pi = 4.0*atan(1.0);
239  // No polydispersion supported in this model.
240  // Get the dispersion points for the radius
241  vector<WeightPoint> weights_radius;
242  radius.get_weights(weights_radius);
243  vector<WeightPoint> weights_edge_separation;
244  edge_separation.get_weights(weights_edge_separation);
245  // Perform the computation, with all weight points
246  double string_vol = 0.0;
247  double pearl_vol = 0.0;
248  double tot_vol = 0.0;
249  // Loop over core weight points
250  for(size_t i=0; i<weights_radius.size(); i++) {
251    dp.radius = weights_radius[i].value;
252    // Loop over thick_inter0 weight points
253    for(size_t j=0; j<weights_edge_separation.size(); j++) {
254      dp.edge_separation = weights_edge_separation[j].value;
255      pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
256      string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
257      tot_vol = (dp.num_pearls - 1.0) * string_vol;
258      tot_vol += dp.num_pearls * pearl_vol;
259      //Find  volume
260      // This may be a too much approximation
261      //Todo: decided whether or not we keep this calculation
262      sum += weights_radius[i].weight * weights_edge_separation[j].weight
263          * pow(3.0*tot_vol/4.0/pi,0.333333);
264      norm += weights_radius[i].weight * weights_edge_separation[j].weight;
265    }
266  }
267
268  if (norm != 0){
269    //return the averaged value
270    rad_out =  sum/norm;}
271  else{
272    //return normal value
273    pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
274    string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
275    tot_vol = (dp.num_pearls - 1.0) * string_vol;
276    tot_vol += dp.num_pearls * pearl_vol;
277
278    rad_out =  pow((3.0*tot_vol/4.0/pi), 0.33333);
279  }
280
281  return rad_out;
282
283}
284double PearlNecklaceModel :: calculate_VR() {
285  return 1.0;
286}
Note: See TracBrowser for help on using the repository browser.