Changeset 0164899a in sasview for sansmodels/src/sans/models/c_extensions
- 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
- Location:
- sansmodels/src/sans/models/c_extensions
- Files:
-
- 6 added
- 3 edited
-
libmultifunc/librefl.c (added)
-
libmultifunc/librefl.h (added)
-
onion.c (modified) (5 diffs)
-
refl.c (modified) (11 diffs)
-
refl.h (modified) (4 diffs)
-
sld_cal.c (added)
-
sld_cal.h (added)
-
spheresld.c (added)
-
spheresld.h (added)
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/sans/models/c_extensions/onion.c
r339ce67 r0164899a 1 1 /** 2 * Scattering model for a sphere2 * Scattering model for a onion 3 3 */ 4 4 … … 7 7 #include <stdio.h> 8 8 #include <stdlib.h> 9 9 // some details can be found in sld_cal.c 10 10 double so_kernel(double dp[], double q) { 11 11 int n = dp[0]; … … 107 107 108 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 }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 112 f += vol * (contr * (fun) + (sld_in[i]-slope[i]) * bes); 113 113 } … … 148 148 } 149 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 }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 153 f += vol * (bes * contr + fun * slope[i]); 154 154 } … … 157 157 158 158 } 159 vol += vol_sub;159 //vol += vol_sub; 160 160 f2 = f * f / vol * 1.0e8; 161 161 f2 *= scale; -
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 -
sansmodels/src/sans/models/c_extensions/refl.h
r339ce67 r0164899a 1 // The original code, of which work was not DANSE funded, 2 // was provided by J. Cho. 1 3 #if !defined(o_h) 2 4 #define refl_h … … 7 9 //[PYTHONCLASS] = ReflModel 8 10 //[DISP_PARAMS] = thick_inter0 9 //[DESCRIPTION] =<text>Form factor of mutishells normalized by the volume. Here each shell is described 10 // by an exponential function; 11 // I) 12 // For A_shell != 0, 13 // f(r) = B*exp(A_shell*(r-r_in)/thick_shell)+C 14 // where 15 // B=(sld_out-sld_in)/(exp(A_shell)-1) 16 // C=sld_in-B. 17 // Note that in the above case, 18 // the function becomes a linear function 19 // as A_shell --> 0+ or 0-. 20 // II) 21 // For the exact point of A_shell == 0, 22 // f(r) = sld_in ,i.e., it crosses over flat function 23 // Note that the 'sld_out' becaomes NULL in this case. 24 // 25 // background:background, 26 // rad_core: radius of sphere(core) 27 // thick_shell#:the thickness of the shell# 28 // sld_core: the SLD of the sphere 29 // sld_solv: the SLD of the solvent 30 // sld_shell: the SLD of the shell# 31 // A_shell#: the coefficient in the exponential function 11 //[DESCRIPTION] =<text>Calculate neutron reflectivity using the Parratt iterative formula 12 // Parameters: 13 // background:background 14 // scale: scale factor 15 // sld_bottom0: the SLD of the substrate 16 // sld_medium: the SLD of the incident medium 17 // or superstrate 18 // sld_flatN: the SLD of the flat region of 19 // the N'th layer 20 // thick_flatN: the thickness of the flat 21 // region of the N'th layer 22 // func_interN: the function used to describe 23 // the interface of the N'th layer 24 // thick_interN: the thickness of the interface 25 // of the N'th layer 26 // Note: the layer number starts to increase 27 // from the bottom (substrate) to the top. 32 28 // </text> 33 29 //[FIXED]= <text></text> … … 48 44 // [DEFAULT]=func_inter0= 0 49 45 double func_inter0; 50 /// sld_ sub0 [1/A^(2)]51 // [DEFAULT]=sld_ sub0= 2.07e-6 [1/A^(2)]52 double sld_ sub0;46 /// sld_bottom0 [1/A^(2)] 47 // [DEFAULT]=sld_bottom0= 2.07e-6 [1/A^(2)] 48 double sld_bottom0; 53 49 /// sld_medium [1/A^(2)] 54 50 // [DEFAULT]=sld_medium= 1.0e-6 [1/A^(2)] … … 142 138 double func_inter10; 143 139 144 // [DEFAULT]=sldIM_flat1=0145 double sldIM_flat1;146 // [DEFAULT]=sldIM_flat2=0147 double sldIM_flat2;148 // [DEFAULT]=sldIM_flat3=0149 double sldIM_flat3;150 // [DEFAULT]=sldIM_flat4=0151 double sldIM_flat4;152 // [DEFAULT]=sldIM_flat5=0153 double sldIM_flat5;154 // [DEFAULT]=sldIM_flat6=0155 double sldIM_flat6;156 // [DEFAULT]=sldIM_flat7=0157 double sldIM_flat7;158 // [DEFAULT]=sldIM_flat8=0159 double sldIM_flat8;160 // [DEFAULT]=sldIM_flat9=0161 double sldIM_flat9;162 // [DEFAULT]=sldIM_flat10=0163 double sldIM_flat10;164 // [DEFAULT]=sldIM_sub0=0165 double sldIM_sub0;166 // [DEFAULT]=sldIM_medium=0167 double sldIM_medium;168 140 169 141 } ReflParameters;
Note: See TracChangeset
for help on using the changeset viewer.
