parallelepiped.cpp @ 530cc84

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 530cc84 was 82c11d3, checked in by Mathieu Doucet <doucetm@…>, 13 years ago
refactored bunch of models
Property mode set to `100644`
File size: 13.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: add 2D function
21	*/
22
23	#include <math.h>
24	#include "parameters.hh"
25	#include <stdio.h>
26	using namespace std;
27
28	extern "C" {
29	#include "libCylinder.h"
30	#include "libStructureFactor.h"
31	}
32	#include "parallelepiped.h"
33
34	// Convenience parameter structure
35	typedef struct {
36	double scale;
37	double short_a;
38	double short_b;
39	double long_c;
40	double sldPipe;
41	double sldSolv;
42	double background;
43	double parallel_theta;
44	double parallel_phi;
45	double parallel_psi;
46	} ParallelepipedParameters;
47
48
49	static double pkernel(double a, double b,double c, double ala, double alb, double alc){
50	// mu passed in is really mu*sqrt(1-sig^2)
51	double argA,argB,argC,tmp1,tmp2,tmp3; //local variables
52
53	//handle arg=0 separately, as sin(t)/t -> 1 as t->0
54	argA = a*ala/2.0;
55	argB = b*alb/2.0;
56	argC = c*alc/2.0;
57	if(argA==0.0) {
58	tmp1 = 1.0;
59	} else {
60	tmp1 = sin(argA)*sin(argA)/argA/argA;
61	}
62	if (argB==0.0) {
63	tmp2 = 1.0;
64	} else {
65	tmp2 = sin(argB)*sin(argB)/argB/argB;
66	}
67
68	if (argC==0.0) {
69	tmp3 = 1.0;
70	} else {
71	tmp3 = sin(argC)*sin(argC)/argC/argC;
72	}
73
74	return (tmp1tmp2tmp3);
75
76	}
77
78	/**
79	* Function to evaluate 2D scattering function
80	* @param pars: parameters of the parallelepiped
81	* @param q: q-value
82	* @param q_x: q_x / q
83	* @param q_y: q_y / q
84	* @return: function value
85	*/
86	static double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) {
87	double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y;
88	double q_z;
89	double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
90
91	double answer;
92	double pi = 4.0*atan(1.0);
93
94	//convert angle degree to radian
95	double theta = pars->parallel_theta * pi/180.0;
96	double phi = pars->parallel_phi * pi/180.0;
97	double psi = pars->parallel_psi * pi/180.0;
98
99	edgeA = pars->short_a;
100	edgeB = pars->short_b;
101	edgeC = pars->long_c;
102
103
104	// parallelepiped c axis orientation
105	cparallel_x = sin(theta) * cos(phi);
106	cparallel_y = sin(theta) * sin(phi);
107	cparallel_z = cos(theta);
108
109	// q vector
110	q_z = 0.0;
111
112	// Compute the angle btw vector q and the
113	// axis of the parallelepiped
114	cos_val_c = cparallel_xq_x + cparallel_yq_y + cparallel_z*q_z;
115	alpha = acos(cos_val_c);
116
117	// parallelepiped a axis orientation
118	parallel_x = sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi);
119	parallel_y = cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi);
120
121	cos_val_a = parallel_xq_x + parallel_yq_y;
122
123
124
125	// parallelepiped b axis orientation
126	bparallel_x = sqrt(1.0-sin(theta)cos(phi))cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi);
127	bparallel_y = sqrt(1.0-sin(theta)cos(phi))sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi);
128	// axis of the parallelepiped
129	cos_val_b = sin(acos(cos_val_a)) ;
130
131
132
133	// The following test should always pass
134	if (fabs(cos_val_c)>1.0) {
135	printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
136	return 0;
137	}
138
139	// Call the IGOR library function to get the kernel
140	answer = pkernel( qedgeA, qedgeB, qedgeC, sin(alpha)cos_val_a,sin(alpha)*cos_val_b,cos_val_c);
141
142	// Multiply by contrast^2
143	answer = (pars->sldPipe - pars->sldSolv) (pars->sldPipe - pars->sldSolv);
144
145	//normalize by cylinder volume
146	//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel
147	vol = edgeA* edgeB * edgeC;
148	answer *= vol;
149
150	//convert to [cm-1]
151	answer *= 1.0e8;
152
153	//Scale
154	answer *= pars->scale;
155
156	// add in the background
157	answer += pars->background;
158
159	return answer;
160	}
161
162	/**
163	* Function to evaluate 2D scattering function
164	* @param pars: parameters of the parallelepiped
165	* @param q: q-value
166	* @return: function value
167	*/
168	static double parallelepiped_analytical_2DXY(ParallelepipedParameters *pars, double qx, double qy) {
169	double q;
170	q = sqrt(qxqx+qyqy);
171	return parallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q);
172	}
173
174	ParallelepipedModel :: ParallelepipedModel() {
175	scale = Parameter(1.0);
176	short_a = Parameter(35.0, true);
177	short_a.set_min(1.0);
178	short_b = Parameter(75.0, true);
179	short_b.set_min(1.0);
180	long_c = Parameter(400.0, true);
181	long_c.set_min(1.0);
182	sldPipe = Parameter(6.3e-6);
183	sldSolv = Parameter(1.0e-6);
184	background = Parameter(0.0);
185	parallel_theta = Parameter(0.0, true);
186	parallel_phi = Parameter(0.0, true);
187	parallel_psi = Parameter(0.0, true);
188	}
189
190	/**
191	* Function to evaluate 1D scattering function
192	* The NIST IGOR library is used for the actual calculation.
193	* @param q: q-value
194	* @return: function value
195	*/
196	double ParallelepipedModel :: operator()(double q) {
197	double dp[7];
198
199	// Fill parameter array for IGOR library
200	// Add the background after averaging
201	dp[0] = scale();
202	dp[1] = short_a();
203	dp[2] = short_b();
204	dp[3] = long_c();
205	dp[4] = sldPipe();
206	dp[5] = sldSolv();
207	dp[6] = 0.0;
208
209	// Get the dispersion points for the short_edgeA
210	vector<WeightPoint> weights_short_a;
211	short_a.get_weights(weights_short_a);
212
213	// Get the dispersion points for the longer_edgeB
214	vector<WeightPoint> weights_short_b;
215	short_b.get_weights(weights_short_b);
216
217	// Get the dispersion points for the longuest_edgeC
218	vector<WeightPoint> weights_long_c;
219	long_c.get_weights(weights_long_c);
220
221
222
223	// Perform the computation, with all weight points
224	double sum = 0.0;
225	double norm = 0.0;
226	double vol = 0.0;
227
228	// Loop over short_edgeA weight points
229	for(int i=0; i< (int)weights_short_a.size(); i++) {
230	dp[1] = weights_short_a[i].value;
231
232	// Loop over longer_edgeB weight points
233	for(int j=0; j< (int)weights_short_b.size(); j++) {
234	dp[2] = weights_short_b[j].value;
235
236	// Loop over longuest_edgeC weight points
237	for(int k=0; k< (int)weights_long_c.size(); k++) {
238	dp[3] = weights_long_c[k].value;
239	//Un-normalize by volume
240	sum += weights_short_a[i].weight * weights_short_b[j].weight
241	* weights_long_c[k].weight * Parallelepiped(dp, q)
242	* weights_short_a[i].value*weights_short_b[j].value
243	* weights_long_c[k].value;
244	//Find average volume
245	vol += weights_short_a[i].weight * weights_short_b[j].weight
246	* weights_long_c[k].weight
247	* weights_short_a[i].value * weights_short_b[j].value
248	* weights_long_c[k].value;
249
250	norm += weights_short_a[i].weight
251	* weights_short_b[j].weight * weights_long_c[k].weight;
252	}
253	}
254	}
255	if (vol != 0.0 && norm != 0.0) {
256	//Re-normalize by avg volume
257	sum = sum/(vol/norm);}
258
259	return sum/norm + background();
260	}
261	/**
262	* Function to evaluate 2D scattering function
263	* @param q_x: value of Q along x
264	* @param q_y: value of Q along y
265	* @return: function value
266	*/
267	double ParallelepipedModel :: operator()(double qx, double qy) {
268	ParallelepipedParameters dp;
269	// Fill parameter array
270	dp.scale = scale();
271	dp.short_a = short_a();
272	dp.short_b = short_b();
273	dp.long_c = long_c();
274	dp.sldPipe = sldPipe();
275	dp.sldSolv = sldSolv();
276	dp.background = 0.0;
277	//dp.background = background();
278	dp.parallel_theta = parallel_theta();
279	dp.parallel_phi = parallel_phi();
280	dp.parallel_psi = parallel_psi();
281
282
283	// Get the dispersion points for the short_edgeA
284	vector<WeightPoint> weights_short_a;
285	short_a.get_weights(weights_short_a);
286
287	// Get the dispersion points for the longer_edgeB
288	vector<WeightPoint> weights_short_b;
289	short_b.get_weights(weights_short_b);
290
291	// Get angular averaging for the longuest_edgeC
292	vector<WeightPoint> weights_long_c;
293	long_c.get_weights(weights_long_c);
294
295	// Get angular averaging for theta
296	vector<WeightPoint> weights_parallel_theta;
297	parallel_theta.get_weights(weights_parallel_theta);
298
299	// Get angular averaging for phi
300	vector<WeightPoint> weights_parallel_phi;
301	parallel_phi.get_weights(weights_parallel_phi);
302
303	// Get angular averaging for psi
304	vector<WeightPoint> weights_parallel_psi;
305	parallel_psi.get_weights(weights_parallel_psi);
306
307	// Perform the computation, with all weight points
308	double sum = 0.0;
309	double norm = 0.0;
310	double norm_vol = 0.0;
311	double vol = 0.0;
312	double pi = 4.0*atan(1.0);
313	// Loop over radius weight points
314	for(int i=0; i< (int)weights_short_a.size(); i++) {
315	dp.short_a = weights_short_a[i].value;
316
317	// Loop over longer_edgeB weight points
318	for(int j=0; j< (int)weights_short_b.size(); j++) {
319	dp.short_b = weights_short_b[j].value;
320
321	// Average over longuest_edgeC distribution
322	for(int k=0; k< (int)weights_long_c.size(); k++) {
323	dp.long_c = weights_long_c[k].value;
324
325	// Average over theta distribution
326	for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
327	dp.parallel_theta = weights_parallel_theta[l].value;
328
329	// Average over phi distribution
330	for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
331	dp.parallel_phi = weights_parallel_phi[m].value;
332
333	// Average over phi distribution
334	for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
335	dp.parallel_psi = weights_parallel_psi[n].value;
336	//Un-normalize by volume
337	double _ptvalue = weights_short_a[i].weight
338	* weights_short_b[j].weight
339	* weights_long_c[k].weight
340	* weights_parallel_theta[l].weight
341	* weights_parallel_phi[m].weight
342	* weights_parallel_psi[n].weight
343	* parallelepiped_analytical_2DXY(&dp, qx, qy)
344	* weights_short_a[i].value*weights_short_b[j].value
345	* weights_long_c[k].value;
346
347	if (weights_parallel_theta.size()>1) {
348	_ptvalue = fabs(sin(weights_parallel_theta[l].valuepi/180.0));
349	}
350	sum += _ptvalue;
351	//Find average volume
352	vol += weights_short_a[i].weight
353	* weights_short_b[j].weight
354	* weights_long_c[k].weight
355	* weights_short_a[i].value*weights_short_b[j].value
356	* weights_long_c[k].value;
357	//Find norm for volume
358	norm_vol += weights_short_a[i].weight
359	* weights_short_b[j].weight
360	* weights_long_c[k].weight;
361
362	norm += weights_short_a[i].weight
363	* weights_short_b[j].weight
364	* weights_long_c[k].weight
365	* weights_parallel_theta[l].weight
366	* weights_parallel_phi[m].weight
367	* weights_parallel_psi[n].weight;
368	}
369	}
370
371	}
372	}
373	}
374	}
375	// Averaging in theta needs an extra normalization
376	// factor to account for the sin(theta) term in the
377	// integration (see documentation).
378	if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
379
380	if (vol != 0.0 && norm_vol != 0.0) {
381	//Re-normalize by avg volume
382	sum = sum/(vol/norm_vol);}
383
384	return sum/norm + background();
385	}
386
387
388	/**
389	* Function to evaluate 2D scattering function
390	* @param pars: parameters of the cylinder
391	* @param q: q-value
392	* @param phi: angle phi
393	* @return: function value
394	*/
395	double ParallelepipedModel :: evaluate_rphi(double q, double phi) {
396	double qx = q*cos(phi);
397	double qy = q*sin(phi);
398	return (*this).operator()(qx, qy);
399	}
400	/**
401	* Function to calculate effective radius
402	* @return: effective radius value
403	*/
404	double ParallelepipedModel :: calculate_ER() {
405	ParallelepipedParameters dp;
406	dp.short_a = short_a();
407	dp.short_b = short_b();
408	dp.long_c = long_c();
409	double rad_out = 0.0;
410	double pi = 4.0*atan(1.0);
411	double suf_rad = sqrt(dp.short_a*dp.short_b/pi);
412
413	// Perform the computation, with all weight points
414	double sum = 0.0;
415	double norm = 0.0;
416
417	// Get the dispersion points for the short_edgeA
418	vector<WeightPoint> weights_short_a;
419	short_a.get_weights(weights_short_a);
420
421	// Get the dispersion points for the longer_edgeB
422	vector<WeightPoint> weights_short_b;
423	short_b.get_weights(weights_short_b);
424
425	// Get angular averaging for the longuest_edgeC
426	vector<WeightPoint> weights_long_c;
427	long_c.get_weights(weights_long_c);
428
429	// Loop over radius weight points
430	for(int i=0; i< (int)weights_short_a.size(); i++) {
431	dp.short_a = weights_short_a[i].value;
432
433	// Loop over longer_edgeB weight points
434	for(int j=0; j< (int)weights_short_b.size(); j++) {
435	dp.short_b = weights_short_b[j].value;
436
437	// Average over longuest_edgeC distribution
438	for(int k=0; k< (int)weights_long_c.size(); k++) {
439	dp.long_c = weights_long_c[k].value;
440	//Calculate surface averaged radius
441	//This is rough approximation.
442	suf_rad = sqrt(dp.short_a*dp.short_b/pi);
443
444	//Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
445	sum +=weights_short_a[i].weight* weights_short_b[j].weight
446	* weights_long_c[k].weight*DiamCyl(dp.long_c, suf_rad)/2.0;
447	norm += weights_short_a[i].weight* weights_short_b[j].weight*weights_long_c[k].weight;
448	}
449	}
450	}
451
452	if (norm != 0){
453	//return the averaged value
454	rad_out = sum/norm;}
455	else{
456	//return normal value
457	//Note: output of "DiamCyl(length,radius)" is DIAMETER.
458	rad_out = DiamCyl(dp.long_c, suf_rad)/2.0;}
459	return rad_out;
460
461	}

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

SasView

source: sasview/sansmodels/src/c_models/parallelepiped.cpp @ 530cc84

Download in other formats: