triaxialellipsoid.cpp @ 4cc8285c

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

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

SasView

source: sasview/sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp @ 4cc8285c

Download in other formats: