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 | // RKH 03Apr2014 a re-parametrised CoreShellEllipsoid with core axial ratio X and shell thickness T |
---|
22 | #include <math.h> |
---|
23 | #include "parameters.hh" |
---|
24 | #include <stdio.h> |
---|
25 | #include <stdlib.h> |
---|
26 | using namespace std; |
---|
27 | #include "spheroidXT.h" |
---|
28 | |
---|
29 | extern "C" { |
---|
30 | #include "libCylinder.h" |
---|
31 | #include "libStructureFactor.h" |
---|
32 | } |
---|
33 | |
---|
34 | typedef struct { |
---|
35 | double scale; |
---|
36 | double equat_core; |
---|
37 | double X_core; |
---|
38 | double T_shell; |
---|
39 | double XpolarShell; |
---|
40 | double sld_core; |
---|
41 | double sld_shell; |
---|
42 | double sld_solvent; |
---|
43 | double background; |
---|
44 | double axis_theta; |
---|
45 | double axis_phi; |
---|
46 | |
---|
47 | } SpheroidXTParameters; |
---|
48 | |
---|
49 | /** |
---|
50 | * Function to evaluate 2D scattering function |
---|
51 | * @param pars: parameters of the prolate |
---|
52 | * @param q: q-value |
---|
53 | * @param q_x: q_x / q |
---|
54 | * @param q_y: q_y / q |
---|
55 | * @return: function value |
---|
56 | */ |
---|
57 | static double spheroidXT_analytical_2D_scaled(SpheroidXTParameters *pars, double q, double q_x, double q_y) { |
---|
58 | |
---|
59 | double cyl_x, cyl_y;//, cyl_z; |
---|
60 | //double q_z; |
---|
61 | double alpha, vol, cos_val; |
---|
62 | double answer; |
---|
63 | double Pi = 4.0*atan(1.0); |
---|
64 | double sldcs,sldss; |
---|
65 | |
---|
66 | //convert angle degree to radian |
---|
67 | double theta = pars->axis_theta * Pi/180.0; |
---|
68 | double phi = pars->axis_phi * Pi/180.0; |
---|
69 | |
---|
70 | |
---|
71 | // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis. |
---|
72 | cyl_x = cos(theta) * cos(phi); |
---|
73 | cyl_y = sin(theta); |
---|
74 | //cyl_z = -cos(theta) * sin(phi); |
---|
75 | //del sld |
---|
76 | sldcs = pars->sld_core - pars->sld_shell; |
---|
77 | sldss = pars->sld_shell- pars->sld_solvent; |
---|
78 | |
---|
79 | // q vector |
---|
80 | //q_z = 0; |
---|
81 | |
---|
82 | // Compute the angle btw vector q and the |
---|
83 | // axis of the cylinder |
---|
84 | cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; |
---|
85 | |
---|
86 | // The following test should always pass |
---|
87 | if (fabs(cos_val)>1.0) { |
---|
88 | printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n"); |
---|
89 | return 0; |
---|
90 | } |
---|
91 | |
---|
92 | // Note: cos(alpha) = 0 and 1 will get an |
---|
93 | // undefined value from CylKernel |
---|
94 | alpha = acos( cos_val ); |
---|
95 | |
---|
96 | // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params. |
---|
97 | // was answer = gfn4(cos_val,pars->equat_core,pars->polar_core,pars->equat_shell,pars->polar_shell,sldcs,sldss,q); |
---|
98 | answer = gfn4(cos_val,pars->equat_core, pars->equat_core*pars->X_core, pars->equat_core + pars->T_shell, pars->equat_core*pars->X_core + pars->T_shell*pars->XpolarShell ,sldcs,sldss,q); |
---|
99 | //It seems that it should be normalized somehow. How??? |
---|
100 | |
---|
101 | //normalize by cylinder volume |
---|
102 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
103 | // was vol = 4.0*Pi/3.0*pars->equat_shell*pars->equat_shell*pars->polar_shell; |
---|
104 | vol = 4.0*Pi/3.0*(pars->equat_core + pars->T_shell)* (pars->equat_core + pars->T_shell) * ( pars->equat_core*pars->X_core + pars->T_shell*pars->XpolarShell); |
---|
105 | answer /= vol; |
---|
106 | |
---|
107 | //convert to [cm-1] |
---|
108 | answer *= 1.0e8; |
---|
109 | |
---|
110 | //Scale |
---|
111 | answer *= pars->scale; |
---|
112 | |
---|
113 | // add in the background |
---|
114 | answer += pars->background; |
---|
115 | |
---|
116 | return answer; |
---|
117 | } |
---|
118 | |
---|
119 | CoreShellEllipsoidXTModel :: CoreShellEllipsoidXTModel() { |
---|
120 | scale = Parameter(1.0); |
---|
121 | equat_core = Parameter(20.0, true); |
---|
122 | equat_core.set_min(0.0); |
---|
123 | X_core = Parameter(3.0, true); |
---|
124 | X_core.set_min(0.001); |
---|
125 | T_shell = Parameter(30.0, true); |
---|
126 | T_shell.set_min(0.0); |
---|
127 | XpolarShell = Parameter(1.0, true); |
---|
128 | XpolarShell.set_min(0.0); |
---|
129 | sld_core = Parameter(2e-6); |
---|
130 | sld_shell = Parameter(1e-6); |
---|
131 | sld_solvent = Parameter(6.3e-6); |
---|
132 | background = Parameter(0.0); |
---|
133 | axis_theta = Parameter(0.0, true); |
---|
134 | axis_phi = Parameter(0.0, true); |
---|
135 | |
---|
136 | } |
---|
137 | |
---|
138 | |
---|
139 | /** |
---|
140 | * Function to evaluate 2D scattering function |
---|
141 | * @param pars: parameters of the prolate |
---|
142 | * @param q: q-value |
---|
143 | * @return: function value |
---|
144 | */ |
---|
145 | |
---|
146 | static double spheroidXT_analytical_2DXY(SpheroidXTParameters *pars, double qx, double qy) { |
---|
147 | double q; |
---|
148 | q = sqrt(qx*qx+qy*qy); |
---|
149 | return spheroidXT_analytical_2D_scaled(pars, q, qx/q, qy/q); |
---|
150 | } |
---|
151 | |
---|
152 | /** |
---|
153 | * Function to evaluate 1D scattering function |
---|
154 | * The NIST IGOR library is used for the actual calculation. |
---|
155 | * @param q: q-value |
---|
156 | * @return: function value |
---|
157 | */ |
---|
158 | double CoreShellEllipsoidXTModel :: operator()(double q) { |
---|
159 | double dp[9]; |
---|
160 | |
---|
161 | // Fill parameter array for IGOR library |
---|
162 | // Add the background after averaging |
---|
163 | dp[0] = scale(); |
---|
164 | dp[1] = equat_core(); |
---|
165 | //dp[2] = polar_core(); |
---|
166 | //dp[3] = equat_shell(); |
---|
167 | //dp[4] = polar_shell(); |
---|
168 | dp[2] = equat_core()*X_core(); |
---|
169 | dp[3] = equat_core() + T_shell(); |
---|
170 | dp[4] = equat_core()*X_core() + T_shell()*XpolarShell(); |
---|
171 | dp[5] = sld_core(); |
---|
172 | dp[6] = sld_shell(); |
---|
173 | dp[7] = sld_solvent(); |
---|
174 | dp[8] = 0.0; |
---|
175 | |
---|
176 | // Get the dispersion points for the major core |
---|
177 | vector<WeightPoint> weights_equat_core; |
---|
178 | equat_core.get_weights(weights_equat_core); |
---|
179 | |
---|
180 | // Get the dispersion points |
---|
181 | vector<WeightPoint> weights_X_core; |
---|
182 | X_core.get_weights(weights_X_core); |
---|
183 | |
---|
184 | // Get the dispersion points |
---|
185 | vector<WeightPoint> weights_T_shell; |
---|
186 | T_shell.get_weights(weights_T_shell); |
---|
187 | |
---|
188 | // Get the dispersion points |
---|
189 | vector<WeightPoint> weights_XpolarShell; |
---|
190 | XpolarShell.get_weights(weights_XpolarShell); |
---|
191 | |
---|
192 | |
---|
193 | // Perform the computation, with all weight points |
---|
194 | double sum = 0.0; |
---|
195 | double norm = 0.0; |
---|
196 | double vol = 0.0; |
---|
197 | double wproduct = 0.0; |
---|
198 | |
---|
199 | // Loop over equat core weight points |
---|
200 | for(int i=0; i<(int)weights_equat_core.size(); i++) { |
---|
201 | dp[1] = weights_equat_core[i].value; |
---|
202 | |
---|
203 | // Loop over polar core weight points |
---|
204 | for(int j=0; j<(int)weights_X_core.size(); j++) { |
---|
205 | dp[2] = dp[1]*weights_X_core[j].value; |
---|
206 | |
---|
207 | // Loop over equat outer weight points |
---|
208 | for(int k=0; k<(int)weights_T_shell.size(); k++) { |
---|
209 | dp[3] = dp[1] + weights_T_shell[k].value; |
---|
210 | |
---|
211 | // Loop over polar outer weight points |
---|
212 | for(int l=0; l<(int)weights_XpolarShell.size(); l++) { |
---|
213 | dp[4] = dp[2] + weights_XpolarShell[l].value*weights_T_shell[k].value; |
---|
214 | |
---|
215 | //Un-normalize by volume |
---|
216 | wproduct =weights_equat_core[i].weight* weights_X_core[j].weight * weights_T_shell[k].weight |
---|
217 | * weights_XpolarShell[l].weight; |
---|
218 | sum += wproduct * OblateForm(dp, q) * dp[3]*dp[3]*dp[4]; |
---|
219 | //Find average volume |
---|
220 | vol += wproduct * dp[3]*dp[3]*dp[4]; |
---|
221 | // was pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value; |
---|
222 | norm += wproduct; |
---|
223 | } |
---|
224 | } |
---|
225 | } |
---|
226 | } |
---|
227 | if (vol != 0.0 && norm != 0.0) { |
---|
228 | //Re-normalize by avg volume |
---|
229 | sum = sum/(vol/norm);} |
---|
230 | return sum/norm + background(); |
---|
231 | } |
---|
232 | |
---|
233 | /** |
---|
234 | * Function to evaluate 2D scattering function |
---|
235 | * @param q_x: value of Q along x |
---|
236 | * @param q_y: value of Q along y |
---|
237 | * @return: function value |
---|
238 | */ |
---|
239 | /* |
---|
240 | double OblateModel :: operator()(double qx, double qy) { |
---|
241 | double q = sqrt(qx*qx + qy*qy); |
---|
242 | |
---|
243 | return (*this).operator()(q); |
---|
244 | } |
---|
245 | */ |
---|
246 | |
---|
247 | /** |
---|
248 | * Function to evaluate 2D scattering function |
---|
249 | * @param pars: parameters of the oblate |
---|
250 | * @param q: q-value |
---|
251 | * @param phi: angle phi |
---|
252 | * @return: function value |
---|
253 | */ |
---|
254 | double CoreShellEllipsoidXTModel :: evaluate_rphi(double q, double phi) { |
---|
255 | double qx = q*cos(phi); |
---|
256 | double qy = q*sin(phi); |
---|
257 | return (*this).operator()(qx, qy); |
---|
258 | } |
---|
259 | |
---|
260 | /** |
---|
261 | * Function to evaluate 2D scattering function |
---|
262 | * @param q_x: value of Q along x |
---|
263 | * @param q_y: value of Q along y |
---|
264 | * @return: function value |
---|
265 | */ |
---|
266 | double CoreShellEllipsoidXTModel :: operator()(double qx, double qy) { |
---|
267 | SpheroidXTParameters dp; |
---|
268 | // Fill parameter array |
---|
269 | dp.scale = scale(); |
---|
270 | dp.equat_core = equat_core(); |
---|
271 | dp.X_core = X_core(); |
---|
272 | dp.T_shell = T_shell(); |
---|
273 | dp.XpolarShell = XpolarShell(); |
---|
274 | dp.sld_core = sld_core(); |
---|
275 | dp.sld_shell = sld_shell(); |
---|
276 | dp.sld_solvent = sld_solvent(); |
---|
277 | dp.background = 0.0; |
---|
278 | dp.axis_theta = axis_theta(); |
---|
279 | dp.axis_phi = axis_phi(); |
---|
280 | |
---|
281 | // Get the dispersion points for the major core |
---|
282 | vector<WeightPoint> weights_equat_core; |
---|
283 | equat_core.get_weights(weights_equat_core); |
---|
284 | |
---|
285 | // Get the dispersion points |
---|
286 | vector<WeightPoint> weights_X_core; |
---|
287 | X_core.get_weights(weights_X_core); |
---|
288 | |
---|
289 | // Get the dispersion points |
---|
290 | vector<WeightPoint> weights_T_shell; |
---|
291 | T_shell.get_weights(weights_T_shell); |
---|
292 | |
---|
293 | // Get the dispersion points |
---|
294 | vector<WeightPoint> weights_XpolarShell; |
---|
295 | XpolarShell.get_weights(weights_XpolarShell); |
---|
296 | |
---|
297 | |
---|
298 | // Get angular averaging for theta |
---|
299 | vector<WeightPoint> weights_theta; |
---|
300 | axis_theta.get_weights(weights_theta); |
---|
301 | |
---|
302 | // Get angular averaging for phi |
---|
303 | vector<WeightPoint> weights_phi; |
---|
304 | axis_phi.get_weights(weights_phi); |
---|
305 | |
---|
306 | // Perform the computation, with all weight points |
---|
307 | double sum = 0.0; |
---|
308 | double norm = 0.0; |
---|
309 | double norm_vol = 0.0; |
---|
310 | double vol = 0.0; |
---|
311 | double pi = 4.0*atan(1.0); |
---|
312 | double equat_outer = 0.0; |
---|
313 | double polar_outer = 0.0; |
---|
314 | |
---|
315 | // Loop over major core weight points |
---|
316 | for(int i=0; i< (int)weights_equat_core.size(); i++) { |
---|
317 | dp.equat_core = weights_equat_core[i].value; |
---|
318 | |
---|
319 | // Loop over minor core weight points |
---|
320 | for(int j=0; j< (int)weights_X_core.size(); j++) { |
---|
321 | dp.X_core = weights_X_core[j].value; |
---|
322 | |
---|
323 | // Loop over equat outer weight points |
---|
324 | for(int k=0; k< (int)weights_T_shell.size(); k++) { |
---|
325 | dp.T_shell = weights_T_shell[k].value; |
---|
326 | equat_outer = weights_equat_core[i].value + weights_T_shell[k].value; |
---|
327 | |
---|
328 | // Loop over polar outer weight points |
---|
329 | for(int l=0; l< (int)weights_XpolarShell.size(); l++) { |
---|
330 | dp.XpolarShell = weights_XpolarShell[l].value; |
---|
331 | polar_outer = weights_equat_core[i].value*weights_X_core[j].value + weights_T_shell[k].value*weights_XpolarShell[l].value; |
---|
332 | |
---|
333 | // Average over theta distribution |
---|
334 | for(int m=0; m< (int)weights_theta.size(); m++) { |
---|
335 | dp.axis_theta = weights_theta[m].value; |
---|
336 | |
---|
337 | // Average over phi distribution |
---|
338 | for(int n=0; n< (int)weights_phi.size(); n++) { |
---|
339 | dp.axis_phi = weights_phi[n].value; |
---|
340 | //Un-normalize by volume |
---|
341 | double _ptvalue = weights_equat_core[i].weight *weights_X_core[j].weight |
---|
342 | * weights_T_shell[k].weight * weights_XpolarShell[l].weight |
---|
343 | * weights_theta[m].weight |
---|
344 | * weights_phi[n].weight |
---|
345 | // rkh NOTE this passes the NEW parameters |
---|
346 | * spheroidXT_analytical_2DXY(&dp, qx, qy) * pow(equat_outer,2)*polar_outer; |
---|
347 | if (weights_theta.size()>1) { |
---|
348 | _ptvalue *= fabs(cos(weights_theta[m].value*pi/180.0)); |
---|
349 | } |
---|
350 | sum += _ptvalue; |
---|
351 | //Find average volume |
---|
352 | // rkh had to change this, original weighted by outer shell volume weights only, see spheroid.cpp, which looks odd, |
---|
353 | // (as has to assume that weights of other loops sume to unity) and here we need all four loops to get the outer size. |
---|
354 | vol += weights_equat_core[i].weight *weights_X_core[j].weight |
---|
355 | * weights_T_shell[k].weight * weights_XpolarShell[l].weight |
---|
356 | * pow(equat_outer,2)*polar_outer; |
---|
357 | //Find norm for volume |
---|
358 | norm_vol += weights_equat_core[i].weight *weights_X_core[j].weight |
---|
359 | * weights_T_shell[k].weight * weights_XpolarShell[l].weight; |
---|
360 | |
---|
361 | norm += weights_equat_core[i].weight *weights_X_core[j].weight |
---|
362 | * weights_T_shell[k].weight * weights_XpolarShell[l].weight |
---|
363 | * weights_theta[m].weight * weights_phi[n].weight; |
---|
364 | } |
---|
365 | } |
---|
366 | } |
---|
367 | } |
---|
368 | } |
---|
369 | } |
---|
370 | // Averaging in theta needs an extra normalization |
---|
371 | // factor to account for the sin(theta) term in the |
---|
372 | // integration (see documentation). |
---|
373 | if (weights_theta.size()>1) norm = norm / asin(1.0); |
---|
374 | |
---|
375 | if (vol != 0.0 && norm_vol != 0.0) { |
---|
376 | //Re-normalize by avg volume |
---|
377 | sum = sum/(vol/norm_vol);} |
---|
378 | |
---|
379 | return sum/norm + background(); |
---|
380 | } |
---|
381 | |
---|
382 | /** |
---|
383 | * Function to calculate effective radius |
---|
384 | * @return: effective radius value |
---|
385 | * rkh This now needs to integrate over all four variables as above not just two |
---|
386 | */ |
---|
387 | double CoreShellEllipsoidXTModel :: calculate_ER() { |
---|
388 | SpheroidXTParameters dp; |
---|
389 | |
---|
390 | dp.equat_core = equat_core(); |
---|
391 | dp.X_core = X_core(); |
---|
392 | dp.T_shell = T_shell(); |
---|
393 | dp.XpolarShell = XpolarShell(); |
---|
394 | |
---|
395 | double rad_out = 0.0; |
---|
396 | double equat_outer = 0.0; |
---|
397 | double polar_outer = 0.0; |
---|
398 | |
---|
399 | // Perform the computation, with all weight points |
---|
400 | double sum = 0.0; |
---|
401 | double norm = 0.0; |
---|
402 | |
---|
403 | // Get the dispersion points for the core |
---|
404 | vector<WeightPoint> weights_equat_core; |
---|
405 | equat_core.get_weights(weights_equat_core); |
---|
406 | |
---|
407 | // Get the dispersion points |
---|
408 | vector<WeightPoint> weights_X_core; |
---|
409 | X_core.get_weights(weights_X_core); |
---|
410 | |
---|
411 | // Get the dispersion points |
---|
412 | vector<WeightPoint> weights_T_shell; |
---|
413 | T_shell.get_weights(weights_T_shell); |
---|
414 | |
---|
415 | // Get the dispersion points |
---|
416 | vector<WeightPoint> weights_XpolarShell; |
---|
417 | XpolarShell.get_weights(weights_XpolarShell); |
---|
418 | |
---|
419 | |
---|
420 | // Loop over core weight points |
---|
421 | for(int i=0; i< (int)weights_equat_core.size(); i++) { |
---|
422 | dp.equat_core = weights_equat_core[i].value; |
---|
423 | // Loop over weight points |
---|
424 | for(int j=0; j< (int)weights_X_core.size(); j++) { |
---|
425 | // Loop over weight points |
---|
426 | for(int k=0; k< (int)weights_T_shell.size(); k++) { |
---|
427 | equat_outer = weights_equat_core[i].value + weights_T_shell[k].value; |
---|
428 | // Loop over polar outer weight points |
---|
429 | for(int l=0; l< (int)weights_XpolarShell.size(); l++) { |
---|
430 | polar_outer = weights_equat_core[i].value*weights_X_core[j].value + weights_T_shell[k].value*weights_XpolarShell[l].value; |
---|
431 | //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER. |
---|
432 | sum +=weights_equat_core[i].weight *weights_X_core[j].weight |
---|
433 | * weights_T_shell[k].weight * weights_XpolarShell[l].weight*DiamEllip(polar_outer,equat_outer)/2.0; |
---|
434 | norm += weights_equat_core[i].weight *weights_X_core[j].weight |
---|
435 | * weights_T_shell[k].weight * weights_XpolarShell[l].weight; |
---|
436 | } |
---|
437 | } |
---|
438 | } |
---|
439 | } |
---|
440 | if (norm != 0){ |
---|
441 | //return the averaged value |
---|
442 | rad_out = sum/norm;} |
---|
443 | else{ |
---|
444 | //return normal value |
---|
445 | //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER. |
---|
446 | rad_out = DiamEllip(dp.equat_core + dp.T_shell, dp.equat_core * dp.X_core + dp.T_shell*dp.XpolarShell)/2.0;} |
---|
447 | |
---|
448 | return rad_out; |
---|
449 | } |
---|
450 | double CoreShellEllipsoidXTModel :: calculate_VR() { |
---|
451 | return 1.0; |
---|
452 | } |
---|