onion.c @ 74b1770

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 74b1770 was 339ce67, checked in by Jae Cho <jhjcho@…>, 14 years ago
added some models and tests
Property mode set to `100644`
File size: 7.7 KB

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

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

SasView

source: sasview/sansmodels/src/sans/models/c_extensions/onion.c @ 74b1770

Download in other formats: