Changeset 0c86b79 in sasmodels


Ignore:
Timestamp:
Mar 31, 2019 7:59:48 AM (7 months ago)
Author:
GitHub <noreply@…>
Parents:
1dd78a2b (diff), 225bf94 (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.
git-author:
Paul Kienzle <pkienzle@…> (03/31/19 07:59:48)
git-committer:
GitHub <noreply@…> (03/31/19 07:59:48)
Message:

Merge 225bf94c769a67d6fdab4091c85fde04b9f7db52 into 1dd78a2b6501dbd20963af88b80e45c6b89beda5

Files:
4 added
144 edited

Legend:

Unmodified
Added
Removed
  • .gitignore

    re9ed2de r6ceca44  
    88*.so 
    99*.obj 
     10*.o 
    1011/doc/_build/ 
    1112/doc/api/ 
     
    1920/.pydevproject 
    2021/.idea 
     22.vscode 
    2123/sasmodels.egg-info/ 
    2224/example/Fit_*/ 
  • sasmodels/modelinfo.py

    ra34b811 r0c86b79  
    1212from os.path import abspath, basename, splitext 
    1313import inspect 
     14import logging 
    1415 
    1516import numpy as np  # type: ignore 
     17 
     18from . import autoc 
    1619 
    1720# Optional typing 
     
    3235    TestCondition = Tuple[ParameterSetUser, TestInput, TestValue] 
    3336# pylint: enable=unused-import 
     37 
     38logger = logging.getLogger(__name__) 
    3439 
    3540# If MAX_PD changes, need to change the loop macros in kernel_iq.c 
     
    860865    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    861866    info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore 
     867    info.random = getattr(kernel_module, 'random', None) 
     868    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
     869 
     870    info.lineno = {} 
     871    _find_source_lines(info, kernel_module) 
     872    if getattr(kernel_module, 'py2c', False): 
     873        try: 
     874            warnings = autoc.convert(info, kernel_module) 
     875        except Exception as exc: 
     876            warnings = [str(exc)] 
     877        if warnings: 
     878            warnings.append("while converting %s from C to python"%name) 
     879            if len(warnings) > 2: 
     880                warnings = "\n".join(warnings) 
     881            else: 
     882                warnings = " ".join(warnings) 
     883            logger.warn(warnings) 
     884 
    862885    # Default single and opencl to True for C models.  Python models have callable Iq. 
     886    # Needs to come after autoc.convert since the Iq symbol may have been 
     887    # converted from python to C 
    863888    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) 
    864889    info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 
    865     info.random = getattr(kernel_module, 'random', None) 
    866     info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
    867  
     890     
    868891    # Set control flag for explicitly set parameters, e.g., in the RPA model. 
    869892    control = getattr(kernel_module, 'control', None) 
     
    874897        raise ValueError("oriented python models not supported") 
    875898 
    876     info.lineno = {} 
    877     _find_source_lines(info, kernel_module) 
    878899    return info 
    879900 
  • doc/guide/index.rst

    rda5536f rbc69321  
    1212   pd/polydispersity.rst 
    1313   resolution.rst 
     14   plugin.rst 
     15   fitting_sq.rst 
    1416   magnetism/magnetism.rst 
    1517   orientation/orientation.rst 
    1618   sesans/sans_to_sesans.rst 
    1719   sesans/sesans_fitting.rst 
    18    plugin.rst 
    1920   scripting.rst 
    2021   refs.rst 
  • doc/guide/plugin.rst

    r9150036 re15a822  
    303303**Note: The order of the parameters in the definition will be the order of the 
    304304parameters in the user interface and the order of the parameters in Fq(), Iq(), 
    305 Iqac(), Iqabc(), form_volume() and shell_volume(). 
     305Iqac(), Iqabc(), radius_effective(), form_volume() and shell_volume(). 
    306306And** *scale* **and** *background* **parameters are implicit to all models, 
    307307so they do not need to be included in the parameter table.** 
     
    387387can take arbitrary values, even for integer parameters, so your model should 
    388388round the incoming parameter value to the nearest integer inside your model 
    389 you should round to the nearest integer.  In C code, you can do this using:: 
     389you should round to the nearest integer.  In C code, you can do this using: 
     390 
     391.. code-block:: c 
    390392 
    391393    static double 
     
    454456............. 
    455457 
     458.. note:: 
     459 
     460   Pure python models do not yet support direct computation of $<F(Q)^2>$ or 
     461   $<F(Q)>^2$. Neither do they support orientational distributions or magnetism 
     462   (use C models if these are required). 
     463 
    456464For pure python models, define the *Iq* function:: 
    457465 
     
    499507Models should define *form_volume(par1, par2, ...)* where the parameter 
    500508list includes the *volume* parameters in order.  This is used for a weighted 
    501 volume normalization so that scattering is on an absolute scale.  If 
    502 *form_volume* is not defined, then the default *form_volume = 1.0* will be 
    503 used. 
     509volume normalization so that scattering is on an absolute scale.  For 
     510solid shapes, the *I(q)* function should use *form_volume* squared 
     511as its scale factor.  If *form_volume* is not defined, then the default 
     512*form_volume = 1.0* will be used. 
    504513 
    505514Hollow shapes, where the volume fraction of particle corresponds to the 
    506515material in the shell rather than the volume enclosed by the shape, must 
    507516also define a *shell_volume(par1, par2, ...)* function.  The parameters 
    508 are the same as for *form_volume*.  The *I(q)* calculation should use 
    509 *shell_volume* squared as its scale factor for the volume normalization. 
    510 The structure factor calculation needs *form_volume* in order to properly 
    511 scale the volume fraction parameter, so both functions are required for 
    512 hollow shapes. 
    513  
    514 Note: Pure python models do not yet support direct computation of the 
    515 average of $F(q)$ and $F^2(q)$. 
     517are the same as for *form_volume*.  Here the *I(q)* function should use 
     518*shell_volume* squared instead of *form_volume* squared so that the scale 
     519parameter corresponds to the volume fraction of material in the sample. 
     520The structure factor calculation needs the volume fraction of the filled 
     521shapes for its calculation, so the volume fraction parameter in the model 
     522is automatically scaled by *form_volume/shell_volume* prior to calling the 
     523structure factor. 
    516524 
    517525Embedded C Models 
     
    524532    """ 
    525533 
    526 This expands into the equivalent C code:: 
     534This expands into the equivalent C code: 
     535 
     536.. code-block:: c 
    527537 
    528538    double Iq(double q, double par1, double par2, ...); 
     
    535545includes only the volume parameters. 
    536546 
    537 *form_volume* defines the volume of the shell for hollow shapes. As in 
     547*shell_volume* defines the volume of the shell for hollow shapes. As in 
    538548python models, it includes only the volume parameters. 
    539549 
     
    567577Rather than returning NAN from Iq, you must define the *INVALID(v)*.  The 
    568578*v* parameter lets you access all the parameters in the model using 
    569 *v.par1*, *v.par2*, etc. For example:: 
     579*v.par1*, *v.par2*, etc. For example: 
     580 
     581.. code-block:: c 
    570582 
    571583    #define INVALID(v) (v.bell_radius < v.radius) 
     
    582594 
    583595Instead of defining the *Iq* function, models can define *Fq* as 
    584 something like:: 
     596something like: 
     597 
     598.. code-block:: c 
    585599 
    586600    double Fq(double q, double *F1, double *F2, double par1, double par2, ...); 
     
    614628laboratory frame and beam travelling along $-z$. 
    615629 
    616 The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
     630The oriented C model (oriented pure Python models are not supported) 
     631is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
    617632*par1*, etc. are the parameters to the model.  If the shape is rotationally 
    618633symmetric about *c* then *psi* is not needed, and the model is called 
     
    642657 
    643658Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the 
    644 numerical integration is then:: 
     659numerical integration is then: 
     660 
     661.. code-block:: c 
    645662 
    646663    double outer_sum = 0.0; 
     
    718735to compute the proper magnetism and orientation, which you can implement 
    719736using *Iqxy(qx, qy, par1, par2, ...)*. 
     737 
     738**Note: Magnetism is not supported in pure Python models.** 
    720739 
    721740Special Functions 
     
    9831002  memory, and wrong answers computed. The conclusion from a very long and 
    9841003  strange debugging session was that any arrays that you declare in your 
    985   model should be a multiple of four. For example:: 
     1004  model should be a multiple of four. For example: 
     1005 
     1006  .. code-block:: c 
    9861007 
    9871008      double Iq(q, p1, p2, ...) 
     
    10151036structure factor is the *hardsphere* interaction, which 
    10161037uses the effective radius of the form factor as an input to the structure 
    1017 factor model.  The effective radius is the average radius of the 
    1018 form averaged over all the polydispersity values. 
    1019  
    1020 :: 
    1021  
    1022     def ER(radius, thickness): 
    1023         """Effective radius of a core-shell sphere.""" 
    1024         return radius + thickness 
    1025  
    1026 Now consider the *core_shell_sphere*, which has a simple effective radius 
    1027 equal to the radius of the core plus the thickness of the shell, as 
    1028 shown above. Given polydispersity over *(r1, r2, ..., rm)* in radius and 
    1029 *(t1, t2, ..., tn)* in thickness, *ER* is called with a mesh 
    1030 grid covering all possible combinations of radius and thickness. 
    1031 That is, *radius* is *(r1, r2, ..., rm, r1, r2, ..., rm, ...)* 
    1032 and *thickness* is *(t1, t1, ... t1, t2, t2, ..., t2, ...)*. 
    1033 The *ER* function returns one effective radius for each combination. 
    1034 The effective radius calculator weights each of these according to 
    1035 the polydispersity distributions and calls the structure factor 
    1036 with the average *ER*. 
    1037  
    1038 :: 
    1039  
    1040     def VR(radius, thickness): 
    1041         """Sphere and shell volumes for a core-shell sphere.""" 
    1042         whole = 4.0/3.0 * pi * (radius + thickness)**3 
    1043         core = 4.0/3.0 * pi * radius**3 
    1044         return whole, whole - core 
    1045  
    1046 Core-shell type models have an additional volume ratio which scales 
    1047 the structure factor.  The *VR* function returns the volume of 
    1048 the whole sphere and the volume of the shell. Like *ER*, there is 
    1049 one return value for each point in the mesh grid. 
    1050  
    1051 *NOTE: we may be removing or modifying this feature soon. As of the 
    1052 time of writing, core-shell sphere returns (1., 1.) for VR, giving a volume 
    1053 ratio of 1.0.* 
     1038factor model.  The effective radius is the weighted average over all 
     1039values of the shape in polydisperse systems. 
     1040 
     1041There can be many notions of effective radius, depending on the shape.  For 
     1042a sphere it is clearly just the radius, but for an ellipsoid of revolution 
     1043we provide average curvature, equivalent sphere radius, minimum radius and 
     1044maximum radius.  These options are listed as *radius_effective_modes* in 
     1045the python model defintion, and must be computed by the *radius_effective* 
     1046function which takes the *radius_effective_mode* parameter as an integer, 
     1047along with the various model parameters.  Unlike normal C/Python arrays, 
     1048the first mode is 1, the second is 2, etc.  Mode 0 indicates that the 
     1049effective radius from the shape is to be ignored in favour of the the 
     1050effective radius parameter in the structure factor model. 
     1051 
     1052 
     1053Consider the core-shell sphere, which defines the following effective radius 
     1054modes in the python model:: 
     1055 
     1056    radius_effective_modes = [ 
     1057        "outer radius", 
     1058        "core radius", 
     1059    ] 
     1060 
     1061and the following function in the C-file for the model: 
     1062 
     1063.. code-block:: c 
     1064 
     1065    static double 
     1066    radius_effective(int mode, double radius, double thickness) 
     1067    { 
     1068        switch (mode) { 
     1069            case 0: return radius + thickness; 
     1070            case 1: return radius; 
     1071            default: return 0.; 
     1072        } 
     1073    } 
     1074 
     1075    static double 
     1076    form_volume(double radius, double thickness) 
     1077    { 
     1078        return M_4PI_3 * cube(radius + thickness); 
     1079    } 
     1080 
     1081Given polydispersity over *(r1, r2, ..., rm)* in radius and *(t1, t2, ..., tn)* 
     1082in thickness, *radius_effective* is called over a mesh grid covering all 
     1083possible combinations of radius and thickness, with a single *(ri, tj)* pair 
     1084in each call. The weights each of these results according to the 
     1085polydispersity distributions and calls the structure factor with the average 
     1086effective radius.  Similarly, for *form_volume*. 
     1087 
     1088Hollow models have an additional volume ratio which is needed to scale the 
     1089structure factor.  The structure factor uses the volume fraction of the filled 
     1090particles as part of its density estimate, but the scale factor for the 
     1091scattering intensity (as non-solvent volume fraction / volume) is determined 
     1092by the shell volume only.  Therefore the *shell_volume* function is 
     1093needed to compute the form:shell volume ratio, which then scales the 
     1094*volfraction* parameter prior to calling the structure factor calculator. 
     1095In the case of a hollow sphere, this would be: 
     1096 
     1097.. code-block:: c 
     1098 
     1099    static double 
     1100    shell_volume(double radius, double thickness) 
     1101    { 
     1102        double whole = M_4PI_3 * cube(radius + thickness); 
     1103        double core = M_4PI_3 * cube(radius); 
     1104        return whole - core; 
     1105    } 
     1106 
     1107If *shell_volume* is not present, then *form_volume* and *shell_volume* are 
     1108assumed to be equal, and the shape is considered solid. 
    10541109 
    10551110Unit Tests 
     
    11081163and a check that the model runs. 
    11091164 
    1110 If you are not using sasmodels from SasView, skip this step. 
    1111  
    11121165Recommended Testing 
    11131166................... 
     1167 
     1168**NB: For now, this more detailed testing is only possible if you have a 
     1169SasView build environment available!** 
    11141170 
    11151171If the model compiles and runs, you can next run the unit tests that 
     
    12481304| 2016-10-25 Steve King 
    12491305| 2017-05-07 Paul Kienzle - Moved from sasview to sasmodels docs 
     1306| 2019-03-28 Paul Kienzle - Update docs for radius_effective and shell_volume 
  • doc/guide/resolution.rst

    r0db85af rdb1d9d5  
    1 .. sm_help.rst 
    2  
    3 .. This is a port of the original SasView html help file to ReSTructured text 
    4 .. by S King, ISIS, during SasView CodeCamp-III in Feb 2015. 
     1.. resolution.rst 
     2 
     3.. This is a port of the original SasView html help file sm_help to ReSTructured 
     4.. text by S King, ISIS, during SasView CodeCamp-III in Feb 2015. 
    55 
    66 
     
    1717resolution contribution into a model calculation/simulation (which by definition 
    1818will be exact) to make it more representative of what has been measured 
    19 experimentally - a process called *smearing*. Sasmodels does the latter. 
     19experimentally - a process called *smearing*. The Sasmodels component of SasView 
     20does the latter. 
    2021 
    2122Both smearing and desmearing rely on functions to describe the resolution 
     
    2930for the instrument and stored with the data file.  If not, they will need to 
    3031be set manually before fitting. 
     32 
     33.. note:: 
     34    Problems may be encountered if the data set loaded by SasView is a 
     35    concatenation of SANS data from several detector distances where, of 
     36    course, the worst Q resolution is next to the beam stop at each detector 
     37    distance. (This will also be noticeable in the residuals plot where 
     38    there will be poor overlap). SasView sensibly orders all the input 
     39    data points by increasing Q for nicer-looking plots, however, the dQ 
     40    data can then vary considerably from point to point. If 'Use dQ data' 
     41    smearing is selected then spikes may appear in the model fits, whereas 
     42    if 'None' or 'Custom Pinhole Smear' are selected the fits look normal. 
     43 
     44    In such instances, possible solutions are to simply remove the data 
     45    with poor Q resolution from the shorter detector distances, or to fit 
     46    the data from different detector distances simultaneously. 
    3147 
    3248 
  • explore/asymint.py

    rc462169 r66dbbfb  
    2727import os, sys 
    2828sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
     29import warnings 
    2930 
    3031import numpy as np 
     
    3637 
    3738import sasmodels.special as sp 
    38  
    39 # Need to parse shape early since it determines the kernel function 
    40 # that will be used for the various integrators 
    41 shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
    42 Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
    4339 
    4440DTYPE = 'd' 
     
    201197    return norm, Fq 
    202198 
    203 if shape == 'sphere': 
    204     RADIUS = 50  # integer for the sake of mpf 
    205     NORM, KERNEL = make_sphere(radius=RADIUS) 
    206     NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 
    207 elif shape == 'cylinder': 
    208     #RADIUS, LENGTH = 10, 100000 
    209     RADIUS, LENGTH = 10, 300  # integer for the sake of mpf 
    210     NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 
    211     NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 
    212 elif shape == 'triaxial_ellipsoid': 
    213     #A, B, C = 4450, 14000, 47 
    214     A, B, C = 445, 140, 47  # integer for the sake of mpf 
    215     NORM, KERNEL = make_triaxial_ellipsoid(A, B, C) 
    216     NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv) 
    217 elif shape == 'parallelepiped': 
    218     #A, B, C = 4450, 14000, 47 
    219     A, B, C = 445, 140, 47  # integer for the sake of mpf 
    220     NORM, KERNEL = make_parallelepiped(A, B, C) 
    221     NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
    222 elif shape == 'core_shell_parallelepiped': 
    223     #A, B, C = 4450, 14000, 47 
    224     #A, B, C = 445, 140, 47  # integer for the sake of mpf 
    225     A, B, C = 114, 1380, 6800 
    226     DA, DB, DC = 21, 58, 2300 
    227     SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 
    228     ## default parameters from sasmodels 
    229     #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
    230     ## swap A-B-C to C-B-A 
    231     #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
    232     #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
    233     #SLD_SOLVENT,CONTRAST = 0, 4 
    234     if 1: # C shortest 
    235         B, C = C, B 
    236         DB, DC = DC, DB 
    237         SLDB, SLDC = SLDC, SLDB 
    238     elif 0: # C longest 
    239         A, C = C, A 
    240         DA, DC = DC, DA 
    241         SLDA, SLDC = SLDC, SLDA 
    242     #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    243     NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    244     NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
    245 elif shape == 'paracrystal': 
    246     LATTICE = 'bcc' 
    247     #LATTICE = 'fcc' 
    248     #LATTICE = 'sc' 
    249     DNN, D_FACTOR = 220, '0.06'  # mpmath needs to initialize floats from string 
    250     RADIUS = 40  # integer for the sake of mpf 
    251     NORM, KERNEL = make_paracrystal( 
    252         radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 
    253     NORM_MP, KERNEL_MP = make_paracrystal( 
    254         radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 
    255 else: 
    256     raise ValueError("Unknown shape %r"%shape) 
     199NORM = 1.0  # type: float 
     200KERNEL = None  # type: CALLABLE[[ndarray, ndarray, ndarray], ndarray] 
     201NORM_MP = 1  # type: mpf 
     202KERNEL = None  # type: CALLABLE[[mpf, mpf, mpf], mpf] 
     203 
     204SHAPES = [ 
     205    'sphere', 
     206    'cylinder', 
     207    'triaxial_ellipsoid', 
     208    'parallelepiped', 
     209    'core_shell_parallelepiped', 
     210    'fcc_paracrystal', 
     211    'bcc_paracrystal', 
     212    'sc_paracrystal', 
     213] 
     214def build_shape(shape, **pars): 
     215    global NORM, KERNEL 
     216    global NORM_MP, KERNEL_MP 
     217 
     218    # Note: using integer or string defaults for the sake of mpf 
     219    if shape == 'sphere': 
     220        RADIUS = pars.get('radius', 50) 
     221        NORM, KERNEL = make_sphere(radius=RADIUS) 
     222        NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 
     223    elif shape == 'cylinder': 
     224        #RADIUS, LENGTH = 10, 100000 
     225        RADIUS = pars.get('radius', 10) 
     226        LENGTH = pars.get('radius', 300) 
     227        NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 
     228        NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 
     229    elif shape == 'triaxial_ellipsoid': 
     230        #A, B, C = 4450, 14000, 47 
     231        A = pars.get('a', 445) 
     232        B = pars.get('b', 140) 
     233        C = pars.get('c', 47) 
     234        NORM, KERNEL = make_triaxial_ellipsoid(A, B, C) 
     235        NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv) 
     236    elif shape == 'parallelepiped': 
     237        #A, B, C = 4450, 14000, 47 
     238        A = pars.get('a', 445) 
     239        B = pars.get('b', 140) 
     240        C = pars.get('c', 47) 
     241        NORM, KERNEL = make_parallelepiped(A, B, C) 
     242        NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
     243    elif shape == 'core_shell_parallelepiped': 
     244        #A, B, C = 4450, 14000, 47 
     245        #A, B, C = 445, 140, 47  # integer for the sake of mpf 
     246        A = pars.get('a', 114) 
     247        B = pars.get('b', 1380) 
     248        C = pars.get('c', 6800) 
     249        DA = pars.get('da', 21) 
     250        DB = pars.get('db', 58) 
     251        DC = pars.get('dc', 2300) 
     252        SLDA = pars.get('slda', "5") 
     253        SLDB = pars.get('sldb', "-0.3") 
     254        SLDC = pars.get('sldc', "11.5") 
     255        ## default parameters from sasmodels 
     256        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
     257        ## swap A-B-C to C-B-A 
     258        #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
     259        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
     260        #SLD_SOLVENT,CONTRAST = 0, 4 
     261        if 1: # C shortest 
     262            B, C = C, B 
     263            DB, DC = DC, DB 
     264            SLDB, SLDC = SLDC, SLDB 
     265        elif 0: # C longest 
     266            A, C = C, A 
     267            DA, DC = DC, DA 
     268            SLDA, SLDC = SLDC, SLDA 
     269        #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
     270        NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
     271        NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
     272    elif shape.endswith('paracrystal'): 
     273        LATTICE, _ = shape.split('_') 
     274        DNN = pars.get('dnn', 220) 
     275        D_FACTOR = pars.get('d_factor', '0.06') 
     276        RADIUS = pars.get('radius', 40) 
     277        NORM, KERNEL = make_paracrystal( 
     278            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 
     279        NORM_MP, KERNEL_MP = make_paracrystal( 
     280            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 
     281    else: 
     282        raise ValueError("Unknown shape %r"%shape) 
    257283 
    258284# Note: hardcoded in mp_quad 
     
    273299 
    274300# 2D integration functions 
    275 def mp_quad_2d(q, shape): 
     301def mp_quad_2d(q): 
    276302    evals = [0] 
    277303    def integrand(theta, phi): 
     
    393419    print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
    394420    print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
     421 
     422def quadpy_method(q, rule): 
     423    """ 
     424    Use *rule*="name:index" where name and index are chosen from below. 
     425 
     426    Available rule names and the corresponding indices:: 
     427 
     428        AlbrechtCollatz: [1-5] 
     429        BazantOh: 9, 11, 13 
     430        HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3], 
     431            29, 31, 33, 35, 37, 39-[1-2] 
     432        FliegeMaier: 4, 9, 16, 25 
     433        Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 
     434            41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131 
     435        McLaren: [1-10] 
     436        Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3], 
     437            U3 11-[1-3], U3 14-1 
     438    """ 
     439    try: 
     440        import quadpy 
     441    except ImportError: 
     442        warnings.warn("use 'pip install quadpy' to enable quadpy.sphere tests") 
     443        return 
     444 
     445    from quadpy.sphere import (AlbrechtCollatz, BazantOh, HeoXu, 
     446        FliegeMaier, Lebedev, McLaren, Stroud, integrate_spherical) 
     447    RULES = { 
     448        'AlbrechtCollatz': AlbrechtCollatz, 
     449        'BazantOh': BazantOh, 
     450        'HeoXu': HeoXu, 
     451        'FliegeMaier': FliegeMaier, 
     452        'Lebedev': Lebedev, 
     453        'McLaren': McLaren, 
     454        'Stroud': Stroud, 
     455    } 
     456    int_index = 'AlbrechtCollatz', 'McLaren' 
     457 
     458    rule_name, rule_index = rule.split(':') 
     459    index = int(rule_index) if rule_name in int_index else rule_index 
     460    rule_obj = RULES[rule_name](index) 
     461    fn = lambda azimuthal, polar: kernel_2d(q=q, theta=polar, phi=azimuthal) 
     462    Iq = integrate_spherical(fn, rule=rule_obj)/(4*pi) 
     463    print("%s degree=%d points=%s => %.15g" 
     464          % (rule, rule_obj.degree, len(rule_obj.points), Iq)) 
    395465 
    396466def plot_2d(q, n=300): 
     
    414484    pylab.show() 
    415485 
    416 def main(Qstr): 
    417     Q = float(Qstr) 
    418     if shape == 'sphere': 
     486def main(): 
     487    import argparse 
     488 
     489    parser = argparse.ArgumentParser( 
     490        description="asymmetric integration explorer", 
     491        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     492        ) 
     493    parser.add_argument('-s', '--shape', choices=SHAPES, 
     494                        default='parallelepiped', 
     495                        help='oriented shape') 
     496    parser.add_argument('-q', '--q_value', type=str, default='0.005', 
     497                        help='Q value to evaluate') 
     498    parser.add_argument('pars', type=str, nargs='*', default=[], 
     499                        help='p=val for p in shape parameters') 
     500    opts = parser.parse_args() 
     501    pars = {k: v for par in opts.pars for k, v in [par.split('=')]} 
     502    build_shape(opts.shape, **pars) 
     503 
     504    Q = float(opts.q_value) 
     505    if opts.shape == 'sphere': 
    419506        print("exact", NORM*sp.sas_3j1x_x(Q*RADIUS)**2) 
    420     print("gauss-20", *gauss_quad_2d(Q, n=20)) 
    421     print("gauss-76", *gauss_quad_2d(Q, n=76)) 
    422     print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    423     print("gauss-500", *gauss_quad_2d(Q, n=500)) 
    424     print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
    425     print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 
    426     print("gauss-20 usub", *gauss_quad_usub(Q, n=20)) 
    427     print("gauss-76 usub", *gauss_quad_usub(Q, n=76)) 
    428     print("gauss-150 usub", *gauss_quad_usub(Q, n=150)) 
     507 
     508    # Methods from quadpy, if quadpy is available 
     509    #  AlbrechtCollatz: [1-5] 
     510    #  BazantOh: 9, 11, 13 
     511    #  HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3], 
     512    #     29, 31, 33, 35, 37, 39-[1-2] 
     513    #  FliegeMaier: 4, 9, 16, 25 
     514    #  Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 
     515    #     41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131 
     516    #  McLaren: [1-10] 
     517    #  Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3], 
     518    #     U3 11-[1-3], U3 14-1 
     519    quadpy_method(Q, "AlbrechtCollatz:5") 
     520    quadpy_method(Q, "HeoXu:39-2") 
     521    quadpy_method(Q, "FliegeMaier:25") 
     522    quadpy_method(Q, "Lebedev:19") 
     523    quadpy_method(Q, "Lebedev:131") 
     524    quadpy_method(Q, "McLaren:10") 
     525    quadpy_method(Q, "Stroud:U3 14-1") 
     526 
     527    print("gauss-20 points=%d => %.15g" % gauss_quad_2d(Q, n=20)) 
     528    print("gauss-76 points=%d => %.15g" % gauss_quad_2d(Q, n=76)) 
     529    print("gauss-150 points=%d => %.15g" % gauss_quad_2d(Q, n=150)) 
     530    print("gauss-500 points=%d => %.15g" % gauss_quad_2d(Q, n=500)) 
     531    print("gauss-1025 points=%d => %.15g" % gauss_quad_2d(Q, n=1025)) 
     532    print("gauss-2049 points=%d => %.15g" % gauss_quad_2d(Q, n=2049)) 
     533    print("gauss-20 usub points=%d => %.15g" % gauss_quad_usub(Q, n=20)) 
     534    print("gauss-76 usub points=%d => %.15g" % gauss_quad_usub(Q, n=76)) 
     535    print("gauss-150 usub points=%d => %.15g" % gauss_quad_usub(Q, n=150)) 
     536 
    429537    #gridded_2d(Q, n=2**8+1) 
    430538    gridded_2d(Q, n=2**10+1) 
    431539    #gridded_2d(Q, n=2**12+1) 
    432540    #gridded_2d(Q, n=2**15+1) 
    433     if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
    434         # adaptive forms on models for which the calculations are fast enough 
     541    # adaptive forms on models for which the calculations are fast enough 
     542    SLOW_SHAPES = { 
     543        'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal', 
     544        'core_shell_parallelepiped', 
     545    } 
     546    if opts.shape not in SLOW_SHAPES: 
    435547        print("dblquad", *scipy_dblquad_2d(Q)) 
    436548        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 
    437549        print("romberg", *scipy_romberg_2d(Q)) 
    438550        with mp.workprec(100): 
    439             print("mpmath", *mp_quad_2d(mp.mpf(Qstr), shape)) 
     551            print("mpmath", *mp_quad_2d(mp.mpf(opts.q_value))) 
    440552    plot_2d(Q, n=200) 
    441553 
    442554if __name__ == "__main__": 
    443     main(Qstr) 
     555    main() 
  • explore/beta/sasfit_compare.py

    r119073a ra34b811  
    408408        #radius_effective=12.59921049894873, 
    409409        ) 
    410     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     410    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    411411    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    412412    title = " ".join(("sasmodels", model, pd_type)) 
  • explore/beta/sasfit_compare_new.py

    r2586ab72 ra34b811  
    2121#Used to calculate F(q) for the cylinder, sphere, ellipsoid models 
    2222# RKH adding vesicle and hollow_cylinder to test sasview special cases of ER and VR 
    23 # BEWARE there are issues here if you run this from python3 (i.e. qt5 version of sasview), so do "source activate sasview" 
     23# There were issues here from python3 (i.e. qt5 version of sasview),fixed by Paul K (else do "source activate sasview") 
    2424def sas_sinx_x(x): 
    2525    with np.errstate(all='ignore'): 
     
    275275    if radius_effective is None: 
    276276        radius_effective = radius 
    277         print("radius_effective  ",radius_effective) 
    278277        average_volume = total_volume/total_weight 
    279278    if norm == 'sasfit': 
     
    285284        F2 *= 1e-12 
    286285        IQD = F2/average_volume*1e8*volfraction 
     286    print("in sphere_r : radius_effective  ",radius_effective,"  volfraction",volfraction) 
    287287    SQ = hardsphere_simple(q, radius_effective, volfraction) 
    288288    beta = F1**2/F2 
     
    332332    if radius_effective is None: 
    333333        radius_effective = radius_eff/total_weight 
    334         print("radius_effective  ",radius_effective) 
     334        print("in vesicle_pe : radius_effective  ",radius_effective) 
    335335    if norm == 'sasfit': 
    336336        IQD = F2 
     
    348348        F2 *= 1e-12 
    349349        IQD = F2/average_volume*1e8*volfraction 
    350     # RKH SORT THIS OUT WHEN HAVE NEW MODEL 
    351     vfff=volfraction*(30./20.)**3 
    352     print("radius_effective, fudged volfaction for S(Q) ",radius_effective,vfff) 
     350    # RKH SORT THIS OUT WHEN HAVE NEW MODEL  Vtot/Vshell = (R+T)^3/ ( (R+T)^3 -R^3 ) 
     351    vfff=volfraction*( 1.0/(1.0 - ( (radius/(radius+thickness))**3 ))) 
     352    print("in vesicle_pe : radius_effective, fudged volfraction for S(Q) ",radius_effective,vfff) 
    353353    SQ = hardsphere_simple(q, radius_effective, vfff) 
    354354    beta = F1**2/F2 
     
    359359# 
    360360#polydispersity for hollow_cylinder 
    361 def hollow_cylinder_pe(q,radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.15, 
    362         radius_pd=0.1, thickness_pd=0.2, length_pd=0.05, radius_pd_type="gaussian", thickness_pd_type="gaussian", 
    363         radius_effective=None, background=0, scale=1, 
     361def hollow_cylinder_pe(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1, volfraction=0.15, 
     362        radius_pd=0.1, thickness_pd=0.2, length_pd=0.05, radius_pd_type="gaussian", length_pd_type="gaussian", 
     363        thickness_pd_type="gaussian", radius_effective=None, background=0, scale=1, 
    364364                 norm='sasview'): 
    365365    if norm not in ['sasview', 'sasfit', 'yun']: 
     
    374374        T_val, T_prob = schulz_distribution(thickness, thickness_pd*thickness, 0, inf) 
    375375    if length_pd_type == 'gaussian': 
    376         T_val, T_prob = gaussian_distribution(length, length_pd*length, 0, inf) 
     376        L_val, L_prob = gaussian_distribution(length, length_pd*length, 0, inf) 
    377377    elif length_pd_type == 'schulz': 
    378378        L_val, L_prob = schulz_distribution(length, length_pd*length, 0, inf) 
     
    402402    if radius_effective is None: 
    403403        radius_effective = radius_eff/total_weight 
    404         print("radius_effective  ",radius_effective) 
     404        print("in hollow_cylinder : radius_effective  ",radius_effective) 
    405405    if norm == 'sasfit': 
    406406        IQD = F2 
     
    420420# RKH SORT THIS OUT WHEN HAVE NEW MODEL 
    421421    vfff=volfraction 
    422     print("radius_effective, volfaction for S(Q) ",radius_effective,vfff) 
     422    print("in hollow_cylinder : radius_effective, volfaction for S(Q) ",radius_effective,vfff) 
    423423    SQ = hardsphere_simple(q, radius_effective, vfff) 
    424424    beta = F1**2/F2 
     
    548548    #Iq = None#= Sq = None 
    549549    r = dict(I._kernel.results()) 
    550     #print(r) 
     550    print(r) 
    551551    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r["S_eff(Q)"], Ibeta=Iq) 
    552552 
     
    591591        radius_pd=.1, radius_pd_type=pd_type, 
    592592        volfraction=0.15, 
    593         radius_effective=12.59921049894873,  # equivalent average sphere radius 
     593        radius_effective=20  # equivalent average sphere radius 
     594#       NOTE sasview computes its own radius_effective in "target" (the print(r) at end of sasmodels_theory will reveal its value), 
     595#       this one is only used locally for "actual" 
    594596        ) 
    595     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     597    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    596598    actual = sphere_r(q, norm='sasview', **pars) 
    597599    title = " ".join(("sasmodels", model, pd_type)) 
     
    604606    q = np.logspace(-5, 0, 250) 
    605607    model = 'vesicle' 
    606     print("F2 seems OK, but big issues with S(Q), wo work in progress") 
    607 # NOTE parameers for vesicle model are soon to be changed to remove "volfraction" 
     608    print("F2 seems OK, but big issues with S(Q), so work in progress") 
     609# NOTE parameters for vesicle model are soon to be changed to remove "volfraction" (though still there in 5.0) 
    608610    pars = dict( 
    609611        radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.03, 
    610612        background=0, 
    611         radius_pd=0.1, thickness_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type 
    612         #,radius_effective=12.59921049894873,  # equivalent average sphere radius 
    613         ) 
    614     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     613        radius_pd=0.01, thickness_pd=0.01, radius_pd_type=pd_type, thickness_pd_type=pd_type, 
     614        radius_effective=30. ) 
     615        # equivalent average sphere radius for local "actual" to match what sasview uses, use None to compute average outer radius here, 
     616 
     617    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    615618    actual = vesicle_pe(q, norm='sasview', **pars) 
    616619    title = " ".join(("sasmodels", model, pd_type)) 
     
    629632        radius=20, thickness=10, length=80, sld=4, sld_solvent=1, 
    630633        background=0, 
    631         radius_pd=0.1, thickness_pd=0.0, length_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type, length_pd_type=pd_type 
    632         #,radius_effective=12.59921049894873,  # equivalent average sphere radius 
    633         ) 
    634     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     634        radius_pd=0.1, thickness_pd=0.0, length_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type, length_pd_type=pd_type, 
     635        radius_effective=40.687) 
     636        # equivalent average sphere radius for local "actual" to match what sasview uses 
     637    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    635638    actual = hollow_cylinder_pe(q, norm='sasview', **pars) 
    636639# RKH monodisp was OK,    actual = hollow_cylinder_theta(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1 ) 
     
    651654        volfraction=0.15, 
    652655        radius_effective=270.7543927018, 
    653         #radius_effective=12.59921049894873, 
     656# if change radius_effective to some other value, the S(Q) from sasview does not agree 
    654657        ) 
    655     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     658    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    656659    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    657660# RKH test       actual = ellipsoid_theta(q, radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, volfraction=0.15, radius_effective=270.) 
     
    751754    } 
    752755 
    753     Q, IQ = load_sasfit(data_file('richard_test.txt')) 
    754     Q, IQSD = load_sasfit(data_file('richard_test2.txt')) 
    755     Q, IQBD = load_sasfit(data_file('richard_test3.txt')) 
     756    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt')) 
     757    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt')) 
     758    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt')) 
    756759    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    757760    actual = sphere_r(Q, norm="sasfit", **pars) 
     
    772775    } 
    773776 
    774     Q, IQ = load_sasfit(data_file('richard_test4.txt')) 
    775     Q, IQSD = load_sasfit(data_file('richard_test5.txt')) 
    776     Q, IQBD = load_sasfit(data_file('richard_test6.txt')) 
     777    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQD.txt')) 
     778    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQSD.txt')) 
     779    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQBD.txt')) 
    777780    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    778781    actual = ellipsoid_pe(Q, norm="sasfit", **pars) 
  • extra/build_linux.sh

    rdb8756e r0d362b7  
    1414# Run tests 
    1515STATUS=0 
    16 python -m sasmodels.model_test opencl_and_dll all || STATUS=$? 
     16python -m sasmodels.model_test -v all || STATUS=$? 
    1717python -m sasmodels.resolution || STATUS=$? 
    1818 
  • extra/pylint.rc

    r6e68289 rdf20715  
    6767    locally-enabled, 
    6868    locally-disabled, 
     69    #no-else-return, 
     70    #invalid-name, 
     71    fixme, 
     72    too-many-lines, 
     73    too-many-locals, 
     74    too-many-branches, 
     75    too-many-statements, 
     76    #too-many-instance-attributes, 
     77    #too-many-return-statements, 
     78    #no-self-use, 
    6979 
    7080[REPORTS] 
     
    109119good-names=_, 
    110120    input, id, 
    111     h, i, j, k, n, p, v, w, x, y, z, 
    112     q, qx, qy, qz, Q, Qx, Qy, Qz, 
    113     dt, dx, dy, dz, dq, 
     121    h, i, j, k, n, p, u, v, w, x, y, z, r, dr, 
     122    q, qx, qy, qz, Q, Qx, Qy, Qz, nx, ny, nz, nq, nr, 
     123    dt, dx, dy, dz, dq, S, P, a, b, c, 
    114124    Iq, dIq, Iqxy, Iq_calc, 
    115125    ER, call_ER, VR, call_VR, 
     
    133143 
    134144# Regular expression matching correct variable names 
    135 variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{2,30}$ 
     145variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{1,30}$ 
    136146 
    137147# Naming hint for variable names 
     
    151161 
    152162# Regular expression matching correct argument names 
    153 argument-rgx=[a-z_][a-z0-9_]{2,30}$ 
     163argument-rgx=[a-z_][a-z0-9_]{1,30}$ 
    154164 
    155165# Naming hint for argument names 
     
    280290# (useful for modules/projects where namespaces are manipulated during runtime 
    281291# and thus existing member attributes cannot be deduced by static analysis 
    282 ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special 
     292ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special,pyopencl,pycuda,matplotlib.cm,ipyvolume,ipywidget 
    283293 
    284294# List of classes names for which member attributes should not be checked 
  • sasmodels/bumps_model.py

    r2c4a190 rb297ba9  
    2626    from .kernel import KernelModel 
    2727    from .modelinfo import ModelInfo 
     28    from .resolution import Resolution 
    2829    Data = Union[Data1D, Data2D] 
    2930except ImportError: 
     
    180181    @property 
    181182    def resolution(self): 
     183        # type: () -> Union[None, Resolution] 
     184        """ 
     185        :class:`sasmodels.Resolution` applied to the data, if any. 
     186        """ 
    182187        return self._resolution 
    183188 
    184189    @resolution.setter 
    185190    def resolution(self, value): 
     191        # type: (Resolution) -> None 
     192        """ 
     193        :class:`sasmodels.Resolution` applied to the data, if any. 
     194        """ 
    186195        self._resolution = value 
    187196 
     
    270279        Save the model parameters and data into a file. 
    271280 
    272         Not Implemented. 
     281        Not Implemented except for sesans fits. 
    273282        """ 
    274283        if self.data_type == "sesans": 
  • sasmodels/compare.py

    rc1799d3 rb297ba9  
    725725        set_integration_size(model_info, ngauss) 
    726726 
    727     if (dtype != "default" and not dtype.endswith('!')  
     727    if (dtype != "default" and not dtype.endswith('!') 
    728728            and not (kernelcl.use_opencl() or kernelcuda.use_cuda())): 
    729729        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) 
     
    772772        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
    773773        if opts['pars'] is None: 
    774             return 
     774            return limits 
    775775        result = run_models(opts, verbose=True) 
    776776        if opts['plot']: 
     
    822822    args = dict((k, v) for k, v in args.items() 
    823823                if "_pd" not in k 
    824                    and ":" not in k 
    825                    and k not in ("background", "scale", "theta", "phi", "psi")) 
     824                and ":" not in k 
     825                and k not in ("background", "scale", "theta", "phi", "psi")) 
    826826    args = args.copy() 
    827827 
     
    896896    # work with trimmed data, not the full set 
    897897    sorted_err = np.sort(abs(err.compressed())) 
    898     if len(sorted_err) == 0: 
     898    if sorted_err.size == 0: 
    899899        print(label + "  no valid values") 
    900900        return 
  • sasmodels/convert.py

    r21c93c3 ref476e6  
    384384 
    385385def revert_name(model_info): 
     386    """Translate model name back to the name used in SasView 3.x""" 
    386387    oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    387388    return oldname 
     
    653654 
    654655def test_backward_forward(): 
     656    """ 
     657    Test conversion of model parameters from 4.x to 3.x and back. 
     658    """ 
    655659    from .core import list_models 
    656660    L = lambda name: _check_one(name, seed=1) 
  • sasmodels/core.py

    rd92182f r9562dd2  
    323323 
    324324def test_composite_order(): 
     325    """ 
     326    Check that mixture models produce the same result independent of ordder. 
     327    """ 
    325328    def test_models(fst, snd): 
    326329        """Confirm that two models produce the same parameters""" 
     
    337340 
    338341    def build_test(first, second): 
     342        """Construct pair model test""" 
    339343        test = lambda description: test_models(first, second) 
    340344        description = first + " vs. " + second 
     
    372376    # type: () -> None 
    373377    """Check that model load works""" 
     378    from .product import RADIUS_ID, VOLFRAC_ID, STRUCTURE_MODE_ID, RADIUS_MODE_ID 
    374379    #Test the the model produces the parameters that we would expect 
    375380    model = load_model("cylinder@hardsphere*sphere") 
    376381    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    377     target = ("sld sld_solvent radius length theta phi" 
    378               " radius_effective volfraction " 
    379               " structure_factor_mode radius_effective_mode" 
    380               " A_sld A_sld_solvent A_radius").split() 
     382    target = ["sld", "sld_solvent", "radius", "length", "theta", "phi", 
     383              RADIUS_ID, VOLFRAC_ID, STRUCTURE_MODE_ID, RADIUS_MODE_ID, 
     384              "A_sld", "A_sld_solvent", "A_radius"] 
    381385    assert target == actual, "%s != %s"%(target, actual) 
    382386 
    383387def list_models_main(): 
    384     # type: () -> None 
     388    # type: () -> int 
    385389    """ 
    386390    Run list_models as a main program.  See :func:`list_models` for the 
     
    391395    try: 
    392396        models = list_models(kind) 
    393     except Exception as exc: 
     397        print("\n".join(models)) 
     398    except Exception: 
    394399        print(list_models.__doc__) 
    395400        return 1 
    396  
    397     print("\n".join(list_models(kind))) 
     401    return 0 
    398402 
    399403if __name__ == "__main__": 
  • sasmodels/data.py

    rc88f983 rb297ba9  
    500500            mdata = masked_array(data.y, data.mask.copy()) 
    501501            mdata[~np.isfinite(mdata)] = masked 
    502             if view is 'log': 
     502            if view == 'log': 
    503503                mdata[mdata <= 0] = masked 
    504504            plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.') 
     
    514514            mtheory = masked_array(theory) 
    515515            mtheory[~np.isfinite(mtheory)] = masked 
    516             if view is 'log': 
     516            if view == 'log': 
    517517                mtheory[mtheory <= 0] = masked 
    518518            plt.plot(theory_x, scale*mtheory, '-') 
  • sasmodels/details.py

    r885753a rb297ba9  
    236236    npars = kernel.info.parameters.npars 
    237237    nvalues = kernel.info.parameters.nvalues 
    238     scalars = [value for value, _dispersity, _weight in mesh] 
     238    scalars = [value for value, dispersity, weight in mesh] 
    239239    # skipping scale and background when building values and weights 
    240     _values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
    241     #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights) 
    242     length = np.array([len(w) for w in weights]) 
     240    value, dispersity, weight = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
     241    #weight = correct_theta_weights(kernel.info.parameters, dispersity, weight) 
     242    length = np.array([len(w) for w in weight]) 
    243243    offset = np.cumsum(np.hstack((0, length))) 
    244244    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
     
    246246    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    247247    extra = (32 - data_len%32)%32 
    248     data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra]) 
     248    data = np.hstack((scalars,) + dispersity + weight + ZEROS[:extra]) 
    249249    data = data.astype(kernel.dtype) 
    250250    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
  • sasmodels/direct_model.py

    r9150036 re220acc  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
    33 from .modelinfo import DEFAULT_BACKGROUND 
    3433from .product import RADIUS_MODE_ID 
    3534 
     
    6059    """ 
    6160    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
    62     #print("pars", list(zip(*mesh))[0]) 
     61    #print("in call_kernel: pars:", list(zip(*mesh))[0]) 
    6362    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
    64     #print("values:", values) 
     63    #print("in call_kernel: values:", values) 
    6564    return calculator(call_details, values, cutoff, is_magnetic) 
    6665 
     
    7372 
    7473    Use parameter *radius_effective_mode* to select the effective radius 
    75     calculation. 
     74    calculation to use amongst the *radius_effective_modes* list given in the 
     75    model. 
    7676    """ 
    7777    R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0)) 
    7878    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
    79     #print("pars", list(zip(*mesh))[0]) 
     79    #print("in call_Fq: pars", list(zip(*mesh))[0]) 
    8080    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
    81     #print("values:", values) 
     81    #print("in call_Fq: values:", values) 
    8282    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) 
    8383 
     
    119119        active = lambda name: True 
    120120 
    121     #print("pars",[p.id for p in parameters.call_parameters]) 
     121    #print("in get_mesh: pars:",[p.id for p in parameters.call_parameters]) 
    122122    mesh = [_get_par_weights(p, values, active(p.name)) 
    123123            for p in parameters.call_parameters] 
  • sasmodels/generate.py

    ra8a1d48 ra34b811  
    2727    which are hollow. 
    2828 
    29     *effective_radius(mode, p1, p2, ...)* returns the effective radius of 
     29    *radius_effective(mode, p1, p2, ...)* returns the effective radius of 
    3030    the form with particular dimensions.  Mode determines the type of 
    3131    effective radius returned, with mode=1 for equivalent volume. 
     
    190190 
    191191def get_data_path(external_dir, target_file): 
     192    """ 
     193    Search for the target file relative in the installed application. 
     194 
     195    Search first in the location of the generate module in case we are 
     196    running directly from the distribution.  Search next to the python 
     197    executable for windows installs.  Search in the ../Resources directory 
     198    next to the executable for Mac OS/X installs. 
     199    """ 
    192200    path = abspath(dirname(__file__)) 
    193201    if exists(joinpath(path, target_file)): 
     
    225233# 
    226234# NOTE: there is an RST_PROLOG at the end of this file which is NOT 
    227 # used for the bundled documentation. Still as long as we are defining the macros 
    228 # in two places any new addition should define the macro in both places. 
     235# used for the bundled documentation. Still as long as we are defining the 
     236# macros in two places any new addition should define the macro in both places. 
    229237RST_UNITS = { 
    230238    "Ang": "|Ang|", 
     
    530538 
    531539def test_tag_float(): 
    532     """check that floating point constants are properly identified and tagged with 'f'""" 
     540    """Check that floating point constants are identified and tagged with 'f'""" 
    533541 
    534542    cases = """ 
     
    707715    return False 
    708716 
    709 _SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) 
     717_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", 
     718                                   flags=re.MULTILINE) 
    710719def contains_shell_volume(source): 
    711720    # type: (List[str]) -> bool 
     
    801810            logger.warn("oriented shapes should define Iqac or Iqabc") 
    802811        else: 
    803             raise ValueError("Expected function Iqac or Iqabc for oriented shape") 
     812            raise ValueError("Expected Iqac or Iqabc for oriented shape") 
    804813 
    805814    # Define the parameter table 
     
    811820                              for p in partable.kernel_parameters)) 
    812821    # Define the function calls 
    813     call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0" 
     822    call_radius_effective = "#define CALL_RADIUS_EFFECTIVE(_mode, _v) 0.0" 
    814823    if partable.form_volume_parameters: 
    815824        refs = _call_pars("_v.", partable.form_volume_parameters) 
    816825        if is_hollow: 
    817             call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2) 
     826            call_volume = ( 
     827                "#define CALL_VOLUME(_form, _shell, _v) " 
     828                "do { _form = form_volume(%s); _shell = shell_volume(%s); } " 
     829                "while (0)") % ((",".join(refs),)*2) 
    818830        else: 
    819             call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs)) 
    820         if model_info.effective_radius_type: 
    821             call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs)) 
     831            call_volume = ( 
     832                "#define CALL_VOLUME(_form, _shell, _v) " 
     833                "do { _form = _shell = form_volume(%s); } " 
     834                "while (0)") % (",".join(refs)) 
     835        if model_info.radius_effective_modes: 
     836            call_radius_effective = ( 
     837                "#define CALL_RADIUS_EFFECTIVE(_mode, _v) " 
     838                "radius_effective(_mode, %s)") % (",".join(refs)) 
    822839    else: 
    823840        # Model doesn't have volume.  We could make the kernel run a little 
    824841        # faster by not using/transferring the volume normalizations, but 
    825842        # the ifdef's reduce readability more than is worthwhile. 
    826         call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)" 
     843        call_volume = ( 
     844            "#define CALL_VOLUME(_form, _shell, _v) " 
     845            "do { _form = _shell = 1.0; } while (0)") 
    827846    source.append(call_volume) 
    828     source.append(call_effective_radius) 
     847    source.append(call_radius_effective) 
    829848    model_refs = _call_pars("_v.", partable.iq_parameters) 
    830849 
     
    879898    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    880899    source.append("#define PROJECTION %d"%PROJECTION) 
    881     # TODO: allow mixed python/opencl kernels? 
    882     ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
    883     dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
    884  
    885     result = { 
    886         'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    887         'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    888     } 
     900    wrappers = _kernels(kernel_code, call_iq, clear_iq, 
     901                        call_iqxy, clear_iqxy, model_info.name) 
     902    code = '\n'.join(source + wrappers[0] + wrappers[1] + wrappers[2]) 
     903 
     904    # Note: Identical code for dll and opencl.  This may not always be the case 
     905    # so leave them as separate entries in the returned value. 
     906    result = {'dll': code, 'opencl': code} 
    889907    return result 
    890908 
  • sasmodels/gengauss.py

    rff31782 rb297ba9  
    11#!/usr/bin/env python 
     2""" 
     3Generate the Gauss-Legendre integration points and save them as a C file. 
     4""" 
    25from __future__ import division, print_function 
    3  
    4 import sys 
    56 
    67import numpy as np 
     
    89 
    910def gengauss(n, path): 
     11    """ 
     12    Save the Gauss-Legendre integration points for length *n* into file *path*. 
     13    """ 
    1014    z, w = leggauss(n) 
    1115 
  • sasmodels/guyou.py

    rb3703f5 rb297ba9  
    6868 
    6969def ellipticJi(u, v, m): 
     70    """Returns [sn, cn, dn](u + iv|m).""" 
    7071    scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m) 
    7172    u, v, m = np.broadcast_arrays(u, v, m) 
     
    100101    ] 
    101102 
    102 # calculate F(phi+ipsi|m). 
    103 # see Abramowitz and Stegun, 17.4.11. 
    104103def ellipticFi(phi, psi, m): 
     104    """Returns F(phi+ipsi|m). See Abramowitz and Stegun, 17.4.11.""" 
    105105    if np.any(phi == 0): 
    106106        scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m) 
     
    207207        plt.plot(x, y, 'g') 
    208208 
    209     for long in range(-limit, limit+1, step): 
    210         x, y = guyou_invert(scale*long, scale*long_line) 
     209    for longitude in range(-limit, limit+1, step): 
     210        x, y = guyou_invert(scale*longitude, scale*long_line) 
    211211        plt.plot(x, y, 'b') 
    212212    plt.xlabel('longitude') 
  • sasmodels/jitter.py

    r7d97437 r4e28511  
    11#!/usr/bin/env python 
     2# -*- coding: utf-8 -*- 
    23""" 
    34Jitter Explorer 
     
    56 
    67Application to explore orientation angle and angular dispersity. 
     8 
     9From the command line:: 
     10 
     11    # Show docs 
     12    python -m sasmodels.jitter --help 
     13 
     14    # Guyou projection jitter, uniform over 20 degree theta and 10 in phi 
     15    python -m sasmodels.jitter --projection=guyou --dist=uniform --jitter=20,10,0 
     16 
     17From a jupyter cell:: 
     18 
     19    import ipyvolume as ipv 
     20    from sasmodels import jitter 
     21    import importlib; importlib.reload(jitter) 
     22    jitter.set_plotter("ipv") 
     23 
     24    size = (10, 40, 100) 
     25    view = (20, 0, 0) 
     26 
     27    #size = (15, 15, 100) 
     28    #view = (60, 60, 0) 
     29 
     30    dview = (0, 0, 0) 
     31    #dview = (5, 5, 0) 
     32    #dview = (15, 180, 0) 
     33    #dview = (180, 15, 0) 
     34 
     35    projection = 'equirectangular' 
     36    #projection = 'azimuthal_equidistance' 
     37    #projection = 'guyou' 
     38    #projection = 'sinusoidal' 
     39    #projection = 'azimuthal_equal_area' 
     40 
     41    dist = 'uniform' 
     42    #dist = 'gaussian' 
     43 
     44    jitter.run(size=size, view=view, jitter=dview, dist=dist, projection=projection) 
     45    #filename = projection+('_theta' if dview[0] == 180 else '_phi' if dview[1] == 180 else '') 
     46    #ipv.savefig(filename+'.png') 
    747""" 
    848from __future__ import division, print_function 
     
    1050import argparse 
    1151 
    12 try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d 
    13     import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot 
    14 except ImportError: 
    15     pass 
    16  
    17 import matplotlib as mpl 
    18 import matplotlib.pyplot as plt 
    19 from matplotlib.widgets import Slider 
    2052import numpy as np 
    21 from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    22  
    23 def draw_beam(axes, view=(0, 0)): 
     53from numpy import pi, cos, sin, sqrt, exp, log, degrees, radians, arccos, arctan2 
     54 
     55# Too many complaints about variable names from pylint: 
     56#    a, b, c, u, v, x, y, z, dx, dy, dz, px, py, pz, R, Rx, Ry, Rz, ... 
     57# pylint: disable=invalid-name 
     58 
     59def draw_beam(axes, view=(0, 0), alpha=0.5, steps=25): 
    2460    """ 
    2561    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
    2662    """ 
    2763    #axes.plot([0,0],[0,0],[1,-1]) 
    28     #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    29  
    30     steps = 25 
    31     u = np.linspace(0, 2 * np.pi, steps) 
    32     v = np.linspace(-1, 1, steps) 
     64    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=alpha) 
     65 
     66    u = np.linspace(0, 2 * pi, steps) 
     67    v = np.linspace(-1, 1, 2) 
    3368 
    3469    r = 0.02 
    35     x = r*np.outer(np.cos(u), np.ones_like(v)) 
    36     y = r*np.outer(np.sin(u), np.ones_like(v)) 
     70    x = r*np.outer(cos(u), np.ones_like(v)) 
     71    y = r*np.outer(sin(u), np.ones_like(v)) 
    3772    z = 1.3*np.outer(np.ones_like(u), v) 
    3873 
     
    4277    points = Rz(phi)*Ry(theta)*points 
    4378    x, y, z = [v.reshape(shape) for v in points] 
    44  
    45     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
     79    axes.plot_surface(x, y, z, color='yellow', alpha=alpha) 
     80 
     81    # TODO: draw endcaps on beam 
     82    ## Drawing tiny balls on the end will work 
     83    #draw_sphere(axes, radius=0.02, center=(0, 0, 1.3), color='yellow', alpha=alpha) 
     84    #draw_sphere(axes, radius=0.02, center=(0, 0, -1.3), color='yellow', alpha=alpha) 
     85    ## The following does not work 
     86    #triangles = [(0, i+1, i+2) for i in range(steps-2)] 
     87    #x_cap, y_cap = x[:, 0], y[:, 0] 
     88    #for z_cap in z[:, 0], z[:, -1]: 
     89    #    axes.plot_trisurf(x_cap, y_cap, z_cap, triangles, 
     90    #                      color='yellow', alpha=alpha) 
     91 
    4692 
    4793def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1): 
    4894    """Draw an ellipsoid.""" 
    4995    a, b, c = size 
    50     u = np.linspace(0, 2 * np.pi, steps) 
    51     v = np.linspace(0, np.pi, steps) 
    52     x = a*np.outer(np.cos(u), np.sin(v)) 
    53     y = b*np.outer(np.sin(u), np.sin(v)) 
    54     z = c*np.outer(np.ones_like(u), np.cos(v)) 
     96    u = np.linspace(0, 2 * pi, steps) 
     97    v = np.linspace(0, pi, steps) 
     98    x = a*np.outer(cos(u), sin(v)) 
     99    y = b*np.outer(sin(u), sin(v)) 
     100    z = c*np.outer(np.ones_like(u), cos(v)) 
    55101    x, y, z = transform_xyz(view, jitter, x, y, z) 
    56102 
    57     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
     103    axes.plot_surface(x, y, z, color='w', alpha=alpha) 
    58104 
    59105    draw_labels(axes, view, jitter, [ 
     
    68114def draw_sc(axes, size, view, jitter, steps=None, alpha=1): 
    69115    """Draw points for simple cubic paracrystal""" 
     116    # pylint: disable=unused-argument 
    70117    atoms = _build_sc() 
    71118    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     
    73120def draw_fcc(axes, size, view, jitter, steps=None, alpha=1): 
    74121    """Draw points for face-centered cubic paracrystal""" 
     122    # pylint: disable=unused-argument 
    75123    # Build the simple cubic crystal 
    76124    atoms = _build_sc() 
     
    88136def draw_bcc(axes, size, view, jitter, steps=None, alpha=1): 
    89137    """Draw points for body-centered cubic paracrystal""" 
     138    # pylint: disable=unused-argument 
    90139    # Build the simple cubic crystal 
    91140    atoms = _build_sc() 
     
    124173    return atoms 
    125174 
    126 def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1): 
    127     """Draw a parallelepiped.""" 
     175def draw_box(axes, size, view): 
     176    """Draw a wireframe box at a particular view.""" 
     177    a, b, c = size 
     178    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     179    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     180    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
     181    x, y, z = transform_xyz(view, None, x, y, z) 
     182    def _draw(i, j): 
     183        axes.plot([x[i], x[j]], [y[i], y[j]], [z[i], z[j]], color='black') 
     184    _draw(0, 1) 
     185    _draw(0, 2) 
     186    _draw(0, 3) 
     187    _draw(7, 4) 
     188    _draw(7, 5) 
     189    _draw(7, 6) 
     190 
     191def draw_parallelepiped(axes, size, view, jitter, steps=None, 
     192                        color=(0.6, 1.0, 0.6), alpha=1): 
     193    """Draw a parallelepiped surface, with view and jitter.""" 
     194    # pylint: disable=unused-argument 
    128195    a, b, c = size 
    129196    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     
    142209 
    143210    x, y, z = transform_xyz(view, jitter, x, y, z) 
    144     axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    145  
    146     # Draw pink face on box. 
     211    axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha, 
     212                      linewidth=0) 
     213 
     214    # Colour the c+ face of the box. 
    147215    # Since I can't control face color, instead draw a thin box situated just 
    148216    # in front of the "c+" face.  Use the c face so that rotations about psi 
    149217    # rotate that face. 
    150     if 1: 
     218    if 0: # pylint: disable=using-constant-test 
     219        color = (1, 0.6, 0.6)  # pink 
    151220        x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
    152221        y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
    153222        z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
    154223        x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) 
    155         axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) 
     224        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha) 
    156225 
    157226    draw_labels(axes, view, jitter, [ 
     
    164233    ]) 
    165234 
    166 def draw_sphere(axes, radius=10., steps=100): 
     235def draw_sphere(axes, radius=1.0, steps=25, 
     236                center=(0, 0, 0), color='w', alpha=1.): 
    167237    """Draw a sphere""" 
    168     u = np.linspace(0, 2 * np.pi, steps) 
    169     v = np.linspace(0, np.pi, steps) 
    170  
    171     x = radius * np.outer(np.cos(u), np.sin(v)) 
    172     y = radius * np.outer(np.sin(u), np.sin(v)) 
    173     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    174     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    175  
    176 def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), 
    177                 draw_shape=draw_parallelepiped): 
     238    u = np.linspace(0, 2 * pi, steps) 
     239    v = np.linspace(0, pi, steps) 
     240 
     241    x = radius * np.outer(cos(u), sin(v)) + center[0] 
     242    y = radius * np.outer(sin(u), sin(v)) + center[1] 
     243    z = radius * np.outer(np.ones(np.size(u)), cos(v)) + center[2] 
     244    axes.plot_surface(x, y, z, color=color, alpha=alpha) 
     245    #axes.plot_wireframe(x, y, z) 
     246 
     247def draw_axes(axes, origin=(-1, -1, -1), length=(2, 2, 2)): 
     248    """Draw wireframe axes lines, with given origin and length""" 
     249    x, y, z = origin 
     250    dx, dy, dz = length 
     251    axes.plot([x, x+dx], [y, y], [z, z], color='black') 
     252    axes.plot([x, x], [y, y+dy], [z, z], color='black') 
     253    axes.plot([x, x], [y, y], [z, z+dz], color='black') 
     254 
     255def draw_person_on_sphere(axes, view, height=0.5, radius=1.0): 
     256    """ 
     257    Draw a person on the surface of a sphere. 
     258 
     259    *view* indicates (latitude, longitude, orientation) 
     260    """ 
     261    limb_offset = height * 0.05 
     262    head_radius = height * 0.10 
     263    head_height = height - head_radius 
     264    neck_length = head_radius * 0.50 
     265    shoulder_height = height - 2*head_radius - neck_length 
     266    torso_length = shoulder_height * 0.55 
     267    torso_radius = torso_length * 0.30 
     268    leg_length = shoulder_height - torso_length 
     269    arm_length = torso_length * 0.90 
     270 
     271    def _draw_part(y, z): 
     272        x = np.zeros_like(y) 
     273        xp, yp, zp = transform_xyz(view, None, x, y, z + radius) 
     274        axes.plot(xp, yp, zp, color='k') 
     275 
     276    # circle for head 
     277    u = np.linspace(0, 2 * pi, 40) 
     278    y = head_radius * cos(u) 
     279    z = head_radius * sin(u) + head_height 
     280    _draw_part(y, z) 
     281 
     282    # rectangle for body 
     283    y = np.array([-torso_radius, torso_radius, torso_radius, -torso_radius, -torso_radius]) 
     284    z = np.array([0., 0, torso_length, torso_length, 0]) + leg_length 
     285    _draw_part(y, z) 
     286 
     287    # arms 
     288    y = np.array([-torso_radius - limb_offset, -torso_radius - limb_offset, -torso_radius]) 
     289    z = np.array([shoulder_height - arm_length, shoulder_height, shoulder_height]) 
     290    _draw_part(y, z) 
     291    _draw_part(-y, z)  # pylint: disable=invalid-unary-operand-type 
     292 
     293    # legs 
     294    y = np.array([-torso_radius + limb_offset, -torso_radius + limb_offset]) 
     295    z = np.array([0, leg_length]) 
     296    _draw_part(y, z) 
     297    _draw_part(-y, z)  # pylint: disable=invalid-unary-operand-type 
     298 
     299    limits = [-radius-height, radius+height] 
     300    axes.set_xlim(limits) 
     301    axes.set_ylim(limits) 
     302    axes.set_zlim(limits) 
     303    axes.set_axis_off() 
     304 
     305def draw_jitter(axes, view, jitter, dist='gaussian', 
     306                size=(0.1, 0.4, 1.0), 
     307                draw_shape=draw_parallelepiped, 
     308                projection='equirectangular', 
     309                alpha=0.8, 
     310                views=None): 
    178311    """ 
    179312    Represent jitter as a set of shapes at different orientations. 
    180313    """ 
     314    project, project_weight = get_projection(projection) 
     315 
    181316    # set max diagonal to 0.95 
    182317    scale = 0.95/sqrt(sum(v**2 for v in size)) 
    183318    size = tuple(scale*v for v in size) 
    184319 
    185     #np.random.seed(10) 
    186     #cloud = np.random.randn(10,3) 
    187     cloud = [ 
    188         [-1, -1, -1], 
    189         [-1, -1, +0], 
    190         [-1, -1, +1], 
    191         [-1, +0, -1], 
    192         [-1, +0, +0], 
    193         [-1, +0, +1], 
    194         [-1, +1, -1], 
    195         [-1, +1, +0], 
    196         [-1, +1, +1], 
    197         [+0, -1, -1], 
    198         [+0, -1, +0], 
    199         [+0, -1, +1], 
    200         [+0, +0, -1], 
    201         [+0, +0, +0], 
    202         [+0, +0, +1], 
    203         [+0, +1, -1], 
    204         [+0, +1, +0], 
    205         [+0, +1, +1], 
    206         [+1, -1, -1], 
    207         [+1, -1, +0], 
    208         [+1, -1, +1], 
    209         [+1, +0, -1], 
    210         [+1, +0, +0], 
    211         [+1, +0, +1], 
    212         [+1, +1, -1], 
    213         [+1, +1, +0], 
    214         [+1, +1, +1], 
    215     ] 
    216320    dtheta, dphi, dpsi = jitter 
    217     if dtheta == 0: 
    218         cloud = [v for v in cloud if v[0] == 0] 
    219     if dphi == 0: 
    220         cloud = [v for v in cloud if v[1] == 0] 
    221     if dpsi == 0: 
    222         cloud = [v for v in cloud if v[2] == 0] 
    223     draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    224     scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    225     for point in cloud: 
    226         delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
    227         draw_shape(axes, size, view, delta, alpha=0.8) 
     321    base = {'gaussian':3, 'rectangle':sqrt(3), 'uniform':1}[dist] 
     322    def _steps(delta): 
     323        if views is None: 
     324            n = max(3, min(25, 2*int(base*delta/5))) 
     325        else: 
     326            n = views 
     327        return base*delta*np.linspace(-1, 1, n) if delta > 0 else [0.] 
     328    for theta in _steps(dtheta): 
     329        for phi in _steps(dphi): 
     330            for psi in _steps(dpsi): 
     331                w = project_weight(theta, phi, 1.0, 1.0) 
     332                if w > 0: 
     333                    dview = project(theta, phi, psi) 
     334                    draw_shape(axes, size, view, dview, alpha=alpha) 
    228335    for v in 'xyz': 
    229336        a, b, c = size 
    230         lim = np.sqrt(a**2 + b**2 + c**2) 
     337        lim = sqrt(a**2 + b**2 + c**2) 
    231338        getattr(axes, 'set_'+v+'lim')([-lim, lim]) 
    232         getattr(axes, v+'axis').label.set_text(v) 
     339        #getattr(axes, v+'axis').label.set_text(v) 
    233340 
    234341PROJECTIONS = [ 
     
    238345    'azimuthal_equal_area', 
    239346] 
    240 def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
    241               projection='equirectangular'): 
    242     """ 
    243     Draw the dispersion mesh showing the theta-phi orientations at which 
    244     the model will be evaluated. 
    245  
     347def get_projection(projection): 
     348 
     349    """ 
    246350    jitter projections 
    247351    <https://en.wikipedia.org/wiki/List_of_map_projections> 
     
    297401        Should allow free movement in theta, but phi is distorted. 
    298402    """ 
     403    # pylint: disable=unused-argument 
    299404    # TODO: try Kent distribution instead of a gaussian warped by projection 
    300405 
    301     dist_x = np.linspace(-1, 1, n) 
    302     weights = np.ones_like(dist_x) 
    303     if dist == 'gaussian': 
    304         dist_x *= 3 
    305         weights = exp(-0.5*dist_x**2) 
    306     elif dist == 'rectangle': 
    307         # Note: uses sasmodels ridiculous definition of rectangle width 
    308         dist_x *= sqrt(3) 
    309     elif dist == 'uniform': 
    310         pass 
    311     else: 
    312         raise ValueError("expected dist to be gaussian, rectangle or uniform") 
    313  
    314406    if projection == 'equirectangular':  #define PROJECTION 1 
    315         def _rotate(theta_i, phi_j): 
    316             return Rx(phi_j)*Ry(theta_i) 
     407        def _project(theta_i, phi_j, psi): 
     408            latitude, longitude = theta_i, phi_j 
     409            return latitude, longitude, psi, 'xyz' 
     410            #return Rx(phi_j)*Ry(theta_i) 
    317411        def _weight(theta_i, phi_j, w_i, w_j): 
    318412            return w_i*w_j*abs(cos(radians(theta_i))) 
    319413    elif projection == 'sinusoidal':  #define PROJECTION 2 
    320         def _rotate(theta_i, phi_j): 
     414        def _project(theta_i, phi_j, psi): 
    321415            latitude = theta_i 
    322416            scale = cos(radians(latitude)) 
    323417            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 
    324418            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    325             return Rx(longitude)*Ry(latitude) 
    326         def _weight(theta_i, phi_j, w_i, w_j): 
     419            return latitude, longitude, psi, 'xyz' 
     420            #return Rx(longitude)*Ry(latitude) 
     421        def _project(theta_i, phi_j, w_i, w_j): 
    327422            latitude = theta_i 
    328423            scale = cos(radians(latitude)) 
     
    330425            return active*w_i*w_j 
    331426    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
    332         def _rotate(theta_i, phi_j): 
     427        def _project(theta_i, phi_j, psi): 
    333428            from .guyou import guyou_invert 
    334429            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
    335430            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
    336             return Rx(longitude[0])*Ry(latitude[0]) 
     431            return latitude, longitude, psi, 'xyz' 
     432            #return Rx(longitude[0])*Ry(latitude[0]) 
    337433        def _weight(theta_i, phi_j, w_i, w_j): 
    338434            return w_i*w_j 
    339     elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
    340         def _rotate(theta_i, phi_j): 
     435    elif projection == 'azimuthal_equidistance': 
     436        # Note that calculates angles for Rz Ry rather than Rx Ry 
     437        def _project(theta_i, phi_j, psi): 
    341438            latitude = sqrt(theta_i**2 + phi_j**2) 
    342             longitude = degrees(np.arctan2(phi_j, theta_i)) 
     439            longitude = degrees(arctan2(phi_j, theta_i)) 
    343440            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    344             return Rz(longitude)*Ry(latitude) 
     441            return latitude, longitude, psi-longitude, 'zyz' 
     442            #R = Rz(longitude)*Ry(latitude)*Rz(psi) 
     443            #return R_to_xyz(R) 
     444            #return Rz(longitude)*Ry(latitude) 
    345445        def _weight(theta_i, phi_j, w_i, w_j): 
    346446            # Weighting for each point comes from the integral: 
     
    376476            return weight*w_i*w_j if latitude < 180 else 0 
    377477    elif projection == 'azimuthal_equal_area': 
    378         def _rotate(theta_i, phi_j): 
     478        # Note that calculates angles for Rz Ry rather than Rx Ry 
     479        def _project(theta_i, phi_j, psi): 
    379480            radius = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
    380             latitude = 180-degrees(2*np.arccos(radius)) 
    381             longitude = degrees(np.arctan2(phi_j, theta_i)) 
     481            latitude = 180-degrees(2*arccos(radius)) 
     482            longitude = degrees(arctan2(phi_j, theta_i)) 
    382483            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    383             return Rz(longitude)*Ry(latitude) 
     484            return latitude, longitude, psi, 'zyz' 
     485            #R = Rz(longitude)*Ry(latitude)*Rz(psi) 
     486            #return R_to_xyz(R) 
     487            #return Rz(longitude)*Ry(latitude) 
    384488        def _weight(theta_i, phi_j, w_i, w_j): 
    385489            latitude = sqrt(theta_i**2 + phi_j**2) 
     
    389493        raise ValueError("unknown projection %r"%projection) 
    390494 
     495    return _project, _weight 
     496 
     497def R_to_xyz(R): 
     498    """ 
     499    Return phi, theta, psi Tait-Bryan angles corresponding to the given rotation matrix. 
     500 
     501    Extracting Euler Angles from a Rotation Matrix 
     502    Mike Day, Insomniac Games 
     503    https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2012/07/euler-angles1.pdf 
     504    Based on: Shoemake’s "Euler Angle Conversion", Graphics Gems IV, pp.  222-229 
     505    """ 
     506    phi = arctan2(R[1, 2], R[2, 2]) 
     507    theta = arctan2(-R[0, 2], sqrt(R[0, 0]**2 + R[0, 1]**2)) 
     508    psi = arctan2(R[0, 1], R[0, 0]) 
     509    return degrees(phi), degrees(theta), degrees(psi) 
     510 
     511def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
     512              projection='equirectangular'): 
     513    """ 
     514    Draw the dispersion mesh showing the theta-phi orientations at which 
     515    the model will be evaluated. 
     516    """ 
     517 
     518    _project, _weight = get_projection(projection) 
     519    def _rotate(theta, phi, z): 
     520        dview = _project(theta, phi, 0.) 
     521        if dview[3] == 'zyz': 
     522            return Rz(dview[1])*Ry(dview[0])*z 
     523        else:  # dview[3] == 'xyz': 
     524            return Rx(dview[1])*Ry(dview[0])*z 
     525 
     526 
     527    dist_x = np.linspace(-1, 1, n) 
     528    weights = np.ones_like(dist_x) 
     529    if dist == 'gaussian': 
     530        dist_x *= 3 
     531        weights = exp(-0.5*dist_x**2) 
     532    elif dist == 'rectangle': 
     533        # Note: uses sasmodels ridiculous definition of rectangle width 
     534        dist_x *= sqrt(3) 
     535    elif dist == 'uniform': 
     536        pass 
     537    else: 
     538        raise ValueError("expected dist to be gaussian, rectangle or uniform") 
     539 
    391540    # mesh in theta, phi formed by rotating z 
    392     dtheta, dphi, dpsi = jitter 
     541    dtheta, dphi, dpsi = jitter  # pylint: disable=unused-variable 
    393542    z = np.matrix([[0], [0], [radius]]) 
    394     points = np.hstack([_rotate(theta_i, phi_j)*z 
     543    points = np.hstack([_rotate(theta_i, phi_j, z) 
    395544                        for theta_i in dtheta*dist_x 
    396545                        for phi_j in dphi*dist_x]) 
     
    470619    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
    471620    """ 
    472     dtheta, dphi, dpsi = jitter 
    473     points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     621    if jitter is None: 
     622        return points 
     623    # Hack to deal with the fact that azimuthal_equidistance uses euler angles 
     624    if len(jitter) == 4: 
     625        dtheta, dphi, dpsi, _ = jitter 
     626        points = Rz(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     627    else: 
     628        dtheta, dphi, dpsi = jitter 
     629        points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
    474630    return points 
    475631 
     
    481637    """ 
    482638    theta, phi, psi = view 
    483     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     639    points = Rz(phi)*Ry(theta)*Rz(psi)*points # viewing angle 
     640    #points = Rz(psi)*Ry(pi/2-theta)*Rz(phi)*points # 1-D integration angles 
     641    #points = Rx(phi)*Ry(theta)*Rz(psi)*points  # angular dispersion angle 
    484642    return points 
     643 
     644def orient_relative_to_beam_quaternion(view, points): 
     645    """ 
     646    Apply the view transform to a set of points. 
     647 
     648    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     649 
     650    This variant uses quaternions rather than rotation matrices for the 
     651    computation.  It works but it is not used because it doesn't solve 
     652    any problems.  The challenge of mapping theta/phi/psi to SO(3) does 
     653    not disappear by calculating the transform differently. 
     654    """ 
     655    theta, phi, psi = view 
     656    x, y, z = [1, 0, 0], [0, 1, 0], [0, 0, 1] 
     657    q = Quaternion(1, [0, 0, 0]) 
     658    ## Compose a rotation about the three axes by rotating 
     659    ## the unit vectors before applying the rotation. 
     660    #q = Quaternion.from_angle_axis(theta, q.rot(x)) * q 
     661    #q = Quaternion.from_angle_axis(phi, q.rot(y)) * q 
     662    #q = Quaternion.from_angle_axis(psi, q.rot(z)) * q 
     663    ## The above turns out to be equivalent to reversing 
     664    ## the order of application, so ignore it and use below. 
     665    q = q * Quaternion.from_angle_axis(theta, x) 
     666    q = q * Quaternion.from_angle_axis(phi, y) 
     667    q = q * Quaternion.from_angle_axis(psi, z) 
     668    ## Reverse the order by post-multiply rather than pre-multiply 
     669    #q = Quaternion.from_angle_axis(theta, x) * q 
     670    #q = Quaternion.from_angle_axis(phi, y) * q 
     671    #q = Quaternion.from_angle_axis(psi, z) * q 
     672    #print("axes psi", q.rot(np.matrix([x, y, z]).T)) 
     673    return q.rot(points) 
     674#orient_relative_to_beam = orient_relative_to_beam_quaternion 
     675 
     676# === Quaterion class definition === BEGIN 
     677# Simple stand-alone quaternion class 
     678 
     679# Note: this code works but isn't unused since quaternions didn't solve the 
     680# representation problem.  Leave it here in case we want to revisit this later. 
     681 
     682#import numpy as np 
     683class Quaternion(object): 
     684    r""" 
     685    Quaternion(w, r) = w + ir[0] + jr[1] + kr[2] 
     686 
     687    Quaternion.from_angle_axis(theta, r) for a rotation of angle theta about 
     688    an axis oriented toward the direction r.  This defines a unit quaternion, 
     689    normalizing $r$ to the unit vector $\hat r$, and setting quaternion 
     690    $Q = \cos \theta + \sin \theta \hat r$ 
     691 
     692    Quaternion objects can be multiplied, which applies a rotation about the 
     693    given axis, allowing composition of rotations without risk of gimbal lock. 
     694    The resulting quaternion is applied to a set of points using *Q.rot(v)*. 
     695    """ 
     696    def __init__(self, w, r): 
     697        self.w = w 
     698        self.r = np.asarray(r, dtype='d') 
     699 
     700    @staticmethod 
     701    def from_angle_axis(theta, r): 
     702        """Build quaternion as rotation theta about axis r""" 
     703        theta = np.radians(theta)/2 
     704        r = np.asarray(r) 
     705        w = np.cos(theta) 
     706        r = np.sin(theta)*r/np.dot(r, r) 
     707        return Quaternion(w, r) 
     708 
     709    def __mul__(self, other): 
     710        """Multiply quaterions""" 
     711        if isinstance(other, Quaternion): 
     712            w = self.w*other.w - np.dot(self.r, other.r) 
     713            r = self.w*other.r + other.w*self.r + np.cross(self.r, other.r) 
     714            return Quaternion(w, r) 
     715        raise NotImplementedError("Quaternion * non-quaternion not implemented") 
     716 
     717    def rot(self, v): 
     718        """Transform point *v* by quaternion""" 
     719        v = np.asarray(v).T 
     720        use_transpose = (v.shape[-1] != 3) 
     721        if use_transpose: 
     722            v = v.T 
     723        v = v + np.cross(2*self.r, np.cross(self.r, v) + self.w*v) 
     724        #v = v + 2*self.w*np.cross(self.r, v) + np.cross(2*self.r, np.cross(self.r, v)) 
     725        if use_transpose: 
     726            v = v.T 
     727        return v.T 
     728 
     729    def conj(self): 
     730        """Conjugate quaternion""" 
     731        return Quaternion(self.w, -self.r) 
     732 
     733    def inv(self): 
     734        """Inverse quaternion""" 
     735        return self.conj()/self.norm()**2 
     736 
     737    def norm(self): 
     738        """Quaternion length""" 
     739        return np.sqrt(self.w**2 + np.sum(self.r**2)) 
     740 
     741    def __str__(self): 
     742        return "%g%+gi%+gj%+gk"%(self.w, self.r[0], self.r[1], self.r[2]) 
     743 
     744def test_qrot(): 
     745    """Quaternion checks""" 
     746    # Define rotation of 60 degrees around an axis in y-z that is 60 degrees 
     747    # from y.  The rotation axis is determined by rotating the point [0, 1, 0] 
     748    # about x. 
     749    ax = Quaternion.from_angle_axis(60, [1, 0, 0]).rot([0, 1, 0]) 
     750    q = Quaternion.from_angle_axis(60, ax) 
     751    # Set the point to be rotated, and its expected rotated position. 
     752    p = [1, -1, 2] 
     753    target = [(10+4*np.sqrt(3))/8, (1+2*np.sqrt(3))/8, (14-3*np.sqrt(3))/8] 
     754    #print(q, q.rot(p) - target) 
     755    assert max(abs(q.rot(p) - target)) < 1e-14 
     756#test_qrot() 
     757#import sys; sys.exit() 
     758# === Quaterion class definition === END 
    485759 
    486760# translate between number of dimension of dispersity and the number of 
     
    504778    or the top of the range, depending on whether *mode* is 'central' or 'top'. 
    505779    """ 
    506     if portion == 1.0: 
    507         return data.min(), data.max() 
    508     elif mode == 'central': 
    509         data = np.sort(data.flatten()) 
    510         offset = int(portion*len(data)/2 + 0.5) 
    511         return data[offset], data[-offset] 
    512     elif mode == 'top': 
    513         data = np.sort(data.flatten()) 
    514         offset = int(portion*len(data) + 0.5) 
    515         return data[offset], data[-1] 
     780    if portion < 1.0: 
     781        if mode == 'central': 
     782            data = np.sort(data.flatten()) 
     783            offset = int(portion*len(data)/2 + 0.5) 
     784            return data[offset], data[-offset] 
     785        if mode == 'top': 
     786            data = np.sort(data.flatten()) 
     787            offset = int(portion*len(data) + 0.5) 
     788            return data[offset], data[-1] 
     789    # Default: full range 
     790    return data.min(), data.max() 
    516791 
    517792def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): 
     
    544819 
    545820    # compute the pattern 
    546     qx, qy = calculator._data.x_bins, calculator._data.y_bins 
     821    qx, qy = calculator.qxqy 
    547822    Iqxy = calculator(**pars).reshape(len(qx), len(qy)) 
    548823 
    549824    # scale it and draw it 
    550     Iqxy = np.log(Iqxy) 
     825    Iqxy = log(Iqxy) 
    551826    if calculator.limits: 
    552827        # use limits from orientation (0,0,0) 
     
    556831        vmin = vmax*10**-7 
    557832        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top') 
     833    #vmin, vmax = Iqxy.min(), Iqxy.max() 
    558834    #print("range",(vmin,vmax)) 
    559835    #qx, qy = np.meshgrid(qx, qy) 
    560     if 0: 
     836    if 0:  # pylint: disable=using-constant-test 
    561837        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
    562838        level[level < 0] = 0 
     839        from matplotlib import pylab as plt 
    563840        colors = plt.get_cmap()(level) 
    564         axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
    565     elif 1: 
     841        #from matplotlib import cm 
     842        #colors = cm.coolwarm(level) 
     843        #colors = cm.gist_yarg(level) 
     844        #colors = cm.Wistia(level) 
     845        colors[level <= 0, 3] = 0.  # set floor to transparent 
     846        x, y = np.meshgrid(qx/qx.max(), qy/qy.max()) 
     847        axes.plot_surface(x, y, -1.1*np.ones_like(x), facecolors=colors) 
     848    elif 1:  # pylint: disable=using-constant-test 
    566849        axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
    567850                      levels=np.linspace(vmin, vmax, 24)) 
     
    594877    calculator = DirectModel(data, model) 
    595878 
     879    # Remember the data axes so we can plot the results 
     880    calculator.qxqy = (q, q) 
     881 
    596882    # stuff the values for non-orientation parameters into the calculator 
    597883    calculator.pars = pars.copy() 
     
    601887    # under rotation or angular dispersion 
    602888    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) 
    603     Iqxy = np.log(Iqxy) 
     889    Iqxy = log(Iqxy) 
    604890    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') 
    605891    calculator.limits = vmin, vmax+1 
     
    692978} 
    693979 
     980 
    694981def run(model_name='parallelepiped', size=(10, 40, 100), 
     982        view=(0, 0, 0), jitter=(0, 0, 0), 
    695983        dist='gaussian', mesh=30, 
    696984        projection='equirectangular'): 
     
    702990 
    703991    *size* gives the dimensions (a, b, c) of the shape. 
     992 
     993    *view* gives the initial view (theta, phi, psi) of the shape. 
     994 
     995    *view* gives the initial jitter (dtheta, dphi, dpsi) of the shape. 
    704996 
    705997    *dist* is the type of dispersition: gaussian, rectangle, or uniform. 
     
    7211013    calculator, size = select_calculator(model_name, n=150, size=size) 
    7221014    draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped) 
     1015    #draw_shape = draw_fcc 
    7231016 
    7241017    ## uncomment to set an independent the colour range for every view 
     
    7261019    calculator.limits = None 
    7271020 
    728     ## initial view 
    729     #theta, dtheta = 70., 10. 
    730     #phi, dphi = -45., 3. 
    731     #psi, dpsi = -45., 3. 
    732     theta, phi, psi = 0, 0, 0 
    733     dtheta, dphi, dpsi = 0, 0, 0 
     1021    PLOT_ENGINE(calculator, draw_shape, size, view, jitter, dist, mesh, projection) 
     1022 
     1023def _mpl_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): 
     1024    # Note: travis-ci does not support mpl_toolkits.mplot3d, but this shouldn't be 
     1025    # an issue since we are lazy-loading the package on a path that isn't tested. 
     1026    # Importing mplot3d adds projection='3d' option to subplot 
     1027    import mpl_toolkits.mplot3d  # pylint: disable=unused-variable 
     1028    import matplotlib as mpl 
     1029    import matplotlib.pyplot as plt 
     1030    from matplotlib.widgets import Slider 
    7341031 
    7351032    ## create the plot window 
     
    7521049 
    7531050    ## add control widgets to plot 
    754     axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props) 
    755     axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props) 
    756     axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props) 
    757     stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
    758     sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
    759     spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
    760  
    761     axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props) 
    762     axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props) 
    763     axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props) 
     1051    axes_theta = plt.axes([0.05, 0.15, 0.50, 0.04], **props) 
     1052    axes_phi = plt.axes([0.05, 0.10, 0.50, 0.04], **props) 
     1053    axes_psi = plt.axes([0.05, 0.05, 0.50, 0.04], **props) 
     1054    stheta = Slider(axes_theta, u'Ξ', -90, 90, valinit=0) 
     1055    sphi = Slider(axes_phi, u'φ', -180, 180, valinit=0) 
     1056    spsi = Slider(axes_psi, u'ψ', -180, 180, valinit=0) 
     1057 
     1058    axes_dtheta = plt.axes([0.70, 0.15, 0.20, 0.04], **props) 
     1059    axes_dphi = plt.axes([0.70, 0.1, 0.20, 0.04], **props) 
     1060    axes_dpsi = plt.axes([0.70, 0.05, 0.20, 0.04], **props) 
     1061 
    7641062    # Note: using ridiculous definition of rectangle distribution, whose width 
    7651063    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
    7661064    # the maximum width to 90. 
    7671065    dlimit = DIST_LIMITS[dist] 
    768     sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
    769     sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
    770     sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
    771  
     1066    sdtheta = Slider(axes_dtheta, u'Δξ', 0, 2*dlimit, valinit=0) 
     1067    sdphi = Slider(axes_dphi, u'Δφ', 0, 2*dlimit, valinit=0) 
     1068    sdpsi = Slider(axes_dpsi, u'Δψ', 0, 2*dlimit, valinit=0) 
     1069 
     1070    ## initial view and jitter 
     1071    theta, phi, psi = view 
     1072    stheta.set_val(theta) 
     1073    sphi.set_val(phi) 
     1074    spsi.set_val(psi) 
     1075    dtheta, dphi, dpsi = jitter 
     1076    sdtheta.set_val(dtheta) 
     1077    sdphi.set_val(dphi) 
     1078    sdpsi.set_val(dpsi) 
    7721079 
    7731080    ## callback to draw the new view 
    774     def update(val, axis=None): 
     1081    def _update(val, axis=None): 
     1082        # pylint: disable=unused-argument 
    7751083        view = stheta.val, sphi.val, spsi.val 
    7761084        jitter = sdtheta.val, sdphi.val, sdpsi.val 
    7771085        # set small jitter as 0 if multiple pd dims 
    7781086        dims = sum(v > 0 for v in jitter) 
    779         limit = [0, 0.5, 5][dims] 
     1087        limit = [0, 0.5, 5, 5][dims] 
    7801088        jitter = [0 if v < limit else v for v in jitter] 
    7811089        axes.cla() 
    782         draw_beam(axes, (0, 0)) 
    783         draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape) 
    784         #draw_jitter(axes, view, (0,0,0)) 
     1090 
     1091        ## Visualize as person on globe 
     1092        #draw_sphere(axes, radius=0.5) 
     1093        #draw_person_on_sphere(axes, view, radius=0.5) 
     1094 
     1095        ## Move beam instead of shape 
     1096        #draw_beam(axes, -view[:2]) 
     1097        #draw_jitter(axes, (0,0,0), (0,0,0), views=3) 
     1098 
     1099        ## Move shape and draw scattering 
     1100        draw_beam(axes, (0, 0), alpha=1.) 
     1101        #draw_person_on_sphere(axes, view, radius=1.2, height=0.5) 
     1102        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1., 
     1103                    draw_shape=draw_shape, projection=projection, views=3) 
    7851104        draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection) 
    7861105        draw_scattering(calculator, axes, view, jitter, dist=dist) 
     1106 
    7871107        plt.gcf().canvas.draw() 
    7881108 
    7891109    ## bind control widgets to view updater 
    790     stheta.on_changed(lambda v: update(v, 'theta')) 
    791     sphi.on_changed(lambda v: update(v, 'phi')) 
    792     spsi.on_changed(lambda v: update(v, 'psi')) 
    793     sdtheta.on_changed(lambda v: update(v, 'dtheta')) 
    794     sdphi.on_changed(lambda v: update(v, 'dphi')) 
    795     sdpsi.on_changed(lambda v: update(v, 'dpsi')) 
     1110    stheta.on_changed(lambda v: _update(v, 'theta')) 
     1111    sphi.on_changed(lambda v: _update(v, 'phi')) 
     1112    spsi.on_changed(lambda v: _update(v, 'psi')) 
     1113    sdtheta.on_changed(lambda v: _update(v, 'dtheta')) 
     1114    sdphi.on_changed(lambda v: _update(v, 'dphi')) 
     1115    sdpsi.on_changed(lambda v: _update(v, 'dpsi')) 
    7961116 
    7971117    ## initialize view 
    798     update(None, 'phi') 
     1118    _update(None, 'phi') 
    7991119 
    8001120    ## go interactive 
    8011121    plt.show() 
    8021122 
     1123 
     1124def map_colors(z, kw): 
     1125    """ 
     1126    Process matplotlib-style colour arguments. 
     1127 
     1128    Pulls 'cmap', 'alpha', 'vmin', and 'vmax' from th *kw* dictionary, setting 
     1129    the *kw['color']* to an RGB array.  These are ignored if 'c' or 'color' are 
     1130    set inside *kw*. 
     1131    """ 
     1132    from matplotlib import cm 
     1133 
     1134    cmap = kw.pop('cmap', cm.coolwarm) 
     1135    alpha = kw.pop('alpha', None) 
     1136    vmin = kw.pop('vmin', z.min()) 
     1137    vmax = kw.pop('vmax', z.max()) 
     1138    c = kw.pop('c', None) 
     1139    color = kw.pop('color', c) 
     1140    if color is None: 
     1141        znorm = ((z - vmin) / (vmax - vmin)).clip(0, 1) 
     1142        color = cmap(znorm) 
     1143    elif isinstance(color, np.ndarray) and color.shape == z.shape: 
     1144        color = cmap(color) 
     1145    if alpha is None: 
     1146        if isinstance(color, np.ndarray): 
     1147            color = color[..., :3] 
     1148    else: 
     1149        color[..., 3] = alpha 
     1150    kw['color'] = color 
     1151 
     1152def make_vec(*args): 
     1153    """Turn all elements of *args* into numpy arrays""" 
     1154    #return [np.asarray(v, 'd').flatten() for v in args] 
     1155    return [np.asarray(v, 'd') for v in args] 
     1156 
     1157def make_image(z, kw): 
     1158    """Convert numpy array *z* into a *PIL* RGB image.""" 
     1159    import PIL.Image 
     1160    from matplotlib import cm 
     1161 
     1162    cmap = kw.pop('cmap', cm.coolwarm) 
     1163 
     1164    znorm = (z-z.min())/z.ptp() 
     1165    c = cmap(znorm) 
     1166    c = c[..., :3] 
     1167    rgb = np.asarray(c*255, 'u1') 
     1168    image = PIL.Image.fromarray(rgb, mode='RGB') 
     1169    return image 
     1170 
     1171 
     1172_IPV_MARKERS = { 
     1173    'o': 'sphere', 
     1174} 
     1175_IPV_COLORS = { 
     1176    'w': 'white', 
     1177    'k': 'black', 
     1178    'c': 'cyan', 
     1179    'm': 'magenta', 
     1180    'y': 'yellow', 
     1181    'r': 'red', 
     1182    'g': 'green', 
     1183    'b': 'blue', 
     1184} 
     1185def _ipv_fix_color(kw): 
     1186    alpha = kw.pop('alpha', None) 
     1187    color = kw.get('color', None) 
     1188    if isinstance(color, str): 
     1189        color = _IPV_COLORS.get(color, color) 
     1190        kw['color'] = color 
     1191    if alpha is not None: 
     1192        color = kw['color'] 
     1193        #TODO: convert color to [r, g, b, a] if not already 
     1194        if isinstance(color, (tuple, list)): 
     1195            if len(color) == 3: 
     1196                color = (color[0], color[1], color[2], alpha) 
     1197            else: 
     1198                color = (color[0], color[1], color[2], alpha*color[3]) 
     1199            color = np.array(color) 
     1200        if isinstance(color, np.ndarray) and color.shape[-1] == 4: 
     1201            color[..., 3] = alpha 
     1202            kw['color'] = color 
     1203 
     1204def _ipv_set_transparency(kw, obj): 
     1205    color = kw.get('color', None) 
     1206    if (isinstance(color, np.ndarray) 
     1207            and color.shape[-1] == 4 
     1208            and (color[..., 3] != 1.0).any()): 
     1209        obj.material.transparent = True 
     1210        obj.material.side = "FrontSide" 
     1211 
     1212def ipv_axes(): 
     1213    """ 
     1214    Build a matplotlib style Axes interface for ipyvolume 
     1215    """ 
     1216    import ipyvolume as ipv 
     1217 
     1218    class Axes(object): 
     1219        """ 
     1220        Matplotlib Axes3D style interface to ipyvolume renderer. 
     1221        """ 
     1222        # pylint: disable=no-self-use,no-init 
     1223        # transparency can be achieved by setting the following: 
     1224        #    mesh.color = [r, g, b, alpha] 
     1225        #    mesh.material.transparent = True 
     1226        #    mesh.material.side = "FrontSide" 
     1227        # smooth(ish) rotation can be achieved by setting: 
     1228        #    slide.continuous_update = True 
     1229        #    figure.animation = 0. 
     1230        #    mesh.material.x = x 
     1231        #    mesh.material.y = y 
     1232        #    mesh.material.z = z 
     1233        # maybe need to synchronize update of x/y/z to avoid shimmy when moving 
     1234        def plot(self, x, y, z, **kw): 
     1235            """mpl style plot interface for ipyvolume""" 
     1236            _ipv_fix_color(kw) 
     1237            x, y, z = make_vec(x, y, z) 
     1238            ipv.plot(x, y, z, **kw) 
     1239        def plot_surface(self, x, y, z, **kw): 
     1240            """mpl style plot_surface interface for ipyvolume""" 
     1241            facecolors = kw.pop('facecolors', None) 
     1242            if facecolors is not None: 
     1243                kw['color'] = facecolors 
     1244            _ipv_fix_color(kw) 
     1245            x, y, z = make_vec(x, y, z) 
     1246            h = ipv.plot_surface(x, y, z, **kw) 
     1247            _ipv_set_transparency(kw, h) 
     1248            #h.material.side = "DoubleSide" 
     1249            return h 
     1250        def plot_trisurf(self, x, y, triangles=None, Z=None, **kw): 
     1251            """mpl style plot_trisurf interface for ipyvolume""" 
     1252            kw.pop('linewidth', None) 
     1253            _ipv_fix_color(kw) 
     1254            x, y, z = make_vec(x, y, Z) 
     1255            if triangles is not None: 
     1256                triangles = np.asarray(triangles) 
     1257            h = ipv.plot_trisurf(x, y, z, triangles=triangles, **kw) 
     1258            _ipv_set_transparency(kw, h) 
     1259            return h 
     1260        def scatter(self, x, y, z, **kw): 
     1261            """mpl style scatter interface for ipyvolume""" 
     1262            x, y, z = make_vec(x, y, z) 
     1263            map_colors(z, kw) 
     1264            marker = kw.get('marker', None) 
     1265            kw['marker'] = _IPV_MARKERS.get(marker, marker) 
     1266            h = ipv.scatter(x, y, z, **kw) 
     1267            _ipv_set_transparency(kw, h) 
     1268            return h 
     1269        def contourf(self, x, y, v, zdir='z', offset=0, levels=None, **kw): 
     1270            """mpl style contourf interface for ipyvolume""" 
     1271            # pylint: disable=unused-argument 
     1272            # Don't use contour for now (although we might want to later) 
     1273            self.pcolor(x, y, v, zdir='z', offset=offset, **kw) 
     1274        def pcolor(self, x, y, v, zdir='z', offset=0, **kw): 
     1275            """mpl style pcolor interface for ipyvolume""" 
     1276            # pylint: disable=unused-argument 
     1277            x, y, v = make_vec(x, y, v) 
     1278            image = make_image(v, kw) 
     1279            xmin, xmax = x.min(), x.max() 
     1280            ymin, ymax = y.min(), y.max() 
     1281            x = np.array([[xmin, xmax], [xmin, xmax]]) 
     1282            y = np.array([[ymin, ymin], [ymax, ymax]]) 
     1283            z = x*0 + offset 
     1284            u = np.array([[0., 1], [0, 1]]) 
     1285            v = np.array([[0., 0], [1, 1]]) 
     1286            h = ipv.plot_mesh(x, y, z, u=u, v=v, texture=image, wireframe=False) 
     1287            _ipv_set_transparency(kw, h) 
     1288            h.material.side = "DoubleSide" 
     1289            return h 
     1290        def text(self, *args, **kw): 
     1291            """mpl style text interface for ipyvolume""" 
     1292            pass 
     1293        def set_xlim(self, limits): 
     1294            """mpl style set_xlim interface for ipyvolume""" 
     1295            ipv.xlim(*limits) 
     1296        def set_ylim(self, limits): 
     1297            """mpl style set_ylim interface for ipyvolume""" 
     1298            ipv.ylim(*limits) 
     1299        def set_zlim(self, limits): 
     1300            """mpl style set_zlim interface for ipyvolume""" 
     1301            ipv.zlim(*limits) 
     1302        def set_axes_on(self): 
     1303            """mpl style set_axes_on interface for ipyvolume""" 
     1304            ipv.style.axis_on() 
     1305        def set_axis_off(self): 
     1306            """mpl style set_axes_off interface for ipyvolume""" 
     1307            ipv.style.axes_off() 
     1308    return Axes() 
     1309 
     1310def _ipv_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): 
     1311    from IPython.display import display 
     1312    import ipywidgets as widgets 
     1313    import ipyvolume as ipv 
     1314 
     1315    axes = ipv_axes() 
     1316 
     1317    def _draw(view, jitter): 
     1318        camera = ipv.gcf().camera 
     1319        #print(ipv.gcf().__dict__.keys()) 
     1320        #print(dir(ipv.gcf())) 
     1321        ipv.figure(animation=0.)  # no animation when updating object mesh 
     1322 
     1323        # set small jitter as 0 if multiple pd dims 
     1324        dims = sum(v > 0 for v in jitter) 
     1325        limit = [0, 0.5, 5, 5][dims] 
     1326        jitter = [0 if v < limit else v for v in jitter] 
     1327 
     1328        ## Visualize as person on globe 
     1329        #draw_beam(axes, (0, 0)) 
     1330        #draw_sphere(axes, radius=0.5) 
     1331        #draw_person_on_sphere(axes, view, radius=0.5) 
     1332 
     1333        ## Move beam instead of shape 
     1334        #draw_beam(axes, view=(-view[0], -view[1])) 
     1335        #draw_jitter(axes, view=(0,0,0), jitter=(0,0,0)) 
     1336 
     1337        ## Move shape and draw scattering 
     1338        draw_beam(axes, (0, 0), steps=25) 
     1339        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.0, 
     1340                    draw_shape=draw_shape, projection=projection) 
     1341        draw_mesh(axes, view, jitter, dist=dist, n=mesh, radius=0.95, 
     1342                  projection=projection) 
     1343        draw_scattering(calculator, axes, view, jitter, dist=dist) 
     1344 
     1345        draw_axes(axes, origin=(-1, -1, -1.1)) 
     1346        ipv.style.box_off() 
     1347        ipv.style.axes_off() 
     1348        ipv.xyzlabel(" ", " ", " ") 
     1349 
     1350        ipv.gcf().camera = camera 
     1351        ipv.show() 
     1352 
     1353 
     1354    trange, prange = (-180., 180., 1.), (-180., 180., 1.) 
     1355    dtrange, dprange = (0., 180., 1.), (0., 180., 1.) 
     1356 
     1357    ## Super simple interfaca, but uses non-ascii variable namese 
     1358    # Ξ φ ψ Δξ Δφ Δψ 
     1359    #def update(**kw): 
     1360    #    view = kw['Ξ'], kw['φ'], kw['ψ'] 
     1361    #    jitter = kw['Δξ'], kw['Δφ'], kw['Δψ'] 
     1362    #    draw(view, jitter) 
     1363    #widgets.interact(update, Ξ=trange, φ=prange, ψ=prange, Δξ=dtrange, Δφ=dprange, Δψ=dprange) 
     1364 
     1365    def _update(theta, phi, psi, dtheta, dphi, dpsi): 
     1366        _draw(view=(theta, phi, psi), jitter=(dtheta, dphi, dpsi)) 
     1367 
     1368    def _slider(name, slice, init=0.): 
     1369        return widgets.FloatSlider( 
     1370            value=init, 
     1371            min=slice[0], 
     1372            max=slice[1], 
     1373            step=slice[2], 
     1374            description=name, 
     1375            disabled=False, 
     1376            #continuous_update=True, 
     1377            continuous_update=False, 
     1378            orientation='horizontal', 
     1379            readout=True, 
     1380            readout_format='.1f', 
     1381            ) 
     1382    theta = _slider(u'Ξ', trange, view[0]) 
     1383    phi = _slider(u'φ', prange, view[1]) 
     1384    psi = _slider(u'ψ', prange, view[2]) 
     1385    dtheta = _slider(u'Δξ', dtrange, jitter[0]) 
     1386    dphi = _slider(u'Δφ', dprange, jitter[1]) 
     1387    dpsi = _slider(u'Δψ', dprange, jitter[2]) 
     1388    fields = { 
     1389        'theta': theta, 'phi': phi, 'psi': psi, 
     1390        'dtheta': dtheta, 'dphi': dphi, 'dpsi': dpsi, 
     1391    } 
     1392    ui = widgets.HBox([ 
     1393        widgets.VBox([theta, phi, psi]), 
     1394        widgets.VBox([dtheta, dphi, dpsi]) 
     1395    ]) 
     1396 
     1397    out = widgets.interactive_output(_update, fields) 
     1398    display(ui, out) 
     1399 
     1400 
     1401_ENGINES = { 
     1402    "matplotlib": _mpl_plot, 
     1403    "mpl": _mpl_plot, 
     1404    #"plotly": _plotly_plot, 
     1405    "ipvolume": _ipv_plot, 
     1406    "ipv": _ipv_plot, 
     1407} 
     1408PLOT_ENGINE = _ENGINES["matplotlib"] 
     1409def set_plotter(name): 
     1410    """ 
     1411    Setting the plotting engine to matplotlib/ipyvolume or equivalently mpl/ipv. 
     1412    """ 
     1413    global PLOT_ENGINE 
     1414    PLOT_ENGINE = _ENGINES[name] 
     1415 
    8031416def main(): 
     1417    """ 
     1418    Command line interface to the jitter viewer. 
     1419    """ 
    8041420    parser = argparse.ArgumentParser( 
    8051421        description="Display jitter", 
     
    8111427    parser.add_argument('-s', '--size', type=str, default='10,40,100', 
    8121428                        help='a,b,c lengths') 
     1429    parser.add_argument('-v', '--view', type=str, default='0,0,0', 
     1430                        help='initial view angles') 
     1431    parser.add_argument('-j', '--jitter', type=str, default='0,0,0', 
     1432                        help='initial angular dispersion') 
    8131433    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, 
    8141434                        default=DISTRIBUTIONS[0], 
     
    8191439                        help='oriented shape') 
    8201440    opts = parser.parse_args() 
    821     size = tuple(int(v) for v in opts.size.split(',')) 
    822     run(opts.shape, size=size, 
     1441    size = tuple(float(v) for v in opts.size.split(',')) 
     1442    view = tuple(float(v) for v in opts.view.split(',')) 
     1443    jitter = tuple(float(v) for v in opts.jitter.split(',')) 
     1444    run(opts.shape, size=size, view=view, jitter=jitter, 
    8231445        mesh=opts.mesh, dist=opts.distribution, 
    8241446        projection=opts.projection) 
  • sasmodels/kernel.py

    r36a2418 ra34b811  
    2525 
    2626class KernelModel(object): 
     27    """ 
     28    Model definition for the compute engine. 
     29    """ 
    2730    info = None  # type: ModelInfo 
    2831    dtype = None # type: np.dtype 
    2932    def make_kernel(self, q_vectors): 
    3033        # type: (List[np.ndarray]) -> "Kernel" 
     34        """ 
     35        Instantiate a kernel for evaluating the model at *q_vectors*. 
     36        """ 
    3137        raise NotImplementedError("need to implement make_kernel") 
    3238 
    3339    def release(self): 
    3440        # type: () -> None 
     41        """ 
     42        Free resources associated with the kernel. 
     43        """ 
    3544        pass 
    3645 
    3746 
    3847class Kernel(object): 
    39     #: kernel dimension, either "1d" or "2d" 
     48    """ 
     49    Instantiated model for the compute engine, applied to a particular *q*. 
     50 
     51    Subclasses should define :meth:`_call_kernel` to evaluate the kernel over 
     52    its inputs. 
     53    """ 
     54    #: Kernel dimension, either "1d" or "2d". 
    4055    dim = None  # type: str 
     56    #: Model info. 
    4157    info = None  # type: ModelInfo 
    42     results = None # type: List[np.ndarray] 
     58    #: Numerical precision for the computation. 
    4359    dtype = None  # type: np.dtype 
     60    #: q values at which the kernel is to be evaluated 
     61    q_input = None  # type: Any 
     62    #: Place to hold result of :meth:`_call_kernel` for subclass. 
     63    result = None # type: np.ndarray 
    4464 
    4565    def Iq(self, call_details, values, cutoff, magnetic): 
     
    5777        built into the model, and do not need an additional scale. 
    5878        """ 
    59         _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, 
    60                               effective_radius_type=0) 
     79        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, 
     80                                            magnetic, radius_effective_mode=0) 
    6181        combined_scale = values[0]/shell_volume 
    6282        background = values[1] 
     
    6484    __call__ = Iq 
    6585 
    66     def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0): 
     86    def Fq(self, call_details, values, cutoff, magnetic, 
     87           radius_effective_mode=0): 
    6788        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    6889        r""" 
     
    81102        .. math:: 
    82103 
    83             I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
     104            I(q) = \text{scale} P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
    84105                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background} 
    85106 
     
    90111        .. math:: 
    91112 
    92             P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
     113            P(q)=\frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
    93114 
    94115        The form factor itself is scaled by volume and contrast to compute the 
     
    122143        volume fraction of the particles.  The model can have several 
    123144        different ways to compute effective radius, with the 
    124         *effective_radius_type* parameter used to select amongst them.  The 
     145        *radius_effective_mode* parameter used to select amongst them.  The 
    125146        volume fraction of particles should be determined from the total 
    126147        volume fraction of the form, not just the shell volume fraction. 
     
    131152        hollow and solid shapes. 
    132153        """ 
    133         self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     154        self._call_kernel(call_details, values, cutoff, magnetic, 
     155                          radius_effective_mode) 
    134156        #print("returned",self.q_input.q, self.result) 
    135157        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     
    143165        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight 
    144166        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight 
    145         effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight 
     167        radius_effective = self.result[nout*self.q_input.nq + 3]/total_weight 
    146168        if shell_volume == 0.: 
    147169            shell_volume = 1. 
    148         F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None 
     170        F1 = (self.result[1:nout*self.q_input.nq:nout]/total_weight 
     171              if nout == 2 else None) 
    149172        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight 
    150         return F1, F2, effective_radius, shell_volume, form_volume/shell_volume 
     173        return F1, F2, radius_effective, shell_volume, form_volume/shell_volume 
    151174 
    152175    def release(self): 
    153176        # type: () -> None 
     177        """ 
     178        Free resources associated with the kernel instance. 
     179        """ 
    154180        pass 
    155181 
    156     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     182    def _call_kernel(self, call_details, values, cutoff, magnetic, 
     183                     radius_effective_mode): 
    157184        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    158185        """ 
  • sasmodels/kernel_iq.c

    r12f4c19 rde032da  
    2727//      parameters in the parameter table. 
    2828//  CALL_VOLUME(form, shell, table) : assign form and shell values 
    29 //  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 
     29//  CALL_RADIUS_EFFECTIVE(mode, table) : call the R_eff function 
    3030//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3131//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
     
    8585static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 
    8686{ 
     87 
     88  double norm; 
    8789  in_spin = clip(in_spin, 0.0, 1.0); 
    8890  out_spin = clip(out_spin, 0.0, 1.0); 
     
    9496  // However, since the weights are applied to the final intensity and 
    9597  // are not interned inside the I(q) function, we want the full 
    96   // weight and not the square root.  Any function using 
    97   // set_spin_weights as part of calculating an amplitude will need to 
    98   // manually take that square root, but there is currently no such 
    99   // function. 
    100   weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd 
    101   weight[1] = (1.0-in_spin) * out_spin;       // du 
    102   weight[2] = in_spin * (1.0-out_spin);       // ud 
    103   weight[3] = in_spin * out_spin;             // uu 
     98  // weight and not the square root.  Anyway no function will ever use 
     99  // set_spin_weights as part of calculating an amplitude, as the weights are 
     100  // related to polarisation efficiency of the instrument. The weights serve to 
     101  // construct various magnet scattering cross sections, which are linear combinations 
     102  // of the spin-resolved cross sections. The polarisation efficiency e_in and e_out 
     103  // are parameters ranging from 0.5 (unpolarised) beam to 1 (perfect optics). 
     104  // For in_spin or out_spin <0.5 one assumes a CS, where the spin is reversed/flipped 
     105  // with respect to the initial supermirror polariser. The actual polarisation efficiency 
     106  // in this case is however e_in/out = 1-in/out_spin. 
     107 
     108  if (out_spin < 0.5){norm=1-out_spin;} 
     109  else{norm=out_spin;} 
     110 
     111 
     112// The norm is needed to make sure that the scattering cross sections are 
     113//correctly weighted, such that the sum of spin-resolved measurements adds up to 
     114// the unpolarised or half-polarised scattering cross section. No intensity weighting 
     115// needed on the incoming polariser side (assuming that a user), has normalised 
     116// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively. 
     117 
     118 
     119  weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd 
     120  weight[1] = (1.0-in_spin) * out_spin / norm;       // du 
     121  weight[2] = in_spin * (1.0-out_spin) / norm;       // ud 
     122  weight[3] = in_spin * out_spin / norm;             // uu 
    104123  weight[4] = weight[1]; // du.imag 
    105124  weight[5] = weight[2]; // ud.imag 
     
    119138    switch (xs) { 
    120139      default: // keep compiler happy; condition ensures xs in [0,1,2,3] 
    121       case 0: // uu => sld - D M_perpx 
     140      case 0: // dd => sld - D M_perpx 
    122141          return sld - px*perp; 
    123       case 1: // ud.real => -D M_perpy 
     142      case 1: // du.real => -D M_perpy 
    124143          return py*perp; 
    125       case 2: // du.real => -D M_perpy 
     144      case 2: // ud.real => -D M_perpy 
    126145          return py*perp; 
    127       case 3: // dd => sld + D M_perpx 
     146      case 3: // uu => sld + D M_perpx 
    128147          return sld + px*perp; 
    129148    } 
     
    285304    pglobal double *result,       // nq+1 return values, again with padding 
    286305    const double cutoff,          // cutoff in the dispersity weight product 
    287     int32_t effective_radius_type // which effective radius to compute 
     306    int32_t radius_effective_mode // which effective radius to compute 
    288307    ) 
    289308{ 
     
    703722      weighted_form += weight * form; 
    704723      weighted_shell += weight * shell; 
    705       if (effective_radius_type != 0) { 
    706         weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 
     724      if (radius_effective_mode != 0) { 
     725        weighted_radius += weight * CALL_RADIUS_EFFECTIVE(radius_effective_mode, local_values.table); 
    707726      } 
    708727      BUILD_ROTATION(); 
  • sasmodels/kernelcl.py

    r0d26e91 r9fac5f5  
    5353from __future__ import print_function 
    5454 
     55import sys 
    5556import os 
    5657import warnings 
     
    6162    from time import perf_counter as clock 
    6263except ImportError: # CRUFT: python < 3.3 
    63     import sys 
    6464    if sys.platform.count("darwin") > 0: 
    6565        from time import time as clock 
     
    9999# CRUFT: pyopencl < 2017.1 (as of June 2016 needs quotes around include path). 
    100100def quote_path(v): 
     101    # type: (str) -> str 
    101102    """ 
    102103    Quote the path if it is not already quoted. 
     
    110111 
    111112def fix_pyopencl_include(): 
     113    # type: (None) -> None 
    112114    """ 
    113115    Monkey patch pyopencl to allow spaces in include file path. 
    114116    """ 
    115     import pyopencl as cl 
    116     if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'): 
    117         cl._DEFAULT_INCLUDE_OPTIONS = [ 
    118             quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS 
     117    # pylint: disable=protected-access 
     118    import pyopencl 
     119    if hasattr(pyopencl, '_DEFAULT_INCLUDE_OPTIONS'): 
     120        pyopencl._DEFAULT_INCLUDE_OPTIONS = [ 
     121            quote_path(v) for v in pyopencl._DEFAULT_INCLUDE_OPTIONS 
    119122            ] 
    120123 
     
    147150 
    148151def use_opencl(): 
     152    # type: () -> bool 
     153    """Return True if OpenCL is the default computational engine""" 
    149154    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower() 
    150155    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda") 
     
    153158ENV = None 
    154159def reset_environment(): 
    155     """ 
    156     Call to create a new OpenCL context, such as after a change to SAS_OPENCL. 
     160    # type: () -> "GpuEnvironment" 
     161    """ 
     162    Return a new OpenCL context, such as after a change to SAS_OPENCL. 
    157163    """ 
    158164    global ENV 
    159165    ENV = GpuEnvironment() if use_opencl() else None 
    160  
     166    return ENV 
    161167 
    162168def environment(): 
     
    323329            return [cl.create_some_context(interactive=False)] 
    324330        except Exception as exc: 
     331            # TODO: Should warnings instead be put into logging.warn? 
    325332            warnings.warn(str(exc)) 
    326             warnings.warn("pyopencl.create_some_context() failed.  The " 
    327                 "environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might " 
    328                 "not be set correctly") 
     333            warnings.warn( 
     334                "pyopencl.create_some_context() failed.  The environment " 
     335                "variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set " 
     336                "correctly") 
    329337 
    330338    return _get_default_context() 
     
    569577 
    570578    def _call_kernel(self, call_details, values, cutoff, magnetic, 
    571                      effective_radius_type): 
     579                     radius_effective_mode): 
    572580        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray 
    573581        env = environment() 
     
    593601            self._result_b,   # Result storage. 
    594602            self._as_dtype(cutoff),  # Probability cutoff. 
    595             np.uint32(effective_radius_type),  # R_eff mode. 
     603            np.uint32(radius_effective_mode),  # R_eff mode. 
    596604        ] 
    597605 
  • sasmodels/kernelcuda.py

    rfa26e78 ra34b811  
    5959 
    6060import os 
    61 import warnings 
    6261import logging 
    6362import time 
     
    110109 
    111110def use_cuda(): 
     111    # type: None -> bool 
     112    """Returns True if CUDA is the default compute engine.""" 
    112113    sas_opencl = os.environ.get("SAS_OPENCL", "CUDA").lower() 
    113114    return HAVE_CUDA and sas_opencl.startswith("cuda") 
     
    136137        if not HAVE_CUDA: 
    137138            raise RuntimeError("CUDA startup failed with ***" 
    138                             + CUDA_ERROR + "***; using C compiler instead") 
     139                               + CUDA_ERROR + "***; using C compiler instead") 
    139140        reset_environment() 
    140141        if ENV is None: 
     
    198199    """ 
    199200    for match in FUNCTION_PATTERN.finditer(source): 
    200         print(match.group('qualifiers').replace('\n',r'\n'), match.group('function'), '(') 
     201        print(match.group('qualifiers').replace('\n', r'\n'), 
     202              match.group('function'), '(') 
    201203    return source 
    202204 
     
    260262 
    261263    def release(self): 
     264        """Free the CUDA device associated with this context.""" 
    262265        if self.context is not None: 
    263266            self.context.pop() 
     
    470473 
    471474    def _call_kernel(self, call_details, values, cutoff, magnetic, 
    472                      effective_radius_type): 
     475                     radius_effective_mode): 
    473476        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray 
    474477 
     
    489492            self._result_b,   # Result storage. 
    490493            self._as_dtype(cutoff),  # Probability cutoff. 
    491             np.uint32(effective_radius_type),  # R_eff mode. 
     494            np.uint32(radius_effective_mode),  # R_eff mode. 
    492495        ] 
    493496        grid = partition(self.q_input.nq) 
  • sasmodels/kerneldll.py

    r3199b17 ra34b811  
    174174 
    175175 
    176 def compile(source, output): 
     176def compile_model(source, output): 
    177177    # type: (str, str) -> None 
    178178    """ 
     
    265265        with os.fdopen(system_fd, "w") as file_handle: 
    266266            file_handle.write(source) 
    267         compile(source=filename, output=dll) 
     267        compile_model(source=filename, output=dll) 
    268268        # Comment the following to keep the generated C file. 
    269269        # Note: If there is a syntax error then compile raises an error 
     
    403403 
    404404    def _call_kernel(self, call_details, values, cutoff, magnetic, 
    405                      effective_radius_type): 
     405                     radius_effective_mode): 
    406406        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray 
    407407 
     
    417417            self.result.ctypes.data,   # Result storage. 
    418418            self._as_dtype(cutoff),  # Probability cutoff. 
    419             effective_radius_type,  # R_eff mode. 
     419            radius_effective_mode,  # R_eff mode. 
    420420        ] 
    421421 
  • sasmodels/kernelpy.py

    r3199b17 ra34b811  
    1616    from numpy import cbrt 
    1717except ImportError: 
    18     def cbrt(x): return x ** (1.0/3.0) 
     18    def cbrt(x): 
     19        """Return cubed root of x.""" 
     20        return x ** (1.0/3.0) 
    1921 
    2022from .generate import F64 
     
    4345        self.info = model_info 
    4446        self.dtype = np.dtype('d') 
    45         logger.info("make python model " + self.info.name) 
     47        logger.info("make python model %s", self.info.name) 
    4648 
    4749    def make_kernel(self, q_vectors): 
     50        """Instantiate the python kernel with input *q_vectors*""" 
    4851        q_input = PyInput(q_vectors, dtype=F64) 
    4952        return PyKernel(self.info, q_input) 
     
    169172        volume = model_info.form_volume 
    170173        shell = model_info.shell_volume 
    171         radius = model_info.effective_radius 
     174        radius = model_info.radius_effective 
    172175        self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume 
    173176                        else (lambda: [volume(*volume_args)]*2) if volume 
     
    177180                        else (lambda mode: 1.0)) 
    178181 
    179     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     182    def _call_kernel(self, call_details, values, cutoff, magnetic, radius_effective_mode): 
    180183        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    181184        if magnetic: 
     
    183186        #print("Calling python kernel") 
    184187        #call_details.show(values) 
    185         radius = ((lambda: 0.0) if effective_radius_type == 0 
    186                   else (lambda: self._radius(effective_radius_type))) 
     188        radius = ((lambda: 0.0) if radius_effective_mode == 0 
     189                  else (lambda: self._radius(radius_effective_mode))) 
    187190        self.result = _loops( 
    188191            self._parameter_vector, self._form, self._volume, radius, 
  • sasmodels/list_pars.py

    r2d81cfe rb297ba9  
    1616from .compare import columnize 
    1717 
    18 def find_pars(type=None): 
     18def find_pars(kind=None): 
    1919    """ 
    2020    Find all parameters in all models. 
     
    2626        model_info = load_model_info(name) 
    2727        for p in model_info.parameters.kernel_parameters: 
    28             if type is None or p.type == type: 
     28            if kind is None or p.type == kind: 
    2929                partable.setdefault(p.name, []) 
    3030                partable[p.name].append(name) 
    3131    return partable 
    3232 
    33 def list_pars(names_only=True, type=None): 
     33def list_pars(names_only=True, kind=None): 
    3434    """ 
    3535    Print all parameters in all models. 
     
    3838    occurs in. 
    3939    """ 
    40     partable = find_pars(type) 
     40    partable = find_pars(kind) 
    4141    if names_only: 
    4242        print(columnize(list(sorted(partable.keys())))) 
     
    5757        help="list models which use this argument") 
    5858    parser.add_argument( 
    59         'type', default="any", nargs='?', 
     59        'kind', default="any", nargs='?', 
    6060        metavar="volume|orientation|sld|none|any", 
    6161        choices=['volume', 'orientation', 'sld', None, 'any'], 
    62         type=lambda v: None if v == 'any' else '' if v == 'none' else v, 
    63         help="only list arguments of the given type") 
     62        type=lambda v: None if v == 'any' else None if v == 'none' else v, 
     63        help="only list arguments of the given kind") 
    6464    args = parser.parse_args() 
    6565 
    66     list_pars(names_only=not args.verbose, type=args.type) 
     66    list_pars(names_only=not args.verbose, kind=args.kind) 
    6767 
    6868if __name__ == "__main__": 
  • sasmodels/mixture.py

    r39a06c9 rb297ba9  
    120120 
    121121    def random(): 
     122        """Random set of model parameters for mixture model""" 
    122123        combined_pars = {} 
    123124        for k, part in enumerate(parts): 
     
    149150 
    150151class MixtureModel(KernelModel): 
     152    """ 
     153    Model definition for mixture of models. 
     154    """ 
    151155    def __init__(self, model_info, parts): 
    152156        # type: (ModelInfo, List[KernelModel]) -> None 
     
    165169        kernels = [part.make_kernel(q_vectors) for part in self.parts] 
    166170        return MixtureKernel(self.info, kernels) 
     171    make_kernel.__doc__ = KernelModel.make_kernel.__doc__ 
    167172 
    168173    def release(self): 
    169174        # type: () -> None 
    170         """ 
    171         Free resources associated with the model. 
    172         """ 
     175        """Free resources associated with the model.""" 
    173176        for part in self.parts: 
    174177            part.release() 
     178    release.__doc__ = KernelModel.release.__doc__ 
    175179 
    176180 
    177181class MixtureKernel(Kernel): 
     182    """ 
     183    Instantiated kernel for mixture of models. 
     184    """ 
    178185    def __init__(self, model_info, kernels): 
    179186        # type: (ModelInfo, List[Kernel]) -> None 
     
    185192        self.results = []  # type: List[np.ndarray] 
    186193 
    187     def __call__(self, call_details, values, cutoff, magnetic): 
     194    def Iq(self, call_details, values, cutoff, magnetic): 
    188195        # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray 
    189196        scale, background = values[0:2] 
     
    191198        # remember the parts for plotting later 
    192199        self.results = []  # type: List[np.ndarray] 
    193         parts = MixtureParts(self.info, self.kernels, call_details, values) 
     200        parts = _MixtureParts(self.info, self.kernels, call_details, values) 
    194201        for kernel, kernel_details, kernel_values in parts: 
    195202            #print("calling kernel", kernel.info.name) 
     
    208215        return scale*total + background 
    209216 
     217    Iq.__doc__ = Kernel.Iq.__doc__ 
     218    __call__ = Iq 
     219 
    210220    def release(self): 
    211221        # type: () -> None 
     222        """Free resources associated with the kernel.""" 
    212223        for k in self.kernels: 
    213224            k.release() 
    214225 
    215226 
    216 class MixtureParts(object): 
     227# Note: _MixtureParts doesn't implement iteration correctly, and only allows 
     228# a single iterator to be active at once.  It doesn't matter in this case 
     229# since _MixtureParts is only used in one place, but it is not clean style. 
     230class _MixtureParts(object): 
     231    """ 
     232    Mixture component iterator. 
     233    """ 
    217234    def __init__(self, model_info, kernels, call_details, values): 
    218235        # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None 
     
    223240        self.values = values 
    224241        self.spin_index = model_info.parameters.npars + 2 
     242        # The following are redefined by __iter__, but set them here so that 
     243        # lint complains a little less. 
     244        self.part_num = -1 
     245        self.par_index = -1 
     246        self.mag_index = -1 
    225247        #call_details.show(values) 
    226248 
  • sasmodels/model_test.py

    rd92182f r8795b6f  
    6161import numpy as np  # type: ignore 
    6262 
    63 from . import core 
    6463from .core import list_models, load_model_info, build_model 
    6564from .direct_model import call_kernel, call_Fq 
     
    137136        test_method_name = "test_%s_python" % model_info.id 
    138137        test = ModelTestCase(test_name, model_info, 
    139                                 test_method_name, 
    140                                 platform="dll",  # so that 
    141                                 dtype="double", 
    142                                 stash=stash) 
     138                             test_method_name, 
     139                             platform="dll",  # so that 
     140                             dtype="double", 
     141                             stash=stash) 
    143142        suite.addTest(test) 
    144143    else:   # kernel implemented in C 
     
    149148            test_method_name = "test_%s_dll" % model_info.id 
    150149            test = ModelTestCase(test_name, model_info, 
    151                                     test_method_name, 
    152                                     platform="dll", 
    153                                     dtype="double", 
    154                                     stash=stash) 
     150                                 test_method_name, 
     151                                 platform="dll", 
     152                                 dtype="double", 
     153                                 stash=stash) 
    155154            suite.addTest(test) 
    156155 
     
    164163            # presence of *single=False* in the model file. 
    165164            test = ModelTestCase(test_name, model_info, 
    166                                     test_method_name, 
    167                                     platform="ocl", dtype=None, 
    168                                     stash=stash) 
     165                                 test_method_name, 
     166                                 platform="ocl", dtype=None, 
     167                                 stash=stash) 
    169168            #print("defining", test_name) 
    170169            suite.addTest(test) 
     
    179178            # presence of *single=False* in the model file. 
    180179            test = ModelTestCase(test_name, model_info, 
    181                                     test_method_name, 
    182                                     platform="cuda", dtype=None, 
    183                                     stash=stash) 
     180                                 test_method_name, 
     181                                 platform="cuda", dtype=None, 
     182                                 stash=stash) 
    184183            #print("defining", test_name) 
    185184            suite.addTest(test) 
     
    241240                    s_name = pars.pop('@S') 
    242241                    ps_test = [pars] + list(test[1:]) 
     242                    #print("PS TEST PARAMS!!!",ps_test) 
    243243                    # build the P@S model 
    244244                    s_info = load_model_info(s_name) 
     
    247247                                           platform=self.platform) 
    248248                    # run the tests 
     249                    #self.info = ps_model.info 
     250                    #print("SELF.INFO PARAMS!!!",[p.id for p in self.info.parameters.call_parameters]) 
     251                    #print("PS MODEL PARAMETERS:",[p.id for p in ps_model.info.parameters.call_parameters]) 
    249252                    results.append(self.run_one(ps_model, ps_test)) 
    250253 
     
    304307            """Run a single test case.""" 
    305308            user_pars, x, y = test[:3] 
    306             pars = expand_pars(self.info.parameters, user_pars) 
    307             invalid = invalid_pars(self.info.parameters, pars) 
     309            #print("PS MODEL PARAMETERS:",[p.id for p in model.info.parameters.call_parameters]) 
     310            pars = expand_pars(model.info.parameters, user_pars) 
     311            invalid = invalid_pars(model.info.parameters, pars) 
    308312            if invalid: 
    309313                raise ValueError("Unknown parameters in test: " + ", ".join(invalid)) 
     
    329333            else: 
    330334                y1 = y 
    331                 y2 = test[3] if not isinstance(test[3], list) else [test[3]] 
    332                 F1, F2, R_eff, volume, volume_ratio = call_Fq(kernel, pars) 
    333                 if F1 is not None:  # F1 is none for models with Iq instead of Fq 
    334                     self._check_vectors(x, y1, F1, 'F') 
    335                 self._check_vectors(x, y2, F2, 'F^2') 
     335                y2 = test[3] if isinstance(test[3], list) else [test[3]] 
     336                F, Fsq, R_eff, volume, volume_ratio = call_Fq(kernel, pars) 
     337                if F is not None:  # F is none for models with Iq instead of Fq 
     338                    self._check_vectors(x, y1, F, 'F') 
     339                self._check_vectors(x, y2, Fsq, 'F^2') 
    336340                self._check_scalar(test[4], R_eff, 'R_eff') 
    337341                self._check_scalar(test[5], volume, 'volume') 
    338342                self._check_scalar(test[6], volume_ratio, 'form:shell ratio') 
    339                 return F2 
     343                return Fsq 
    340344 
    341345        def _check_scalar(self, target, actual, name): 
    342             if target is None: 
    343                 # smoke test --- make sure it runs and produces a value 
    344                 self.assertTrue(not np.isnan(actual), 
    345                                 'invalid %s: %s' % (name, actual)) 
    346             elif np.isnan(target): 
    347                 # make sure nans match 
    348                 self.assertTrue(np.isnan(actual), 
    349                                 '%s: expected:%s; actual:%s' 
    350                                 % (name, target, actual)) 
    351             else: 
    352                 # is_near does not work for infinite values, so also test 
    353                 # for exact values. 
    354                 self.assertTrue(target == actual or is_near(target, actual, 5), 
    355                                 '%s: expected:%s; actual:%s' 
    356                                 % (name, target, actual)) 
     346            self.assertTrue(is_near(target, actual, 5), 
     347                            '%s: expected:%s; actual:%s' 
     348                            % (name, target, actual)) 
    357349 
    358350        def _check_vectors(self, x, target, actual, name='I'): 
     
    364356                             '%s(...) returned wrong length'%name) 
    365357            for xi, yi, actual_yi in zip(x, target, actual): 
    366                 if yi is None: 
    367                     # smoke test --- make sure it runs and produces a value 
    368                     self.assertTrue(not np.isnan(actual_yi), 
    369                                     'invalid %s(%s): %s' % (name, xi, actual_yi)) 
    370                 elif np.isnan(yi): 
    371                     # make sure nans match 
    372                     self.assertTrue(np.isnan(actual_yi), 
    373                                     '%s(%s): expected:%s; actual:%s' 
    374                                     % (name, xi, yi, actual_yi)) 
    375                 else: 
    376                     # is_near does not work for infinite values, so also test 
    377                     # for exact values. 
    378                     self.assertTrue(yi == actual_yi or is_near(yi, actual_yi, 5), 
    379                                     '%s(%s); expected:%s; actual:%s' 
    380                                     % (name, xi, yi, actual_yi)) 
     358                self.assertTrue(is_near(yi, actual_yi, 5), 
     359                                '%s(%s): expected:%s; actual:%s' 
     360                                % (name, xi, target, actual)) 
    381361 
    382362    return ModelTestCase 
     
    390370    invalid = [] 
    391371    for par in sorted(pars.keys()): 
    392         # special handling of R_eff mode, which is not a usual parameter 
     372        # Ignore the R_eff mode parameter when checking for valid parameters. 
     373        # It is an allowed parameter for a model even though it does not exist 
     374        # in the parameter table.  The call_Fq() function pops it from the 
     375        # parameter list and sends it directly to kernel.Fq(). 
    393376        if par == product.RADIUS_MODE_ID: 
    394377            continue 
     
    406389    """ 
    407390    Returns true if *actual* is within *digits* significant digits of *target*. 
    408     """ 
    409     import math 
    410     shift = 10**math.ceil(math.log10(abs(target))) 
    411     return abs(target-actual)/shift < 1.5*10**-digits 
     391 
     392    *taget* zero and inf should match *actual* zero and inf.  If you want to 
     393    accept eps for zero, choose a value such as 1e-10, which must match up to 
     394    +/- 1e-15 when *digits* is the default value of 5. 
     395 
     396    If *target* is None, then just make sure that *actual* is not NaN. 
     397 
     398    If *target* is NaN, make sure *actual* is NaN. 
     399    """ 
     400    if target is None: 
     401        # target is None => actual cannot be NaN 
     402        return not np.isnan(actual) 
     403    elif target == 0.: 
     404        # target is 0. => actual must be 0. 
     405        # Note: if small values are allowed, then use maybe test zero against eps instead? 
     406        return actual == 0. 
     407    elif np.isfinite(target): 
     408        shift = np.ceil(np.log10(abs(target))) 
     409        return abs(target-actual) < 1.5*10**(shift-digits) 
     410    elif target == actual: 
     411        # target is inf => actual must be inf of same sign 
     412        return True 
     413    else: 
     414        # target is NaN => actual must be NaN 
     415        return np.isnan(target) == np.isnan(actual) 
    412416 
    413417# CRUFT: old interface; should be deprecated and removed 
    414418def run_one(model_name): 
     419    # type: (str) -> str 
     420    """ 
     421    [Deprecated] Run the tests associated with *model_name*. 
     422 
     423    Use the following instead:: 
     424 
     425        succss, output = check_model(load_model_info(model_name)) 
     426    """ 
    415427    # msg = "use check_model(model_info) rather than run_one(model_name)" 
    416428    # warnings.warn(msg, category=DeprecationWarning, stacklevel=2) 
     
    421433        return output 
    422434 
    423     success, output = check_model(model_info) 
     435    _, output = check_model(model_info) 
    424436    return output 
    425437 
    426438def check_model(model_info): 
    427     # type: (ModelInfo) -> str 
     439    # type: (ModelInfo) -> Tuple[bool, str] 
    428440    """ 
    429441    Run the tests for a single model, capturing the output. 
     
    492504        loaders.append('cuda') 
    493505    tests = make_suite(loaders, ['all']) 
    494     def build_test(test): 
     506    def _build_test(test): 
    495507        # In order for nosetest to show the test name, wrap the test.run_all 
    496508        # instance in function that takes the test name as a parameter which 
     
    510522 
    511523    for test in tests: 
    512         yield build_test(test) 
     524        yield _build_test(test) 
    513525 
    514526 
     
    539551                        "models will be tested.  See core.list_models() for " 
    540552                        "names of other groups, such as 'py' or 'single'.") 
    541     args, models = parser.parse_known_args() 
    542  
    543     if args.engine == "opencl": 
     553    opts = parser.parse_args() 
     554 
     555    if opts.engine == "opencl": 
    544556        if not use_opencl(): 
    545557            print("opencl is not available") 
    546558            return 1 
    547559        loaders = ['opencl'] 
    548     elif args.engine == "dll": 
     560    elif opts.engine == "dll": 
    549561        loaders = ["dll"] 
    550     elif args.engine == "cuda": 
     562    elif opts.engine == "cuda": 
    551563        if not use_cuda(): 
    552564            print("cuda is not available") 
    553565            return 1 
    554566        loaders = ['cuda'] 
    555     elif args.engine == "all": 
     567    elif opts.engine == "all": 
    556568        loaders = ['dll'] 
    557569        if use_opencl(): 
     
    560572            loaders.append('cuda') 
    561573    else: 
    562         print("unknown engine " + args.engine) 
     574        print("unknown engine " + opts.engine) 
    563575        return 1 
    564576 
    565     runner = TestRunner(verbosity=args.verbose, **test_args) 
    566     result = runner.run(make_suite(loaders, args.models)) 
     577    runner = TestRunner(verbosity=opts.verbose, **test_args) 
     578    result = runner.run(make_suite(loaders, opts.models)) 
    567579    return 1 if result.failures or result.errors else 0 
    568580 
  • sasmodels/models/_spherepy.py

    r304c775 ra34b811  
    3333to the output of the software provided by the NIST (Kline, 2006). 
    3434 
    35  
    3635References 
    3736---------- 
     
    4039John Wiley and Sons, New York, (1955) 
    4140 
     41Source 
     42------ 
     43 
     44`_spherepy.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/_spherepy.py>`_ 
     45`sphere.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/sphere.c>`_ 
     46 
     47Authorship and Verification 
     48---------------------------- 
     49 
     50* **Author: P Kienzle** 
     51* **Last Modified by:** 
    4252* **Last Reviewed by:** S King and P Parker **Date:** 2013/09/09 and 2014/01/06 
     53* **Source added by :** Steve King **Date:** March 25, 2019 
    4354""" 
    4455 
     
    6980 
    7081def form_volume(radius): 
     82    """Calculate volume for sphere""" 
    7183    return 1.333333333333333 * pi * radius ** 3 
    7284 
    73 def effective_radius(mode, radius): 
    74     return radius 
     85def radius_effective(mode, radius): 
     86    """Calculate R_eff for sphere""" 
     87    return radius if mode else 0. 
    7588 
    7689def Iq(q, sld, sld_solvent, radius): 
     90    """Calculate I(q) for sphere""" 
    7791    #print "q",q 
    7892    #print "sld,r",sld,sld_solvent,radius 
  • sasmodels/models/adsorbed_layer.py

    r2d81cfe r0507e09  
    4949   Layers*, *Macromol. Symp.*, 190 (2002) 33-42. 
    5050 
     51Source 
     52------ 
     53 
     54`adsorbed_layer.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/adsorbed_layer.py>`_ 
     55 
    5156Authorship and Verification 
    5257---------------------------- 
     
    5560* **Last Modified by:** Paul Kienzle **Date:** April 14, 2016 
    5661* **Last Reviewed by:** Steve King **Date:** March 18, 2016 
     62* **Source added by :** Steve King **Date:** March 25, 2019 
    5763""" 
    5864 
     
    8793def Iq(q, second_moment, adsorbed_amount, density_shell, radius, 
    8894       volfraction, sld_shell, sld_solvent): 
     95    """Return I(q) for adsorbed layer model.""" 
    8996    with errstate(divide='ignore'): 
    9097        aa = ((sld_shell - sld_solvent)/density_shell * adsorbed_amount) / q 
     
    96103 
    97104def random(): 
     105    """Return a random parameter set for the model.""" 
    98106    # only care about the value of second_moment: 
    99107    #    curve = scale * e**(-second_moment^2 q^2)/q^2 
  • sasmodels/models/barbell.c

    r99658f6 ra34b811  
    8585 
    8686static double 
    87 effective_radius(int mode, double radius_bell, double radius, double length) 
     87radius_effective(int mode, double radius_bell, double radius, double length) 
    8888{ 
    8989    switch (mode) { 
  • sasmodels/models/barbell.py

    r99658f6 ra34b811  
    7676---------- 
    7777 
    78 .. [#] H Kaya, *J. Appl. Cryst.*, 37 (2004) 37 223-230 
     78.. [#] H Kaya, *J. Appl. Cryst.*, 37 (2004) 223-230 
    7979.. [#] H Kaya and N R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda 
    8080   and errata) 
    81 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     81.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     82 
     83Source 
     84------ 
     85 
     86`barbell.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/barbell.py>`_ 
     87 
     88`barbell.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/barbell.c>`_ 
    8289 
    8390Authorship and Verification 
     
    8794* **Last Modified by:** Paul Butler **Date:** March 20, 2016 
    8895* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
     96* **Source added by :** Steve King **Date:** March 25, 2019 
    8997""" 
    9098 
     
    117125 
    118126source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
    119 have_Fq = True  
    120 effective_radius_type = [ 
    121     "equivalent cylinder excluded volume","equivalent volume sphere", "radius", "half length", "half total length", 
     127have_Fq = True 
     128radius_effective_modes = [ 
     129    "equivalent cylinder excluded volume", "equivalent volume sphere", 
     130    "radius", "half length", "half total length", 
    122131    ] 
    123132 
    124133def random(): 
     134    """Return a random parameter set for the model.""" 
    125135    # TODO: increase volume range once problem with bell radius is fixed 
    126136    # The issue is that bell radii of more than about 200 fail at high q 
  • sasmodels/models/bcc_paracrystal.py

    rda7b26b r0507e09  
    11r""" 
    2 .. warning:: This model and this model description are under review following  
    3              concerns raised by SasView users. If you need to use this model,  
    4              please email help@sasview.org for the latest situation. *The  
     2.. warning:: This model and this model description are under review following 
     3             concerns raised by SasView users. If you need to use this model, 
     4             please email help@sasview.org for the latest situation. *The 
    55             SasView Developers. September 2018.* 
    66 
     
    100100   (Corrections to FCC and BCC lattice structure calculation) 
    101101 
     102Source 
     103------ 
     104 
     105`bcc_paracrystal.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/bcc_paracrystal.py>`_ 
     106 
     107`bcc_paracrystal.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/bcc_paracrystal.c>`_ 
     108 
    102109Authorship and Verification 
    103110--------------------------- 
     
    106113* **Last Modified by:** Paul Butler **Date:** September 29, 2016 
    107114* **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 
     115* **Source added by :** Steve King **Date:** March 25, 2019 
    108116""" 
    109117 
     
    141149 
    142150def random(): 
     151    """Return a random parameter set for the model.""" 
    143152    # Define lattice spacing as a multiple of the particle radius 
    144153    # using the formulat a = 4 r/sqrt(3).  Systems which are ordered 
  • sasmodels/models/be_polyelectrolyte.py

    rca77fc1 r0507e09  
    3636  constituting the polymer monomer and the solvent molecules, respectively. 
    3737 
    38 - $v_p$ and $v_s$ are the partial molar volume of the polymer and the  
     38- $v_p$ and $v_s$ are the partial molar volume of the polymer and the 
    3939  solvent, respectively. 
    4040 
     
    5050- $C_s$ is the concentration of monovalent salt(1/|Ang^3| - internally converted from mol/L). 
    5151 
    52 - $\alpha$ is the degree of ionization (the ratio of charged monomers to the total  
     52- $\alpha$ is the degree of ionization (the ratio of charged monomers to the total 
    5353  number of monomers) 
    5454 
     
    7474dimensionally useful units of 1/|Ang^3|, only the converted version of the 
    7575polymer concentration was actually being used in the calculation while the 
    76 unconverted salt concentration (still in apparent units of mol/L) was being  
    77 used.  This was carried through to Sasmodels as used for SasView 4.1 (though  
    78 the line of code converting the salt concentration to the new units was removed  
    79 somewhere along the line). Simple dimensional analysis of the calculation shows  
    80 that the converted salt concentration should be used, which the original code  
    81 suggests was the intention, so this has now been corrected (for SasView 4.2).  
     76unconverted salt concentration (still in apparent units of mol/L) was being 
     77used.  This was carried through to Sasmodels as used for SasView 4.1 (though 
     78the line of code converting the salt concentration to the new units was removed 
     79somewhere along the line). Simple dimensional analysis of the calculation shows 
     80that the converted salt concentration should be used, which the original code 
     81suggests was the intention, so this has now been corrected (for SasView 4.2). 
    8282Once better validation has been performed this note will be removed. 
    8383 
     
    9191.. [#] E Raphael, J F Joanny, *Europhysics Letters*, 11 (1990) 179 
    9292 
     93Source 
     94------ 
     95 
     96`be_polyelectrolyte.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/be_polyelectrolyte.py>`_ 
     97 
    9398Authorship and Verification 
    9499---------------------------- 
     
    97102* **Last Modified by:** Paul Butler **Date:** September 25, 2018 
    98103* **Last Reviewed by:** Paul Butler **Date:** September 25, 2018 
     104* **Source added by :** Steve King **Date:** March 25, 2019 
    99105""" 
    100106 
     
    141147    :params: see parameter table 
    142148    :return: 1-D form factor for polyelectrolytes in low salt 
    143      
     149 
    144150    parameter names, units, default values, and behavior (volume, sld etc) are 
    145151    defined in the parameter table.  The concentrations are converted from 
     
    167173 
    168174def random(): 
     175    """Return a random parameter set for the model.""" 
    169176    # TODO: review random be_polyelectrolyte model generation 
    170177    pars = dict( 
  • sasmodels/models/binary_hard_sphere.py

    r2d81cfe r0507e09  
    6565.. [#] S R Kline, *J Appl. Cryst.*, 39 (2006) 895 
    6666 
     67Source 
     68------ 
     69 
     70`binary_hard_sphere.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/binary_hard_sphere.py>`_ 
     71 
     72`binary_hard_sphere.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/binary_hard_sphere.c>`_ 
     73 
    6774Authorship and Verification 
    6875---------------------------- 
     
    7178* **Last Modified by:** Paul Butler **Date:** March 20, 2016 
    7279* **Last Reviewed by:** Paul Butler **Date:** March 20, 2016 
     80* **Source added by :** Steve King **Date:** March 25, 2019 
    7381""" 
    7482 
     
    112120 
    113121def random(): 
     122    """Return a random parameter set for the model.""" 
    114123    # TODO: binary_hard_sphere fails at low qr 
    115124    radius_lg = 10**np.random.uniform(2, 4.7) 
  • sasmodels/models/broad_peak.py

    r2d81cfe r0507e09  
    3333None. 
    3434 
     35Source 
     36------ 
     37 
     38`broad_peak.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/broad_peak.py>`_ 
     39 
    3540Authorship and Verification 
    3641---------------------------- 
     
    3944* **Last Modified by:** Paul kienle **Date:** July 24, 2016 
    4045* **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 
     46* **Source added by :** Steve King **Date:** March 25, 2019 
    4147""" 
    4248 
     
    96102 
    97103def random(): 
     104    """Return a random parameter set for the model.""" 
    98105    pars = dict( 
    99106        scale=1, 
  • sasmodels/models/capped_cylinder.c

    r99658f6 ra34b811  
    107107 
    108108static double 
    109 effective_radius(int mode, double radius, double radius_cap, double length) 
     109radius_effective(int mode, double radius, double radius_cap, double length) 
    110110{ 
    111111    switch (mode) { 
  • sasmodels/models/capped_cylinder.py

    r99658f6 ra34b811  
    8282.. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda 
    8383   and errata) 
    84 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     84.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     85 
     86Source 
     87------ 
     88 
     89`capped_cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/capped_cylinder.py>`_ 
     90 
     91`capped_cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/capped_cylinder.c>`_ 
    8592 
    8693Authorship and Verification 
     
    9097* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    9198* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
     99* **Source added by :** Steve King **Date:** March 25, 2019 
    92100""" 
    93101 
     
    138146source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
    139147have_Fq = True 
    140 effective_radius_type = [ 
    141     "equivalent cylinder excluded volume", "equivalent volume sphere", "radius", "half length", "half total length", 
     148radius_effective_modes = [ 
     149    "equivalent cylinder excluded volume", "equivalent volume sphere", 
     150    "radius", "half length", "half total length", 
    142151    ] 
    143152 
    144153def random(): 
     154    """Return a random parameter set for the model.""" 
    145155    # TODO: increase volume range once problem with bell radius is fixed 
    146156    # The issue is that bell radii of more than about 200 fail at high q 
  • sasmodels/models/core_multi_shell.c

    rd42dd4a ra34b811  
    2626 
    2727static double 
    28 effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 
     28radius_effective(int mode, double core_radius, double fp_n, double thickness[]) 
    2929{ 
    3030  switch (mode) { 
  • sasmodels/models/core_multi_shell.py

    ree60aa7 ra34b811  
    3939   Neutron Scattering*, Plenum Press, New York, 1987. 
    4040 
     41Source 
     42------ 
     43 
     44`core_multi_shell.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_multi_shell.py>`_ 
     45 
     46`core_multi_shell.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_multi_shell.c>`_ 
     47 
    4148Authorship and Verification 
    4249---------------------------- 
     
    4552* **Last Modified by:** Paul Kienzle **Date:** September 12, 2016 
    4653* **Last Reviewed by:** Paul Kienzle **Date:** September 12, 2016 
     54* **Source added by :** Steve King **Date:** March 25, 2019 
    4755""" 
    4856from __future__ import division 
     
    100108source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
    101109have_Fq = True 
    102 effective_radius_type = ["outer radius", "core radius"] 
     110radius_effective_modes = ["outer radius", "core radius"] 
    103111 
    104112def random(): 
     113    """Return a random parameter set for the model.""" 
    105114    num_shells = np.minimum(np.random.poisson(3)+1, 10) 
    106115    total_radius = 10**np.random.uniform(1.7, 4) 
  • sasmodels/models/core_shell_bicelle.c

    r99658f6 ra34b811  
    6060 
    6161static double 
    62 effective_radius(int mode, double radius, double thick_rim, double thick_face, double length) 
     62radius_effective(int mode, double radius, double thick_rim, double thick_face, double length) 
    6363{ 
    6464    switch (mode) { 
  • sasmodels/models/core_shell_bicelle.py

    r99658f6 ra34b811  
    8989   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    9090   =26379>`_ 
    91     
    92    L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     91 
     92.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     93 
     94Source 
     95------ 
     96 
     97`core_shell_bicelle.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_bicelle.py>`_ 
     98 
     99`core_shell_bicelle.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_bicelle.c>`_ 
    93100 
    94101Authorship and Verification 
     
    98105* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    99106* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
     107* **Source added by :** Steve King **Date:** March 25, 2019 
    100108""" 
    101109 
     
    157165          "core_shell_bicelle.c"] 
    158166have_Fq = True 
    159 effective_radius_type = [ 
    160     "excluded volume","equivalent volume sphere", "outer rim radius", 
     167radius_effective_modes = [ 
     168    "excluded volume", "equivalent volume sphere", "outer rim radius", 
    161169    "half outer thickness", "half diagonal", 
    162170    ] 
    163171 
    164172def random(): 
     173    """Return a random parameter set for the model.""" 
    165174    pars = dict( 
    166175        radius=10**np.random.uniform(1.3, 3), 
  • sasmodels/models/core_shell_bicelle_elliptical.c

    r99658f6 ra34b811  
    3535 
    3636static double 
    37 effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     37radius_effective(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
    3838{ 
    3939    switch (mode) { 
  • sasmodels/models/core_shell_bicelle_elliptical.py

    r99658f6 ra34b811  
    9999---------- 
    100100 
    101 .. [#] 
    102 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     101.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     102 
     103Source 
     104------ 
     105 
     106`core_shell_bicelle_elliptical.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_bicelle_elliptical.py>`_ 
     107 
     108`core_shell_bicelle_elliptical.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_bicelle_elliptical.c>`_ 
    103109 
    104110Authorship and Verification 
     
    108114* **Last Modified by:**  Richard Heenan **Date:** December 14, 2016 
    109115* **Last Reviewed by:**  Paul Kienzle **Date:** Feb 28, 2018 
     116* **Source added by :**  Steve King **Date:** March 25, 2019 
    110117""" 
    111118 
     
    148155          "core_shell_bicelle_elliptical.c"] 
    149156have_Fq = True 
    150 effective_radius_type = [ 
    151     "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", 
     157radius_effective_modes = [ 
     158    "equivalent cylinder excluded volume", "equivalent volume sphere", 
     159    "outer rim average radius", "outer rim min radius", 
    152160    "outer max radius", "half outer thickness", "half diagonal", 
    153161    ] 
    154162 
    155163def random(): 
     164    """Return a random parameter set for the model.""" 
    156165    outer_major = 10**np.random.uniform(1, 4.7) 
    157166    outer_minor = 10**np.random.uniform(1, 4.7) 
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r99658f6 ra34b811  
    3636 
    3737static double 
    38 effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     38radius_effective(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
    3939{ 
    4040    switch (mode) { 
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r99658f6 ra34b811  
    111111---------- 
    112112 
    113 .. [#] 
    114 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     113.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     114 
     115Source 
     116------ 
     117 
     118`core_shell_bicelle_elliptical_belt_rough.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py>`_ 
     119 
     120`core_shell_bicelle_elliptical_belt_rough.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c>`_ 
    115121 
    116122Authorship and Verification 
     
    120126* **Last Modified by:**  Richard Heenan new 2d orientation **Date:** October 5, 2017 
    121127* **Last Reviewed by:**  Richard Heenan 2d calc seems agree with 1d **Date:** Nov 2, 2017 
     128* **Source added by :**  Steve King **Date:** March 25, 2019 
    122129""" 
    123130 
     
    161168          "core_shell_bicelle_elliptical_belt_rough.c"] 
    162169have_Fq = True 
    163 effective_radius_type = [ 
    164     "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", 
     170radius_effective_modes = [ 
     171    "equivalent cylinder excluded volume", "equivalent volume sphere", 
     172    "outer rim average radius", "outer rim min radius", 
    165173    "outer max radius", "half outer thickness", "half diagonal", 
    166174    ] 
     175 
     176# TODO: No random() for core-shell bicelle elliptical belt rough 
    167177 
    168178demo = dict(scale=1, background=0, 
     
    190200    #[{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'VR', 1], 
    191201 
    192     [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0, 
    193       'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0}, 
     202    [{'radius': 30.0, 'x_core': 3.0, 'thick_rim': 8.0, 'thick_face': 14.0, 
     203      'length': 50.0, 'sld_core': 4.0, 'sld_face': 7.0, 'sld_rim': 1.0, 
     204      'sld_solvent': 6.0, 'background': 0.0}, 
    194205     0.015, 189.328], 
    195206    #[{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
  • sasmodels/models/core_shell_cylinder.c

    r99658f6 ra34b811  
    3737 
    3838static double 
    39 effective_radius(int mode, double radius, double thickness, double length) 
     39radius_effective(int mode, double radius, double thickness, double length) 
    4040{ 
    4141    switch (mode) { 
  • sasmodels/models/core_shell_cylinder.py

    r99658f6 rdb1d9d5  
    4646density of the solvent, and *background* is the background level.  The outer 
    4747radius of the shell is given by $R+T$ and the total length of the outer 
    48 shell is given by $L+2T$. $J1$ is the first order Bessel function. 
     48shell is given by $L+2T$. $J_1$ is the first order Bessel function. 
    4949 
    5050.. _core-shell-cylinder-geometry: 
     
    7070   1445-1452 
    7171.. [#kline] S R Kline, *J Appl. Cryst.*, 39 (2006) 895 
    72 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     72.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     73 
     74Source 
     75------ 
     76 
     77`core_shell_cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_cylinder.py>`_ 
     78 
     79`core_shell_cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_cylinder.c>`_ 
    7380 
    7481Authorship and Verification 
     
    7885* **Last Modified by:** Paul Kienzle **Date:** Aug 8, 2016 
    7986* **Last Reviewed by:** Richard Heenan **Date:** March 18, 2016 
     87* **Source added by :** Steve King **Date:** March 25, 2019 
    8088""" 
    8189 
     
    133141source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
    134142have_Fq = True 
    135 effective_radius_type = [ 
     143radius_effective_modes = [ 
    136144    "excluded volume", "equivalent volume sphere", "outer radius", "half outer length", 
    137145    "half min outer dimension", "half max outer dimension", "half outer diagonal", 
     
    139147 
    140148def random(): 
     149    """Return a random parameter set for the model.""" 
    141150    outer_radius = 10**np.random.uniform(1, 4.7) 
    142151    # Use a distribution with a preference for thin shell or thin core 
  • sasmodels/models/core_shell_ellipsoid.c

    r99658f6 ra34b811  
    6868 
    6969static double 
    70 effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     70radius_effective(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
    7171{ 
    7272    const double radius_equat_tot = radius_equat_core + thick_shell; 
  • sasmodels/models/core_shell_ellipsoid.py

    r99658f6 r830cf6b  
    33---------- 
    44 
    5 Parameters for this model are the core axial ratio X and a shell thickness, 
    6 which are more often what we would like to determine and makes the model 
    7 better behaved, particularly when polydispersity is applied than the four 
    8 independent radii used in the original parameterization of this model. 
     5Parameters for this model are the core axial ratio $X_{core}$ and a shell 
     6thickness $t_{shell}$, which are more often what we would like to determine 
     7and make the model better behaved, particularly when polydispersity is 
     8applied, than the four independent radii used in the original parameterization 
     9of this model. 
    910 
    1011 
     
    1516the poles, of a prolate ellipsoid. 
    1617 
    17 When *X_core < 1* the core is oblate; when *X_core > 1* it is prolate. 
    18 *X_core = 1* is a spherical core. 
    19  
    20 For a fixed shell thickness *XpolarShell = 1*, to scale the shell thickness 
    21 pro-rata with the radius set or constrain *XpolarShell = X_core*. 
    22  
    23 When including an $S(q)$, the radius in $S(q)$ is calculated to be that of 
    24 a sphere with the same 2nd virial coefficient of the outer surface of the 
    25 ellipsoid. This may have some undesirable effects if the aspect ratio of the 
    26 ellipsoid is large (ie, if $X << 1$ or $X >> 1$ ), when the $S(q)$ 
    27 - which assumes spheres - will not in any case be valid.  Generating a 
    28 custom product model will enable separate effective volume fraction and 
    29 effective radius in the $S(q)$. 
     18When $X_{core}$ < 1 the core is oblate; when $X_{core}$ > 1 it is prolate. 
     19$X_{core}$ = 1 is a spherical core. 
     20 
     21For a fixed shell thickness $X_{polar shell}$ = 1, to scale $t_{shell}$ 
     22pro-rata with the radius set or constrain $X_{polar shell}$ = $X_{core}$. 
     23 
     24.. note:: 
     25 
     26   When including an $S(q)$, the radius in $S(q)$ is calculated to be that of 
     27   a sphere with the same 2nd virial coefficient of the outer surface of the 
     28   ellipsoid. This may have some undesirable effects if the aspect ratio of the 
     29   ellipsoid is large (ie, if $X << 1$ or $X >> 1$), when the $S(q)$ 
     30   - which assumes spheres - will not in any case be valid.  Generating a 
     31   custom product model will enable separate effective volume fraction and 
     32   effective radius in the $S(q)$. 
    3033 
    3134If SAS data are in absolute units, and the SLDs are correct, then scale should 
     
    4346where 
    4447 
     48.. In following equation SK changed radius\_equat\_core to R_e 
    4549.. math:: 
    4650    :nowrap: 
    4751 
    4852    \begin{align*} 
    49     F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\ 
    50     &+ f(q,radius\_equat\_core + thick\_shell, 
    51          radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha) 
     53    F(q,\alpha) = &f(q,R_e,R_e.x_{core},\alpha) \\ 
     54    &+ f(q,R_e + t_{shell}, 
     55         R_e.x_{core} + t_{shell}.x_{polar shell},\alpha) 
    5256    \end{align*} 
    5357 
     
    7175$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the 
    7276polar radius along the rotational axis of the ellipsoid, $R_e$ is the 
    73 equatorial radius perpendicular to the rotational axis of the ellipsoid 
    74 and $\Delta \rho$ (contrast) is the scattering length density difference, 
    75 either $(sld\_core - sld\_shell)$ or $(sld\_shell - sld\_solvent)$. 
     77equatorial radius perpendicular to the rotational axis of the ellipsoid, 
     78$t_{shell}$ is the thickness of the shell at the equator, 
     79and $\Delta \rho$ (the contrast) is the scattering length density difference, 
     80either $(\rho_{core} - \rho_{shell})$ or $(\rho_{shell} - \rho_{solvent})$. 
    7681 
    7782For randomly oriented particles: 
     
    8893---------- 
    8994see for example: 
    90 Kotlarchyk, M.; Chen, S.-H. J. Chem. Phys., 1983, 79, 2461. 
    91 Berr, S.  J. Phys. Chem., 1987, 91, 4760. 
     95 
     96.. [#] Kotlarchyk, M.; Chen, S.-H. *J. Chem. Phys.*, 1983, 79, 2461 
     97.. [#] Berr, S. *J. Phys. Chem.*, 1987, 91, 4760 
     98 
     99Source 
     100------ 
     101 
     102`core_shell_ellipsoid.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_ellipsoid.py>`_ 
     103 
     104`core_shell_ellipsoid.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_ellipsoid.c>`_ 
    92105 
    93106Authorship and Verification 
     
    96109* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    97110* **Last Modified by:** Richard Heenan (reparametrised model) **Date:** 2015 
    98 * **Last Reviewed by:** Richard Heenan **Date:** October 6, 2016 
     111* **Last Reviewed by:** Steve King **Date:** March 27, 2019 
     112* **Source added by :** Steve King **Date:** March 25, 2019 
    99113""" 
    100114 
     
    146160source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    147161have_Fq = True 
    148 effective_radius_type = [ 
    149     "average outer curvature", "equivalent volume sphere",      
     162radius_effective_modes = [ 
     163    "average outer curvature", "equivalent volume sphere", 
    150164    "min outer radius", "max outer radius", 
    151165    ] 
    152166 
    153167def random(): 
     168    """Return a random parameter set for the model.""" 
    154169    volume = 10**np.random.uniform(5, 12) 
    155170    outer_polar = 10**np.random.uniform(1.3, 4) 
  • sasmodels/models/core_shell_parallelepiped.c

    r99658f6 ra34b811  
    7575 
    7676static double 
    77 effective_radius(int mode, double length_a, double length_b, double length_c, 
     77radius_effective(int mode, double length_a, double length_b, double length_c, 
    7878                 double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    7979{ 
  • sasmodels/models/core_shell_parallelepiped.py

    r99658f6 ra34b811  
    173173   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    174174   =26379>`_ 
    175 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     175.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     176 
     177Source 
     178------ 
     179 
     180`core_shell_parallelepiped.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_parallelepiped.py>`_ 
     181 
     182`core_shell_parallelepiped.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_parallelepiped.c>`_ 
    176183 
    177184Authorship and Verification 
     
    183190* **Last Reviewed by:** Paul Butler **Date:** May 24, 2018 - documentation 
    184191  updated 
     192* **Source added by :** Steve King **Date:** March 25, 2019 
    185193""" 
    186194 
    187195import numpy as np 
    188 from numpy import pi, inf, sqrt, cos, sin 
     196from numpy import inf 
    189197 
    190198name = "core_shell_parallelepiped" 
     
    228236source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 
    229237have_Fq = True 
    230 effective_radius_type = [ 
    231     "equivalent cylinder excluded volume",  
     238radius_effective_modes = [ 
     239    "equivalent cylinder excluded volume", 
    232240    "equivalent volume sphere", 
    233241    "half outer length_a", "half outer length_b", "half outer length_c", 
     
    237245 
    238246def random(): 
     247    """Return a random parameter set for the model.""" 
    239248    outer = 10**np.random.uniform(1, 4.7, size=3) 
    240249    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1 
     
    269278# 2d values not tested against other codes or models 
    270279if 0:  # pak: model rewrite; need to update tests 
    271     qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
     280    qx, qy = 0.2 * np.cos(np.pi/6.), 0.2 * np.sin(np.pi/6.) 
    272281    tests = [[{}, 0.2, 0.533149288477], 
    273282             [{}, [0.2], [0.533149288477]], 
  • sasmodels/models/core_shell_sphere.c

    rd42dd4a ra34b811  
    66 
    77static double 
    8 effective_radius(int mode, double radius, double thickness) 
     8radius_effective(int mode, double radius, double thickness) 
    99{ 
    1010    switch (mode) { 
  • sasmodels/models/core_shell_sphere.py

    r304c775 ra34b811  
    3838effective radius for $S(Q)$ when $P(Q) \cdot S(Q)$ is applied. 
    3939 
    40 References 
    41 ---------- 
    42  
    43 A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*, 
    44 John Wiley and Sons, New York, (1955) 
    45  
    4640Validation 
    4741---------- 
     
    5044the output of the software provided by NIST (Kline, 2006). Figure 1 shows a 
    5145comparison of the output of our model and the output of the NIST software. 
     46 
     47References 
     48---------- 
     49 
     50.. [#] A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*, John Wiley and Sons, New York, (1955) 
     51 
     52Source 
     53------ 
     54 
     55`core_shell_sphere.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_sphere.py>`_ 
     56 
     57`core_shell_sphere.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/core_shell_sphere.c>`_ 
     58 
     59Authorship and Verification 
     60---------------------------- 
     61 
     62* **Author:** 
     63* **Last Modified by:** 
     64* **Last Reviewed by:** 
     65* **Source added by :** Steve King **Date:** March 25, 2019 
    5266""" 
    5367 
     
    7892source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 
    7993have_Fq = True 
    80 effective_radius_type = ["outer radius", "core radius"] 
     94radius_effective_modes = ["outer radius", "core radius"] 
    8195 
    8296demo = dict(scale=1, background=0, radius=60, thickness=10, 
     
    8498 
    8599def random(): 
     100    """Return a random parameter set for the model.""" 
    86101    outer_radius = 10**np.random.uniform(1.3, 4.3) 
    87102    # Use a distribution with a preference for thin shell or thin core 
  • sasmodels/models/correlation_length.py

    ra807206 r0507e09  
    3030---------- 
    3131 
    32 B Hammouda, D L Ho and S R Kline, Insight into Clustering in 
    33 Poly(ethylene oxide) Solutions, Macromolecules, 37 (2004) 6932-6937 
     32.. [#] B Hammouda, D L Ho and S R Kline, Insight into Clustering in Poly(ethylene oxide) Solutions, Macromolecules, 37 (2004) 6932-6937 
     33 
     34Source 
     35------ 
     36 
     37`correlation_length.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/correlation_length.py>`_ 
     38 
     39Authorship and Verification 
     40---------------------------- 
     41 
     42* **Author:**  
     43* **Last Modified by:**  
     44* **Last Reviewed by:**  
     45* **Source added by :** Steve King **Date:** March 25, 2019 
    3446""" 
    3547 
  • sasmodels/models/cylinder.c

    r99658f6 ra34b811  
    3232 
    3333static double 
    34 effective_radius(int mode, double radius, double length) 
     34radius_effective(int mode, double radius, double length) 
    3535{ 
    3636    switch (mode) { 
  • sasmodels/models/cylinder.py

    r99658f6 ra34b811  
    9696---------- 
    9797 
    98 J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997). 
    99 G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951). 
    100 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     98.. [#] J. Pedersen, *Adv. Colloid Interface Sci.*, 70 (1997) 171-210 
     99.. [#] G. Fournet, *Bull. Soc. Fr. Mineral. Cristallogr.*, 74 (1951) 39-113 
     100.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     101 
     102Source 
     103------ 
     104 
     105`cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/cylinder.py>`_ 
     106 
     107`cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/cylinder.c>`_ 
     108 
     109Authorship and Verification 
     110---------------------------- 
     111 
     112* **Author:** 
     113* **Last Modified by:** 
     114* **Last Reviewed by:** 
     115* **Source added by :** Steve King **Date:** March 25, 2019 
    101116""" 
    102117 
     
    140155source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
    141156have_Fq = True 
    142 effective_radius_type = [ 
     157radius_effective_modes = [ 
    143158    "excluded volume", "equivalent volume sphere", "radius", 
    144159    "half length", "half min dimension", "half max dimension", "half diagonal", 
     
    146161 
    147162def random(): 
     163    """Return a random parameter set for the model.""" 
    148164    volume = 10**np.random.uniform(5, 12) 
    149165    length = 10**np.random.uniform(-2, 2)*volume**0.333 
     
    168184            phi_pd=10, phi_pd_n=5) 
    169185 
    170 # pylint: disable=bad-whitespace, line-too-long 
     186# Test 1-D and 2-D models 
    171187qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    172 # After redefinition of angles, find new tests values.  Was 10 10 in old coords 
     188theta, phi = 80.1534480601659, 10.1510817110481  # (10, 10) in sasview 3.x 
    173189tests = [ 
    174190    [{}, 0.2, 0.042761386790780453], 
    175191    [{}, [0.2], [0.042761386790780453]], 
    176     #  new coords 
    177     [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 
    178     [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 
    179     # old coords 
    180     #[{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
    181     #[{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
     192    [{'theta': theta, 'phi': phi}, (qx, qy), 0.03514647218513852], 
     193    [{'theta': theta, 'phi': phi}, [(qx, qy)], [0.03514647218513852]], 
    182194] 
    183 del qx, qy  # not necessary to delete, but cleaner 
    184  
    185 # Default radius and length 
    186 radius, length = parameters[2][2], parameters[3][2] 
    187 tests.extend([ 
    188     ({'radius_effective_mode': 0}, 0.1, None, None, 0., pi*radius**2*length, 1.0),    
    189     ({'radius_effective_mode': 1}, 0.1, None, None, 0.5*(0.75*radius*(2.0*radius*length + (radius + length)*(pi*radius + length)))**(1./3.), None, None),     
    190     ({'radius_effective_mode': 2}, 0.1, None, None, (0.75*radius**2*length)**(1./3.), None, None), 
    191     ({'radius_effective_mode': 3}, 0.1, None, None, radius, None, None), 
    192     ({'radius_effective_mode': 4}, 0.1, None, None, length/2., None, None), 
    193     ({'radius_effective_mode': 5}, 0.1, None, None, min(radius, length/2.), None, None), 
    194     ({'radius_effective_mode': 6}, 0.1, None, None, max(radius, length/2.), None, None), 
    195     ({'radius_effective_mode': 7}, 0.1, None, None, np.sqrt(4*radius**2 + length**2)/2., None, None), 
    196 ]) 
    197 del radius, length 
    198 # pylint: enable=bad-whitespace, line-too-long 
     195del qx, qy, theta, phi  # not necessary to delete, but cleaner 
     196 
     197def _extend_with_reff_tests(radius, length): 
     198    """Test R_eff and form volume calculations""" 
     199    # V and Vr are the same for each R_eff mode 
     200    V = pi*radius**2*length  # shell volume = form volume for solid objects 
     201    Vr = 1.0  # form:shell volume ratio 
     202    # Use test value for I(0.2) from above to check Fsq value.  Need to 
     203    # remove scale and background before testing. 
     204    q = 0.2 
     205    scale, background = V, 0.001 
     206    Fsq = (0.042761386790780453 - background)*scale 
     207    F = None  # Need target value for <F> 
     208    # Various values for R_eff, depending on mode 
     209    r_effs = [ 
     210        0., 
     211        0.5*(0.75*radius*(2.0*radius*length 
     212                          + (radius + length)*(pi*radius + length)))**(1./3.), 
     213        (0.75*radius**2*length)**(1./3.), 
     214        radius, 
     215        length/2., 
     216        min(radius, length/2.), 
     217        max(radius, length/2.), 
     218        np.sqrt(4*radius**2 + length**2)/2., 
     219    ] 
     220    tests.extend([ 
     221        ({'radius_effective_mode': 0}, q, F, Fsq, r_effs[0], V, Vr), 
     222        ({'radius_effective_mode': 1}, q, F, Fsq, r_effs[1], V, Vr), 
     223        ({'radius_effective_mode': 2}, q, F, Fsq, r_effs[2], V, Vr), 
     224        ({'radius_effective_mode': 3}, q, F, Fsq, r_effs[3], V, Vr), 
     225        ({'radius_effective_mode': 4}, q, F, Fsq, r_effs[4], V, Vr), 
     226        ({'radius_effective_mode': 5}, q, F, Fsq, r_effs[5], V, Vr), 
     227        ({'radius_effective_mode': 6}, q, F, Fsq, r_effs[6], V, Vr), 
     228        ({'radius_effective_mode': 7}, q, F, Fsq, r_effs[7], V, Vr), 
     229    ]) 
     230 
     231# Test Reff and volume with default model parameters 
     232_extend_with_reff_tests(parameters[2][2], parameters[3][2]) 
     233del _extend_with_reff_tests 
    199234 
    200235# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc 
  • sasmodels/models/dab.py

    r304c775 r0507e09  
    3131---------- 
    3232 
    33 P Debye, H R Anderson, H Brumberger, *Scattering by an Inhomogeneous Solid. II. 
    34 The Correlation Function and its Application*, *J. Appl. Phys.*, 28(6) (1957) 679 
     33.. [#] P Debye, H R Anderson, H Brumberger, *Scattering by an Inhomogeneous Solid. II. The Correlation Function and its Application*, *J. Appl. Phys.*, 28(6) (1957) 679 
     34.. [#] P Debye, A M Bueche, *Scattering by an Inhomogeneous Solid*, *J. Appl. Phys.*, 20 (1949) 518 
    3535 
    36 P Debye, A M Bueche, *Scattering by an Inhomogeneous Solid*, *J. Appl. Phys.*, 
    37 20 (1949) 518 
     36Source 
     37------ 
    3838 
    39 *2013/09/09 - Description reviewed by King, S and Parker, P.* 
     39`dab.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/dab.py>`_ 
     40 
     41Authorship and Verification 
     42---------------------------- 
     43 
     44* **Author:**  
     45* **Last Modified by:**  
     46* **Last Reviewed by:** Steve King & Peter Parker **Date:** September 09, 2013 
     47* **Source added by :** Steve King **Date:** March 25, 2019 
    4048""" 
    4149 
     
    6674 
    6775def random(): 
     76    """Return a random parameter set for the model.""" 
    6877    pars = dict( 
    6978        scale=10**np.random.uniform(1, 4), 
    7079        cor_length=10**np.random.uniform(0.3, 3), 
    71         #background = 0, 
     80#        background = 0, 
    7281    ) 
    7382    pars['scale'] /= pars['cor_length']**3 
  • sasmodels/models/ellipsoid.c

    r99658f6 ra34b811  
    3232 
    3333static double 
    34 effective_radius(int mode, double radius_polar, double radius_equatorial) 
     34radius_effective(int mode, double radius_polar, double radius_equatorial) 
    3535{ 
    3636    switch (mode) { 
  • sasmodels/models/ellipsoid.py

    r99658f6 ra34b811  
    6767 
    6868The $\theta$ and $\phi$ parameters are not used for the 1D output. 
    69  
    70  
    7169 
    7270Validation 
     
    108106---------- 
    109107 
    110 L A Feigin and D I Svergun. 
    111 *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 
    112 Plenum Press, New York, 1987. 
    113  
    114 A. Isihara. J. Chem. Phys. 18(1950) 1446-1449 
     108.. [#] L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum Press, New York, 1987 
     109.. [#] A. Isihara. *J. Chem. Phys.*, 18 (1950) 1446-1449 
     110 
     111Source 
     112------ 
     113 
     114`ellipsoid.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/ellipsoid.py>`_ 
     115 
     116`ellipsoid.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/ellipsoid.c>`_ 
    115117 
    116118Authorship and Verification 
     
    120122* **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014 
    121123* **Last Modified by:** Paul Kienzle **Date:** March 22, 2017 
     124* **Source added by :** Steve King **Date:** March 25, 2019 
    122125""" 
    123126from __future__ import division 
     
    125128import numpy as np 
    126129from numpy import inf, sin, cos, pi 
    127  
    128 try: 
    129     from numpy import cbrt 
    130 except ImportError: 
    131     def cbrt(x): return x ** (1.0/3.0) 
    132130 
    133131name = "ellipsoid" 
     
    169167source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    170168have_Fq = True 
    171 effective_radius_type = [ 
     169radius_effective_modes = [ 
    172170    "average curvature", "equivalent volume sphere", "min radius", "max radius", 
    173171    ] 
    174172 
    175173def random(): 
     174    """Return a random parameter set for the model.""" 
    176175    volume = 10**np.random.uniform(5, 12) 
    177176    radius_polar = 10**np.random.uniform(1.3, 4) 
  • sasmodels/models/elliptical_cylinder.c

    r99658f6 ra34b811  
    4141 
    4242static double 
    43 effective_radius(int mode, double radius_minor, double r_ratio, double length) 
     43radius_effective(int mode, double radius_minor, double r_ratio, double length) 
    4444{ 
    4545    switch (mode) { 
  • sasmodels/models/elliptical_cylinder.py

    r99658f6 ra34b811  
    8686---------- 
    8787 
    88 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    89 Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] 
    90 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
     88.. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] 
     89.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 
     90 
     91Source 
     92------ 
     93 
     94`elliptical_cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/elliptical_cylinder.py>`_ 
     95 
     96`elliptical_cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/elliptical_cylinder.c>`_ 
    9197 
    9298Authorship and Verification 
     
    96102* **Last Modified by:** 
    97103* **Last Reviewed by:**  Richard Heenan - corrected equation in docs **Date:** December 21, 2016 
     104* **Source added by :** Steve King **Date:** March 25, 2019 
    98105""" 
    99106 
     
    124131source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    125132have_Fq = True 
    126 effective_radius_type = [ 
    127     "equivalent cylinder excluded volume", "equivalent volume sphere", "average radius", "min radius", "max radius", 
     133radius_effective_modes = [ 
     134    "equivalent cylinder excluded volume", 
     135    "equivalent volume sphere", "average radius", "min radius", "max radius", 
    128136    "equivalent circular cross-section", 
    129137    "half length", "half min dimension", "half max dimension", "half diagonal", 
     
    135143 
    136144def random(): 
     145    """Return a random parameter set for the model.""" 
    137146    # V = pi * radius_major * radius_minor * length; 
    138147    volume = 10**np.random.uniform(3, 9) 
    139148    length = 10**np.random.uniform(1, 3) 
    140149    axis_ratio = 10**np.random.uniform(0, 2) 
    141     radius_minor = np.sqrt(volume/length/axis_ratio) 
     150    radius_minor = sqrt(volume/length/axis_ratio) 
    142151    volfrac = 10**np.random.uniform(-4, -2) 
    143152    pars = dict( 
     
    156165 
    157166tests = [ 
    158 #    [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 
    159 #    [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 
     167    #[{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 
     168    #[{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 
    160169 
    161170    # The SasView test result was 0.00169, with a background of 0.001 
  • sasmodels/models/fcc_paracrystal.py

    rda7b26b r0507e09  
    33#note - calculation requires double precision 
    44r""" 
    5 .. warning:: This model and this model description are under review following  
    6              concerns raised by SasView users. If you need to use this model,  
    7              please email help@sasview.org for the latest situation. *The  
     5.. warning:: This model and this model description are under review following 
     6             concerns raised by SasView users. If you need to use this model, 
     7             please email help@sasview.org for the latest situation. *The 
    88             SasView Developers. September 2018.* 
    99 
     
    9999   (Corrections to FCC and BCC lattice structure calculation) 
    100100 
     101Source 
     102------ 
     103 
     104`fcc_paracrystal.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/fcc_paracrystal.py>`_ 
     105 
     106`fcc_paracrystal.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/fcc_paracrystal.c>`_ 
     107 
    101108Authorship and Verification 
    102109--------------------------- 
     
    105112* **Last Modified by:** Paul Butler **Date:** September 29, 2016 
    106113* **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 
     114* **Source added by :** Steve King **Date:** March 25, 2019 
    107115""" 
    108116 
     
    137145 
    138146def random(): 
     147    """Return a random parameter set for the model.""" 
    139148    # copied from bcc_paracrystal 
    140149    radius = 10**np.random.uniform(1.3, 4) 
  • sasmodels/models/flexible_cylinder.py

    r2d81cfe r830cf6b  
    3030The Kuhn length $(b = 2*l_p)$ is also used to describe the stiffness of a chain. 
    3131 
    32 The returned value is in units of $cm^{-1}$, on absolute scale. 
    33  
    3432In the parameters, the sld and sld\_solvent represent the SLD of the cylinder 
    3533and solvent respectively. 
    3634 
    37 Our model uses the form factor calculations implemented in a c-library provided 
    38 by the NIST Center for Neutron Research (Kline, 2006). 
    39  
    40  
    41 From the reference: 
     35Our model uses the form factor calculations in reference [1] as implemented in a 
     36c-library provided by the NIST Center for Neutron Research (Kline, 2006). This states: 
    4237 
    4338    'Method 3 With Excluded Volume' is used. 
     
    4742    See equations (13,26-27) in the original reference for the details. 
    4843 
     44.. note:: 
     45 
     46    There are several typos in the original reference that have been corrected 
     47    by WRC [2]. Details of the corrections are in the reference below. Most notably 
     48 
     49    - Equation (13): the term $(1 - w(QR))$ should swap position with $w(QR)$ 
     50 
     51    - Equations (23) and (24) are incorrect; WRC has entered these into 
     52      Mathematica and solved analytically. The results were then converted to 
     53      code. 
     54 
     55    - Equation (27) should be $q0 = max(a3/(Rg^2)^{1/2},3)$ instead of 
     56      $max(a3*b(Rg^2)^{1/2},3)$ 
     57 
     58    - The scattering function is negative for a range of parameter values and 
     59      q-values that are experimentally accessible. A correction function has been 
     60      added to give the proper behavior. 
     61 
     62 
     63**This is a model with complex behaviour depending on the ratio of** $L/b$ **and the 
     64reader is strongly encouraged to read reference [1] before use.** 
     65 
     66.. note:: 
     67 
     68    There are several typos in the original reference that have been corrected 
     69    by WRC [2]. Details of the corrections are in the reference below. Most notably 
     70 
     71    - Equation (13): the term $(1 - w(QR))$ should swap position with $w(QR)$ 
     72 
     73    - Equations (23) and (24) are incorrect; WRC has entered these into 
     74      Mathematica and solved analytically. The results were then converted to 
     75      code. 
     76 
     77    - Equation (27) should be $q0 = max(a3/(Rg^2)^{1/2},3)$ instead of 
     78      $max(a3*b(Rg^2)^{1/2},3)$ 
     79 
     80    - The scattering function is negative for a range of parameter values and 
     81      q-values that are experimentally accessible. A correction function has been 
     82      added to give the proper behavior. 
     83 
     84 
     85**This is a model with complex behaviour depending on the ratio of** $L/b$ **and the 
     86reader is strongly encouraged to read reference [1] before use.** 
     87 
    4988References 
    5089---------- 
    5190 
    52 J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible 
    53 polymers with and without excluded volume effects.* Macromolecules, 
    54 29 (1996) 7602-7612 
     91.. [#] J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612 
    5592 
    5693Correction of the formula can be found in 
    5794 
    58 W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions 
    59 in the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 
    60 22(15) 2006 6539-6548 
     95.. [#] W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 22(15) 2006 6539-6548 
     96 
     97Source 
     98------ 
     99 
     100`flexible_cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/flexible_cylinder.py>`_ 
     101 
     102`flexible_cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/flexible_cylinder.c>`_ 
     103 
     104`wrc_cyl.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/lib/wrc_cyl.c>`_ 
     105 
     106Authorship and Verification 
     107---------------------------- 
     108 
     109* **Author:** 
     110* **Last Modified by:** 
     111* **Last Reviewed by:** Steve King **Date:** March 26, 2019 
     112* **Source added by :** Steve King **Date:** March 25, 2019 
    61113""" 
    62114 
     
    65117 
    66118name = "flexible_cylinder" 
    67 title = "Flexible cylinder where the form factor is normalized by the volume" \ 
     119title = "Flexible cylinder where the form factor is normalized by the volume " \ 
    68120        "of the cylinder." 
    69121description = """Note : scale and contrast = (sld - sld_solvent) are both 
     
    89141 
    90142def random(): 
     143    """Return a random parameter set for the model.""" 
    91144    length = 10**np.random.uniform(2, 6) 
    92145    radius = 10**np.random.uniform(1, 3) 
  • sasmodels/models/flexible_cylinder_elliptical.py

    r2d81cfe rdb1d9d5  
    44The non-negligible diameter of the cylinder is included by accounting 
    55for excluded volume interactions within the walk of a single cylinder. 
     6**Inter-cylinder interactions are NOT provided for.** 
     7 
    68The form factor is normalized by the particle volume such that 
    79 
     
    2426----------- 
    2527 
    26 The function calculated in a similar way to that for the flexible_cylinder model 
    27 from the reference given below using the author's "Method 3 With Excluded Volume". 
     28The function is calculated in a similar way to that for the 
     29:ref:`flexible-cylinder` model in reference [1] below using the author's 
     30"Method 3 With Excluded Volume". 
     31 
    2832The model is a parameterization of simulations of a discrete representation of 
    2933the worm-like chain model of Kratky and Porod applied in the pseudo-continuous 
     
    3337 
    3438    There are several typos in the original reference that have been corrected 
    35     by WRC. Details of the corrections are in the reference below. Most notably 
     39    by WRC [2]. Details of the corrections are in the reference below. Most notably 
    3640 
    3741    - Equation (13): the term $(1 - w(QR))$ should swap position with $w(QR)$ 
     
    4145      code. 
    4246 
    43     - Equation (27) should be $q0 = max(a3/sqrt(RgSquare),3)$ instead of 
    44       $max(a3*b/sqrt(RgSquare),3)$ 
     47    - Equation (27) should be $q0 = max(a3/(Rg^2)^{1/2},3)$ instead of 
     48      $max(a3*b(Rg^2)^{1/2},3)$ 
    4549 
    4650    - The scattering function is negative for a range of parameter values and 
     
    5862 
    5963The cross section of the cylinder is elliptical, with minor radius $a$ . 
    60 The major radius is larger, so of course, **the axis ratio (parameter 5) must be 
     64The major radius is larger, so of course, **the axis_ratio must be 
    6165greater than one.** Simple constraints should be applied during curve fitting to 
    6266maintain this inequality. 
    63  
    64 The returned value is in units of $cm^{-1}$, on absolute scale. 
    6567 
    6668In the parameters, the $sld$ and $sld\_solvent$ represent the SLD of the 
     
    6971these parameters must be held fixed during model fitting. 
    7072 
    71 **No inter-cylinder interference effects are included in this calculation.** 
     73**This is a model with complex behaviour depending on the ratio of** $L/b$ **and the 
     74reader is strongly encouraged to read reference [1] before use.** 
    7275 
    7376References 
    7477---------- 
    7578 
    76 J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers 
    77 with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612 
     79.. [#] J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612 
    7880 
    7981Correction of the formula can be found in 
    8082 
    81 W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in 
    82 the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 
    83 22(15) 2006 6539-6548 
     83.. [#] W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 22(15) 2006 6539-6548 
     84 
     85Source 
     86------ 
     87 
     88`flexible_cylinder_elliptical.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/flexible_cylinder_elliptical.py>`_ 
     89 
     90`flexible_cylinder_elliptical.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/flexible_cylinder_elliptical.c>`_ 
     91 
     92`wrc_cyl.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/lib/wrc_cyl.c>`_ 
     93 
     94Authorship and Verification 
     95---------------------------- 
     96 
     97* **Author:** 
     98* **Last Modified by:** Richard Heenan **Date:** December, 2016 
     99* **Last Reviewed by:** Steve King **Date:** March 26, 2019 
     100* **Source added by :** Steve King **Date:** March 25, 2019 
    84101""" 
    85102 
     
    115132 
    116133def random(): 
     134    """Return a random parameter set for the model.""" 
    117135    length = 10**np.random.uniform(2, 6) 
    118136    radius = 10**np.random.uniform(1, 3) 
  • sasmodels/models/fractal.py

    r2d81cfe r0507e09  
    4646.. [#] J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 
    4747 
     48Source 
     49------ 
     50 
     51`fractal.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/fractal.py>`_ 
     52 
     53`fractal.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/fractal.c>`_ 
     54 
    4855Authorship and Verification 
    4956---------------------------- 
     
    5360* **Last Modified by:** Paul Butler **Date:** March 12, 2017 
    5461* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 
     62* **Source added by :** Steve King **Date:** March 25, 2019 
    5563""" 
    5664from __future__ import division 
     
    100108 
    101109def random(): 
     110    """Return a random parameter set for the model.""" 
    102111    radius = 10**np.random.uniform(0.7, 4) 
    103112    #radius = 5 
  • sasmodels/models/fractal_core_shell.py

    r2cc8aa2 r0507e09  
    5151.. [#Kline]  S R Kline, *J Appl. Cryst.*, 39 (2006) 895 
    5252 
     53Source 
     54------ 
     55 
     56`fractal_core_shell.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/fractal_core_shell.py>`_ 
     57 
     58`fractal_core_shell.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/fractal_core_shell.c>`_ 
     59 
    5360Authorship and Verification 
    5461---------------------------- 
     
    5764* **Last Modified by:** Paul Butler and Paul Kienzle **Date:** November 27, 2016 
    5865* **Last Reviewed by:** Paul Butler and Paul Kienzle **Date:** November 27, 2016 
     66* **Source added by :** Steve King **Date:** March 25, 2019 
    5967""" 
    6068 
    6169import numpy as np 
    62 from numpy import pi, inf 
     70from numpy import inf 
    6371 
    6472name = "fractal_core_shell" 
     
    98106 
    99107def random(): 
     108    """Return a random parameter set for the model.""" 
    100109    outer_radius = 10**np.random.uniform(0.7, 4) 
    101110    # Use a distribution with a preference for thin shell or thin core 
     
    126135            cor_length=100.0) 
    127136 
     137# TODO: why is there an ER function here? 
    128138def ER(radius, thickness): 
    129139    """ 
     
    136146#tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    137147tests = [ 
    138 #         # At some point the SasView 3.x test result was deemed incorrect. The 
    139           #following tests were verified against NIST IGOR macros ver 7.850. 
    140           #NOTE: NIST macros do only provide for a polydispers core (no option 
    141           #for a poly shell or for a monodisperse core.  The results seemed 
    142           #extremely sensitive to the core PD, varying non monotonically all 
    143           #the way to a PD of 1e-6. From 1e-6 to 1e-9 no changes in the 
    144           #results were observed and the values below were taken using PD=1e-9. 
    145           #Non-monotonically = I(0.001)=188 to 140 to 177 back to 160 etc. 
    146     &n