refl_adv.c @ 8bac371

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

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

SasView

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

Download in other formats: