ellipticalcylinder.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 34c2649, checked in by Mathieu Doucet <doucetm@…>, 13 years ago
Re #4 Fixed warnings
Property mode set to `100644`
File size: 8.9 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	*/
22
23	#include <math.h>
24	#include "models.hh"
25	#include "parameters.hh"
26	#include <stdio.h>
27	using namespace std;
28
29	extern "C" {
30	#include "libCylinder.h"
31	#include "libStructureFactor.h"
32	#include "elliptical_cylinder.h"
33	}
34
35	EllipticalCylinderModel :: EllipticalCylinderModel() {
36	scale = Parameter(1.0);
37	r_minor = Parameter(20.0, true);
38	r_minor.set_min(0.0);
39	r_ratio = Parameter(1.5, true);
40	r_ratio.set_min(0.0);
41	length = Parameter(400.0, true);
42	length.set_min(0.0);
43	sldCyl = Parameter(4.e-6);
44	sldSolv = Parameter(1.e-6);
45	background = Parameter(0.0);
46	cyl_theta = Parameter(57.325, true);
47	cyl_phi = Parameter(0.0, true);
48	cyl_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 EllipticalCylinderModel :: operator()(double q) {
58	double dp[7];
59
60	dp[0] = scale();
61	dp[1] = r_minor();
62	dp[2] = r_ratio();
63	dp[3] = length();
64	dp[4] = sldCyl();
65	dp[5] = sldSolv();
66	dp[6] = 0.0;
67
68	// Get the dispersion points for the r_minor
69	vector<WeightPoint> weights_rad;
70	r_minor.get_weights(weights_rad);
71
72	// Get the dispersion points for the r_ratio
73	vector<WeightPoint> weights_rat;
74	r_ratio.get_weights(weights_rat);
75
76	// Get the dispersion points for the length
77	vector<WeightPoint> weights_len;
78	length.get_weights(weights_len);
79
80	// Perform the computation, with all weight points
81	double sum = 0.0;
82	double norm = 0.0;
83	double vol = 0.0;
84
85	// Loop over r_minor weight points
86	for(size_t i=0; i<weights_rad.size(); i++) {
87	dp[1] = weights_rad[i].value;
88
89	// Loop over r_ratio weight points
90	for(size_t j=0; j<weights_rat.size(); j++) {
91	dp[2] = weights_rat[j].value;
92
93	// Loop over length weight points
94	for(size_t k=0; k<weights_len.size(); k++) {
95	dp[3] = weights_len[k].value;
96	//Un-normalize by volume
97	sum += weights_rad[i].weight
98	* weights_len[k].weight
99	* weights_rat[j].weight
100	* EllipCyl20(dp, q)
101	* pow(weights_rad[i].value,2) * weights_rat[j].value
102	* weights_len[k].value;
103	//Find average volume
104	vol += weights_rad[i].weight
105	* weights_len[k].weight
106	* weights_rat[j].weight
107	* pow(weights_rad[i].value,2) * weights_rat[j].value
108	* weights_len[k].value;
109	norm += weights_rad[i].weight
110	* weights_len[k].weight
111	* weights_rat[j].weight;
112	}
113	}
114	}
115
116	if (vol != 0.0 && norm != 0.0) {
117	//Re-normalize by avg volume
118	sum = sum/(vol/norm);}
119
120	return sum/norm + background();
121	}
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 EllipticalCylinderModel :: operator()(double qx, double qy) {
130	EllipticalCylinderParameters dp;
131	// Fill parameter array
132	dp.scale = scale();
133	dp.r_minor = r_minor();
134	dp.r_ratio = r_ratio();
135	dp.length = length();
136	dp.sldCyl = sldCyl();
137	dp.sldSolv = sldSolv();
138	dp.background = 0.0;
139	dp.cyl_theta = cyl_theta();
140	dp.cyl_phi = cyl_phi();
141	dp.cyl_psi = cyl_psi();
142
143	// Get the dispersion points for the r_minor
144	vector<WeightPoint> weights_rad;
145	r_minor.get_weights(weights_rad);
146
147	// Get the dispersion points for the r_ratio
148	vector<WeightPoint> weights_rat;
149	r_ratio.get_weights(weights_rat);
150
151	// Get the dispersion points for the length
152	vector<WeightPoint> weights_len;
153	length.get_weights(weights_len);
154
155	// Get angular averaging for theta
156	vector<WeightPoint> weights_theta;
157	cyl_theta.get_weights(weights_theta);
158
159	// Get angular averaging for phi
160	vector<WeightPoint> weights_phi;
161	cyl_phi.get_weights(weights_phi);
162
163	// Get angular averaging for psi
164	vector<WeightPoint> weights_psi;
165	cyl_psi.get_weights(weights_psi);
166
167	// Perform the computation, with all weight points
168	double sum = 0.0;
169	double norm = 0.0;
170	double norm_vol = 0.0;
171	double vol = 0.0;
172	double pi = 4.0*atan(1.0);
173	// Loop over minor radius weight points
174	for(size_t i=0; i<weights_rad.size(); i++) {
175	dp.r_minor = weights_rad[i].value;
176
177
178	// Loop over length weight points
179	for(size_t j=0; j<weights_len.size(); j++) {
180	dp.length = weights_len[j].value;
181
182	// Loop over r_ration weight points
183	for(size_t m=0; m<weights_rat.size(); m++) {
184	dp.r_ratio = weights_rat[m].value;
185
186	// Average over theta distribution
187	for(size_t k=0; k<weights_theta.size(); k++) {
188	dp.cyl_theta = weights_theta[k].value;
189
190	// Average over phi distribution
191	for(size_t l=0; l<weights_phi.size(); l++) {
192	dp.cyl_phi = weights_phi[l].value;
193
194	// Average over phi distribution
195	for(size_t o=0; o<weights_psi.size(); o++) {
196	dp.cyl_psi = weights_psi[o].value;
197	//Un-normalize by volume
198	double _ptvalue = weights_rad[i].weight
199	* weights_len[j].weight
200	* weights_rat[m].weight
201	* weights_theta[k].weight
202	* weights_phi[l].weight
203	* weights_psi[o].weight
204	* elliptical_cylinder_analytical_2DXY(&dp, qx, qy)
205	* pow(weights_rad[i].value,2)
206	* weights_len[j].value
207	* weights_rat[m].value;
208	if (weights_theta.size()>1) {
209	_ptvalue = fabs(sin(weights_theta[k].valuepi/180.0));
210	}
211	sum += _ptvalue;
212	//Find average volume
213	vol += weights_rad[i].weight
214	* weights_len[j].weight
215	* weights_rat[m].weight
216	* pow(weights_rad[i].value,2)
217	* weights_len[j].value
218	* weights_rat[m].value;
219	//Find norm for volume
220	norm_vol += weights_rad[i].weight
221	* weights_len[j].weight
222	* weights_rat[m].weight;
223
224	norm += weights_rad[i].weight
225	* weights_len[j].weight
226	* weights_rat[m].weight
227	* weights_theta[k].weight
228	* weights_phi[l].weight
229	* weights_psi[o].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_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 EllipticalCylinderModel :: 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 EllipticalCylinderModel :: calculate_ER() {
267	EllipticalCylinderParameters dp;
268	dp.r_minor = r_minor();
269	dp.r_ratio = r_ratio();
270	dp.length = length();
271	double rad_out = 0.0;
272	double suf_rad = sqrt(dp.r_minordp.r_minordp.r_ratio);
273
274	// Perform the computation, with all weight points
275	double sum = 0.0;
276	double norm = 0.0;
277
278	// Get the dispersion points for the r_minor
279	vector<WeightPoint> weights_rad;
280	r_minor.get_weights(weights_rad);
281
282	// Get the dispersion points for the r_ratio
283	vector<WeightPoint> weights_rat;
284	r_ratio.get_weights(weights_rat);
285
286	// Get the dispersion points for the length
287	vector<WeightPoint> weights_len;
288	length.get_weights(weights_len);
289
290	// Loop over minor radius weight points
291	for(size_t i=0; i<weights_rad.size(); i++) {
292	dp.r_minor = weights_rad[i].value;
293
294	// Loop over length weight points
295	for(size_t j=0; j<weights_len.size(); j++) {
296	dp.length = weights_len[j].value;
297
298	// Loop over r_ration weight points
299	for(size_t m=0; m<weights_rat.size(); m++) {
300	dp.r_ratio = weights_rat[m].value;
301	//Calculate surface averaged radius
302	suf_rad = sqrt(dp.r_minor * dp.r_minor * dp.r_ratio);
303
304	//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
305	sum +=weights_rad[i].weight *weights_len[j].weight
306	* weights_rat[m].weight*DiamCyl(dp.length, suf_rad)/2.0;
307	norm += weights_rad[i].weight weights_len[j].weight weights_rat[m].weight;
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(dp.length,dp.radius)" is DIAMETER.
317	rad_out = DiamCyl(dp.length, suf_rad)/2.0;}
318
319	return rad_out;
320	}

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

SasView

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

Download in other formats: