csparallelepiped.cpp @ acd0fd10

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 acd0fd10 was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago
Reorganizing models in preparation of cpp cleanup
Property mode set to `100644`
File size: 10.2 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 2010, 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
21	#include <math.h>
22	#include "models.hh"
23	#include "parameters.hh"
24	#include <stdio.h>
25	using namespace std;
26
27	extern "C" {
28	#include "libCylinder.h"
29	#include "libStructureFactor.h"
30	#include "csparallelepiped.h"
31	}
32
33	CSParallelepipedModel :: CSParallelepipedModel() {
34	scale = Parameter(1.0);
35	shortA = Parameter(35.0, true);
36	shortA.set_min(1.0);
37	midB = Parameter(75.0, true);
38	midB.set_min(1.0);
39	longC = Parameter(400.0, true);
40	longC.set_min(1.0);
41	rimA = Parameter(10.0, true);
42	rimB = Parameter(10.0, true);
43	rimC = Parameter(10.0, true);
44	sld_rimA = Parameter(2.0e-6, true);
45	sld_rimB = Parameter(4.0e-6, true);
46	sld_rimC = Parameter(2.0e-6, true);
47	sld_pcore = Parameter(1.0e-6);
48	sld_solv = Parameter(6.0e-6);
49	background = Parameter(0.06);
50	parallel_theta = Parameter(0.0, true);
51	parallel_phi = Parameter(0.0, true);
52	parallel_psi = Parameter(0.0, true);
53	}
54
55	/**
56	* Function to evaluate 1D scattering function
57	* The NIST IGOR library is used for the actual calculation.
58	* @param q: q-value
59	* @return: function value
60	*/
61	double CSParallelepipedModel :: operator()(double q) {
62	double dp[13];
63
64	// Fill parameter array for IGOR library
65	// Add the background after averaging
66	dp[0] = scale();
67	dp[1] = shortA();
68	dp[2] = midB();
69	dp[3] = longC();
70	dp[4] = rimA();
71	dp[5] = rimB();
72	dp[6] = rimC();
73	dp[7] = sld_rimA();
74	dp[8] = sld_rimB();
75	dp[9] = sld_rimC();
76	dp[10] = sld_pcore();
77	dp[11] = sld_solv();
78	dp[12] = 0.0;
79
80	// Get the dispersion points for the short_edgeA
81	vector<WeightPoint> weights_shortA;
82	shortA.get_weights(weights_shortA);
83
84	// Get the dispersion points for the longer_edgeB
85	vector<WeightPoint> weights_midB;
86	midB.get_weights(weights_midB);
87
88	// Get the dispersion points for the longuest_edgeC
89	vector<WeightPoint> weights_longC;
90	longC.get_weights(weights_longC);
91
92
93
94	// Perform the computation, with all weight points
95	double sum = 0.0;
96	double norm = 0.0;
97	double vol = 0.0;
98
99	// Loop over short_edgeA weight points
100	for(int i=0; i< (int)weights_shortA.size(); i++) {
101	dp[1] = weights_shortA[i].value;
102
103	// Loop over longer_edgeB weight points
104	for(int j=0; j< (int)weights_midB.size(); j++) {
105	dp[2] = weights_midB[j].value;
106
107	// Loop over longuest_edgeC weight points
108	for(int k=0; k< (int)weights_longC.size(); k++) {
109	dp[3] = weights_longC[k].value;
110	//Un-normalize by volume
111	sum += weights_shortA[i].weight * weights_midB[j].weight
112	* weights_longC[k].weight * CSParallelepiped(dp, q)
113	* weights_shortA[i].value*weights_midB[j].value
114	* weights_longC[k].value;
115	//Find average volume
116	vol += weights_shortA[i].weight * weights_midB[j].weight
117	* weights_longC[k].weight
118	* weights_shortA[i].value * weights_midB[j].value
119	* weights_longC[k].value;
120
121	norm += weights_shortA[i].weight
122	* weights_midB[j].weight * weights_longC[k].weight;
123	}
124	}
125	}
126	if (vol != 0.0 && norm != 0.0) {
127	//Re-normalize by avg volume
128	sum = sum/(vol/norm);}
129
130	return sum/norm + background();
131	}
132	/**
133	* Function to evaluate 2D scattering function
134	* @param q_x: value of Q along x
135	* @param q_y: value of Q along y
136	* @return: function value
137	*/
138	double CSParallelepipedModel :: operator()(double qx, double qy) {
139	CSParallelepipedParameters dp;
140	// Fill parameter array
141	dp.scale = scale();
142	dp.shortA = shortA();
143	dp.midB = midB();
144	dp.longC = longC();
145	dp.rimA = rimA();
146	dp.rimB = rimB();
147	dp.rimC = rimC();
148	dp.sld_rimA = sld_rimA();
149	dp.sld_rimB = sld_rimB();
150	dp.sld_rimC = sld_rimC();
151	dp.sld_pcore = sld_pcore();
152	dp.sld_solv = sld_solv();
153	dp.background = 0.0;
154	//dp.background = background();
155	dp.parallel_theta = parallel_theta();
156	dp.parallel_phi = parallel_phi();
157	dp.parallel_psi = parallel_psi();
158
159
160
161	// Get the dispersion points for the short_edgeA
162	vector<WeightPoint> weights_shortA;
163	shortA.get_weights(weights_shortA);
164
165	// Get the dispersion points for the longer_edgeB
166	vector<WeightPoint> weights_midB;
167	midB.get_weights(weights_midB);
168
169	// Get the dispersion points for the longuest_edgeC
170	vector<WeightPoint> weights_longC;
171	longC.get_weights(weights_longC);
172
173	// Get angular averaging for theta
174	vector<WeightPoint> weights_parallel_theta;
175	parallel_theta.get_weights(weights_parallel_theta);
176
177	// Get angular averaging for phi
178	vector<WeightPoint> weights_parallel_phi;
179	parallel_phi.get_weights(weights_parallel_phi);
180
181	// Get angular averaging for psi
182	vector<WeightPoint> weights_parallel_psi;
183	parallel_psi.get_weights(weights_parallel_psi);
184
185	// Perform the computation, with all weight points
186	double sum = 0.0;
187	double norm = 0.0;
188	double norm_vol = 0.0;
189	double vol = 0.0;
190	double pi = 4.0*atan(1.0);
191
192	// Loop over radius weight points
193	for(int i=0; i< (int)weights_shortA.size(); i++) {
194	dp.shortA = weights_shortA[i].value;
195
196	// Loop over longer_edgeB weight points
197	for(int j=0; j< (int)weights_midB.size(); j++) {
198	dp.midB = weights_midB[j].value;
199
200	// Average over longuest_edgeC distribution
201	for(int k=0; k< (int)weights_longC.size(); k++) {
202	dp.longC = weights_longC[k].value;
203
204	// Average over theta distribution
205	for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
206	dp.parallel_theta = weights_parallel_theta[l].value;
207
208	// Average over phi distribution
209	for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
210	dp.parallel_phi = weights_parallel_phi[m].value;
211
212	// Average over phi distribution
213	for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
214	dp.parallel_psi = weights_parallel_psi[n].value;
215	//Un-normalize by volume
216	double _ptvalue = weights_shortA[i].weight
217	* weights_midB[j].weight
218	* weights_longC[k].weight
219	* weights_parallel_theta[l].weight
220	* weights_parallel_phi[m].weight
221	* weights_parallel_psi[n].weight
222	* csparallelepiped_analytical_2DXY(&dp, qx, qy)
223	* weights_shortA[i].value*weights_midB[j].value
224	* weights_longC[k].value;
225
226	if (weights_parallel_theta.size()>1) {
227	_ptvalue = fabs(sin(weights_parallel_theta[l].valuepi/180.0));
228	}
229	sum += _ptvalue;
230	//Find average volume
231	vol += weights_shortA[i].weight
232	* weights_midB[j].weight
233	* weights_longC[k].weight
234	* weights_shortA[i].value*weights_midB[j].value
235	* weights_longC[k].value;
236	//Find norm for volume
237	norm_vol += weights_shortA[i].weight
238	* weights_midB[j].weight
239	* weights_longC[k].weight;
240
241	norm += weights_shortA[i].weight
242	* weights_midB[j].weight
243	* weights_longC[k].weight
244	* weights_parallel_theta[l].weight
245	* weights_parallel_phi[m].weight
246	* weights_parallel_psi[n].weight;
247	}
248	}
249
250	}
251	}
252	}
253	}
254	// Averaging in theta needs an extra normalization
255	// factor to account for the sin(theta) term in the
256	// integration (see documentation).
257	if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
258
259	if (vol != 0.0 && norm_vol != 0.0) {
260	//Re-normalize by avg volume
261	sum = sum/(vol/norm_vol);}
262
263	return sum/norm + background();
264	}
265
266
267	/**
268	* Function to evaluate 2D scattering function
269	* @param pars: parameters of the cylinder
270	* @param q: q-value
271	* @param phi: angle phi
272	* @return: function value
273	*/
274	double CSParallelepipedModel :: evaluate_rphi(double q, double phi) {
275	double qx = q*cos(phi);
276	double qy = q*sin(phi);
277	return (*this).operator()(qx, qy);
278	}
279	/**
280	* Function to calculate effective radius
281	* @return: effective radius value
282	*/
283	double CSParallelepipedModel :: calculate_ER() {
284	CSParallelepipedParameters dp;
285	dp.shortA = shortA();
286	dp.midB = midB();
287	dp.longC = longC();
288	dp.rimA = rimA();
289	dp.rimB = rimB();
290	dp.rimC = rimC();
291
292	double rad_out = 0.0;
293	double pi = 4.0*atan(1.0);
294	double suf_rad = sqrt((dp.shortAdp.midB+2.0dp.rimAdp.midB+2.0dp.rimA*dp.shortA)/pi);
295	double height =(dp.longC + 2.0*dp.rimC);
296	// Perform the computation, with all weight points
297	double sum = 0.0;
298	double norm = 0.0;
299
300	// Get the dispersion points for the short_edgeA
301	vector<WeightPoint> weights_shortA;
302	shortA.get_weights(weights_shortA);
303
304	// Get the dispersion points for the longer_edgeB
305	vector<WeightPoint> weights_midB;
306	midB.get_weights(weights_midB);
307
308	// Get angular averaging for the longuest_edgeC
309	vector<WeightPoint> weights_longC;
310	longC.get_weights(weights_longC);
311
312	// Loop over radius weight points
313	for(int i=0; i< (int)weights_shortA.size(); i++) {
314	dp.shortA = weights_shortA[i].value;
315
316	// Loop over longer_edgeB weight points
317	for(int j=0; j< (int)weights_midB.size(); j++) {
318	dp.midB = weights_midB[j].value;
319
320	// Average over longuest_edgeC distribution
321	for(int k=0; k< (int)weights_longC.size(); k++) {
322	dp.longC = weights_longC[k].value;
323	//Calculate surface averaged radius
324	//This is rough approximation.
325	suf_rad = sqrt((dp.shortAdp.midB+2.0dp.rimAdp.midB+2.0dp.rimA*dp.shortA)/pi);
326	height =(dp.longC + 2.0*dp.rimC);
327	//Note: output of "DiamCyl(dp.length,dp.radius)" is a DIAMETER.
328	sum +=weights_shortA[i].weight* weights_midB[j].weight
329	* weights_longC[k].weight*DiamCyl(height, suf_rad)/2.0;
330	norm += weights_shortA[i].weight* weights_midB[j].weight*weights_longC[k].weight;
331	}
332	}
333	}
334
335	if (norm != 0){
336	//return the averaged value
337	rad_out = sum/norm;}
338	else{
339	//return normal value
340	//Note: output of "DiamCyl(length,radius)" is DIAMETER.
341	rad_out = DiamCyl(dp.longC, suf_rad)/2.0;}
342	return rad_out;
343
344	}

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

SasView

source: sasview/sansmodels/src/c_models/csparallelepiped.cpp @ acd0fd10

Download in other formats: