Changeset 0ce5710 in sasmodels for sasmodels/models
- Timestamp:
- Apr 14, 2016 8:32:31 AM (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:
- b151003
- Parents:
- a5b8477 (diff), 62cf915 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Location:
- sasmodels/models
- Files:
-
- 23 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/adsorbed_layer.py
rec45c4f r62cf915 3 3 #converted by Steve King, Mar 2016 4 4 5 r""" 6 This model describes the scattering from a layer of surfactant or polymer 7 adsorbed on large, smooth, notionally spherical particles under the conditions 8 that (i) the particles (cores) are contrast-matched to the dispersion medium, 9 (ii) $S(Q) \sim 1$ (ie, the particle volume fraction is dilute), (iii) the 10 particle radius is >> layer thickness (ie, the interface is locally flat), 11 and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is 12 absent or has been corrected for. 5 13 6 7 r""" 8 This model describes the scattering from a layer of surfactant or polymer adsorbed on large, smooth, notionally spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for. 9 10 Unlike many other core-shell models, this model does not assume any form for the density distribution of the adsorbed species normal to the interface (cf, a core-shell model normally assumes the density distribution to be a homogeneous step-function). For comparison, if the thickness of a (traditional core-shell like) step function distribution is *t*, the second moment about the mean of the density distribution (ie, the distance of the centre-of-mass of the distribution from the interface), |sigma| = sqrt((*t* :sup:`2` )/12). 14 Unlike many other core-shell models, this model does not assume any form 15 for the density distribution of the adsorbed species normal to the interface 16 (cf, a core-shell model normally assumes the density distribution to be a 17 homogeneous step-function). For comparison, if the thickness of a (traditional 18 core-shell like) step function distribution is $t$, the second moment about 19 the mean of the density distribution (ie, the distance of the centre-of-mass 20 of the distribution from the interface), $\sigma = \sqrt{t^2/12}$. 11 21 12 22 Definition … … 15 25 .. math:: 16 26 17 I(q) = \text{scale} \cdot(\rho_\text{poly}-\rho_\text{solvent})^2 \left[\frac{6\pi\phi_\text{core}}{Q^2}\frac{\Gamma^2}{\delta_\text{poly}^2R_\text{core}} \exp(-Q^2\sigma^2)\right] + \text{background} 27 I(q) = \text{scale} \cdot (\rho_\text{poly}-\rho_\text{solvent})^2 28 \left[ 29 \frac{6\pi\phi_\text{core}}{Q^2} 30 \frac{\Gamma^2}{\delta_\text{poly}^2R_\text{core}} 31 \exp(-Q^2\sigma^2) 32 \right] + \text{background} 18 33 19 where *scale* is a scale factor, |rho|\ :sub:`poly` is the sld of the polymer (or surfactant) layer, |rho|\ :sub:`solv` is the sld of the solvent/medium and cores, |phi|\ :sub:`core` is the volume fraction of the core particles, |delta|\ :sub:`poly` is the bulk density of the polymer, |biggamma| is the adsorbed amount, and |sigma| is the second moment of the thickness distribution. 34 where *scale* is a scale factor, $\rho_\text{poly}$ is the sld of the 35 polymer (or surfactant) layer, $\rho_\text{solv}$ is the sld of the 36 solvent/medium and cores, $\phi_\text{core}$ is the volume fraction of 37 the core particles, $\delta_\text{poly}$ is the bulk density of the 38 polymer, $\Gamma$ is the adsorbed amount, and $\sigma$ is the second 39 moment of the thickness distribution. 20 40 21 Note that all parameters except the |sigma| are correlated so fitting more than one of these parameters will generally fail. Also note that unlike other shape models, no volume normalization is applied to this model (the calculation is exact). 41 Note that all parameters except $\sigma$ are correlated so fitting more 42 than one of these parameters will generally fail. Also note that unlike 43 other shape models, no volume normalization is applied to this model (the 44 calculation is exact). 22 45 23 46 References 24 47 ---------- 25 48 26 S King, P Griffiths, J Hone, and T Cosgrove, *SANS from Adsorbed Polymer Layers*,27 * Macromol. Symp.*, 190 (2002) 33-42.49 S King, P Griffiths, J Hone, and T Cosgrove, 50 *SANS from Adsorbed Polymer Layers*, *Macromol. Symp.*, 190 (2002) 33-42. 28 51 """ 29 52 30 53 from numpy import inf, sqrt, pi, exp 31 54 32 name = 33 title = 55 name = "adsorbed_layer" 56 title = "Scattering from an adsorbed layer on particles" 34 57 35 description = 58 description = """ 36 59 Evaluates the scattering from large particles 37 60 with an adsorbed layer of surfactant or … … 39 62 density distribution. 40 63 """ 41 category = 64 category = "shape:sphere" 42 65 66 # pylint: disable=bad-whitespace, line-too-long 43 67 # ["name", "units", default, [lower, upper], "type", "description"], 44 parameters = [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment of polymer distribution"], 45 ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount of polymer"], 46 ["density_shell", "g/cm3", 0.7, [0.0, inf], "", "Bulk density of polymer in the shell"], 47 ["radius", "Ang", 500.0, [0.0, inf], "", "Core particle radius"], 48 ["volfraction", "None", 0.14, [0.0, inf], "", "Core particle volume fraction"], 49 ["sld_shell", "1e-6/Ang^2", 1.5, [-inf, inf], "sld", "Polymer shell SLD"], 50 ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld", "Solvent SLD"]] 68 parameters = [ 69 ["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment of polymer distribution"], 70 ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount of polymer"], 71 ["density_shell", "g/cm3", 0.7, [0.0, inf], "", "Bulk density of polymer in the shell"], 72 ["radius", "Ang", 500.0, [0.0, inf], "", "Core particle radius"], 73 ["volfraction", "None", 0.14, [0.0, inf], "", "Core particle volume fraction"], 74 ["sld_shell", "1e-6/Ang^2", 1.5, [-inf, inf], "sld", "Polymer shell SLD"], 75 ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld", "Solvent SLD"], 76 ] 77 # pylint: disable=bad-whitespace, line-too-long 51 78 52 79 # NB: Scale and Background are implicit parameters on every model 53 def Iq(q, second_moment, adsorbed_amount, density_shell, radius, 54 80 def Iq(q, second_moment, adsorbed_amount, density_shell, radius, 81 volfraction, sld_shell, sld_solvent): 55 82 # pylint: disable = missing-docstring 56 #deltarhosqrd = (sld_shell - sld_solvent) * (sld_shell - sld_solvent)57 #numerator = 6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount)58 #denominator = (q * q) * (density_shell * density_shell) * radius59 #eterm = exp(-1.0 * (q * q) * (second_moment * second_moment))60 ##scale by 10^-2 for units conversion to cm^-161 #inten = 1.0e-02 * deltarhosqrd * ((numerator / denominator) * eterm)62 aa = (sld_shell - sld_solvent) * adsorbed_amount / q / density_shell83 #deltarhosqrd = (sld_shell - sld_solvent) * (sld_shell - sld_solvent) 84 #numerator = 6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount) 85 #denominator = (q * q) * (density_shell * density_shell) * radius 86 #eterm = exp(-1.0 * (q * q) * (second_moment * second_moment)) 87 ##scale by 10^-2 for units conversion to cm^-1 88 #inten = 1.0e-02 * deltarhosqrd * ((numerator / denominator) * eterm) 89 aa = (sld_shell - sld_solvent) * adsorbed_amount / q / density_shell 63 90 bb = q * second_moment 64 91 #scale by 10^-2 for units conversion to cm^-1 65 inten = 92 inten = 6.0e-02 * pi * volfraction * aa * aa * exp(-bb * bb) / radius 66 93 return inten 67 94 Iq.vectorized = True # Iq accepts an array of q values … … 70 97 # pylint: disable = missing-docstring 71 98 return Iq(sqrt(qx ** 2 + qy ** 2), *args) 72 Iqxy.vectorized = 99 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values 73 100 74 demo = dict(scale = 1.0, 75 second_moment = 23.0, 76 adsorbed_amount = 1.9, 77 density_shell = 0.7, 78 radius = 500.0, 79 volfraction = 0.14, 80 sld_shell = 1.5, 81 sld_solvent = 6.3, 82 background = 0.0) 101 # unit test values taken from SasView 3.1.2 102 tests = [ 103 [{'scale': 1.0, 'second_moment': 23.0, 'adsorbed_amount': 1.9, 104 'density_shell': 0.7, 'radius': 500.0, 'volfraction': 0.14, 105 'sld_shell': 1.5, 'sld_solvent': 6.3, 'background': 0.0}, 106 [0.0106939, 0.1], [73.741, 4.51684e-3]], 107 ] 83 108 84 # these unit test values taken from SasView 3.1.2 85 tests = [ 86 [{'scale': 1.0, 'second_moment': 23.0, 'adsorbed_amount': 1.9, 87 'density_shell': 0.7, 'radius': 500.0, 'volfraction': 0.14, 88 'sld_shell': 1.5, 'sld_solvent': 6.3, 'background': 0.0}, 89 [0.0106939, 0.1], [73.741, 4.51684e-3]], 90 ] 91 # ADDED by: SMK ON: 16Mar2016 convert from sasview, check vs SANDRA, 18Mar2016 RKH some edits & renaming 109 # 2016-03-16 SMK converted from sasview, checked vs SANDRA 110 # 2016-03-18 RKH some edits & renaming 111 # 2016-04-14 PAK reformatting -
sasmodels/models/cylinder.c
r26141cb re9b1663d 3 3 double Iqxy(double qx, double qy, double sld, double solvent_sld, 4 4 double radius, double length, double theta, double phi); 5 6 #define INVALID(v) (v.radius<0 || v.length<0) 5 7 6 8 double form_volume(double radius, double length) … … 15 17 double length) 16 18 { 17 // TODO: return NaN if radius<0 or length<0?18 19 // precompute qr and qh to save time in the loop 19 20 const double qr = q*radius; … … 47 48 double phi) 48 49 { 49 // TODO: return NaN if radius<0 or length<0?50 50 double sn, cn; // slots to hold sincos function output 51 51 -
sasmodels/models/cylinder.py
rf247314 r7ae2b7f 82 82 """ 83 83 84 import numpy as np 85 from numpy import pi, inf 84 import numpy as np # type: ignore 85 from numpy import pi, inf # type: ignore 86 86 87 87 name = "cylinder" -
sasmodels/models/flexible_cylinder.c
re6408d0 r4937980 1 double form_volume(double length, double kuhn_length, double radius); 2 double Iq(double q, double length, double kuhn_length, double radius, 3 double sld, double solvent_sld); 4 double Iqxy(double qx, double qy, double length, double kuhn_length, 5 double radius, double sld, double solvent_sld); 6 double flexible_cylinder_kernel(double q, double length, double kuhn_length, 7 double radius, double sld, double solvent_sld); 8 9 10 double form_volume(double length, double kuhn_length, double radius) 1 static double 2 form_volume(length, kuhn_length, radius) 11 3 { 12 4 return 1.0; 13 5 } 14 6 15 double flexible_cylinder_kernel(double q, 16 double length, 17 double kuhn_length, 18 double radius, 19 double sld, 20 double solvent_sld) 7 static double 8 Iq(double q, 9 double length, 10 double kuhn_length, 11 double radius, 12 double sld, 13 double solvent_sld) 21 14 { 22 23 const double cont = sld-solvent_sld; 24 const double qr = q*radius; 25 //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 26 const double crossSect = sas_J1c(qr); 27 double flex = Sk_WR(q,length,kuhn_length); 28 flex *= crossSect*crossSect; 29 flex *= M_PI*radius*radius*length; 30 flex *= cont*cont; 31 flex *= 1.0e-4; 32 return flex; 15 const double contrast = sld - solvent_sld; 16 const double cross_section = sas_J1c(q*radius); 17 const double volume = M_PI*radius*radius*length; 18 const double flex = Sk_WR(q, length, kuhn_length); 19 return 1.0e-4 * volume * square(contrast*cross_section) * flex; 33 20 } 34 35 double Iq(double q,36 double length,37 double kuhn_length,38 double radius,39 double sld,40 double solvent_sld)41 {42 43 double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);44 return result;45 }46 47 double Iqxy(double qx, double qy,48 double length,49 double kuhn_length,50 double radius,51 double sld,52 double solvent_sld)53 {54 double q;55 q = sqrt(qx*qx+qy*qy);56 double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);57 58 return result;59 } -
sasmodels/models/gel_fit.c
r30b4ddf r03cac08 1 double form_volume(void); 2 3 double Iq(double q, 4 double guinier_scale, 5 double lorentzian_scale, 6 double gyration_radius, 7 double fractal_exp, 8 double cor_length); 9 10 double Iqxy(double qx, double qy, 11 double guinier_scale, 12 double lorentzian_scale, 13 double gyration_radius, 14 double fractal_exp, 15 double cor_length); 16 17 static double _gel_fit_kernel(double q, 1 static double Iq(double q, 18 2 double guinier_scale, 19 3 double lorentzian_scale, … … 24 8 // Lorentzian Term 25 9 ////////////////////////double a(x[i]*x[i]*zeta*zeta); 26 double lorentzian_term = q*q*cor_length*cor_length;10 double lorentzian_term = square(q*cor_length); 27 11 lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term; 28 12 lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) ); … … 30 14 // Exponential Term 31 15 ////////////////////////double d(x[i]*x[i]*rg*rg); 32 double exp_term = q*q*gyration_radius*gyration_radius;16 double exp_term = square(q*gyration_radius); 33 17 exp_term = exp(-1.0*(exp_term/3.0)); 34 18 … … 37 21 return result; 38 22 } 39 double form_volume(void){40 // Unused, so free to return garbage.41 return NAN;42 }43 44 double Iq(double q,45 double guinier_scale,46 double lorentzian_scale,47 double gyration_radius,48 double fractal_exp,49 double cor_length)50 {51 return _gel_fit_kernel(q,52 guinier_scale,53 lorentzian_scale,54 gyration_radius,55 fractal_exp,56 cor_length);57 }58 59 // Iqxy is never called since no orientation or magnetic parameters.60 double Iqxy(double qx, double qy,61 double guinier_scale,62 double lorentzian_scale,63 double gyration_radius,64 double fractal_exp,65 double cor_length)66 {67 double q = sqrt(qx*qx + qy*qy);68 return _gel_fit_kernel(q,69 guinier_scale,70 lorentzian_scale,71 gyration_radius,72 fractal_exp,73 cor_length);74 }75 -
sasmodels/models/hardsphere.py
rec45c4f rd2bb604 149 149 """ 150 150 151 Iqxy = """152 // never called since no orientation or magnetic parameters.153 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);154 """155 156 151 # ER defaults to 0.0 157 152 # VR defaults to 1.0 -
sasmodels/models/hayter_msa.py
rec45c4f rd2bb604 87 87 return 1.0; 88 88 """ 89 Iqxy = """90 // never called since no orientation or magnetic parameters.91 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);92 """93 89 # ER defaults to 0.0 94 90 # VR defaults to 1.0 -
sasmodels/models/lamellar.py
rec45c4f rd2bb604 82 82 """ 83 83 84 Iqxy = """85 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);86 """87 88 84 # ER defaults to 0.0 89 85 # VR defaults to 1.0 -
sasmodels/models/lamellar_hg.py
rec45c4f rd2bb604 101 101 """ 102 102 103 Iqxy = """104 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);105 """106 107 103 # ER defaults to 0.0 108 104 # VR defaults to 1.0 -
sasmodels/models/lamellar_hg_stack_caille.py
rec45c4f rd2bb604 120 120 """ 121 121 122 Iqxy = """123 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);124 """125 126 122 # ER defaults to 0.0 127 123 # VR defaults to 1.0 -
sasmodels/models/lamellar_stack_caille.py
rec45c4f rd2bb604 104 104 """ 105 105 106 Iqxy = """107 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);108 """109 110 106 # ER defaults to 0.0 111 107 # VR defaults to 1.0 -
sasmodels/models/lamellar_stack_paracrystal.py
rec45c4f rd2bb604 132 132 """ 133 133 134 Iqxy = """135 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);136 """137 138 134 # ER defaults to 0.0 139 135 # VR defaults to 1.0 -
sasmodels/models/lib/sas_JN.c
re6408d0 r4937980 48 48 */ 49 49 50 static double 51 sas_JN( int n, double x ) { 50 double sas_JN( int n, double x ); 51 52 double sas_JN( int n, double x ) { 52 53 53 54 const double MACHEP = 1.11022302462515654042E-16; -
sasmodels/models/lib/sph_j1c.c
re6f1410 rba32cdd 7 7 * using double precision that are the source. 8 8 */ 9 double sph_j1c(double q); 9 10 10 11 // The choice of the number of terms in the series and the cutoff value for … … 43 44 #endif 44 45 45 double sph_j1c(double q);46 46 double sph_j1c(double q) 47 47 { -
sasmodels/models/lib/sphere_form.c
rad90df9 rba32cdd 1 inline double 2 sphere_volume(double radius) 1 double sphere_volume(double radius); 2 double sphere_form(double q, double radius, double sld, double solvent_sld); 3 4 double sphere_volume(double radius) 3 5 { 4 6 return M_4PI_3*cube(radius); 5 7 } 6 8 7 inline double 8 sphere_form(double q, double radius, double sld, double solvent_sld) 9 double sphere_form(double q, double radius, double sld, double solvent_sld) 9 10 { 10 11 const double fq = sphere_volume(radius) * sph_j1c(q*radius); -
sasmodels/models/lib/wrc_cyl.c
re7678b2 rba32cdd 2 2 Functions for WRC implementation of flexible cylinders 3 3 */ 4 double Sk_WR(double q, double L, double b); 5 6 4 7 static double 5 8 AlphaSquare(double x) … … 363 366 } 364 367 365 double Sk_WR(double q, double L, double b);366 368 double Sk_WR(double q, double L, double b) 367 369 { -
sasmodels/models/onion.c
rfdb1487 rce896fd 4 4 double thickness, double A) 5 5 { 6 const double vol = 4.0/3.0 * M_PI * r * r * r;6 const double vol = M_4PI_3 * cube(r); 7 7 const double qr = q * r; 8 8 const double alpha = A * r/thickness; … … 19 19 double sinqr, cosqr; 20 20 SINCOS(qr, sinqr, cosqr); 21 fun = -3.0 (21 fun = -3.0*( 22 22 ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 23 23 - (alpha*sinqr/qr - cosqr) … … 32 32 f_linear(double q, double r, double sld, double slope) 33 33 { 34 const double vol = 4.0/3.0 * M_PI * r * r * r;34 const double vol = M_4PI_3 * cube(r); 35 35 const double qr = q * r; 36 36 const double bes = sph_j1c(qr); … … 52 52 { 53 53 const double bes = sph_j1c(q * r); 54 const double vol = 4.0/3.0 * M_PI * r * r * r;54 const double vol = M_4PI_3 * cube(r); 55 55 return sld * vol * bes; 56 56 } … … 64 64 r += thickness[i]; 65 65 } 66 return 4.0/3.0 * M_PI * r * r * r;66 return M_4PI_3*cube(r); 67 67 } 68 68 69 69 static double 70 Iq(double q, double core_sld, double core_radius, double solvent_sld,71 double n, double in_sld[], double out_sld[], double thickness[],70 Iq(double q, double sld_core, double core_radius, double sld_solvent, 71 double n, double sld_in[], double sld_out[], double thickness[], 72 72 double A[]) 73 73 { 74 74 int i; 75 r = core_radius;76 f = f_constant(q, r, core_sld);75 double r = core_radius; 76 double f = f_constant(q, r, sld_core); 77 77 for (i=0; i<n; i++){ 78 78 const double r0 = r; … … 92 92 } 93 93 } 94 f -= f_constant(q, r, s olvent_sld);95 f2 = f * f * 1.0e-4;94 f -= f_constant(q, r, sld_solvent); 95 const double f2 = f * f * 1.0e-4; 96 96 97 97 return f2; -
sasmodels/models/onion.py
r416609b rfa5fd8d 293 293 294 294 # ["name", "units", default, [lower, upper], "type","description"], 295 parameters = [[" core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",295 parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 296 296 "Core scattering length density"], 297 297 ["core_radius", "Ang", 200., [0, inf], "volume", 298 298 "Radius of the core"], 299 ["s olvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "",299 ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", 300 300 "Solvent scattering length density"], 301 ["n ", "", 1, [0, 10], "volume",301 ["n_shells", "", 1, [0, 10], "volume", 302 302 "number of shells"], 303 [" in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",303 ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 304 304 "scattering length density at the inner radius of shell k"], 305 [" out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "",305 ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 306 306 "scattering length density at the outer radius of shell k"], 307 ["thickness[n ]", "Ang", 40., [0, inf], "volume",307 ["thickness[n_shells]", "Ang", 40., [0, inf], "volume", 308 308 "Thickness of shell k"], 309 ["A[n ]", "", 1.0, [-inf, inf], "",309 ["A[n_shells]", "", 1.0, [-inf, inf], "", 310 310 "Decay rate of shell k"], 311 311 ] 312 312 313 #source = ["lib/sph_j1c.c", "onion.c"] 314 315 def Iq(q, *args, **kw): 316 return q 317 318 def Iqxy(qx, *args, **kw): 319 return qx 320 321 322 def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 313 source = ["lib/sph_j1c.c", "onion.c"] 314 315 #def Iq(q, *args, **kw): 316 # return q 317 318 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 319 def profile(core_sld, core_radius, solvent_sld, n_shells, 320 in_sld, out_sld, thickness, A): 323 321 """ 324 SLD profile322 Returns shape profile with x=radius, y=SLD. 325 323 """ 326 324 327 total_radius = 1.25*(sum(thickness[:n ]) + core_radius + 1)325 total_radius = 1.25*(sum(thickness[:n_shells]) + core_radius + 1) 328 326 dr = total_radius/400 # 400 points for a smooth plot 329 327 … … 338 336 339 337 # add in the shells 340 for k in range(n ):338 for k in range(n_shells): 341 339 # Left side of each shells 342 340 r0 = r[-1] … … 365 363 beta.append(solvent_sld) 366 364 367 return np.asarray(r), np.asarray(beta) 365 return np.asarray(r), np.asarray(beta)*1e-6 368 366 369 367 def ER(core_radius, n, thickness): … … 374 372 375 373 demo = { 376 "s olvent_sld": 2.2,377 " core_sld": 1.0,374 "sld_solvent": 2.2, 375 "sld_core": 1.0, 378 376 "core_radius": 100, 379 377 "n": 4, 380 " in_sld": [0.5, 1.5, 0.9, 2.0],381 " out_sld": [nan, 0.9, 1.2, 1.6],378 "sld_in": [0.5, 1.5, 0.9, 2.0], 379 "sld_out": [nan, 0.9, 1.2, 1.6], 382 380 "thickness": [50, 75, 150, 75], 383 381 "A": [0, -1, 1e-4, 1], -
sasmodels/models/rpa.c
r13ed84c rd2bb604 1 1 double Iq(double q, double case_num, 2 double Na, double Phia, double va, double a_sld, double ba, 3 double Nb, double Phib, double vb, double b_sld, double bb, 4 double Nc, double Phic, double vc, double c_sld, double bc, 5 double Nd, double Phid, double vd, double d_sld, double bd, 2 double N[], double Phi[], double v[], double L[], double b[], 6 3 double Kab, double Kac, double Kad, 7 4 double Kbc, double Kbd, double Kcd 8 5 ); 9 6 10 double Iqxy(double qx, double qy, double case_num,11 double Na, double Phia, double va, double a_sld, double ba,12 double Nb, double Phib, double vb, double b_sld, double bb,13 double Nc, double Phic, double vc, double c_sld, double bc,14 double Nd, double Phid, double vd, double d_sld, double bd,15 double Kab, double Kac, double Kad,16 double Kbc, double Kbd, double Kcd17 );18 19 double form_volume(void);20 21 double form_volume(void)22 {23 return 1.0;24 }25 26 7 double Iq(double q, double case_num, 27 double Na, double Phia, double va, double La, double ba, 28 double Nb, double Phib, double vb, double Lb, double bb, 29 double Nc, double Phic, double vc, double Lc, double bc, 30 double Nd, double Phid, double vd, double Ld, double bd, 8 double N[], double Phi[], double v[], double L[], double b[], 31 9 double Kab, double Kac, double Kad, 32 10 double Kbc, double Kbd, double Kcd … … 36 14 #if 0 // Sasview defaults 37 15 if (icase <= 1) { 38 N a=Nb=1000.0;39 Phi a=Phib=0.0000001;16 N[0]=N[1]=1000.0; 17 Phi[0]=Phi[1]=0.0000001; 40 18 Kab=Kac=Kad=Kbc=Kbd=-0.0004; 41 L a=Lb=1e-12;42 v a=vb=100.0;43 b a=bb=5.0;19 L[0]=L[1]=1e-12; 20 v[0]=v[1]=100.0; 21 b[0]=b[1]=5.0; 44 22 } else if (icase <= 4) { 45 Phi a=0.0000001;23 Phi[0]=0.0000001; 46 24 Kab=Kac=Kad=-0.0004; 47 L a=1e-12;48 v a=100.0;49 b a=5.0;25 L[0]=1e-12; 26 v[0]=100.0; 27 b[0]=5.0; 50 28 } 51 29 #else 52 30 if (icase <= 1) { 53 N a=Nb=0.0;54 Phi a=Phib=0.0;31 N[0]=N[1]=0.0; 32 Phi[0]=Phi[1]=0.0; 55 33 Kab=Kac=Kad=Kbc=Kbd=0.0; 56 L a=Lb=Ld;57 v a=vb=vd;58 b a=bb=0.0;34 L[0]=L[1]=L[3]; 35 v[0]=v[1]=v[3]; 36 b[0]=b[1]=0.0; 59 37 } else if (icase <= 4) { 60 N a= 0.0;61 Phi a=0.0;38 N[0] = 0.0; 39 Phi[0]=0.0; 62 40 Kab=Kac=Kad=0.0; 63 L a=Ld;64 v a=vd;65 b a=0.0;41 L[0]=L[3]; 42 v[0]=v[3]; 43 b[0]=0.0; 66 44 } 67 45 #endif 68 46 69 const double Xa = q*q*b a*ba*Na/6.0;70 const double Xb = q*q*b b*bb*Nb/6.0;71 const double Xc = q*q*b c*bc*Nc/6.0;72 const double Xd = q*q*b d*bd*Nd/6.0;47 const double Xa = q*q*b[0]*b[0]*N[0]/6.0; 48 const double Xb = q*q*b[1]*b[1]*N[1]/6.0; 49 const double Xc = q*q*b[2]*b[2]*N[2]/6.0; 50 const double Xd = q*q*b[3]*b[3]*N[3]/6.0; 73 51 74 52 // limit as Xa goes to 0 is 1 … … 98 76 #if 0 99 77 const double S0aa = icase<5 100 ? 1.0 : N a*Phia*va*Paa;78 ? 1.0 : N[0]*Phi[0]*v[0]*Paa; 101 79 const double S0bb = icase<2 102 ? 1.0 : N b*Phib*vb*Pbb;103 const double S0cc = N c*Phic*vc*Pcc;104 const double S0dd = N d*Phid*vd*Pdd;80 ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; 81 const double S0cc = N[2]*Phi[2]*v[2]*Pcc; 82 const double S0dd = N[3]*Phi[3]*v[3]*Pdd; 105 83 const double S0ab = icase<8 106 ? 0.0 : sqrt(N a*va*Phia*Nb*vb*Phib)*Pa*Pb;84 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 107 85 const double S0ac = icase<9 108 ? 0.0 : sqrt(N a*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);86 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 109 87 const double S0ad = icase<9 110 ? 0.0 : sqrt(N a*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);88 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 111 89 const double S0bc = (icase!=4 && icase!=7 && icase!= 9) 112 ? 0.0 : sqrt(N b*vb*Phib*Nc*vc*Phic)*Pb*Pc;90 ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 113 91 const double S0bd = (icase!=4 && icase!=7 && icase!= 9) 114 ? 0.0 : sqrt(N b*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);92 ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 115 93 const double S0cd = (icase==0 || icase==2 || icase==5) 116 ? 0.0 : sqrt(N c*vc*Phic*Nd*vd*Phid)*Pc*Pd;94 ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 117 95 #else // sasview equivalent 118 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N c,Phic,vc,Pcc);119 double S0aa = N a*Phia*va*Paa;120 double S0bb = N b*Phib*vb*Pbb;121 double S0cc = N c*Phic*vc*Pcc;122 double S0dd = N d*Phid*vd*Pdd;123 double S0ab = sqrt(N a*va*Phia*Nb*vb*Phib)*Pa*Pb;124 double S0ac = sqrt(N a*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);125 double S0ad = sqrt(N a*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);126 double S0bc = sqrt(N b*vb*Phib*Nc*vc*Phic)*Pb*Pc;127 double S0bd = sqrt(N b*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);128 double S0cd = sqrt(N c*vc*Phic*Nd*vd*Phid)*Pc*Pd;96 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); 97 double S0aa = N[0]*Phi[0]*v[0]*Paa; 98 double S0bb = N[1]*Phi[1]*v[1]*Pbb; 99 double S0cc = N[2]*Phi[2]*v[2]*Pcc; 100 double S0dd = N[3]*Phi[3]*v[3]*Pdd; 101 double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 102 double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 103 double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 104 double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 105 double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 106 double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 129 107 switch(icase){ 130 108 case 0: … … 311 289 // Note: 1e-13 to convert from fm to cm for scattering length 312 290 const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; 313 const double Lad = icase<5 ? 0.0 : (L a/va - Ld/vd)*sqrt_Nav;314 const double Lbd = icase<2 ? 0.0 : (L b/vb - Ld/vd)*sqrt_Nav;315 const double Lcd = (L c/vc - Ld/vd)*sqrt_Nav;291 const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; 292 const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; 293 const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; 316 294 317 295 const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 … … 321 299 322 300 } 323 324 double Iqxy(double qx, double qy,325 double case_num,326 double Na, double Phia, double va, double a_sld, double ba,327 double Nb, double Phib, double vb, double b_sld, double bb,328 double Nc, double Phic, double vc, double c_sld, double bc,329 double Nd, double Phid, double vd, double d_sld, double bd,330 double Kab, double Kac, double Kad,331 double Kbc, double Kbd, double Kcd332 )333 {334 double q = sqrt(qx*qx + qy*qy);335 return Iq(q,336 case_num,337 Na, Phia, va, a_sld, ba,338 Nb, Phib, vb, b_sld, bb,339 Nc, Phic, vc, c_sld, bc,340 Nd, Phid, vd, d_sld, bd,341 Kab, Kac, Kad,342 Kbc, Kbd, Kcd);343 } -
sasmodels/models/rpa.py
rec45c4f ra5b8477 86 86 # ["name", "units", default, [lower, upper], "type","description"], 87 87 parameters = [ 88 ["case_num", CASES, 0, [0, 10], "", "Component organization"],88 ["case_num", "", 1, [CASES], "", "Component organization"], 89 89 90 ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"], 91 ["Phia", "", 0.25, [0, 1], "", "volume fraction"], 92 ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"], 93 ["La", "fm", 10.0, [-inf, inf], "", "scattering length"], 94 ["ba", "Ang", 5.0, [0, inf], "", "segment length"], 95 96 ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"], 97 ["Phib", "", 0.25, [0, 1], "", "volume fraction"], 98 ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"], 99 ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"], 100 ["bb", "Ang", 5.0, [0, inf], "", "segment length"], 101 102 ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"], 103 ["Phic", "", 0.25, [0, 1], "", "volume fraction"], 104 ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"], 105 ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"], 106 ["bc", "Ang", 5.0, [0, inf], "", "segment length"], 107 108 ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"], 109 ["Phid", "", 0.25, [0, 1], "", "volume fraction"], 110 ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"], 111 ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"], 112 ["bd", "Ang", 5.0, [0, inf], "", "segment length"], 90 ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 91 ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 92 ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 93 ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 94 ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 113 95 114 96 ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"], -
sasmodels/models/spherical_sld.py
rec45c4f rd2bb604 170 170 # pylint: disable=bad-whitespace, line-too-long 171 171 # ["name", "units", default, [lower, upper], "type", "description"], 172 parameters = [["n _shells","", 1, [0, 9], "", "number of shells"],172 parameters = [["n", "", 1, [0, 9], "", "number of shells"], 173 173 ["radius_core", "Ang", 50.0, [0, inf], "", "intern layer thickness"], 174 174 ["sld_core", "1e-6/Ang^2", 2.07, [-inf, inf], "", "sld function flat"], … … 192 192 193 193 demo = dict( 194 n _shells=4,194 n=4, 195 195 scale=1.0, 196 196 solvent_sld=1.0, -
sasmodels/models/squarewell.py
rec45c4f rd2bb604 123 123 """ 124 124 125 Iqxy = """126 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);127 """128 129 125 # ER defaults to 0.0 130 126 # VR defaults to 1.0 -
sasmodels/models/stickyhardsphere.py
rec45c4f rd2bb604 171 171 """ 172 172 173 Iqxy = """174 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);175 """176 177 173 # ER defaults to 0.0 178 174 # VR defaults to 1.0
Note: See TracChangeset
for help on using the changeset viewer.