triaxialellipsoid.cpp @ 8a6d4af

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 8a6d4af was f9bf661, checked in by Jae Cho <jhjcho@…>, 15 years ago
updated documents
Property mode set to `100644`
File size: 8.5 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
84	// Loop over semi axis A weight points
85	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
86	dp[1] = weights_semi_axisA[i].value;
87
88	// Loop over semi axis B weight points
89	for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
90	dp[2] = weights_semi_axisB[j].value;
91
92	// Loop over semi axis C weight points
93	for(int k=0; k< (int)weights_semi_axisC.size(); k++) {
94	dp[3] = weights_semi_axisC[k].value;
95
96	sum += weights_semi_axisA[i].weight
97	* weights_semi_axisB[j].weight * weights_semi_axisC[k].weight* TriaxialEllipsoid(dp, q);
98	norm += weights_semi_axisA[i].weight
99	* weights_semi_axisB[j].weight * weights_semi_axisC[k].weight;
100	}
101	}
102	}
103	return sum/norm + background();
104	}
105
106	/**
107	* Function to evaluate 2D scattering function
108	* @param q_x: value of Q along x
109	* @param q_y: value of Q along y
110	* @return: function value
111	*/
112	double TriaxialEllipsoidModel :: operator()(double qx, double qy) {
113	TriaxialEllipsoidParameters dp;
114	// Fill parameter array
115	dp.scale = scale();
116	dp.semi_axisA = semi_axisA();
117	dp.semi_axisB = semi_axisB();
118	dp.semi_axisC = semi_axisC();
119	dp.contrast = contrast();
120	dp.background = 0.0;
121	dp.axis_theta = axis_theta();
122	dp.axis_phi = axis_phi();
123	dp.axis_psi = axis_psi();
124
125	// Get the dispersion points for the semi_axis A
126	vector<WeightPoint> weights_semi_axisA;
127	semi_axisA.get_weights(weights_semi_axisA);
128
129	// Get the dispersion points for the semi_axis B
130	vector<WeightPoint> weights_semi_axisB;
131	semi_axisB.get_weights(weights_semi_axisB);
132
133	// Get the dispersion points for the semi_axis C
134	vector<WeightPoint> weights_semi_axisC;
135	semi_axisC.get_weights(weights_semi_axisC);
136
137	// Get angular averaging for theta
138	vector<WeightPoint> weights_theta;
139	axis_theta.get_weights(weights_theta);
140
141	// Get angular averaging for phi
142	vector<WeightPoint> weights_phi;
143	axis_phi.get_weights(weights_phi);
144
145	// Get angular averaging for psi
146	vector<WeightPoint> weights_psi;
147	axis_psi.get_weights(weights_psi);
148
149	// Perform the computation, with all weight points
150	double sum = 0.0;
151	double norm = 0.0;
152
153	// Loop over semi axis A weight points
154	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
155	dp.semi_axisA = weights_semi_axisA[i].value;
156
157	// Loop over semi axis B weight points
158	for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
159	dp.semi_axisB = weights_semi_axisB[j].value;
160
161	// Loop over semi axis C weight points
162	for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
163	dp.semi_axisC = weights_semi_axisC[k].value;
164
165	// Average over theta distribution
166	for(int l=0; l< (int)weights_theta.size(); l++) {
167	dp.axis_theta = weights_theta[l].value;
168
169	// Average over phi distribution
170	for(int m=0; m <(int)weights_phi.size(); m++) {
171	dp.axis_phi = weights_phi[m].value;
172	// Average over psi distribution
173	for(int n=0; n <(int)weights_psi.size(); n++) {
174	dp.axis_psi = weights_psi[n].value;
175
176	double _ptvalue = weights_semi_axisA[i].weight
177	* weights_semi_axisB[j].weight
178	* weights_semi_axisC[k].weight
179	* weights_theta[l].weight
180	* weights_phi[m].weight
181	* weights_psi[n].weight
182	* triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy);
183	if (weights_theta.size()>1) {
184	_ptvalue *= sin(weights_theta[k].value);
185	}
186	sum += _ptvalue;
187
188	norm += weights_semi_axisA[i].weight
189	* weights_semi_axisB[j].weight
190	* weights_semi_axisC[k].weight
191	* weights_theta[l].weight
192	* weights_phi[m].weight
193	* weights_psi[n].weight;
194	}
195	}
196
197	}
198	}
199	}
200	}
201	// Averaging in theta needs an extra normalization
202	// factor to account for the sin(theta) term in the
203	// integration (see documentation).
204	if (weights_theta.size()>1) norm = norm / asin(1.0);
205	return sum/norm + background();
206	}
207
208	/**
209	* Function to evaluate 2D scattering function
210	* @param pars: parameters of the triaxial ellipsoid
211	* @param q: q-value
212	* @param phi: angle phi
213	* @return: function value
214	*/
215	double TriaxialEllipsoidModel :: evaluate_rphi(double q, double phi) {
216	double qx = q*cos(phi);
217	double qy = q*sin(phi);
218	return (*this).operator()(qx, qy);
219	}
220	/**
221	* Function to calculate effective radius
222	* @return: effective radius value
223	*/
224	double TriaxialEllipsoidModel :: calculate_ER() {
225	TriaxialEllipsoidParameters dp;
226
227	dp.semi_axisA = semi_axisA();
228	dp.semi_axisB = semi_axisB();
229	//polar axis C
230	dp.semi_axisC = semi_axisC();
231
232	double rad_out = 0.0;
233	//Surface average radius at the equat. cross section.
234	double suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
235
236	// Perform the computation, with all weight points
237	double sum = 0.0;
238	double norm = 0.0;
239
240	// Get the dispersion points for the semi_axis A
241	vector<WeightPoint> weights_semi_axisA;
242	semi_axisA.get_weights(weights_semi_axisA);
243
244	// Get the dispersion points for the semi_axis B
245	vector<WeightPoint> weights_semi_axisB;
246	semi_axisB.get_weights(weights_semi_axisB);
247
248	// Get the dispersion points for the semi_axis C
249	vector<WeightPoint> weights_semi_axisC;
250	semi_axisC.get_weights(weights_semi_axisC);
251
252	// Loop over semi axis A weight points
253	for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
254	dp.semi_axisA = weights_semi_axisA[i].value;
255
256	// Loop over semi axis B weight points
257	for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
258	dp.semi_axisB = weights_semi_axisB[j].value;
259
260	// Loop over semi axis C weight points
261	for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
262	dp.semi_axisC = weights_semi_axisC[k].value;
263
264	//Calculate surface averaged radius
265	suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
266
267	//Sum
268	sum += weights_semi_axisA[i].weight
269	* weights_semi_axisB[j].weight
270	* weights_semi_axisC[k].weight * DiamEllip(dp.semi_axisC, suf_rad)/2.0;
271	//Norm
272	norm += weights_semi_axisA[i].weight* weights_semi_axisB[j].weight
273	* weights_semi_axisC[k].weight;
274	}
275	}
276	}
277	if (norm != 0){
278	//return the averaged value
279	rad_out = sum/norm;}
280	else{
281	//return normal value
282	rad_out = DiamEllip(dp.semi_axisC, suf_rad)/2.0;}
283
284	return rad_out;
285	}

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

SasView

source: sasview/sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp @ 8a6d4af

Download in other formats: