parallelepiped.cpp @ 97603c0

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 97603c0 was c451be9, checked in by Jae Cho <jhjcho@…>, 15 years ago
corrections on the definition of polydispersity as suggested by steve K: should be normalized by average volume
Property mode set to `100644`
File size: 9.6 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	* TODO: refactor so that we pull in the old sansmodels.c_extensions
21	* TODO: add 2D function
22	*/
23
24	#include <math.h>
25	#include "models.hh"
26	#include "parameters.hh"
27	#include <stdio.h>
28	using namespace std;
29
30	extern "C" {
31	#include "libCylinder.h"
32	#include "libStructureFactor.h"
33	#include "parallelepiped.h"
34	}
35
36	ParallelepipedModel :: ParallelepipedModel() {
37	scale = Parameter(1.0);
38	short_a = Parameter(35.0, true);
39	short_a.set_min(1.0);
40	short_b = Parameter(75.0, true);
41	short_b.set_min(1.0);
42	long_c = Parameter(400.0, true);
43	long_c.set_min(1.0);
44	contrast = Parameter(53.e-7);
45	background = Parameter(0.0);
46	parallel_theta = Parameter(0.0, true);
47	parallel_phi = Parameter(0.0, true);
48	parallel_psi = Parameter(0.0, true);
49	}
50
51	/**
52	* Function to evaluate 1D scattering function
53	* The NIST IGOR library is used for the actual calculation.
54	* @param q: q-value
55	* @return: function value
56	*/
57	double ParallelepipedModel :: operator()(double q) {
58	double dp[6];
59
60	// Fill parameter array for IGOR library
61	// Add the background after averaging
62	dp[0] = scale();
63	dp[1] = short_a();
64	dp[2] = short_b();
65	dp[3] = long_c();
66	dp[4] = contrast();
67	dp[5] = 0.0;
68
69	// Get the dispersion points for the short_edgeA
70	vector<WeightPoint> weights_short_a;
71	short_a.get_weights(weights_short_a);
72
73	// Get the dispersion points for the longer_edgeB
74	vector<WeightPoint> weights_short_b;
75	short_b.get_weights(weights_short_b);
76
77	// Get the dispersion points for the longuest_edgeC
78	vector<WeightPoint> weights_long_c;
79	long_c.get_weights(weights_long_c);
80
81
82
83	// Perform the computation, with all weight points
84	double sum = 0.0;
85	double norm = 0.0;
86	double vol = 0.0;
87
88	// Loop over short_edgeA weight points
89	for(int i=0; i< (int)weights_short_a.size(); i++) {
90	dp[1] = weights_short_a[i].value;
91
92	// Loop over longer_edgeB weight points
93	for(int j=0; j< (int)weights_short_b.size(); j++) {
94	dp[2] = weights_short_b[j].value;
95
96	// Loop over longuest_edgeC weight points
97	for(int k=0; k< (int)weights_long_c.size(); k++) {
98	dp[3] = weights_long_c[k].value;
99	//Un-normalize by volume
100	sum += weights_short_a[i].weight * weights_short_b[j].weight
101	* weights_long_c[k].weight * Parallelepiped(dp, q)
102	* weights_short_a[i].value*weights_short_b[j].value
103	* weights_long_c[k].value;
104	//Find average volume
105	vol += weights_short_a[i].weight * weights_short_b[j].weight
106	* weights_long_c[k].weight
107	* weights_short_a[i].value * weights_short_b[j].value
108	* weights_long_c[k].value;
109
110	norm += weights_short_a[i].weight
111	* weights_short_b[j].weight * weights_long_c[k].weight;
112	}
113	}
114	}
115	if (vol != 0.0 && norm != 0.0) {
116	//Re-normalize by avg volume
117	sum = sum/(vol/norm);}
118
119	return sum/norm + background();
120	}
121	/**
122	* Function to evaluate 2D scattering function
123	* @param q_x: value of Q along x
124	* @param q_y: value of Q along y
125	* @return: function value
126	*/
127	double ParallelepipedModel :: operator()(double qx, double qy) {
128	ParallelepipedParameters dp;
129	// Fill parameter array
130	dp.scale = scale();
131	dp.short_a = short_a();
132	dp.short_b = short_b();
133	dp.long_c = long_c();
134	dp.contrast = contrast();
135	dp.background = 0.0;
136	//dp.background = background();
137	dp.parallel_theta = parallel_theta();
138	dp.parallel_phi = parallel_phi();
139	dp.parallel_psi = parallel_psi();
140
141
142	// Get the dispersion points for the short_edgeA
143	vector<WeightPoint> weights_short_a;
144	short_a.get_weights(weights_short_a);
145
146	// Get the dispersion points for the longer_edgeB
147	vector<WeightPoint> weights_short_b;
148	short_b.get_weights(weights_short_b);
149
150	// Get angular averaging for the longuest_edgeC
151	vector<WeightPoint> weights_long_c;
152	long_c.get_weights(weights_long_c);
153
154	// Get angular averaging for theta
155	vector<WeightPoint> weights_parallel_theta;
156	parallel_theta.get_weights(weights_parallel_theta);
157
158	// Get angular averaging for phi
159	vector<WeightPoint> weights_parallel_phi;
160	parallel_phi.get_weights(weights_parallel_phi);
161
162	// Get angular averaging for psi
163	vector<WeightPoint> weights_parallel_psi;
164	parallel_psi.get_weights(weights_parallel_psi);
165
166	// Perform the computation, with all weight points
167	double sum = 0.0;
168	double norm = 0.0;
169	double norm_vol = 0.0;
170	double vol = 0.0;
171
172	// Loop over radius weight points
173	for(int i=0; i< (int)weights_short_a.size(); i++) {
174	dp.short_a = weights_short_a[i].value;
175
176	// Loop over longer_edgeB weight points
177	for(int j=0; j< (int)weights_short_b.size(); j++) {
178	dp.short_b = weights_short_b[j].value;
179
180	// Average over longuest_edgeC distribution
181	for(int k=0; k< (int)weights_long_c.size(); k++) {
182	dp.long_c = weights_long_c[k].value;
183
184	// Average over theta distribution
185	for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
186	dp.parallel_theta = weights_parallel_theta[l].value;
187
188	// Average over phi distribution
189	for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
190	dp.parallel_phi = weights_parallel_phi[m].value;
191
192	// Average over phi distribution
193	for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
194	dp.parallel_psi = weights_parallel_psi[n].value;
195	//Un-normalize by volume
196	double _ptvalue = weights_short_a[i].weight
197	* weights_short_b[j].weight
198	* weights_long_c[k].weight
199	* weights_parallel_theta[l].weight
200	* weights_parallel_phi[m].weight
201	* weights_parallel_psi[n].weight
202	* parallelepiped_analytical_2DXY(&dp, qx, qy)
203	* weights_short_a[i].value*weights_short_b[j].value
204	* weights_long_c[k].value;
205
206	if (weights_parallel_theta.size()>1) {
207	_ptvalue *= sin(weights_parallel_theta[l].value);
208	}
209	sum += _ptvalue;
210	//Find average volume
211	vol += weights_short_a[i].weight
212	* weights_short_b[j].weight
213	* weights_long_c[k].weight
214	* weights_short_a[i].value*weights_short_b[j].value
215	* weights_long_c[k].value;
216	//Find norm for volume
217	norm_vol += weights_short_a[i].weight
218	* weights_short_b[j].weight
219	* weights_long_c[k].weight;
220
221	norm += weights_short_a[i].weight
222	* weights_short_b[j].weight
223	* weights_long_c[k].weight
224	* weights_parallel_theta[l].weight
225	* weights_parallel_phi[m].weight
226	* weights_parallel_psi[n].weight;
227	}
228	}
229
230	}
231	}
232	}
233	}
234	// Averaging in theta needs an extra normalization
235	// factor to account for the sin(theta) term in the
236	// integration (see documentation).
237	if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
238
239	if (vol != 0.0 && norm_vol != 0.0) {
240	//Re-normalize by avg volume
241	sum = sum/(vol/norm_vol);}
242
243	return sum/norm + background();
244	}
245
246
247	/**
248	* Function to evaluate 2D scattering function
249	* @param pars: parameters of the cylinder
250	* @param q: q-value
251	* @param phi: angle phi
252	* @return: function value
253	*/
254	double ParallelepipedModel :: 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	* Function to calculate effective radius
261	* @return: effective radius value
262	*/
263	double ParallelepipedModel :: calculate_ER() {
264	ParallelepipedParameters dp;
265	dp.short_a = short_a();
266	dp.short_b = short_b();
267	dp.long_c = long_c();
268	double rad_out = 0.0;
269	double pi = 4.0*atan(1.0);
270	double suf_rad = sqrt(dp.short_a*dp.short_b/pi);
271
272	// Perform the computation, with all weight points
273	double sum = 0.0;
274	double norm = 0.0;
275
276	// Get the dispersion points for the short_edgeA
277	vector<WeightPoint> weights_short_a;
278	short_a.get_weights(weights_short_a);
279
280	// Get the dispersion points for the longer_edgeB
281	vector<WeightPoint> weights_short_b;
282	short_b.get_weights(weights_short_b);
283
284	// Get angular averaging for the longuest_edgeC
285	vector<WeightPoint> weights_long_c;
286	long_c.get_weights(weights_long_c);
287
288	// Loop over radius weight points
289	for(int i=0; i< (int)weights_short_a.size(); i++) {
290	dp.short_a = weights_short_a[i].value;
291
292	// Loop over longer_edgeB weight points
293	for(int j=0; j< (int)weights_short_b.size(); j++) {
294	dp.short_b = weights_short_b[j].value;
295
296	// Average over longuest_edgeC distribution
297	for(int k=0; k< (int)weights_long_c.size(); k++) {
298	dp.long_c = weights_long_c[k].value;
299	//Calculate surface averaged radius
300	//This is rough approximation.
301	suf_rad = sqrt(dp.short_a*dp.short_b/pi);
302
303	//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
304	sum +=weights_short_a[i].weight* weights_short_b[j].weight
305	* weights_long_c[k].weight*DiamCyl(dp.long_c, suf_rad)/2.0;
306	norm += weights_short_a[i].weight* weights_short_b[j].weight*weights_long_c[k].weight;
307	}
308	}
309	}
310
311	if (norm != 0){
312	//return the averaged value
313	rad_out = sum/norm;}
314	else{
315	//return normal value
316	//Note: output of "DiamCyl(length,radius)" is DIAMETER.
317	rad_out = DiamCyl(dp.long_c, suf_rad)/2.0;}
318	return rad_out;
319
320	}

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/sansmodels/src/sans/models/c_models/parallelepiped.cpp @ 97603c0

Download in other formats: