Changeset 0164899a in sasview for sansmodels/src/sans/models/c_extensions/refl.c
- Timestamp:
- Nov 1, 2010 4:22:12 PM (14 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- 4b3d25b
- Parents:
- 6cda91f
- File:
-
- 1 edited
-
sansmodels/src/sans/models/c_extensions/refl.c (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/sans/models/c_extensions/refl.c
r339ce67 r0164899a 1 1 /** 2 * Scattering model for a sphere2 * model for NR 3 3 */ 4 4 // The original code, of which work was not DANSE funded, 5 // was provided by J. Cho. 5 6 #include <math.h> 6 7 #include "refl.h" 8 #include "libmultifunc/librefl.h" 7 9 #include <stdio.h> 8 10 #include <stdlib.h> … … 10 12 #define lamda 4.62 11 13 12 typedef struct {13 double re;14 double im;15 } complex;16 17 typedef struct {18 complex a;19 complex b;20 complex c;21 complex d;22 } matrix;23 24 complex cassign(real, imag)25 double real, imag;26 {27 complex x;28 x.re = real;29 x.im = imag;30 return x;31 }32 33 34 complex cadd(x,y)35 complex x,y;36 {37 complex z;38 z.re = x.re + y.re;39 z.im = x.im + y.im;40 return z;41 }42 43 complex rcmult(x,y)44 double x;45 complex y;46 {47 complex z;48 z.re = x*y.re;49 z.im = x*y.im;50 return z;51 }52 53 complex csub(x,y)54 complex x,y;55 {56 complex z;57 z.re = x.re - y.re;58 z.im = x.im - y.im;59 return z;60 }61 62 63 complex cmult(x,y)64 complex x,y;65 {66 complex z;67 z.re = x.re*y.re - x.im*y.im;68 z.im = x.re*y.im + x.im*y.re;69 return z;70 }71 72 complex cdiv(x,y)73 complex x,y;74 {75 complex z;76 z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);77 z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);78 return z;79 }80 81 complex cexp(b)82 complex b;83 {84 complex z;85 double br,bi;86 br=b.re;87 bi=b.im;88 z.re = exp(br)*cos(bi);89 z.im = exp(br)*sin(bi);90 return z;91 }92 93 94 complex csqrt(z) /* see Schaum`s Math Handbook p. 22, 6.6 and 6.10 */95 complex z;96 {97 complex c;98 double zr,zi,x,y,r,w;99 100 zr=z.re;101 zi=z.im;102 103 if (zr==0.0 && zi==0.0)104 {105 c.re=0.0;106 c.im=0.0;107 return c;108 }109 else110 {111 x=fabs(zr);112 y=fabs(zi);113 if (x>y)114 {115 r=y/x;116 w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));117 }118 else119 {120 r=x/y;121 w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));122 }123 if (zr >=0.0)124 {125 c.re=w;126 c.im=zi/(2.0*w);127 }128 else129 {130 c.im=(zi >= 0) ? w : -w;131 c.re=zi/(2.0*c.im);132 }133 return c;134 }135 }136 137 complex ccos(b)138 complex b;139 {140 complex neg,negb,zero,two,z,i,bi,negbi;141 zero = cassign(0.0,0.0);142 two = cassign(2.0,0.0);143 i = cassign(0.0,1.0);144 bi = cmult(b,i);145 negbi = csub(zero,bi);146 z = cdiv(cadd(cexp(bi),cexp(negbi)),two);147 return z;148 }149 150 double errfunc(n_sub, i)151 double n_sub;152 int i;153 {154 double bin_size, ind, func;155 ind = i;156 // i range = [ -4..4], x range = [ -2.5..2.5]157 bin_size = n_sub/2.0/2.5; //size of each sub-layer158 // rescale erf so that 0 < erf < 1 in -2.5 <= x <= 2.5159 func = (erf(ind/bin_size)/2.0+0.5);160 return func;161 }162 163 double linefunc(n_sub, i)164 double n_sub;165 int i;166 {167 double bin_size, ind, func;168 ind = i + 0.5;169 // i range = [ -4..4], x range = [ -2.5..2.5]170 bin_size = 1.0/n_sub; //size of each sub-layer171 // rescale erf so that 0 < erf < 1 in -2.5 <= x <= 2.5172 func = ((ind + floor(n_sub/2.0))*bin_size);173 return func;174 }175 176 double parabolic_r(n_sub, i)177 double n_sub;178 int i;179 {180 double bin_size, ind, func;181 ind = i + 0.5;182 // i range = [ -4..4], x range = [ 0..1]183 bin_size = 1.0/n_sub; //size of each sub-layer; n_sub = 0 is a singular point (error)184 func = ((ind + floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size);185 return func;186 }187 188 double parabolic_l(n_sub, i)189 double n_sub;190 int i;191 {192 double bin_size,ind, func;193 ind = i + 0.5;194 bin_size = 1.0/n_sub; //size of each sub-layer; n_sub = 0 is a singular point (error)195 func =1.0-(((ind + floor(n_sub/2.0))*bin_size) - 1.0) *(((ind + floor(n_sub/2.0))*bin_size) - 1.0);196 return func;197 }198 199 double cubic_r(n_sub, i)200 double n_sub;201 int i;202 {203 double bin_size,ind, func;204 ind = i + 0.5;205 // i range = [ -4..4], x range = [ 0..1]206 bin_size = 1.0/n_sub; //size of each sub-layer; n_sub = 0 is a singular point (error)207 func = ((ind+ floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size);208 return func;209 }210 211 double cubic_l(n_sub, i)212 double n_sub;213 int i;214 {215 double bin_size,ind, func;216 ind = i + 0.5;217 bin_size = 1.0/n_sub; //size of each sub-layer; n_sub = 0 is a singular point (error)218 func = 1.0+(((ind + floor(n_sub/2.0)))*bin_size - 1.0)*(((ind + floor(n_sub/2.0)))*bin_size - 1.0)*(((ind + floor(n_sub/2.0)))*bin_size - 1.0);219 return func;220 }221 222 double interfunc(fun_type, n_sub, i, sld_l, sld_r)223 double n_sub, sld_l, sld_r;224 int fun_type, i;225 {226 double sld_i, func;227 switch(fun_type){228 case 1 :229 func = linefunc(n_sub, i);230 break;231 case 2 :232 func = parabolic_r(n_sub, i);233 break;234 case 3 :235 func = parabolic_l(n_sub, i);236 break;237 case 4 :238 func = cubic_r(n_sub, i);239 break;240 case 5 :241 func = cubic_l(n_sub, i);242 break;243 default:244 func = errfunc(n_sub, i);245 break;246 }247 if (sld_r>sld_l){248 sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);249 }250 else if (sld_r<sld_l){251 func = 1.0-func;252 sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);253 }254 else{255 sld_i = sld_r;256 }257 return sld_i;258 }259 14 260 15 double re_kernel(double dp[], double q) { … … 268 23 double background = dp[6]; 269 24 270 double sld[n+2], sld_im[n+2],thick_inter[n+2],thick[n+2],total_thick;25 double sld[n+2],thick_inter[n+2],thick[n+2],total_thick; 271 26 fun_type[0] = dp[3]; 272 27 for (i =1; i<=n; i++){ … … 275 30 thick[i] = dp[i+26]; 276 31 fun_type[i] = dp[i+36]; 277 sld_im[i] = dp[i+46];278 279 32 total_thick += thick[i] + thick_inter[i]; 280 33 } 281 34 sld[0] = sld_sub; 282 35 sld[n+1] = sld_super; 283 sld_im[0] = fabs(dp[0+56]); 284 sld_im[n+1] = fabs(dp[1+56]); 36 285 37 thick[0] = total_thick/5.0; 286 38 thick[n+1] = total_thick/5.0; … … 288 40 thick_inter[n+1] = 0.0; 289 41 290 double nsl=21.0; //nsl = Num_sub_layer: MUST ODD number in double //no other number works now291 int n_s , floor_nsl;42 double nsl=21.0; //nsl = Num_sub_layer: 43 int n_s; 292 44 double sld_i,sldim_i,dz,phi,R,ko2; 293 double sign,erfunc ;45 double sign,erfunc, fun; 294 46 double pi; 295 47 complex inv_n,phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,zero,two,n_sub,n_sup,knp1,Xnp1; … … 299 51 two= cassign(0.0,-2.0); 300 52 301 floor_nsl = floor(nsl/2.0);302 303 53 //Checking if floor is available. 304 54 //no imaginary sld inputs in this function yet 305 n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi), pow(lamda,2.0)/(2.0*pi)*sld_im[0]);306 n_sup=cassign(1.0-sld_super*pow(lamda,2.0)/(2.0*pi), pow(lamda,2.0)/(2.0*pi)*sld_im[n+1]);55 n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi),0.0); 56 n_sup=cassign(1.0-sld_super*pow(lamda,2.0)/(2.0*pi),0.0); 307 57 ko2 = pow(2.0*pi/lamda,2.0); 308 58 … … 318 68 // iteration for # of layers +sub from the top 319 69 for (i=1;i<=n+1; i++){ 70 if (fun_type[i-1]==1) 71 fun = 5; 72 else 73 fun = 0; 320 74 //iteration for 9 sub-layers 321 75 for (j=0;j<2;j++){ 322 for (n_s= -floor_nsl;n_s<=floor_nsl; n_s++){76 for (n_s=0;n_s<nsl; n_s++){ 323 77 if (j==1){ 324 78 if (i==n+1) … … 326 80 dz = thick[i]; 327 81 sld_i = sld[i]; 328 sldim_i = sld_im[i];329 82 } 330 83 else{ 331 dz = thick_inter[i-1]/nsl; 84 dz = thick_inter[i-1]/nsl;//nsl; 332 85 if (sld[i-1] == sld[i]){ 333 86 sld_i = sld[i]; 334 87 } 335 88 else{ 336 sld_i = interfunc(fun_type[i-1],nsl, n_s, sld[i-1], sld[i]); 337 } 338 if (sld_im[i-1] == sld_im[i]){ 339 sldim_i = sld_im[i]; 340 } 341 else{ 342 sldim_i = interfunc(fun_type[i-1],nsl, n_s, sld_im[i-1], sld_im[i]); 89 sld_i = intersldfunc(fun,nsl, n_s+0.5, 2.5, sld[i-1], sld[i]); 343 90 } 344 91 } 345 nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi), pow(lamda,2.0)/(2.0*pi)*sldim_i);92 nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),0.0); 346 93 nn2=cmult(nn,nn); 347 94 … … 375 122 } 376 123 /** 377 * Function to evaluate 1D scatteringfunction378 * @param pars: parameters of the sphere124 * Function to evaluate NR function 125 * @param pars: parameters 379 126 * @param q: q-value 380 127 * @return: function value 381 128 */ 382 129 double refl_analytical_1D(ReflParameters *pars, double q) { 383 double dp[ 59];130 double dp[47]; 384 131 385 132 dp[0] = pars->n_layers; … … 387 134 dp[2] = pars->thick_inter0; 388 135 dp[3] = pars->func_inter0; 389 dp[4] = pars->sld_ sub0;136 dp[4] = pars->sld_bottom0; 390 137 dp[5] = pars->sld_medium; 391 138 dp[6] = pars->background; … … 435 182 dp[46] = pars->func_inter10; 436 183 437 dp[47] = pars->sldIM_flat1;438 dp[48] = pars->sldIM_flat2;439 dp[49] = pars->sldIM_flat3;440 dp[50] = pars->sldIM_flat4;441 dp[51] = pars->sldIM_flat5;442 dp[52] = pars->sldIM_flat6;443 dp[53] = pars->sldIM_flat7;444 dp[54] = pars->sldIM_flat8;445 dp[55] = pars->sldIM_flat9;446 dp[56] = pars->sldIM_flat10;447 448 dp[57] = pars->sldIM_sub0;449 dp[58] = pars->sldIM_medium;450 451 184 return re_kernel(dp, q); 452 185 } 453 186 454 187 /** 455 * Function to evaluate 2D scatteringfunction456 * @param pars: parameters of the sphere188 * Function to evaluate NR function 189 * @param pars: parameters of NR 457 190 * @param q: q-value 458 191 * @return: function value
Note: See TracChangeset
for help on using the changeset viewer.
