Changeset d38f244 in sasmodels


Ignore:
Timestamp:
Mar 21, 2017 9:04:50 AM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
3a45c2c
Parents:
b32caab (diff), 4050e6a (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.
Message:

Merge branch 'master' into ticket-835

Files:
10 edited

Legend:

Unmodified
Added
Removed

  • doc/ref/index.rst

    rc34a31f r9f12fbe  
    1010   refs.rst 
    1111   gpu/gpu_computations.rst 
     12   gpu/opencl_installation.rst 
    1213   magnetism/magnetism.rst 
    1314   sesans/sans_to_sesans.rst 

  • sasmodels/conversion_table.py

    r790b16bb r0c2da4b  
    549549            "radius": "core_radius", 
    550550            "sld_solvent": "core_sld", 
    551             "n_pairs": "n_pairs", 
     551            "n_shells": "n_pairs", 
    552552            "thick_shell": "s_thickness", 
    553553            "sld": "shell_sld", 

  • sasmodels/modelinfo.py

    r5124c969 r9c44b7b  
    230230    defined as a sublist with the following elements: 
    231231 
    232     *name* is the name that will be used in the call to the kernel 
    233     function and the name that will be displayed to the user.  Names 
     232    *name* is the name that will be displayed to the user.  Names 
    234233    should be lower case, with words separated by underscore.  If 
    235     acronyms are used, the whole acronym should be upper case. 
     234    acronyms are used, the whole acronym should be upper case. For vector 
     235    parameters, the name will be followed by *[len]* where *len* is an 
     236    integer length of the vector, or the name of the parameter which 
     237    controls the length.  The attribute *id* will be created from name 
     238    without the length. 
    236239 
    237240    *units* should be one of *degrees* for angles, *Ang* for lengths, 
    … …  
    603606        # Using the call_parameters table, we already have expanded forms 
    604607        # for each of the vector parameters; put them in a lookup table 
    605         expanded_pars = dict((p.name, p) for p in self.call_parameters) 
     608        # Note: p.id and p.name are currently identical for the call parameters 
     609        expanded_pars = dict((p.id, p) for p in self.call_parameters) 
    606610 
    607611        def append_group(name): 

  • sasmodels/models/fractal.c

    r925ad6e r4788822  
    11#define INVALID(p) (p.fractal_dim < 0.0) 
     2 
     3 static double 
     4 form_volume(double radius) 
     5 { 
     6     return M_4PI_3 * cube(radius); 
     7 } 
    28 
    39static double 
    … …  
    1319 
    1420    //calculate P(q) for the spherical subunits 
    15     const double V = M_4PI_3*cube(radius); 
    16     const double pq = V * square((sld_block-sld_solvent)*sas_3j1x_x(q*radius)); 
     21    const double pq = square(form_volume(radius) * (sld_block-sld_solvent) 
     22                      *sas_3j1x_x(q*radius)); 
    1723 
    1824    // scale to units cm-1 sr-1 (assuming data on absolute scale) 

  • sasmodels/models/fractal.py

    r925ad6e rdf89d77  
    2020.. math:: 
    2121 
    22     P(q)&= F(qR_0)^2 
    23  
    24     F(q)&= \frac{3 (\sin x - x \cos x)}{x^3} 
    25  
    26     V_\text{particle} &= \frac{4}{3}\ \pi R_0 
    27  
     22    P(q)&= F(qR_0)^2 \\ 
     23    F(q)&= \frac{3 (\sin x - x \cos x)}{x^3} \\ 
     24    V_\text{particle} &= \frac{4}{3}\ \pi R_0 \\ 
    2825    S(q) &= 1 + \frac{D_f\  \Gamma\!(D_f-1)}{[1+1/(q \xi)^2\  ]^{(D_f -1)/2}} 
    2926    \frac{\sin[(D_f-1) \tan^{-1}(q \xi) ]}{(q R_0)^{D_f}} 
    … …  
    3229is the fractal dimension, representing the self similarity of the structure.  
    3330Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function  
    34 diverges at $D_f$=0 and $D_f$=1.   
     31diverges at $D_f=0$ and $D_f=1$. 
    3532 
    3633**Polydispersity on the radius is provided for.** 
    … …  
    4744---------- 
    4845 
    49 J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 
     46.. [#] J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 
    5047 
    51 **Author:** NIST IGOR/DANSE **on:** pre 2010 
     48Authorship and Verification 
     49---------------------------- 
    5250 
    53 **Last Modified by:** Paul Butler **on:** March 20, 2016 
    54  
    55 **Last Reviewed by:** Paul Butler **on:** March 20, 2016 
     51* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     52* **Converted to sasmodels by:** Paul Butler **Date:** March 19, 2016 
     53* **Last Modified by:** Paul Butler **Date:** March 12, 2017 
     54* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 
    5655 
    5756""" 
    … …  
    8483parameters = [["volfraction", "", 0.05, [0.0, 1], "", 
    8584               "volume fraction of blocks"], 
    86               ["radius",    "Ang",  5.0, [0.0, inf], "", 
     85              ["radius",    "Ang",  5.0, [0.0, inf], "volume", 
    8786               "radius of particles"], 
    8887              ["fractal_dim",      "",  2.0, [0.0, 6.0], "", 

  • sasmodels/models/multilayer_vesicle.c

    rc3ccaec rec1d4bc  
    1 static 
    2 double multilayer_vesicle_kernel(double q, 
     1static double 
     2form_volume(double radius, 
     3          double thick_shell, 
     4          double thick_solvent, 
     5          double fp_n_shells) 
     6{ 
     7    int n_shells = (int)(fp_n_shells + 0.5); 
     8    double R_N = radius + n_shells*(thick_shell+thick_solvent) - thick_solvent; 
     9    return M_4PI_3*cube(R_N); 
     10} 
     11 
     12static double 
     13multilayer_vesicle_kernel(double q, 
    314          double volfraction, 
    415          double radius, 
    … …  
    718          double sld_solvent, 
    819          double sld, 
    9           int n_pairs) 
     20          int n_shells) 
    1021{ 
    1122    //calculate with a loop, two shells at a time 
    … …  
    2940 
    3041        //do 2 layers at a time 
    31         ii += 1; 
     42        ii++; 
    3243 
    33     } while(ii <= n_pairs-1);  //change to make 0 < n_pairs < 2 correspond to 
     44    } while(ii <= n_shells-1);  //change to make 0 < n_shells < 2 correspond to 
    3445                               //unilamellar vesicles (C. Glinka, 11/24/03) 
    3546 
    36     fval *= volfraction*1.0e-4*fval/voli; 
    37  
    38     return(fval); 
     47    return 1.0e-4*volfraction*fval*fval;  // Volume normalization happens in caller 
    3948} 
    4049 
    41 static 
    42 double Iq(double q, 
     50static double 
     51Iq(double q, 
    4352          double volfraction, 
    4453          double radius, 
    … …  
    4756          double sld_solvent, 
    4857          double sld, 
    49           double fp_n_pairs) 
     58          double fp_n_shells) 
    5059{ 
    51     int n_pairs = (int)(fp_n_pairs + 0.5); 
     60    int n_shells = (int)(fp_n_shells + 0.5); 
    5261    return multilayer_vesicle_kernel(q, 
    5362           volfraction, 
    … …  
    5766           sld_solvent, 
    5867           sld, 
    59            n_pairs); 
     68           n_shells); 
    6069} 
    6170 

  • sasmodels/models/multilayer_vesicle.py

    rc3ccaec r5d23de2  
    33---------- 
    44 
    5 This model is a trivial extension of the core_shell_sphere function to include 
    6 *N* shells where the core is filled with solvent and the shells are interleaved 
    7 with layers of solvent. For *N = 1*, this returns the same as the vesicle model, 
    8 except for the normalisation, which here is to outermost volume. 
    9 The shell thicknessess and SLD are constant for all shells as expected for 
    10 a multilayer vesicle. 
     5This model is a trivial extension of the core_shell_sphere function where the 
     6core is filled with solvent and is surrounded by $N$ shells of material 
     7(such as lipids) interleaved with $N - 1$ layers of solvent. For $N = 1$, this 
     8returns the same as the vesicle model, except for the normalisation, which here 
     9is to outermost volume. The shell thicknesses and SLD are constant for all 
     10shells as expected for a multilayer vesicle. 
    1111 
    1212.. figure:: img/multi_shell_geometry.jpg 
    … …  
    1919 
    2020.. math:: 
    21     P(q) = \text{scale} \cdot \frac{V_f}{V_t} F^2(q) + \text{background} 
     21    P(q) = \text{scale} \cdot \frac{\phi}{V(R_N)} F^2(q) + \text{background} 
     22 
     23where 
     24 
     25.. math:: 
     26     F(q) = (\rho_\text{shell}-\rho_\text{solv}) \sum_{i=1}^{N} \left[ 
     27     3V(r_i)\frac{\sin(qr_i) - qr_i\cos(qr_i)}{(qr_i)^3} 
     28     - 3V(R_i)\frac{\sin(qR_i) - qR_i\cos(qR_i)}{(qR_i)^3} 
     29     \right] 
    2230 
    2331for 
    2432 
    2533.. math:: 
    26     F(q) = (\rho_\text{shell}-\rho_\text{solv}) \sum_{i=1}^{n_\text{pairs}} 
    27         \left[ 
    28           3V(R_i)\frac{\sin(qR_i)-qR_i\cos(qR_i)}{(qR_i)^3} \\ 
    29           - 3V(R_i+t_s)\frac{\sin(q(R_i+t_s))-q(R_i+t_s)\cos(q(R_i+t_s))}{(q(R_i+t_s))^3} 
    30         \right] 
    3134 
    32 and 
     35     r_i &= r_c + (i-1)(t_s + t_w) && \text{ solvent radius before shell } i \\ 
     36     R_i &= r_i + t_s && \text{ shell radius for shell } i 
    3337 
    34 .. math:: 
    35      R_i = r_c + (i-1)(t_s + t_w) 
     38$\phi$ is the volume fraction of particles, $V(r)$ is the volume of a sphere 
     39of radius $r$, $r_c$ is the radius of the core, $t_s$ is the thickness of 
     40the shell, $t_w$ is the thickness of the solvent layer between the shells, 
     41$\rho_\text{shell}$ is the scattering length density of a shell, and 
     42$\rho_\text{solv}$ is the scattering length density of the solvent. 
    3643 
    37 where $V_f$ is the volume fraction of particles, $V_t$ is the volume of the 
    38 whole particle, $V(r)$ is the volume of a sphere of radius $r$, $r_c$ is the 
    39 radius of the core, $\rho_\text{shell}$ is the scattering length density of a 
    40 shell, $\rho_\text{solv}$ is the scattering length density of the solvent. 
     44USAGE NOTES 
    4145 
    42 The outer most radius, $r_o = R_n + t_s$, is used for both the volume fraction 
    43 normalization and for the effective radius for *S(Q)* when $P(Q) * S(Q)$ 
    44 is applied. 
     46* The outer-most shell radius $R_N$ is used as the effective radius 
     47  for $P(Q)$ when $P(Q) * S(Q)$ is applied. 
     48  calculations rather slow. 
     49* The number of shells is always rounded to an integer value as a non interger 
     50  number of layers is not physical. 
     51* Thus Polydispersity should only be applied to number of shells **VERY 
     52  CAREFULLY**.  A possible legitimate use would be for mixed systems in which 
     53  some vesicles have 1 shell, some have 2, etc. A polydispersity on $N$ can be 
     54  used to model the data by using the "array distriubtion" feature. First 
     55  create a file such as *shell_dist.txt* containing the relative portion 
     56  of each vesicle size:: 
     57 
     58    1 20 
     59    2  4 
     60    3  1 
     61 
     62  Turn on polydispersity and select an array distribution for the *n_shells* 
     63  parameter.  Choose the above *shell_dist.txt* file, and the model will be 
     64  computed with 80% 1-shell vesicles, 16% 2-shell vesicles and 4% 
     65  3-shell vesicles. 
     66* This is a highly non-linear, highly oscillatory (especially around the 
     67  q-values that correspond to the repeat distance of the layers), model 
     68  function complicated by the fact that the number of water/shell pairs must 
     69  physically be an integer value, although the optimization treats it as a 
     70  floating point value. Thus it may be that the resolution interpolation is not 
     71  sufficiently fine grained in certain cases. Please report any such occurences 
     72  to the SasView team. Generally, for the best possible experience: 
     73 * Start with the best possible guess 
     74 * Using a priori knowledge, hold as many parameters fixed as possible 
     75 * if N=1, tw (water thickness) must by definition be zero. Both N and tw should 
     76   be fixed during fitting. 
     77 * If N>1, use constraints to keep N > 1 
     78 * Because N only really moves in integer steps, it may get "stuck" if the 
     79   optimizer step size is too small so care should be taken 
     80   If you experience problems with this please contact the SasView team and let 
     81   them know the issue preferably with example data and model which fail to 
     82   converge. 
    4583 
    4684The 2D scattering intensity is the same as 1D, regardless of the orientation 
    … …  
    5492the :ref:`magnetism` documentation. 
    5593 
    56 This code is based on the form factor calculations implemented in the NIST 
    57 Center for Neutron Research provided c-library (Kline, 2006). 
    58  
    5994References 
    6095---------- 
    6196 
    62 B Cabane, *Small Angle Scattering Methods*, 
    63 in *Surfactant Solutions: New Methods of Investigation*, 
    64 Ch.2, Surfactant Science Series Vol. 22, Ed. R Zana and M Dekker, 
    65 New York, (1987). 
     97.. [#] B Cabane, *Small Angle Scattering Methods*, in *Surfactant Solutions: 
     98   New Methods of Investigation*, Ch.2, Surfactant Science Series Vol. 22, Ed. 
     99   R Zana and M Dekker, New York, (1987). 
    66100 
    67 **Author:** NIST IGOR/DANSE **on:** pre 2010 
     101Authorship and Verification 
     102---------------------------- 
    68103 
    69 **Last Modified by:** Piotr Rozyczko **on:** Feb 24, 2016 
    70  
    71 **Last Reviewed by:** Paul Butler **on:** March 20, 2016 
     104* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     105* **Converted to sasmodels by:** Piotr Rozyczko **Date:** Feb 24, 2016 
     106* **Last Modified by:** Paul Kienzle **Date:** Feb 7, 2017 
     107* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 
    72108 
    73109""" 
    … …  
    86122    sld_solvent: solvent scattering length density 
    87123    sld: shell scattering length density 
    88     n_pairs:number of "shell plus solvent" layer pairs 
     124    n_shells:number of "shell plus solvent" layer pairs 
    89125    background: incoherent background 
    90126        """ 
    … …  
    95131parameters = [ 
    96132    ["volfraction", "",  0.05, [0.0, 1],  "", "volume fraction of vesicles"], 
    97     ["radius", "Ang", 60.0, [0.0, inf],  "", "radius of solvent filled core"], 
    98     ["thick_shell", "Ang",        10.0, [0.0, inf],  "", "thickness of one shell"], 
    99     ["thick_solvent", "Ang",        10.0, [0.0, inf],  "", "solvent thickness between shells"], 
     133    ["radius", "Ang", 60.0, [0.0, inf],  "volume", "radius of solvent filled core"], 
     134    ["thick_shell", "Ang",        10.0, [0.0, inf],  "volume", "thickness of one shell"], 
     135    ["thick_solvent", "Ang",        10.0, [0.0, inf],  "volume", "solvent thickness between shells"], 
    100136    ["sld_solvent",    "1e-6/Ang^2",  6.4, [-inf, inf], "sld", "solvent scattering length density"], 
    101137    ["sld",   "1e-6/Ang^2",  0.4, [-inf, inf], "sld", "Shell scattering length density"], 
    102     ["n_pairs",     "",            2.0, [1.0, inf],  "", "Number of shell plus solvent layer pairs"], 
     138    ["n_shells",     "",            2.0, [1.0, inf],  "volume", "Number of shell plus solvent layer pairs"], 
    103139    ] 
    104140# pylint: enable=bad-whitespace, line-too-long 
    105141 
     142# TODO: proposed syntax for specifying which parameters can be polydisperse 
     143#polydispersity = ["radius", "thick_shell"] 
     144 
    106145source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    107146 
    108 # TODO: the following line does nothing 
    109 polydispersity = ["radius", "n_pairs"] 
     147def ER(radius, thick_shell, thick_solvent, n_shells): 
     148    n_shells = int(n_shells+0.5) 
     149    return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 
    110150 
    111151demo = dict(scale=1, background=0, 
    … …  
    116156            sld_solvent=6.4, 
    117157            sld=0.4, 
    118             n_pairs=2.0) 
     158            n_shells=2.0) 
    119159 
    120160tests = [ 
    … …  
    125165      'sld_solvent': 6.4, 
    126166      'sld': 0.4, 
    127       'n_pairs': 2.0, 
     167      'n_shells': 2.0, 
    128168      'scale': 1.0, 
    129169      'background': 0.001, 
    … …  
    136176      'sld_solvent': 6.4, 
    137177      'sld': 0.4, 
    138       'n_pairs': 2.0, 
     178      'n_shells': 2.0, 
    139179      'scale': 1.0, 
    140180      'background': 0.001, 

  • sasmodels/models/rpa.py

    reb63414 r4f9e288  
    3030    These case numbers are different from those in the NIST SANS package! 
    3131 
     32The models are based on the papers by Akcasu et al. [#Akcasu]_ and by 
     33Hammouda [#Hammouda]_ assuming the polymer follows Gaussian statistics such 
     34that $R_g^2 = n b^2/6$ where $b$ is the statistical segment length and $n$ is 
     35the number of statistical segment lengths. A nice tutorial on how these are 
     36constructed and implemented can be found in chapters 28 and 39 of Boualem 
     37Hammouda's 'SANS Toolbox'[#toolbox]_. 
     38 
     39In brief the macroscopic cross sections are derived from the general forms 
     40for homopolymer scattering and the multiblock cross-terms while the inter 
     41polymer cross terms are described in the usual way by the $\chi$ parameter. 
     42 
    3243USAGE NOTES: 
    3344 
    … …  
    3950  for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
    4051* Depending on which case is being used, the number of fitting parameters can  
    41   vary.  Note that in general the degrees of polymerization, the volume 
    42   fractions, the molar volumes, and the neutron scattering lengths for each 
    43   component are obtained from other methods and held fixed while the segment 
    44   lengths (b\ :sub:`a`, b\ :sub:`b`, etc) and $\chi$ parameters (K\ :sub:`ab`, 
    45   K\ :sub:`ac`, etc). The *scale* parameter should be held equal to unity. 
     52  vary. 
     53 
     54  .. Note:: 
     55    * In general the degrees of polymerization, the volume 
     56      fractions, the molar volumes, and the neutron scattering lengths for each 
     57      component are obtained from other methods and held fixed while The *scale* 
     58      parameter should be held equal to unity. 
     59    * The variables are normally the segment lengths (b\ :sub:`a`, b\ :sub:`b`, 
     60      etc) and $\chi$ parameters (K\ :sub:`ab`, K\ :sub:`ac`, etc). 
    4661 
    4762 
    … …  
    4964---------- 
    5065 
    51 .. [#] A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
     66.. [#Akcasu] A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 
     67   4136. 
     68.. [#Hammouda] B. Hammouda, *Advances in Polymer Science* 106 (1993) 87. 
     69.. [#toolbox] https://www.ncnr.nist.gov/staff/hammouda/the_sans_toolbox.pdf 
     70 
     71Authorship and Verification 
     72---------------------------- 
     73 
     74* **Author:** Boualem Hammouda - NIST IGOR/DANSE **Date:** pre 2010 
     75* **Converted to sasmodels by:** Paul Kienzle **Date:** July 18, 2016 
     76* **Last Modified by:** Paul Butler **Date:** March 12, 2017 
     77* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 
    5278""" 
    5379 

  • sasmodels/product.py

    r9951a86 rf88e248  
    4545    # structure factor calculator.  Structure factors should not 
    4646    # have any magnetic parameters 
    47     assert(s_info.parameters.kernel_parameters[0].id == ER_ID) 
    48     assert(s_info.parameters.kernel_parameters[1].id == VF_ID) 
    49     assert(s_info.parameters.magnetism_index == []) 
     47    if not s_info.parameters.kernel_parameters[0].id == ER_ID: 
     48        raise TypeError("S needs %s as first parameter"%ER_ID) 
     49    if not s_info.parameters.kernel_parameters[1].id == VF_ID: 
     50        raise TypeError("S needs %s as second parameter"%VF_ID) 
     51    if not s_info.parameters.magnetism_index == []: 
     52        raise TypeError("S should not have SLD parameters") 
    5053    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 
    5154    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    52     p_set = set(p.id for p in p_pars.call_parameters) 
    53     s_set = set(p.id for p in s_pars.call_parameters) 
    54  
    55     if p_set & s_set: 
    56         # there is some overlap between the parameter names; tag the 
    57         # overlapping S parameters with name_S. 
    58         # Skip the first parameter of s, which is effective radius 
    59         s_list = [(suffix_parameter(par) if par.id in p_set else par) 
    60                   for par in s_pars.kernel_parameters[1:]] 
    61     else: 
    62         # Skip the first parameter of s, which is effective radius 
    63         s_list = s_pars.kernel_parameters[1:] 
     55 
     56    # Create list of parameters for the combined model.  Skip the first 
     57    # parameter of S, which we verified above is effective radius.  If there 
     58    # are any names in P that overlap with those in S, modify the name in S 
     59    # to distinguish it. 
     60    p_set = set(p.id for p in p_pars.kernel_parameters) 
     61    s_list = [(_tag_parameter(par) if par.id in p_set else par) 
     62              for par in s_pars.kernel_parameters[1:]] 
     63    # Check if still a collision after renaming.  This could happen if for 
     64    # example S has volfrac and P has both volfrac and volfrac_S. 
     65    if any(p.id in p_set for p in s_list): 
     66        raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 
     67 
    6468    translate_name = dict((old.id, new.id) for old, new 
    6569                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    6670    demo = {} 
    67     demo.update((k, v) for k, v in p_info.demo.items() 
    68                 if k not in ("background", "scale")) 
     71    demo.update(p_info.demo.items()) 
    6972    demo.update((translate_name[k], v) for k, v in s_info.demo.items() 
    7073                if k not in ("background", "scale") and not k.startswith(ER_ID)) 
    … …  
    9093    # Remember the component info blocks so we can build the model 
    9194    model_info.composition = ('product', [p_info, s_info]) 
    92     model_info.demo = {} 
     95    model_info.demo = demo 
     96 
     97    ## Show the parameter table with the demo values 
     98    #from .compare import get_pars, parlist 
     99    #print("==== %s ====="%model_info.name) 
     100    #values = get_pars(model_info, use_demo=True) 
     101    #print(parlist(model_info, values, is2d=True)) 
    93102    return model_info 
    94103 
    95 def suffix_parameter(par, suffix): 
     104def _tag_parameter(par): 
     105    """ 
     106    Tag the parameter name with _S to indicate that the parameter comes from 
     107    the structure factor parameter set.  This is only necessary if the 
     108    form factor model includes a parameter of the same name as a parameter 
     109    in the structure factor. 
     110    """ 
    96111    par = copy(par) 
    97     par.name = par.name + " S" 
     112    # Protect against a vector parameter in S by appending the vector length 
     113    # to the renamed parameter.  Note: haven't tested this since no existing 
     114    # structure factor models contain vector parameters. 
     115    vector_length = par.name[len(par.id):] 
    98116    par.id = par.id + "_S" 
     117    par.name = par.id + vector_length 
    99118    return par 
    100119 

  • sasmodels/resolution.py

    rb397165 rb32caab  
    1717 
    1818MINIMUM_RESOLUTION = 1e-8 
    19  
    20  
    21 # When extrapolating to -q, what is the minimum positive q relative to q_min 
    22 # that we wish to calculate? 
    23 MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION = 0.01 
     19MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
    2420 
    2521class Resolution(object): 
    … …  
    8278        self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma) 
    8379                       if q_calc is None else np.sort(q_calc)) 
     80 
     81        # Protect against models which are not defined for very low q.  Limit 
     82        # the smallest q value evaluated (in absolute) to 0.02*min 
     83        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
     84        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
     85 
     86        # Build weight matrix from calculated q values 
    8487        self.weight_matrix = pinhole_resolution(self.q_calc, self.q, 
    8588                                np.maximum(q_width, MINIMUM_RESOLUTION)) 
     89        self.q_calc = abs(self.q_calc) 
    8690 
    8791    def apply(self, theory): 
    … …  
    123127        self.q_calc = slit_extend_q(q, qx_width, qy_width) \ 
    124128            if q_calc is None else np.sort(q_calc) 
     129 
     130        # Protect against models which are not defined for very low q.  Limit 
     131        # the smallest q value evaluated (in absolute) to 0.02*min 
     132        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
     133        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
     134 
     135        # Build weight matrix from calculated q values 
    125136        self.weight_matrix = \ 
    126137            slit_resolution(self.q_calc, self.q, qx_width, qy_width) 
     138        self.q_calc = abs(self.q_calc) 
    127139 
    128140    def apply(self, theory): 
    … …  
    153165    # neither trapezoid nor Simpson's rule improved the accuracy. 
    154166    edges = bin_edges(q_calc) 
    155     edges[edges < 0.0] = 0.0 # clip edges below zero 
     167    #edges[edges < 0.0] = 0.0 # clip edges below zero 
    156168    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    157169    weights = cdf[1:] - cdf[:-1] 
    … …  
    286298    # The current algorithm is a midpoint rectangle rule. 
    287299    q_edges = bin_edges(q_calc) # Note: requires q > 0 
    288     q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
     300    #q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
    289301    weights = np.zeros((len(q), len(q_calc)), 'd') 
    290302 
    … …  
    392404    interval. 
    393405 
    394     if *q_min* is zero or less then *q[0]/10* is used instead. 
     406    Note that extrapolated values may be negative. 
    395407    """ 
    396408    q = np.sort(q) 
    397409    if q_min + 2*MINIMUM_RESOLUTION < q[0]: 
    398         if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    399410        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
    400411        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    … …  
    448459        log_delta_q = log(10.) / points_per_decade 
    449460    if q_min < q[0]: 
    450         if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
     461        if q_min < 0: q_min = q[0]*MINIMUM_ABSOLUTE_Q 
    451462        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    452463        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
Note: See TracChangeset for help on using the changeset viewer.