onion.c @ 01a870d

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 01a870d was c25f1fa, checked in by Mathieu Doucet <doucetm@…>, 13 years ago
Re #5 fixing onion compilation on MSVC
Property mode set to `100644`
File size: 8.0 KB

Line
1	/**
2	* Scattering model for a onion
3	*/
4
5	#include <math.h>
6	#include "onion.h"
7	#include <stdio.h>
8	#include <stdlib.h>
9	// some details can be found in sld_cal.c
10	double so_kernel(double dp[], double q) {
11	int n = dp[0];
12	double scale = dp[1];
13	double rad_core0 = dp[2];
14	double sld_core0 = dp[3];
15	double sld_solv = dp[4];
16	double background = dp[5];
17	int i,j;
18	double bes,fun,alpha,f,vol,vol_pre,vol_sub,qr,r,contr,f2;
19	double sign;
20	double pi;
21	double r0 = 0.0;
22
23	double *sld_out;
24	double *slope;
25	double *sld_in;
26	double *thick;
27	double *A;
28	int *fun_type;
29
30	sld_out = (double*) malloc(n+2);
31	slope = (double*) malloc(n+2);
32	sld_in = (double*) malloc(n+2);
33	thick = (double*) malloc(n+2);
34	A = (double*) malloc(n+2);
35	fun_type = (int*) malloc(n+2);
36
37	for (i =1; i<=n; i++){
38	sld_out[i] = dp[i+5];
39	sld_in[i] = dp[i+15];
40	A[i] = dp[i+25];
41	thick[i] = dp[i+35];
42	fun_type[i] = dp[i+45];
43	}
44	sld_out[0] = sld_core0;
45	sld_out[n+1] = sld_solv;
46	sld_in[0] = sld_core0;
47	sld_in[n+1] = sld_solv;
48	thick[0] = rad_core0;
49	thick[n+1] = 1e+10;
50	A[0] = 0.0;
51	A[n+1] = 0.0;
52	fun_type[0] = 0;
53	fun_type[n+1] = 0;
54
55
56	pi = 4.0*atan(1.0);
57	f = 0.0;
58	r = 0.0;
59	vol = 0.0;
60	vol_pre = 0.0;
61	vol_sub = 0.0;
62
63	for (i =0; i<= n+1; i++){
64	if (thick[i] == 0.0){
65	continue;
66	}
67	if (fun_type[i]== 0 ){
68	slope[i] = 0.0;
69	A[i] = 0.0;
70	}
71	vol_pre = vol;
72	switch(fun_type[i]){
73	case 2 :
74	r0 = r;
75	if (A[i] == 0.0){
76	slope[i] = 0.0;
77	}
78	else{
79	slope[i] = (sld_out[i]-sld_in[i])/(exp(A[i])-1.0);
80	}
81	for (j=0; j<2; j++){
82	if ( i == 0 && j == 0){
83	continue;
84	}
85	if (i == n+1 && j == 1){
86	continue;
87	}
88	if ( j == 1){
89	sign = 1.0;
90	r += thick[i];
91	}
92	else{
93	sign = -1.0;
94	}
95	qr = q * r;
96	alpha = A[i] * r/thick[i];
97	fun = 0.0;
98	if(qr == 0.0){
99	fun = sign * 1.0;
100	bes = sign * 1.0;
101	}
102	else{
103	if (fabs(A[i]) > 0.0 ){
104	fun = 3.0 * ((alphaalpha - qr qr) * sin(qr) - 2.0 * alpha * qr * cos(qr))/ ((alpha * alpha + qr * qr) * (alpha * alpha + qr * qr) * qr);
105	fun = fun - 3.0 * (alpha * sin(qr) - qr * cos(qr)) / ((alpha * alpha + qr * qr) * qr);
106	fun = - sign *fun;
107	bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
108	}
109	else {
110	fun = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
111	bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
112	}
113	}
114	contr = slope[i]exp(A[i](r-r0)/thick[i]);
115
116	vol = 4.0 * pi / 3.0 * r * r * r;
117	//if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && A[i]==0.0){
118	// vol_sub += (vol_pre - vol);
119	//}
120	f += vol * (contr * (fun) + (sld_in[i]-slope[i]) * bes);
121	}
122	break;
123	default :
124	if (fun_type[i]==0){
125	slope[i] = 0.0;
126	}
127	else{
128	slope[i]= (sld_out[i] -sld_in[i])/thick[i];
129	}
130	contr = sld_in[i]-slope[i]*r;
131	for (j=0; j<2; j++){
132	if ( i == 0 && j == 0){
133	continue;
134	}
135	if (i == n+1 && j == 1){
136	continue;
137	}
138	if ( j == 1){
139	sign = 1.0;
140	r += thick[i];
141	}
142	else{
143	sign = -1.0;
144	}
145
146	qr = q * r;
147	fun = 0.0;
148	if(qr == 0.0){
149	bes = sign * 1.0;
150	}
151	else{
152	bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
153	if (fabs(slope[i]) > 0.0 ){
154	fun = sign * 3.0 * r * (2.0qrsin(qr)-((qrqr)-2.0)cos(qr))/(qr * qr * qr * qr);
155	}
156	}
157	vol = 4.0 * pi / 3.0 * r * r * r;
158	//if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
159	// vol_sub += (vol_pre - vol);
160	//}
161	f += vol * (bes * contr + fun * slope[i]);
162	}
163	break;
164	}
165
166	}
167	//vol += vol_sub;
168	f2 = f * f / vol * 1.0e8;
169	f2 *= scale;
170	f2 += background;
171
172	free(sld_out);
173	free(slope);
174	free(sld_in);
175	free(thick);
176	free(A);
177	free(fun_type);
178
179	return (f2);
180	}
181	/**
182	* Function to evaluate 1D scattering function
183	* @param pars: parameters of the sphere
184	* @param q: q-value
185	* @return: function value
186	*/
187	double onion_analytical_1D(OnionParameters *pars, double q) {
188	double dp[56];
189
190	dp[0] = pars->n_shells;
191	dp[1] = pars->scale;
192	dp[2] = pars->rad_core0;
193	dp[3] = pars->sld_core0;
194	dp[4] = pars->sld_solv;
195	dp[5] = pars->background;
196
197	dp[6] = pars->sld_out_shell1;
198	dp[7] = pars->sld_out_shell2;
199	dp[8] = pars->sld_out_shell3;
200	dp[9] = pars->sld_out_shell4;
201	dp[10] = pars->sld_out_shell5;
202	dp[11] = pars->sld_out_shell6;
203	dp[12] = pars->sld_out_shell7;
204	dp[13] = pars->sld_out_shell8;
205	dp[14] = pars->sld_out_shell9;
206	dp[15] = pars->sld_out_shell10;
207
208	dp[16] = pars->sld_in_shell1;
209	dp[17] = pars->sld_in_shell2;
210	dp[18] = pars->sld_in_shell3;
211	dp[19] = pars->sld_in_shell4;
212	dp[20] = pars->sld_in_shell5;
213	dp[21] = pars->sld_in_shell6;
214	dp[22] = pars->sld_in_shell7;
215	dp[23] = pars->sld_in_shell8;
216	dp[24] = pars->sld_in_shell9;
217	dp[25] = pars->sld_in_shell10;
218
219	dp[26] = pars->A_shell1;
220	dp[27] = pars->A_shell2;
221	dp[28] = pars->A_shell3;
222	dp[29] = pars->A_shell4;
223	dp[30] = pars->A_shell5;
224	dp[31] = pars->A_shell6;
225	dp[32] = pars->A_shell7;
226	dp[33] = pars->A_shell8;
227	dp[34] = pars->A_shell9;
228	dp[35] = pars->A_shell10;
229
230	dp[36] = pars->thick_shell1;
231	dp[37] = pars->thick_shell2;
232	dp[38] = pars->thick_shell3;
233	dp[39] = pars->thick_shell4;
234	dp[40] = pars->thick_shell5;
235	dp[41] = pars->thick_shell6;
236	dp[42] = pars->thick_shell7;
237	dp[43] = pars->thick_shell8;
238	dp[44] = pars->thick_shell9;
239	dp[45] = pars->thick_shell10;
240
241	dp[46] = pars->func_shell1;
242	dp[47] = pars->func_shell2;
243	dp[48] = pars->func_shell3;
244	dp[49] = pars->func_shell4;
245	dp[50] = pars->func_shell5;
246	dp[51] = pars->func_shell6;
247	dp[52] = pars->func_shell7;
248	dp[53] = pars->func_shell8;
249	dp[54] = pars->func_shell9;
250	dp[55] = pars->func_shell10;
251
252	return so_kernel(dp, q);
253	}
254
255	/**
256	* Function to evaluate 2D scattering function
257	* @param pars: parameters of the sphere
258	* @param q: q-value
259	* @return: function value
260	*/
261	double onion_analytical_2D(OnionParameters *pars, double q, double phi) {
262	return onion_analytical_1D(pars,q);
263	}
264
265	double onion_analytical_2DXY(OnionParameters *pars, double qx, double qy){
266	return onion_analytical_1D(pars,sqrt(qxqx+qyqy));
267	}

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

SasView

source: sasview/sansmodels/src/sans/models/c_extensions/onion.c @ 01a870d

Download in other formats: