triaxialellipsoid.cpp @ 2db1d66

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 2db1d66 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	*/
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 "triaxial_ellipsoid.h"
33	}
34
35	TriaxialEllipsoidModel :: TriaxialEllipsoidModel() {
36	scale = Parameter(1.0);
37	semi_axisA = Parameter(35.0, true);
38	semi_axisA.set_min(0.0);
39	semi_axisB = Parameter(100.0, true);
40	semi_axisB.set_min(0.0);
41	semi_axisC = Parameter(400.0, true);
42	semi_axisC.set_min(0.0);
43	contrast = Parameter(5.3e-6);
44	background = Parameter(0.0);
45	axis_theta = Parameter(1.0, true);
46	axis_phi = Parameter(1.0, true);
47	axis_psi = Parameter(0.0, true);
48	}
49
50	/**
51	* Function to evaluate 1D scattering function
52	* The NIST IGOR library is used for the actual calculation.
53	* @param q: q-value
54	* @return: function value
55	*/
56	double TriaxialEllipsoidModel :: operator()(double q) {
57	double dp[6];
58
59	// Fill parameter array for IGOR library
60	// Add the background after averaging
61	dp[0] = scale();
62	dp[1] = semi_axisA();
63	dp[2] = semi_axisB();
64	dp[3] = semi_axisC();
65	dp[4] = contrast();
66	dp[5] = 0.0;
67
68	// Get the dispersion points for the semi axis A
69	vector<WeightPoint> weights_semi_axisA;
70	semi_axisA.get_weights(weights_semi_axisA);
71
72	// Get the dispersion points for the semi axis B
73	vector<WeightPoint> weights_semi_axisB;
74	semi_axisB.get_weights(weights_semi_axisB);
75
76	// Get the dispersion points for the semi axis C
77	vector<WeightPoint> weights_semi_axisC;
78	semi_axisC.get_weights(weights_semi_axisC);
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 semi axis A weight points
86	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
87	dp[1] = weights_semi_axisA[i].value;
88
89	// Loop over semi axis B weight points
90	for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
91	dp[2] = weights_semi_axisB[j].value;
92
93	// Loop over semi axis C weight points
94	for(int k=0; k< (int)weights_semi_axisC.size(); k++) {
95	dp[3] = weights_semi_axisC[k].value;
96	//Un-normalize by volume
97	sum += weights_semi_axisA[i].weight
98	* weights_semi_axisB[j].weight * weights_semi_axisC[k].weight* TriaxialEllipsoid(dp, q)
99	* weights_semi_axisA[i].valueweights_semi_axisB[j].valueweights_semi_axisC[k].value;
100	//Find average volume
101	vol += weights_semi_axisA[i].weight
102	* weights_semi_axisB[j].weight * weights_semi_axisC[k].weight
103	* weights_semi_axisA[i].valueweights_semi_axisB[j].valueweights_semi_axisC[k].value;
104
105	norm += weights_semi_axisA[i].weight
106	* weights_semi_axisB[j].weight * weights_semi_axisC[k].weight;
107	}
108	}
109	}
110	if (vol != 0.0 && norm != 0.0) {
111	//Re-normalize by avg volume
112	sum = sum/(vol/norm);}
113
114	return sum/norm + background();
115	}
116
117	/**
118	* Function to evaluate 2D scattering function
119	* @param q_x: value of Q along x
120	* @param q_y: value of Q along y
121	* @return: function value
122	*/
123	double TriaxialEllipsoidModel :: operator()(double qx, double qy) {
124	TriaxialEllipsoidParameters dp;
125	// Fill parameter array
126	dp.scale = scale();
127	dp.semi_axisA = semi_axisA();
128	dp.semi_axisB = semi_axisB();
129	dp.semi_axisC = semi_axisC();
130	dp.contrast = contrast();
131	dp.background = 0.0;
132	dp.axis_theta = axis_theta();
133	dp.axis_phi = axis_phi();
134	dp.axis_psi = axis_psi();
135
136	// Get the dispersion points for the semi_axis A
137	vector<WeightPoint> weights_semi_axisA;
138	semi_axisA.get_weights(weights_semi_axisA);
139
140	// Get the dispersion points for the semi_axis B
141	vector<WeightPoint> weights_semi_axisB;
142	semi_axisB.get_weights(weights_semi_axisB);
143
144	// Get the dispersion points for the semi_axis C
145	vector<WeightPoint> weights_semi_axisC;
146	semi_axisC.get_weights(weights_semi_axisC);
147
148	// Get angular averaging for theta
149	vector<WeightPoint> weights_theta;
150	axis_theta.get_weights(weights_theta);
151
152	// Get angular averaging for phi
153	vector<WeightPoint> weights_phi;
154	axis_phi.get_weights(weights_phi);
155
156	// Get angular averaging for psi
157	vector<WeightPoint> weights_psi;
158	axis_psi.get_weights(weights_psi);
159
160	// Perform the computation, with all weight points
161	double sum = 0.0;
162	double norm = 0.0;
163	double norm_vol = 0.0;
164	double vol = 0.0;
165
166	// Loop over semi axis A weight points
167	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
168	dp.semi_axisA = weights_semi_axisA[i].value;
169
170	// Loop over semi axis B weight points
171	for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
172	dp.semi_axisB = weights_semi_axisB[j].value;
173
174	// Loop over semi axis C weight points
175	for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
176	dp.semi_axisC = weights_semi_axisC[k].value;
177
178	// Average over theta distribution
179	for(int l=0; l< (int)weights_theta.size(); l++) {
180	dp.axis_theta = weights_theta[l].value;
181
182	// Average over phi distribution
183	for(int m=0; m <(int)weights_phi.size(); m++) {
184	dp.axis_phi = weights_phi[m].value;
185	// Average over psi distribution
186	for(int n=0; n <(int)weights_psi.size(); n++) {
187	dp.axis_psi = weights_psi[n].value;
188	//Un-normalize by volume
189	double _ptvalue = weights_semi_axisA[i].weight
190	* weights_semi_axisB[j].weight
191	* weights_semi_axisC[k].weight
192	* weights_theta[l].weight
193	* weights_phi[m].weight
194	* weights_psi[n].weight
195	* triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy)
196	* weights_semi_axisA[i].valueweights_semi_axisB[j].valueweights_semi_axisC[k].value;
197	if (weights_theta.size()>1) {
198	_ptvalue *= sin(weights_theta[k].value);
199	}
200	sum += _ptvalue;
201	//Find average volume
202	vol += weights_semi_axisA[i].weight
203	* weights_semi_axisB[j].weight
204	* weights_semi_axisC[k].weight
205	* weights_semi_axisA[i].valueweights_semi_axisB[j].valueweights_semi_axisC[k].value;
206	//Find norm for volume
207	norm_vol += weights_semi_axisA[i].weight
208	* weights_semi_axisB[j].weight
209	* weights_semi_axisC[k].weight;
210
211	norm += weights_semi_axisA[i].weight
212	* weights_semi_axisB[j].weight
213	* weights_semi_axisC[k].weight
214	* weights_theta[l].weight
215	* weights_phi[m].weight
216	* weights_psi[n].weight;
217	}
218	}
219
220	}
221	}
222	}
223	}
224	// Averaging in theta needs an extra normalization
225	// factor to account for the sin(theta) term in the
226	// integration (see documentation).
227	if (weights_theta.size()>1) norm = norm / asin(1.0);
228
229	if (vol != 0.0 && norm_vol != 0.0) {
230	//Re-normalize by avg volume
231	sum = sum/(vol/norm_vol);}
232
233	return sum/norm + background();
234	}
235
236	/**
237	* Function to evaluate 2D scattering function
238	* @param pars: parameters of the triaxial ellipsoid
239	* @param q: q-value
240	* @param phi: angle phi
241	* @return: function value
242	*/
243	double TriaxialEllipsoidModel :: evaluate_rphi(double q, double phi) {
244	double qx = q*cos(phi);
245	double qy = q*sin(phi);
246	return (*this).operator()(qx, qy);
247	}
248	/**
249	* Function to calculate effective radius
250	* @return: effective radius value
251	*/
252	double TriaxialEllipsoidModel :: calculate_ER() {
253	TriaxialEllipsoidParameters dp;
254
255	dp.semi_axisA = semi_axisA();
256	dp.semi_axisB = semi_axisB();
257	//polar axis C
258	dp.semi_axisC = semi_axisC();
259
260	double rad_out = 0.0;
261	//Surface average radius at the equat. cross section.
262	double suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
263
264	// Perform the computation, with all weight points
265	double sum = 0.0;
266	double norm = 0.0;
267
268	// Get the dispersion points for the semi_axis A
269	vector<WeightPoint> weights_semi_axisA;
270	semi_axisA.get_weights(weights_semi_axisA);
271
272	// Get the dispersion points for the semi_axis B
273	vector<WeightPoint> weights_semi_axisB;
274	semi_axisB.get_weights(weights_semi_axisB);
275
276	// Get the dispersion points for the semi_axis C
277	vector<WeightPoint> weights_semi_axisC;
278	semi_axisC.get_weights(weights_semi_axisC);
279
280	// Loop over semi axis A weight points
281	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
282	dp.semi_axisA = weights_semi_axisA[i].value;
283
284	// Loop over semi axis B weight points
285	for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
286	dp.semi_axisB = weights_semi_axisB[j].value;
287
288	// Loop over semi axis C weight points
289	for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
290	dp.semi_axisC = weights_semi_axisC[k].value;
291
292	//Calculate surface averaged radius
293	suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
294
295	//Sum
296	sum += weights_semi_axisA[i].weight
297	* weights_semi_axisB[j].weight
298	* weights_semi_axisC[k].weight * DiamEllip(dp.semi_axisC, suf_rad)/2.0;
299	//Norm
300	norm += weights_semi_axisA[i].weight* weights_semi_axisB[j].weight
301	* weights_semi_axisC[k].weight;
302	}
303	}
304	}
305	if (norm != 0){
306	//return the averaged value
307	rad_out = sum/norm;}
308	else{
309	//return normal value
310	rad_out = DiamEllip(dp.semi_axisC, suf_rad)/2.0;}
311
312	return rad_out;
313	}

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

SasView

source: sasview/sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp @ 2db1d66

Download in other formats: