Changeset 54bcd4a in sasmodels
- Timestamp:
- Aug 4, 2016 4:48:37 PM (8 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 745b7bb
- Parents:
- bd49c79
- Location:
- sasmodels
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified sasmodels/conversion_table.py ¶
rf67f26c r54bcd4a 688 688 } 689 689 ], 690 "magsphere": [691 "SphereModel",692 {693 "sld": "sldSph",694 "radius": "radius",695 "sld_solvent": "sldSolv",696 "sld_M0": "M0_sld_sph",697 "sphere_theta": "M_theta_sph",698 "sphere_phi": "M_phi_sph",699 "sld_solvent_M0": "M0_sld_solv",700 "solvent_theta": "M_theta_solv",701 "solvent_phi": "M_phi_solv",702 "in_spin": "Up_frac_i",703 "out_spin": "Up_frac_f",704 "spin_theta": "Up_theta"705 }706 ],707 "_spherepy": [708 "SphereModel",709 {710 "sld": "sldSph",711 "radius": "radius",712 "sld_solvent": "sldSolv"713 }714 ],715 690 "spherical_sld": [ 716 691 "SphericalSLDModel", 717 { 718 "radius_core": "rad_core0", 719 "sld_core": "sld_core0", 720 "sld_solvent": "sld_solv" 721 } 692 # Be lazy and use a generator expression to define 693 # sld1: sld_flat0, ... 694 # thickness1: thick_flat0, ... 695 # interface1: thick_inter0, ... 696 # shape1: func_inter0, ... 697 # nu1: nu_inter0, ... 698 # but override thickness1 => rad_cor0 and sld1 => sld_core0. 699 # Note: explicit key,value pairs given by **{...} override the 700 # keys from the gnerator expression ((k,v) for k,v in seq) when 701 # used as dict((generator), **{...}) 702 dict(((field_new+str(index+1),field_old+str(index)) 703 for field_new, field_old in [("sld","sld_flat"), 704 ("thickness","thick_flat"), 705 ("interface","thick_inter"), 706 ("shape","func_inter"), 707 ("nu","nu_inter"),] 708 for index in range(11)), 709 **{ 710 "n_shells": "n_shells", 711 "n_steps": "npts_inter", 712 "sld_solvent": "sld_solv", 713 "thickness1": "rad_core0", 714 "sld1": "sld_core0", 715 }) 722 716 ], 723 717 "squarewell": [ -
TabularUnified sasmodels/convert.py ¶
rbd49c79 r54bcd4a 229 229 oldpars['ndensity'] *= 1e15 230 230 elif name == 'spherical_sld': 231 for k in range(1, int(pars['n_shells'])+1): 232 _remove_pd(oldpars, 'thick_flat'+str(k), 'thick_flat') 233 _remove_pd(oldpars, 'thick_inter'+str(k), 'thick_inter') 231 oldpars["CONTROL"] -= 1 232 # remove polydispersity from shells 233 for k in range(1, 11): 234 _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness') 235 _remove_pd(oldpars, 'thick_inter'+str(k), 'interface') 236 # remove extra shells 237 for k in range(int(pars['n_shells']), 11): 238 oldpars.pop('sld_flat'+str(k), 0) 239 oldpars.pop('thick_flat'+str(k), 0) 240 oldpars.pop('thick_inter'+str(k), 0) 241 oldpars.pop('func_inter'+str(k), 0) 242 oldpars.pop('nu_inter'+str(k), 0) 234 243 elif name == 'core_multi_shell': 235 244 # kill extra shells … … 265 274 oldpars.pop(k, None) 266 275 276 #print("convert from",list(sorted(pars))) 277 #print("convert to",list(sorted(oldpars.items()))) 267 278 return oldpars 268 279 … … 323 334 elif name == 'spherical_sld': 324 335 pars['n_shells'] = math.ceil(pars['n_shells']) 325 pars['n pts_inter'] = math.ceil(pars['npts_inter'])326 pars['func_inter0'] = math.trunc(pars['func_inter0']+0.5)327 for k in range(1, 11):328 pars['func_inter%d'%k] = math.trunc(pars['func_inter%d'%k]+0.5)329 pars['thick _flat%d_pd_n'%k] = 0330 pars[' thick_inter%d_pd_n'%k] = 0331 336 pars['n_steps'] = math.ceil(pars['n_steps']) 337 for k in range(1, 12): 338 pars['shape%d'%k] = math.trunc(pars['shape%d'%k]+0.5) 339 for k in range(2, 12): 340 pars['thickness%d_pd_n'%k] = 0 341 pars['interface%d_pd_n'%k] = 0 342 -
TabularUnified sasmodels/models/spherical_sld.c ¶
rdd71228 r54bcd4a 1 1 static double form_volume( 2 2 int n_shells, 3 double radius_core, 4 double thick_inter0, 5 double thick_flat[], 6 double thick_inter[]) 3 double thickness[], 4 double interface[]) 7 5 { 8 int i; 9 double r = radius_core; 10 r += thick_inter0; 11 for (i=0; i < n_shells; i++) { 12 r += thick_inter[i]; 13 r += thick_flat[i]; 6 double r = 0.0; 7 for (int i=0; i < n_shells; i++) { 8 r += thickness[i] + interface[i]; 14 9 } 15 10 return M_4PI_3*cube(r); 16 11 } 17 12 13 static double blend(int shape, double nu, double z) 14 { 15 if (shape==0) { 16 const double num = sas_erf(nu * M_SQRT1_2 * (2.0*z - 1.0)); 17 const double denom = 2.0 * sas_erf(nu * M_SQRT1_2); 18 return num/denom + 0.5; 19 } else if (shape==1) { 20 return pow(z, nu); 21 } else if (shape==2) { 22 return 1.0 - pow(1. - z, nu); 23 } else if (shape==3) { 24 return expm1(-nu*z)/expm1(-nu); 25 } else if (shape==4) { 26 return expm1(nu*z)/expm1(nu); 27 } else { 28 return NAN; 29 } 30 } 18 31 19 static double sphere_sld_kernel( 32 static double f_linear(double q, double r, double contrast, double slope) 33 { 34 const double qr = q * r; 35 const double qrsq = qr * qr; 36 const double bes = sph_j1c(qr); 37 double sinqr, cosqr; 38 SINCOS(qr, sinqr, cosqr); 39 const double fun = 3.0*r*(2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq); 40 const double vol = M_4PI_3 * cube(r); 41 return vol*(bes*contrast + fun*slope); 42 } 43 44 static double Iq( 20 45 double q, 21 46 int n_shells, 22 int npts_inter,23 double radius_core,24 double sld_core,25 47 double sld_solvent, 26 double func_inter_core, 27 double thick_inter_core, 28 double nu_inter_core, 29 double sld_flat[], 30 double thick_flat[], 31 double func_inter[], 32 double thick_inter[], 33 double nu_inter[] ) { 48 double sld[], 49 double thickness[], 50 double interface[], 51 double shape[], 52 double nu[], 53 int n_steps) 54 { 55 // iteration for # of shells + core + solvent 56 double f=0.0; 57 double r=0.0; 58 for (int shell=0; shell<n_shells; shell++){ 59 const double sld_l = sld[shell]; 34 60 35 int i,j,k; 36 int n_s; 61 // uniform shell; r=0 => r^3=0 => f=0, so works for core as well. 62 f -= M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r); 63 r += thickness[shell]; 64 f += M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r); 37 65 38 double sld_i,sld_f,dz,bes,fun,f,vol,qr,r,contr,f2; 39 double sign,slope=0.0; 40 double pi; 66 // iterate over sub_shells in the interface 67 const double dr = interface[shell]/n_steps; 68 const double delta = (shell==n_shells-1 ? sld_solvent : sld[shell+1]) - sld_l; 69 const double nu_shell = fmax(fabs(nu[shell]), 1.e-14); 70 const int shape_shell = (int)(shape[shell]); 41 71 42 double total_thick=0.0; 72 // if there is no interface the equations don't work 73 if (dr == 0.) continue; 43 74 44 //TODO: This part can be further simplified 45 int fun_type[12]; 46 double sld[12]; 47 double thick_internal[12]; 48 double thick[12]; 49 double fun_coef[12]; 75 double sld_in = sld_l; 76 for (int step=1; step <= n_steps; step++) { 77 // find sld_i at the outer boundary of sub-shell step 78 //const double z = (double)step/(double)n_steps; 79 const double z = (double)step/(double)n_steps; 80 const double fraction = blend(shape_shell, nu_shell, z); 81 const double sld_out = fraction*delta + sld_l; 82 // calculate slope 83 const double slope = (sld_out - sld_in)/dr; 84 const double contrast = sld_in - slope*r; 50 85 51 fun_type[0] = func_inter_core; 52 fun_coef[0] = fabs(nu_inter_core); 53 sld[0] = sld_core; 54 thick[0] = radius_core; 55 thick_internal[0] = thick_inter_core; 56 57 for (i =1; i<=n_shells; i++){ 58 sld[i] = sld_flat[i-1]; 59 thick_internal[i]= thick_inter[i-1]; 60 thick[i] = thick_flat[i-1]; 61 fun_type[i] = func_inter[i-1]; 62 fun_coef[i] = fabs(nu_inter[i-1]); 63 total_thick += thick[i]; 64 total_thick += thick_internal[i]; //doesn't account for core layer 65 } 66 67 sld[n_shells+1] = sld_solvent; 68 thick[n_shells+1] = total_thick/5.0; 69 thick_internal[n_shells+1] = 0.0; 70 fun_coef[n_shells+1] = 0.0; 71 fun_type[n_shells+1] = 0; 72 73 pi = 4.0*atan(1.0); 74 f = 0.0; 75 r = 0.0; 76 vol = 0.0; 77 sld_f = sld_core; 78 79 dz = 0.0; 80 // iteration for # of shells + core + solvent 81 for (i=0;i<=n_shells+1; i++){ 82 // iteration for flat and interface 83 for (j=0;j<2;j++){ 84 // iteration for sub_shells in the interface 85 // starts from #1 sub-layer 86 for (n_s=1;n_s<=npts_inter; n_s++){ 87 // for solvent, it doesn't have an interface 88 if (i==n_shells+1 && j==1) 89 break; 90 // for flat layers 91 if (j==0){ 92 dz = thick[i]; 93 sld_i = sld[i]; 94 slope = 0.0; 95 } 96 // for interfacial sub_shells 97 else{ 98 dz = thick_internal[i]/npts_inter; 99 // find sld_i at the outer boundary of sub-layer #n_s 100 sld_i = intersldfunc(fun_type[i], npts_inter, n_s, 101 fun_coef[i], sld[i], sld[i+1]); 102 // calculate slope 103 slope= (sld_i -sld_f)/dz; 104 } 105 contr = sld_f-slope*r; 106 // iteration for the left and right boundary of the shells 107 for (k=0; k<2; k++){ 108 // At r=0, the contribution to I is zero, so skip it. 109 if ( i == 0 && j == 0 && k == 0){ 110 continue; 111 } 112 // On the top of slovent there is no interface; skip it. 113 if (i == n_shells+1 && k == 1){ 114 continue; 115 } 116 // At the right side (outer) boundary 117 if ( k == 1){ 118 sign = 1.0; 119 r += dz; 120 } 121 // At the left side (inner) boundary 122 else{ 123 sign = -1.0; 124 } 125 126 qr = q * r; 127 fun = 0.0; 128 129 if(qr == 0.0){ 130 bes = sign * 1.0; 131 } 132 else{ 133 // for flat sub-layer 134 //TODO: Single precision calculation fails here 135 bes = sign * sph_j1c(qr); 136 if (fabs(slope) > 0.0 ){ 137 const double qrsq = qr*qr; 138 double sinqr, cosqr; 139 SINCOS(qr, sinqr, cosqr); 140 fun = sign * 3.0 * r * 141 (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq * qrsq); 142 // In the onion model Jae-He's formula is rederived 143 // and gives following: 144 //fun = 3.0 * sign * r * 145 //(2.0*cosqr + qr*sinqr)/(qrsq*qrsq); 146 //But this seems not to be working in this case... 147 } 148 } 149 150 // update total volume 151 vol = M_4PI_3 * cube(r); 152 f += vol * (bes * contr + fun * slope); 153 } 154 sld_f = sld_i; 155 // no sub-layer iteration (n_s loop) for the flat layer 156 if (j==0) 157 break; 158 } 86 // iteration for the left and right boundary of the shells 87 f -= f_linear(q, r, contrast, slope); 88 r += dr; 89 f += f_linear(q, r, contrast, slope); 90 sld_in = sld_out; 159 91 } 160 92 } 93 // add in solvent effect 94 f -= M_4PI_3 * cube(r) * sld_solvent * sph_j1c(q*r); 161 95 162 f2 = f * f * 1.0e-4;163 return (f2);96 const double f2 = f * f * 1.0e-4; 97 return f2; 164 98 } 165 99 166 167 /**168 * Function to evaluate 1D SphereSLD function169 * @param q: q-value170 * @return: function value171 */172 static double Iq(double q,173 int n_shells,174 int npts_inter,175 double radius_core,176 double sld_core,177 double sld_solvent,178 double func_inter0,179 double thick_inter0,180 double nu_inter0,181 double sld_flat[],182 double thick_flat[],183 double func_inter[],184 double thick_inter[],185 double nu_inter[] ) {186 187 double intensity = sphere_sld_kernel(q, n_shells, npts_inter, radius_core,188 sld_core, sld_solvent, func_inter0, thick_inter0, nu_inter0,189 sld_flat, thick_flat, func_inter, thick_inter, nu_inter);190 191 return intensity;192 }193 194 /**195 * Function to evaluate 2D SphereSLD function196 * @param q_x: value of Q along x197 * @param q_y: value of Q along y198 * @return: function value199 */200 201 /*static double Iqxy(double qx, double qy,202 int n_shells,203 int npts_inter,204 double radius_core205 double sld_core,206 double sld_solvent,207 double sld_flat[],208 double thick_flat[],209 double func_inter[],210 double thick_inter[],211 double nu_inter[],212 ) {213 214 double q = sqrt(qx*qx + qy*qy);215 return Iq(q, n_shells, npts_inter, radius_core, sld_core, sld_solvent,216 sld_flat[], thick_flat[], func_inter[], thick_inter[], nu_inter[])217 }*/218 -
TabularUnified sasmodels/models/spherical_sld.py ¶
r4fd2c63 r54bcd4a 190 190 191 191 import numpy as np 192 from numpy import inf 192 from numpy import inf, expm1, sqrt 193 from scipy.special import erf 193 194 194 195 name = "spherical_sld" … … 200 201 category = "shape:sphere" 201 202 203 SHAPES = ["erf(|nu|*z)", "Rpow(z^|nu|)", "Lpow(z^|nu|)", 204 "Rexp(-|nu|z)", "Lexp(-|nu|z)"], 205 202 206 # pylint: disable=bad-whitespace, line-too-long 203 207 # ["name", "units", default, [lower, upper], "type", "description"], 204 parameters = [["n_shells", "", 1, [0, 10], "volume", "number of shells"], 205 ["npts_inter", "", 35, [0, inf], "", "number of points in each sublayer Must be odd number"], 206 ["radius_core", "Ang", 50.0, [0, inf], "volume", "intern layer thickness"], 207 ["sld_core", "1e-6/Ang^2", 2.07, [-inf, inf], "sld", "sld function flat"], 208 ["sld_solvent", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "sld function solvent"], 209 ["func_inter0", "", 0, [0, 5], "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4, Linear:5"], 210 ["thick_inter0", "Ang", 50.0, [0, inf], "volume", "intern layer thickness for core layer"], 211 ["nu_inter0", "", 2.5, [-inf, inf], "", "steepness parameter for core layer"], 212 ["sld_flat[n_shells]", "1e-6/Ang^2", 4.06, [-inf, inf], "sld", "sld function flat"], 213 ["thick_flat[n_shells]", "Ang", 100.0, [0, inf], "volume", "flat layer_thickness"], 214 ["func_inter[n_shells]", "", 0, [0, 5], "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4, Linear:5"], 215 ["thick_inter[n_shells]", "Ang", 50.0, [0, inf], "volume", "intern layer thickness"], 216 ["nu_inter[n_shells]", "", 2.5, [-inf, inf], "", "steepness parameter"], 208 parameters = [["n_shells", "", 1, [1, 11], "volume", "number of shells"], 209 ["sld_solvent", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "solvent sld"], 210 ["sld[n_shells]", "1e-6/Ang^2", 4.06, [-inf, inf], "sld", "sld of the shell"], 211 ["thickness[n_shells]", "Ang", 100.0, [0, inf], "volume", "thickness shell"], 212 ["interface[n_shells]", "Ang", 50.0, [0, inf], "volume", "thickness of the interface"], 213 ["shape[n_shells]", "", 0, SHAPES, "", "interface shape"], 214 ["nu[n_shells]", "", 2.5, [0, inf], "", "interface shape exponent"], 215 ["n_steps", "", 35, [0, inf], "", "number of steps in each interface (must be an odd integer)"], 217 216 ] 218 217 # pylint: enable=bad-whitespace, line-too-long 219 source = ["lib/polevl.c", "lib/sas_erf.c", "lib/ librefl.c", "lib/sph_j1c.c", "spherical_sld.c"]218 source = ["lib/polevl.c", "lib/sas_erf.c", "lib/sph_j1c.c", "spherical_sld.c"] 220 219 single = False # TODO: fix low q behaviour 221 220 222 221 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 223 def profile(n_shells, radius_core, sld_core, sld_solvent, sld_flat, 224 thick_flat, func_inter, thick_inter, nu_inter, npts_inter): 222 223 SHAPE_FUNCTIONS = [ 224 lambda z, nu: erf(nu/sqrt(2)*(2*z-1))/(2*erf(nu/sqrt(2))) + 0.5, # erf 225 lambda z, nu: z**nu, # Rpow 226 lambda z, nu: 1 - (1-z)**nu, # Lpow 227 lambda z, nu: expm1(-nu*z)/expm1(-nu), # Rexp 228 lambda z, nu: expm1(nu*z)/expm1(nu), # Lexp 229 ] 230 231 def profile(n_shells, sld_solvent, sld, thickness, 232 interface, shape, nu, n_steps): 225 233 """ 226 234 Returns shape profile with x=radius, y=SLD. … … 228 236 229 237 z = [] 230 beta= []238 rho = [] 231 239 z0 = 0 232 240 # two sld points for core 233 241 z.append(0) 234 beta.append(sld_core) 235 z.append(radius_core) 236 beta.append(sld_core) 237 z0 += radius_core 238 239 for i in range(1, n_shells+2): 240 dz = thick_inter[i-1]/npts_inter 241 # j=0 for interface, j=1 for flat layer 242 for j in range(0, 2): 243 # interation for sub-layers 244 for n_s in range(0, npts_inter+1): 245 if j == 1: 246 if i == n_shells+1: 247 break 248 # shift half sub thickness for the first point 249 z0 -= dz#/2.0 250 z.append(z0) 251 #z0 -= dz/2.0 252 z0 += thick_flat[i] 253 sld_i = sld_flat[i] 254 beta.append(sld_flat[i]) 255 dz = 0 256 else: 257 nu = nu_inter[i-1] 258 # decide which sld is which, sld_r or sld_l 259 if i == 1: 260 sld_l = sld_core 261 else: 262 sld_l = sld_flat[i-1] 263 if i == n_shells+1: 264 sld_r = sld_solvent 265 else: 266 sld_r = sld_flat[i] 267 # get function type 268 func_idx = func_inter[i-1] 269 # calculate the sld 270 sld_i = intersldfunc(func_idx, npts_inter, n_s, nu, 271 sld_l, sld_r) 272 # append to the list 273 z.append(z0) 274 beta.append(sld_i) 275 z0 += dz 276 if j == 1: 277 break 278 z.append(z0) 279 beta.append(sld_solvent) 280 z_ext = z0/5.0 281 z.append(z0+z_ext) 282 beta.append(sld_solvent) 242 rho.append(sld[0]) 243 244 for i in range(0, n_shells): 245 z0 += thickness[i] 246 z.append(z0) 247 rho.append(sld[i]) 248 dz = interface[i]/n_steps 249 sld_l = sld[i] 250 sld_r = sld[i+1] if i < n_shells-1 else sld_solvent 251 interface = SHAPE_FUNCTIONS[int(np.clip(shape[i], 0, len(SHAPES)-1))] 252 for step in range(1, n_steps+1): 253 portion = interface(float(step)/n_steps, max(abs(nu[i]), 1e-14)) 254 z0 += dz 255 z.append(z0) 256 rho.append((sld_r - sld_l)*portion + sld_l) 257 z.append(z0*1.2) 258 rho.append(sld_solvent) 283 259 # return sld profile (r, beta) 284 return np.asarray(z), np.asarray(beta)*1e-6 285 286 def ER(n_shells, radius_core, thick_inter0, thick_inter, thick_flat): 260 return np.asarray(z), np.asarray(rho)*1e-6 261 262 263 def ER(n_shells, thickness, interface): 287 264 n_shells = int(n_shells) 288 total_thickness = thick_inter0 289 total_thickness += np.sum(thick_inter[:n_shells], axis=0) 290 total_thickness += np.sum(thick_flat[:n_shells], axis=0) 291 return total_thickness + radius_core 265 total = (np.sum(thickness[:n_shells], axis=1) 266 + np.sum(interface[:n_shells], axis=1)) 267 return total 292 268 293 269 294 270 demo = { 295 "n_shells": 4, 296 "npts_inter": 35.0, 297 "radius_core": 50.0, 298 "sld_core": 2.07, 271 "n_shells": 5, 272 "n_steps": 35.0, 299 273 "sld_solvent": 1.0, 300 "thick_inter0": 50.0, 301 "func_inter0": 0, 302 "nu_inter0": 2.5, 303 "sld_flat":[4.0,3.5,4.0,3.5], 304 "thick_flat":[100.0,100.0,100.0,100.0], 305 "func_inter":[0,0,0,0], 306 "thick_inter":[50.0,50.0,50.0,50.0], 307 "nu_inter":[2.5,2.5,2.5,2.5], 274 "sld":[2.07,4.0,3.5,4.0,3.5], 275 "thickness":[50.0,100.0,100.0,100.0,100.0], 276 "interface":[50.0,50.0,50.0,50.0], 277 "shape": [0,0,0,0,0], 278 "nu":[2.5,2.5,2.5,2.5,2.5], 308 279 } 309 280 … … 312 283 tests = [ 313 284 # Accuracy tests based on content in test/utest_extra_models.py 314 [{"n_shells":4, 315 'npts_inter':35, 316 "radius_core":50.0, 317 "sld_core":2.07, 285 [{"n_shells": 5, 286 "n_steps": 35, 318 287 "sld_solvent": 1.0, 319 "sld _flat":[4.0,3.5,4.0,3.5],320 "thick _flat":[100.0,100.0,100.0,100.0],321 " func_inter":[0,0,0,0],322 " thick_inter":[50.0,50.0,50.0,50.0],323 "nu _inter":[2.5,2.5,2.5,2.5]288 "sld": [2.07, 4.0, 3.5, 4.0, 3.5], 289 "thickness": [50.0, 100.0, 100.0, 100.0, 100.0], 290 "interface": [50]*5, 291 "shape": [0]*5, 292 "nu": [2.5]*5, 324 293 }, 0.001, 0.001], 325 294 ]
Note: See TracChangeset
for help on using the changeset viewer.