refl_adv.c @ a24f530

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 a24f530 was a24f530, checked in by Mathieu Doucet <doucetm@…>, 12 years ago
Re #4 A few more warnings
Property mode set to `100644`
File size: 7.0 KB

Line
1	// The original code, of which work was not DANSE funded,
2	// was provided by J. Cho.
3	/**
4	* NR model Parratt method
5	*/
6	#include <math.h>
7	#include "refl_adv.h"
8	#include "libmultifunc/librefl.h"
9	#include <stdio.h>
10	#include <stdlib.h>
11
12	#define lamda 4.62
13
14
15	double re_adv_kernel(double dp[], double q) {
16	int n = dp[0];
17	int i,j;
18	double nsl;
19
20	double scale = dp[1];
21	double thick_inter_sub = dp[2];
22	double sld_sub = dp[4];
23	double sld_super = dp[5];
24	double background = dp[6];
25	double npts = dp[69]; //number of sub_layers in each interface
26
27	double total_thick=0.0;
28
29	int n_s;
30	double sld_i,sldim_i,dz,phi,R,ko2;
31	double pi;
32
33	int* fun_type;
34	double* sld;
35	double* sld_im;
36	double* thick_inter;
37	double* thick;
38	double* fun_coef;
39	complex phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,two,n_sub,n_sup,knp1,Xnp1;
40
41	fun_type = (int)malloc((n+2)sizeof(int));
42	sld = (double)malloc((n+2)sizeof(double));
43	sld_im = (double)malloc((n+2)sizeof(double));
44	thick_inter = (double)malloc((n+2)sizeof(double));
45	thick = (double)malloc((n+2)sizeof(double));
46	fun_coef = (double)malloc((n+2)sizeof(double));
47
48	fun_type[0] = dp[3];
49	fun_coef[0] = fabs(dp[70]);
50	for (i =1; i<=n; i++){
51	sld[i] = dp[i+6];
52	thick_inter[i]= dp[i+16];
53	thick[i] = dp[i+26];
54	fun_type[i] = dp[i+36];
55	sld_im[i] = dp[i+46];
56	fun_coef[i] = fabs(dp[i+56]);
57	//printf("type_func2 =%g\n",fun_coef[i]);
58
59	total_thick += thick[i] + thick_inter[i];
60	}
61	sld[0] = sld_sub;
62	sld[n+1] = sld_super;
63	sld_im[0] = fabs(dp[1+66]);
64	sld_im[n+1] = fabs(dp[2+66]);
65	thick[0] = total_thick/5.0;
66	thick[n+1] = total_thick/5.0;
67	thick_inter[0] = thick_inter_sub;
68	thick_inter[n+1] = 0.0;
69	fun_coef[n+1] = 0.0;
70
71	nsl=npts;//21.0; //nsl = Num_sub_layer: MUST ODD number in double //no other number works now
72
73	pi = 4.0*atan(1.0);
74	one = cassign(1.0,0.0);
75	Xn = cassign(0.0,0.0);
76	two = cassign(0.0,-2.0);
77
78	//Checking if floor is available.
79	//no imaginary sld inputs in this function yet
80	n_sub=cassign(1.0-sld_subpow(lamda,2.0)/(2.0pi),pow(lamda,2.0)/(2.0pi)sld_im[0]);
81	n_sup=cassign(1.0-sld_superpow(lamda,2.0)/(2.0pi),pow(lamda,2.0)/(2.0pi)sld_im[n+1]);
82	ko2 = pow(2.0*pi/lamda,2.0);
83
84	phi = asin(lamdaq/(4.0pi));
85	phi1 = cdiv(rcmult(phi,one),n_sup);
86	alpha = cmult(n_sup,ccos(phi1));
87	alpha2 = cmult(alpha,alpha);
88
89	nnp1=n_sub;
90	knp1=csqrt(rcmult(ko2,csub(cmult(nnp1,nnp1),alpha2))); //nnp1kosin(phinp1)
91	Xnp1=cassign(0.0,0.0);
92	dz = 0.0;
93	// iteration for # of layers +sub from the top
94	for (i=1;i<=n+1; i++){
95	//if (fun_coef[i]==0.0)
96	// // this condition protects an error in numerical multiplication
97	// fun_coef[i] = 1e-14;
98	//iteration for 9 sub-layers
99	for (j=0;j<2;j++){
100	for (n_s=0;n_s<nsl; n_s++){
101	// for flat layer
102	if (j==1){
103	if (i==n+1)
104	break;
105	dz = thick[i];
106	sld_i = sld[i];
107	sldim_i = sld_im[i];
108	}
109	// for interface
110	else{
111	dz = thick_inter[i-1]/nsl;
112	if (sld[i-1] == sld[i]){
113	sld_i = sld[i];
114	}
115	else{
116	sld_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld[i-1], sld[i]);
117	}
118	if (sld_im[i-1] == sld_im[i]){
119	sldim_i = sld_im[i];
120	}
121	else{
122	sldim_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld_im[i-1], sld_im[i]);
123	}
124	}
125	nn = cassign(1.0-sld_ipow(lamda,2.0)/(2.0pi),pow(lamda,2.0)/(2.0pi)sldim_i);
126	nn2=cmult(nn,nn);
127
128	kn=csqrt(rcmult(ko2,csub(nn2,alpha2))); //nnkosin(phin)
129	an=cexp(rcmult(dz,cmult(two,kn)));
130
131	fnm=csub(kn,knp1);
132	fnp=cadd(kn,knp1);
133	rn=cdiv(fnm,fnp);
134	Xn=cmult(an,cdiv(cadd(rn,Xnp1),cadd(one,cmult(rn,Xnp1)))); //Xn=an((rn+Xnp1anp1)/(1+rnXnp1anp1))
135
136	Xnp1=Xn;
137	knp1=kn;
138	// no for-loop for flat layer
139	if (j==1)
140	break;
141	}
142	}
143	}
144	R=pow(Xn.re,2.0)+pow(Xn.im,2.0);
145	// This temperarily fixes the total reflection for Rfunction and linear.
146	// ToDo: Show why it happens that Xn.re=0 and Xn.im >1!
147	if (Xn.im == 0.0 \|\| R > 1){
148	R=1.0;
149	}
150	R *= scale;
151	R += background;
152
153	free(fun_type);
154	free(sld);
155	free(sld_im);
156	free(thick_inter);
157	free(thick);
158	free(fun_coef);
159
160	return R;
161
162	}
163
164	/**
165	* Function to evaluate NR function
166	* @param pars: parameters of refl
167	* @param q: q-value
168	* @return: function value
169	*/
170
171	double refl_adv_analytical_1D(ReflAdvParameters *pars, double q) {
172	double dp[71];
173
174	dp[0] = pars->n_layers;
175	dp[1] = pars->scale;
176	dp[2] = pars->thick_inter0;
177	dp[3] = pars->func_inter0;
178	dp[4] = pars->sld_bottom0;
179	dp[5] = pars->sld_medium;
180	dp[6] = pars->background;
181
182	dp[7] = pars->sld_flat1;
183	dp[8] = pars->sld_flat2;
184	dp[9] = pars->sld_flat3;
185	dp[10] = pars->sld_flat4;
186	dp[11] = pars->sld_flat5;
187	dp[12] = pars->sld_flat6;
188	dp[13] = pars->sld_flat7;
189	dp[14] = pars->sld_flat8;
190	dp[15] = pars->sld_flat9;
191	dp[16] = pars->sld_flat10;
192
193	dp[17] = pars->thick_inter1;
194	dp[18] = pars->thick_inter2;
195	dp[19] = pars->thick_inter3;
196	dp[20] = pars->thick_inter4;
197	dp[21] = pars->thick_inter5;
198	dp[22] = pars->thick_inter6;
199	dp[23] = pars->thick_inter7;
200	dp[24] = pars->thick_inter8;
201	dp[25] = pars->thick_inter9;
202	dp[26] = pars->thick_inter10;
203
204	dp[27] = pars->thick_flat1;
205	dp[28] = pars->thick_flat2;
206	dp[29] = pars->thick_flat3;
207	dp[30] = pars->thick_flat4;
208	dp[31] = pars->thick_flat5;
209	dp[32] = pars->thick_flat6;
210	dp[33] = pars->thick_flat7;
211	dp[34] = pars->thick_flat8;
212	dp[35] = pars->thick_flat9;
213	dp[36] = pars->thick_flat10;
214
215	dp[37] = pars->func_inter1;
216	dp[38] = pars->func_inter2;
217	dp[39] = pars->func_inter3;
218	dp[40] = pars->func_inter4;
219	dp[41] = pars->func_inter5;
220	dp[42] = pars->func_inter6;
221	dp[43] = pars->func_inter7;
222	dp[44] = pars->func_inter8;
223	dp[45] = pars->func_inter9;
224	dp[46] = pars->func_inter10;
225
226	dp[47] = pars->sldIM_flat1;
227	dp[48] = pars->sldIM_flat2;
228	dp[49] = pars->sldIM_flat3;
229	dp[50] = pars->sldIM_flat4;
230	dp[51] = pars->sldIM_flat5;
231	dp[52] = pars->sldIM_flat6;
232	dp[53] = pars->sldIM_flat7;
233	dp[54] = pars->sldIM_flat8;
234	dp[55] = pars->sldIM_flat9;
235	dp[56] = pars->sldIM_flat10;
236
237	dp[57] = pars->nu_inter1;
238	dp[58] = pars->nu_inter2;
239	dp[59] = pars->nu_inter3;
240	dp[60] = pars->nu_inter4;
241	dp[61] = pars->nu_inter5;
242	dp[62] = pars->nu_inter6;
243	dp[63] = pars->nu_inter7;
244	dp[64] = pars->nu_inter8;
245	dp[65] = pars->nu_inter9;
246	dp[66] = pars->nu_inter10;
247
248	dp[67] = pars->sldIM_sub0;
249	dp[68] = pars->sldIM_medium;
250	dp[69] = pars->npts_inter;
251	dp[70] = pars->nu_inter0;
252
253	return re_adv_kernel(dp, q);
254	}
255
256	/**
257	* Function to evaluate NR function
258	* @param pars: parameters of NR
259	* @param q: q-value
260	* @return: function value
261	*/
262	double refl_adv_analytical_2D(ReflAdvParameters *pars, double q, double phi) {
263	return refl_adv_analytical_1D(pars,q);
264	}
265
266	double refl_adv_analytical_2DXY(ReflAdvParameters *pars, double qx, double qy){
267	return refl_adv_analytical_1D(pars,sqrt(qxqx+qyqy));
268	}

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

SasView

source: sasview/sansmodels/src/sans/models/c_extensions/refl_adv.c @ a24f530

Download in other formats: