parallelepiped.cpp @ 3eac3816

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

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

SasView

source: sasview/sansmodels/src/sans/models/c_models/parallelepiped.cpp @ 3eac3816

Download in other formats: