source: sasview/sansmodels/src/c_models/refl.cpp @ 9e9e9a5

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 9e9e9a5 was e08bd5b, checked in by Jae Cho <jhjcho@…>, 13 years ago

c models fix: scale fix for P*S

  • Property mode set to 100644
File size: 7.5 KB
Line 
1
2#include <math.h>
3#include "parameters.hh"
4#include <stdio.h>
5#include <stdlib.h>
6#include "refl.h"
7using namespace std;
8
9extern "C" {
10#include "libmultifunc/librefl.h"
11}
12
13#define lamda 4.62
14
15static double re_kernel(double dp[], double q) {
16  int n = dp[0];
17  int i,j;
18
19  double scale = dp[1];
20  double thick_inter_sub = dp[2];
21  double sld_sub = dp[4];
22  double sld_super = dp[5];
23  double background = dp[6];
24
25  double total_thick=0.0;
26  double nsl=21.0; //nsl = Num_sub_layer:
27  int n_s;
28  double sld_i,dz,phi,R,ko2;
29  double fun;
30  double pi;
31
32  double* sld;
33  double* thick_inter;
34  double* thick;
35  int*fun_type;
36  complex  phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,two,n_sub,n_sup,knp1,Xnp1;
37
38  sld = (double*)malloc((n+2)*sizeof(double));
39  thick_inter = (double*)malloc((n+2)*sizeof(double));
40  thick = (double*)malloc((n+2)*sizeof(double));
41  fun_type = (int*)malloc((n+2)*sizeof(int));
42
43  fun_type[0] = dp[3];
44  for (i =1; i<=n; i++){
45    sld[i] = dp[i+6];
46    thick_inter[i]= dp[i+16];
47    thick[i] = dp[i+26];
48    fun_type[i] = dp[i+36];
49    total_thick += thick[i] + thick_inter[i];
50  }
51  sld[0] = sld_sub;
52  sld[n+1] = sld_super;
53
54  thick[0] = total_thick/5.0;
55  thick[n+1] = total_thick/5.0;
56  thick_inter[0] = thick_inter_sub;
57  thick_inter[n+1] = 0.0;
58
59  pi = 4.0*atan(1.0);
60  Xn = cassign(0.0,0.0);
61  one = cassign(1.0,0.0);
62  two = cassign(0.0,-2.0);
63
64  //Checking if floor is available.
65  //no imaginary sld inputs in this function yet
66  n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi),0.0);
67  n_sup=cassign(1.0-sld_super*pow(lamda,2.0)/(2.0*pi),0.0);
68  ko2 = pow(2.0*pi/lamda,2.0);
69
70  phi = asin(lamda*q/(4.0*pi));
71  phi1 = cplx_div(rcmult(phi,one),n_sup);
72  alpha = cplx_mult(n_sup,cplx_cos(phi1));
73  alpha2 = cplx_mult(alpha,alpha);
74
75  nnp1=n_sub;
76  knp1=cplx_sqrt(rcmult(ko2,cplx_sub(cplx_mult(nnp1,nnp1),alpha2)));  //nnp1*ko*sin(phinp1)
77  Xnp1=cassign(0.0,0.0);
78  dz = 0.0;
79  // iteration for # of layers +sub from the top
80  for (i=1;i<=n+1; i++){
81    if (fun_type[i-1]==1)
82      fun = 5;
83    else
84      fun = 0;
85    //iteration for 9 sub-layers
86    for (j=0;j<2;j++){
87      for (n_s=0;n_s<nsl; n_s++){
88        if (j==1){
89          if (i==n+1)
90            break;
91          dz = thick[i];
92          sld_i = sld[i];
93        }
94        else{
95          dz = thick_inter[i-1]/nsl;//nsl;
96          if (sld[i-1] == sld[i]){
97            sld_i = sld[i];
98          }
99          else{
100            sld_i = intersldfunc(fun,nsl, n_s+0.5, 2.5, sld[i-1], sld[i]);
101          }
102        }
103        nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),0.0);
104        nn2=cplx_mult(nn,nn);
105
106        kn=cplx_sqrt(rcmult(ko2,cplx_sub(nn2,alpha2)));        //nn*ko*sin(phin)
107        an=cplx_exp(rcmult(dz,cplx_mult(two,kn)));
108
109        fnm=cplx_sub(kn,knp1);
110        fnp=cplx_add(kn,knp1);
111        rn=cplx_div(fnm,fnp);
112        Xn=cplx_mult(an,cplx_div(cplx_add(rn,Xnp1),cplx_add(one,cplx_mult(rn,Xnp1))));    //Xn=an*((rn+Xnp1*anp1)/(1+rn*Xnp1*anp1))
113
114        Xnp1=Xn;
115        knp1=kn;
116
117        if (j==1)
118          break;
119      }
120    }
121  }
122  R=pow(Xn.re,2.0)+pow(Xn.im,2.0);
123  // This temperarily fixes the total reflection for Rfunction and linear.
124  // ToDo: Show why it happens that Xn.re=0 and Xn.im >1!
125  if (Xn.im == 0.0){
126    R=1.0;
127  }
128  R *= scale;
129  R += background;
130
131  free(sld);
132  free(thick_inter);
133  free(thick);
134  free(fun_type);
135
136  return R;
137
138}
139ReflModel :: ReflModel() {
140  n_layers = Parameter(1.0);
141  scale = Parameter(1.0);
142  thick_inter0 = Parameter(1.0);
143  func_inter0 = Parameter(0);
144  sld_bottom0 = Parameter(2.07e-06);
145  sld_medium = Parameter(1.0e-06);
146  background = Parameter(0.0);
147
148
149  sld_flat1 = Parameter(3.0e-06);
150  sld_flat2 = Parameter(3.5e-06);
151  sld_flat3 = Parameter(4.0e-06);
152  sld_flat4 = Parameter(3.5e-06);
153  sld_flat5 = Parameter(4.0e-06);
154  sld_flat6 = Parameter(3.5e-06);
155  sld_flat7 = Parameter(4.0e-06);
156  sld_flat8 = Parameter(3.5e-06);
157  sld_flat9 = Parameter(4.0e-06);
158  sld_flat10 = Parameter(3.5e-06);
159
160
161  thick_inter1 = Parameter(1);
162  thick_inter2 = Parameter(1);
163  thick_inter3 = Parameter(1);
164  thick_inter4 = Parameter(1);
165  thick_inter5 = Parameter(1);
166  thick_inter6 = Parameter(1);
167  thick_inter7 = Parameter(1);
168  thick_inter8 = Parameter(1);
169  thick_inter9 = Parameter(1);
170  thick_inter10 = Parameter(1);
171
172
173  thick_flat1 = Parameter(15);
174  thick_flat2 = Parameter(100);
175  thick_flat3 = Parameter(100);
176  thick_flat4 = Parameter(100);
177  thick_flat5 = Parameter(100);
178  thick_flat6 = Parameter(100);
179  thick_flat7 = Parameter(100);
180  thick_flat8 = Parameter(100);
181  thick_flat9 = Parameter(100);
182  thick_flat10 = Parameter(100);
183
184
185  func_inter1 = Parameter(0);
186  func_inter2 = Parameter(0);
187  func_inter3 = Parameter(0);
188  func_inter4 = Parameter(0);
189  func_inter5 = Parameter(0);
190  func_inter6 = Parameter(0);
191  func_inter7 = Parameter(0);
192  func_inter8 = Parameter(0);
193  func_inter9 = Parameter(0);
194  func_inter10 = Parameter(0);
195
196}
197
198/**
199 * Function to evaluate 1D NR function
200 * @param q: q-value
201 * @return: function value
202 */
203double ReflModel :: operator()(double q) {
204  double dp[47];
205  // Fill parameter array for IGOR library
206  // Add the background after averaging
207  dp[0] = n_layers();
208  dp[1] = scale();
209  dp[2] = thick_inter0();
210  dp[3] = func_inter0();
211  dp[4] = sld_bottom0();
212  dp[5] = sld_medium();
213  dp[6] = background();
214
215  dp[7] = sld_flat1();
216  dp[8] = sld_flat2();
217  dp[9] = sld_flat3();
218  dp[10] = sld_flat4();
219  dp[11] = sld_flat5();
220  dp[12] = sld_flat6();
221  dp[13] = sld_flat7();
222  dp[14] = sld_flat8();
223  dp[15] = sld_flat9();
224  dp[16] = sld_flat10();
225
226  dp[17] = thick_inter1();
227  dp[18] = thick_inter2();
228  dp[19] = thick_inter3();
229  dp[20] = thick_inter4();
230  dp[21] = thick_inter5();
231  dp[22] = thick_inter6();
232  dp[23] = thick_inter7();
233  dp[24] = thick_inter8();
234  dp[25] = thick_inter9();
235  dp[26] = thick_inter10();
236
237  dp[27] = thick_flat1();
238  dp[28] = thick_flat2();
239  dp[29] = thick_flat3();
240  dp[30] = thick_flat4();
241  dp[31] = thick_flat5();
242  dp[32] = thick_flat6();
243  dp[33] = thick_flat7();
244  dp[34] = thick_flat8();
245  dp[35] = thick_flat9();
246  dp[36] = thick_flat10();
247
248  dp[37] = func_inter1();
249  dp[38] = func_inter2();
250  dp[39] = func_inter3();
251  dp[40] = func_inter4();
252  dp[41] = func_inter5();
253  dp[42] = func_inter6();
254  dp[43] = func_inter7();
255  dp[44] = func_inter8();
256  dp[45] = func_inter9();
257  dp[46] = func_inter10();
258
259  // Get the dispersion points for the radius
260  //vector<WeightPoint> weights_thick;
261  //thick_inter0.get_weights(weights_thick_inter0);
262
263
264  return re_kernel(dp,q);
265}
266
267/**
268 * Function to evaluate 2D NR function
269 * @param q_x: value of Q along x
270 * @param q_y: value of Q along y
271 * @return: function value
272 */
273double ReflModel :: operator()(double qx, double qy) {
274  // For 2D set qy as q, ignoring qx.
275  double q = qy;//sqrt(qx*qx + qy*qy);
276  if (q < 0.0){
277    return 0.0;
278  }
279  return (*this).operator()(q);
280}
281
282/**
283 * Function to evaluate 2D NR function
284 * @param pars: parameters of the sphere
285 * @param q: q-value
286 * @param phi: angle phi
287 * @return: function value
288 */
289double ReflModel :: evaluate_rphi(double q, double phi) {
290  return (*this).operator()(q);
291}
292
293/**
294 * Function to calculate effective radius
295 * @return: effective radius value
296 */
297double ReflModel :: calculate_ER() {
298  //NOT implemented yet!!!
299  return 0.0;
300}
301double ReflModel :: calculate_VR() {
302  return 1.0;
303}
Note: See TracBrowser for help on using the repository browser.