Changeset 1662ebe in sasmodels


Ignore:
Timestamp:
May 2, 2018 5:02:29 PM (2 years ago)
Author:
Paul Kienzle <pkienzle@…>
Children:
765d025
Parents:
6f91c91 (diff), 33969b6 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge remote-tracking branch 'upstream/master'

Files:
9 added
40 edited
2 moved

Legend:

Unmodified
Added
Removed
  • .travis.yml

    r2d09df1 r335271e  
    3737- conda update --yes conda 
    3838- conda info -a 
    39 - conda install --yes python=$PY numpy scipy cython mako cffi 
     39- conda install --yes python=$PY numpy scipy matplotlib docutils setuptools pytest 
    4040install: 
    4141- pip install bumps 
    42 - pip install unittest-xml-reporting 
     42- pip install unittest-xml-reporting tinycc 
    4343script: 
    4444- python --version 
    45 - python -m sasmodels.model_test -v dll all 
     45#- python -m sasmodels.model_test -v dll all 
     46#- python -m pytest -v --cache-clear 
     47- python setup.py test --pytest-args -v 
    4648before_deploy: 
    4749- echo -e "Host danse.chem.utk.edu\n\tStrictHostKeyChecking no\n" >> ~/.ssh/config 
  • appveyor.yml

    r1258e32 r335271e  
    44  # /E:ON and /V:ON options are not enabled in the batch script interpreter 
    55  # See: http://stackoverflow.com/a/13751649/163740 
    6   CMD_IN_ENV: "cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd" 
     6  # 2018-01-19 PAK: probably irrelevant now that we are using tinycc rather than msvc 
     7  #CMD_IN_ENV: "cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd" 
    78 
    89  # Workaround for https://github.com/conda/conda-build/issues/636 
     
    2223    - TARGET_ARCH: "x86" 
    2324      CONDA_PY: "35" 
    24       PY_CONDITION: "python >=3.5" 
     25      PY_CONDITION: "python >=3.5,<3.6" 
    2526      CONDA_INSTALL_LOCN: "C:\\Miniconda35" 
    2627    - TARGET_ARCH: "x64" 
     
    3435    - TARGET_ARCH: "x64" 
    3536      CONDA_PY: "35" 
    36       PY_CONDITION: "python >=3.5" 
     37      PY_CONDITION: "python >=3.5,<3.6" 
    3738      CONDA_INSTALL_LOCN: "C:\\Miniconda35-x64" 
    3839 
     
    4344 
    4445install: 
    45     # Set the CONDA_NPY, although it has no impact on the actual build. We need this because of a test within conda-build. 
     46    # Set the CONDA_NPY, although it has no impact on the actual build. 
     47    # We need this because of a test within conda-build. 
    4648    - cmd: set CONDA_NPY=19 
    4749 
    4850    # Remove cygwin (and therefore the git that comes with it). 
    49     - cmd: rmdir C:\cygwin /s /q 
     51    # 2018-01-19 PAK: irrelevant since we already pulled the repo 
     52    #- cmd: rmdir C:\cygwin /s /q 
     53 
     54    # Set the conda path; would be nice to do this 
     55    - cmd: path %CONDA_INSTALL_LOCN%\Scripts;%PATH% 
    5056 
    5157    # Use the pre-installed Miniconda for the desired arch 
     
    5662    # so that we can update it. Then we remove it so that 
    5763    # we can do a proper activation. 
    58     - cmd: set "OLDPATH=%PATH%" 
    59     - cmd: set "PATH=%CONDA_INSTALL_LOCN%\\Scripts;%CONDA_INSTALL_LOCN%\\Library\\bin;%PATH%" 
    6064    - cmd: conda update --yes --quiet conda python 
    61     - cmd: set "PATH=%OLDPATH%" 
    6265    - cmd: call %CONDA_INSTALL_LOCN%\Scripts\activate.bat 
    63  
    64     - cmd: conda config --add channels conda-forge 
    65     - cmd: conda config --set show_channel_urls true 
    66     - cmd: conda update --yes --quiet conda 
    67     - cmd: conda install --yes --quiet obvious-ci 
    68     - cmd: conda install --yes --quiet numpy toolchain scipy cython cffi 
     66    #- cmd: conda config --add channels conda-forge 
     67    #- cmd: conda config --set show_channel_urls true 
     68    #- cmd: conda install --yes --quiet obvious-ci 
     69    # 2018-01-19 PAK: skipping toolchain cython and cffi (these were for pyopencl?) 
     70    - cmd: conda install --yes --quiet numpy scipy matplotlib docutils setuptools pytest 
     71    # 2018-01-19 PAK: skipping pyopencl; this would be needed for deploy but maybe not for test 
    6972    #- cmd: conda install --yes --channel conda-forge pyopencl 
     73    # 2018-01-19 PAK: 3rd party packages might need msvc, so %CMD_IN_ENV% may be needed for pip 
    7074    - cmd: pip install bumps unittest-xml-reporting tinycc 
    7175 
    7276build_script: 
    7377    # Build the project 
    74     - "%CMD_IN_ENV% python setup.py build" 
     78    # 2018-01-19 PAK: maybe need one of this if using msvc? 
     79    #- "%CMD_IN_ENV% python setup.py build" 
     80    #- cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd python setup.py build 
     81    - python setup.py build 
    7582 
    7683test_script: 
    7784    # Run the project tests 
    78     - "%CMD_IN_ENV% python -m sasmodels.model_test dll all" 
     85    # 2018-01-19 PAK: maybe need one of this if using msvc? 
     86    #- "%CMD_IN_ENV% python -m sasmodels.model_test dll all" 
     87    #- cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd python -m sasmoels.model_test dll all 
     88    #- python -m sasmodels.model_test dll all 
     89    #- nosetests -v sasmodels/*.py 
     90    #- python -m pytest -v --cache-clear 
     91    - python setup.py test --pytest-args -v 
  • doc/developer/overview.rst

    r2a7e20e r0d5a655  
    190190*sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for 
    191191equirectangular and *PROJECTION=2* for sinusoidal.  The more complicated 
    192 Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh 
     192Guyou and Postel projections are not implemented. See jitter.draw_mesh 
    193193for details. 
    194194 
     
    274274code. 
    275275 
    276 *explore/jitter.py* is for exploring different options for handling 
    277 orientation and orientation dispersity.  It uses *explore/guyou.py* to 
     276*sasmodels/jitter.py* is for exploring different options for handling 
     277orientation and orientation dispersity.  It uses *sasmodels/guyou.py* to 
    278278generate the Guyou projection. 
    279279 
  • doc/genapi.py

    r706f466 r0d5a655  
    6969    ('exception', 'Annotate exceptions'), 
    7070    ('generate', 'Model parser'), 
     71    ('jitter', 'Orientation explorer'), 
     72    ('guyou', 'Guyou map projection'), 
    7173    ('kernel', 'Evaluator type definitions'), 
    7274    ('kernelcl', 'OpenCL model evaluator'), 
  • doc/guide/orientation/orientation.rst

    r5fb0634 r0d5a655  
    5252yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the 
    5353$a$-axis. 
     54 
     55You can explore the view and jitter angles interactively using 
     56:func:`sasmodels.jitter.run`.  Enter the following into the python 
     57interpreter:: 
     58 
     59    from sasmodels import jitter 
     60    jitter.run() 
    5461 
    5562More formally, starting with axes $a$-$b$-$c$ of the particle aligned 
  • doc/guide/pd/polydispersity.rst

    reda8b30 r29afc50  
    2323average over the size distribution. 
    2424 
    25 Each distribution is characterized by its center $\bar x$, its width $\sigma$, 
    26 the number of sigmas $N_\sigma$ to include from the tails, and the number of 
    27 points used to compute the average. The center of the distribution is set by the 
    28 value of the model parameter.  Volume parameters have polydispersity *PD* 
    29 (not to be confused with a molecular weight distributions in polymer science) 
    30 leading to a size distribution of width $\text{PD} = \sigma / \bar x$, but 
    31 orientation parameters use an angular distributions of width $\sigma$. 
    32 $N_\sigma$ determines how far into the tails to evaluate the distribution, with 
    33 larger values of $N_\sigma$ required for heavier tailed distributions. 
     25Each distribution is characterized by a center value $\bar x$ or 
     26$x_\text{med}$, a width parameter $\sigma$ (note this is *not necessarily* 
     27the standard deviation, so read the description carefully), the number of 
     28sigmas $N_\sigma$ to include from the tails of the distribution, and the 
     29number of points used to compute the average. The center of the distribution 
     30is set by the value of the model parameter. The meaning of a polydispersity  
     31parameter *PD* (not to be confused with a molecular weight distributions  
     32in polymer science) in a model depends on the type of parameter it is being  
     33applied too. 
     34 
     35The distribution width applied to *volume* (ie, shape-describing) parameters  
     36is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$.  
     37However, the distribution width applied to *orientation* (ie, angle-describing)  
     38parameters is just $\sigma = \mathrm{PD}$. 
     39 
     40$N_\sigma$ determines how far into the tails to evaluate the distribution, 
     41with larger values of $N_\sigma$ required for heavier tailed distributions. 
    3442The scattering in general falls rapidly with $qr$ so the usual assumption 
    3543that $G(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)G(r - 3\sigma_r)$ 
     
    4250calculations are generally more robust with more data points or more angles. 
    4351 
    44 The following five distribution functions are provided: 
    45  
     52The following distribution functions are provided: 
     53 
     54*  *Uniform Distribution* 
    4655*  *Rectangular Distribution* 
    4756*  *Gaussian Distribution* 
     57*  *Boltzmann Distribution* 
    4858*  *Lognormal Distribution* 
    4959*  *Schulz Distribution* 
     
    5262These are all implemented as *number-average* distributions. 
    5363 
     64Additional distributions are under consideration. 
     65 
     66Suggested Applications 
     67^^^^^^^^^^^^^^^^^^^^^^ 
     68 
     69If applying polydispersion to parameters describing particle sizes, use 
     70the Lognormal or Schulz distributions. 
     71 
     72If applying polydispersion to parameters describing interfacial thicknesses 
     73or angular orientations, use the Gaussian or Boltzmann distributions. 
     74 
     75If applying polydispersion to parameters describing angles, use the Uniform  
     76distribution. Beware of using distributions that are always positive (eg, the  
     77Lognormal) because angles can be negative! 
     78 
     79The array distribution allows a user-defined distribution to be applied. 
     80 
     81.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     82 
     83Uniform Distribution 
     84^^^^^^^^^^^^^^^^^^^^ 
     85 
     86The Uniform Distribution is defined as 
     87 
     88.. math:: 
     89 
     90    f(x) = \frac{1}{\text{Norm}} 
     91    \begin{cases} 
     92        1 & \text{for } |x - \bar x| \leq \sigma \\ 
     93        0 & \text{for } |x - \bar x| > \sigma 
     94    \end{cases} 
     95 
     96where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 
     97distribution, $\sigma$ is the half-width, and *Norm* is a normalization 
     98factor which is determined during the numerical calculation. 
     99 
     100The polydispersity in sasmodels is given by 
     101 
     102.. math:: \text{PD} = \sigma / \bar x 
     103 
     104.. figure:: pd_uniform.jpg 
     105 
     106    Uniform distribution. 
     107 
     108The value $N_\sigma$ is ignored for this distribution. 
     109 
    54110.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    55111 
     
    63119    f(x) = \frac{1}{\text{Norm}} 
    64120    \begin{cases} 
    65       1 & \text{for } |x - \bar x| \leq w \\ 
    66       0 & \text{for } |x - \bar x| > w 
     121        1 & \text{for } |x - \bar x| \leq w \\ 
     122        0 & \text{for } |x - \bar x| > w 
    67123    \end{cases} 
    68124 
    69 where $\bar x$ is the mean of the distribution, $w$ is the half-width, and 
    70 *Norm* is a normalization factor which is determined during the numerical 
    71 calculation. 
     125where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 
     126distribution, $w$ is the half-width, and *Norm* is a normalization 
     127factor which is determined during the numerical calculation. 
    72128 
    73129Note that the standard deviation and the half width $w$ are different! 
     
    77133.. math:: \sigma = w / \sqrt{3} 
    78134 
    79 whilst the polydispersity is 
     135whilst the polydispersity in sasmodels is given by 
    80136 
    81137.. math:: \text{PD} = \sigma / \bar x 
     
    84140 
    85141    Rectangular distribution. 
     142 
     143.. note:: The Rectangular Distribution is deprecated in favour of the 
     144            Uniform Distribution above and is described here for backwards 
     145            compatibility with earlier versions of SasView only. 
    86146 
    87147.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     
    95155 
    96156    f(x) = \frac{1}{\text{Norm}} 
    97            \exp\left(-\frac{(x - \bar x)^2}{2\sigma^2}\right) 
    98  
    99 where $\bar x$ is the mean of the distribution and *Norm* is a normalization 
    100 factor which is determined during the numerical calculation. 
    101  
    102 The polydispersity is 
     157            \exp\left(-\frac{(x - \bar x)^2}{2\sigma^2}\right) 
     158 
     159where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 
     160distribution and *Norm* is a normalization factor which is determined 
     161during the numerical calculation. 
     162 
     163The polydispersity in sasmodels is given by 
    103164 
    104165.. math:: \text{PD} = \sigma / \bar x 
     
    107168 
    108169    Normal distribution. 
     170 
     171.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     172 
     173Boltzmann Distribution 
     174^^^^^^^^^^^^^^^^^^^^^^ 
     175 
     176The Boltzmann Distribution is defined as 
     177 
     178.. math:: 
     179 
     180    f(x) = \frac{1}{\text{Norm}} 
     181            \exp\left(-\frac{ | x - \bar x | }{\sigma}\right) 
     182 
     183where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 
     184distribution and *Norm* is a normalization factor which is determined 
     185during the numerical calculation. 
     186 
     187The width is defined as 
     188 
     189.. math:: \sigma=\frac{k T}{E} 
     190 
     191which is the inverse Boltzmann factor, where $k$ is the Boltzmann constant, 
     192$T$ the temperature in Kelvin and $E$ a characteristic energy per particle. 
     193 
     194.. figure:: pd_boltzmann.jpg 
     195 
     196    Boltzmann distribution. 
    109197 
    110198.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     
    113201^^^^^^^^^^^^^^^^^^^^^^ 
    114202 
     203The Lognormal Distribution describes a function of $x$ where $\ln (x)$ has 
     204a normal distribution. The result is a distribution that is skewed towards 
     205larger values of $x$. 
     206 
    115207The Lognormal Distribution is defined as 
    116208 
    117209.. math:: 
    118210 
    119     f(x) = \frac{1}{\text{Norm}} 
    120            \frac{1}{xp}\exp\left(-\frac{(\ln(x) - \mu)^2}{2p^2}\right) 
    121  
    122 where $\mu=\ln(x_\text{med})$ when $x_\text{med}$ is the median value of the 
    123 distribution, and *Norm* is a normalization factor which will be determined 
    124 during the numerical calculation. 
    125  
    126 The median value for the distribution will be the value given for the 
    127 respective size parameter, for example, *radius=60*. 
    128  
    129 The polydispersity is given by $\sigma$ 
    130  
    131 .. math:: \text{PD} = p 
    132  
    133 For the angular distribution 
    134  
    135 .. math:: p = \sigma / x_\text{med} 
    136  
    137 The mean value is given by $\bar x = \exp(\mu+ p^2/2)$. The peak value 
    138 is given by $\max x = \exp(\mu - p^2)$. 
     211    f(x) = \frac{1}{\text{Norm}}\frac{1}{x\sigma} 
     212            \exp\left(-\frac{1}{2} 
     213                        \bigg(\frac{\ln(x) - \mu}{\sigma}\bigg)^2\right) 
     214 
     215where *Norm* is a normalization factor which will be determined during 
     216the numerical calculation, $\mu=\ln(x_\text{med})$ and $x_\text{med}$ 
     217is the *median* value of the *lognormal* distribution, but $\sigma$ is 
     218a parameter describing the width of the underlying *normal* distribution. 
     219 
     220$x_\text{med}$ will be the value given for the respective size parameter 
     221in sasmodels, for example, *radius=60*. 
     222 
     223The polydispersity in sasmodels is given by 
     224 
     225.. math:: \text{PD} = \sigma = p / x_\text{med} 
     226 
     227The mean value of the distribution is given by $\bar x = \exp(\mu+ \sigma^2/2)$ 
     228and the peak value by $\max x = \exp(\mu - \sigma^2)$. 
     229 
     230The variance (the square of the standard deviation) of the *lognormal* 
     231distribution is given by 
     232 
     233.. math:: 
     234 
     235    \nu = [\exp({\sigma}^2) - 1] \exp({2\mu + \sigma^2}) 
     236 
     237Note that larger values of PD might need a larger number of points 
     238and $N_\sigma$. 
    139239 
    140240.. figure:: pd_lognormal.jpg 
    141241 
    142     Lognormal distribution. 
    143  
    144 This distribution function spreads more, and the peak shifts to the left, as 
    145 $p$ increases, so it requires higher values of $N_\sigma$ and more points 
    146 in the distribution. 
     242    Lognormal distribution for PD=0.1. 
     243 
     244For further information on the Lognormal distribution see: 
     245http://en.wikipedia.org/wiki/Log-normal_distribution and 
     246http://mathworld.wolfram.com/LogNormalDistribution.html 
    147247 
    148248.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     
    151251^^^^^^^^^^^^^^^^^^^ 
    152252 
     253The Schulz (sometimes written Schultz) distribution is similar to the 
     254Lognormal distribution, in that it is also skewed towards larger values of 
     255$x$, but which has computational advantages over the Lognormal distribution. 
     256 
    153257The Schulz distribution is defined as 
    154258 
    155259.. math:: 
    156260 
    157     f(x) = \frac{1}{\text{Norm}} 
    158            (z+1)^{z+1}(x/\bar x)^z\frac{\exp[-(z+1)x/\bar x]}{\bar x\Gamma(z+1)} 
    159  
    160 where $\bar x$ is the mean of the distribution and *Norm* is a normalization 
    161 factor which is determined during the numerical calculation, and $z$ is a 
    162 measure of the width of the distribution such that 
     261    f(x) = \frac{1}{\text{Norm}} (z+1)^{z+1}(x/\bar x)^z 
     262            \frac{\exp[-(z+1)x/\bar x]}{\bar x\Gamma(z+1)} 
     263 
     264where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 
     265distribution, *Norm* is a normalization factor which is determined 
     266during the numerical calculation, and $z$ is a measure of the width 
     267of the distribution such that 
    163268 
    164269.. math:: z = (1-p^2) / p^2 
    165270 
    166 The polydispersity is 
    167  
    168 .. math:: p = \sigma / \bar x 
    169  
    170 Note that larger values of PD might need larger number of points and $N_\sigma$. 
    171 For example, at PD=0.7 and radius=60 |Ang|, Npts>=160 and Nsigmas>=15 at least. 
     271where $p$ is the polydispersity in sasmodels given by 
     272 
     273.. math:: PD = p = \sigma / \bar x 
     274 
     275and $\sigma$ is the RMS deviation from $\bar x$. 
     276 
     277Note that larger values of PD might need a larger number of points 
     278and $N_\sigma$. For example, for PD=0.7 with radius=60 |Ang|, at least 
     279Npts>=160 and Nsigmas>=15 are required. 
    172280 
    173281.. figure:: pd_schulz.jpg 
     
    176284 
    177285For further information on the Schulz distribution see: 
    178 M Kotlarchyk & S-H Chen, *J Chem Phys*, (1983), 79, 2461. 
     286M Kotlarchyk & S-H Chen, *J Chem Phys*, (1983), 79, 2461 and 
     287M Kotlarchyk, RB Stephens, and JS Huang, *J Phys Chem*, (1988), 92, 1533 
    179288 
    180289.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     
    183292^^^^^^^^^^^^^^^^^^ 
    184293 
    185 This user-definable distribution should be given as as a simple ASCII text 
     294This user-definable distribution should be given as a simple ASCII text 
    186295file where the array is defined by two columns of numbers: $x$ and $f(x)$. 
    187296The $f(x)$ will be normalized to 1 during the computation. 
     
    202311given for the model will have no affect, and will be ignored when computing 
    203312the average.  This means that any parameter with an array distribution will 
    204 not be fittable. 
     313not be fitable. 
    205314 
    206315.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     
    216325related) except when the DLS polydispersity parameter is <0.13. 
    217326 
     327.. math:: 
     328 
     329    p_{DLS} = \sqrt(\nu / \bar x^2) 
     330 
     331where $\nu$ is the variance of the distribution and $\bar x$ is the mean 
     332value of $x$. 
     333 
    218334For more information see: 
    219335S King, C Washington & R Heenan, *Phys Chem Chem Phys*, (2005), 7, 143 
     
    225341| 2015-05-01 Steve King 
    226342| 2017-05-08 Paul Kienzle 
     343| 2018-03-20 Steve King 
     344| 2018-04-04 Steve King 
  • example/SiO2_100pc_H2O_0pc_D2O.ses

    r31c64d9 rb314f89  
    1 name    SiO2 stober                                              
    2 sample  SiO2 100pcH2O                                            
    3 collimation     unknown                                          
    4 ID      CPD                                              
    5 Date    di 27 jan 2015 18:17:50                                          
    6 Method  Offspec                                          
    7 Thickness [mm]  1                                                
    8 z-acceptance [nm^-1]    0.3                                              
    9 y-acceptance [nm^-1]    0.3                                              
     1FileFormatVersion       1.0 
     2DataFileTitle           SiO2 stober 
     3Sample                  SiO2 100pcH2O 
     4collimation             unknown 
     5ID                      CPD 
     6Date                    di 27 jan 2015 18:17:50 
     7Method                  Offspec 
     8Thickness               1 
     9Thickness_unit          mm 
     10Theta_zmax              0.0100 
     11Theta_zmax_unit         radians 
     12Theta_ymax              0.0100 
     13Theta_ymax_unit         radians 
     14Orientation             Z 
     15SpinEchoLength_unit     A 
     16Depolarisation_unit     A-2 cm-1 
     17Wavelength_unit         A 
    1018 
    11 spin echo length [nm]   error spin-echo length  wavelength [nm] error wavelength        Sample P        Sample P -error          
    12 244.098 5       0.220907401     0.006627222     0.96321 0.0113507       0.029882023     29.882023 
    13 340.136 5       0.260768096     0.007823043     0.973305        0.0034846       0.085257531     85.257531 
    14 452.181 5       0.300665994     0.00901998      0.951588        0.00287498      0.140798634     140.798634 
    15 580.232 5       0.340587727     0.010217632     0.911935        0.00310965      0.196560502     196.560502 
    16 724.29  5       0.380526057     0.011415782     0.856129        0.00356166      0.252528125     252.528125 
    17 884.354 5       0.420476016     0.01261428      0.78817 0.0036651       0.308729274     308.729274 
    18 1060.42 5       0.460433709     0.013813011     0.719608        0.00386121      0.365037117     365.037117 
    19 1252.5  5       0.50039968      0.01501199      0.653424        0.00418035      0.421684677     421.684677 
    20 1460.58 5       0.540369504     0.016211085     0.603456        0.00472218      0.478564852     478.564852 
    21 1684.67 5       0.580344105     0.017410323     0.558174        0.00537231      0.535715521     535.715521 
    22 1924.77 5       0.620322561     0.018609677     0.532672        0.00617699      0.59324408      593.24408 
    23 2180.87 5       0.660302658     0.01980908      0.505018        0.00707792      0.650985498     650.985498 
    24 2452.98 5       0.700285542     0.021008566     0.475536        0.00829962      0.709203958     709.203958 
    25 2741.1  5       0.740270761     0.022208123     0.440819        0.00957178      0.767809943     767.809943 
    26 3045.22 5       0.780256676     0.0234077       0.397642        0.0111225       0.826762925     826.762925 
    27 3365.35 5       0.820244402     0.024607332     0.363717        0.0132055       0.886302296     886.302296 
    28 3701.48 5       0.860232527     0.025806976     0.302206        0.0158338       0.946236366     946.236366 
    29 4053.62 5       0.900222106     0.027006663     0.302188        0.0204351       1.006784115     1006.784115 
    30 4421.77 5       0.940212955     0.028206389     0.263252        0.0264158       1.067891213     1067.891213 
    31 4805.92 5       0.980203897     0.029406117     0.230149        0.0339972       1.129598722     1129.598722 
    32 5206.08 5       1.020195903     0.030605877     0.268749        0.0480167       1.192022694     1192.022694 
    33 5622.25 5       1.060188851     0.031805666     0.340439        0.0701945       1.255271282     1255.271282 
    34 6054.42 5       1.10018173      0.033005452     0.238197        0.0986631       1.319307134     1319.307134 
    35 6502.6  5       1.140175425     0.034205263     0.0985987       0.170273        1.384197821     1384.197821 
    36 6966.79 5       1.180169852     0.035405096     0.0814573       0.256765        1.448518618     1448.518618 
     19BEGIN_DATA 
     20 SpinEchoLength  Depolarisation  Depolarisation_error   Wavelength  Wavelength_error  SpinEchoLength_error  Polarisation  Polarisation_error  
     21 2440.98   -7.6810986e-3          2.4147998e-3   2.20907401        0.06627222                    50       0.96321           0.0113507  
     22 3401.36   -3.9790857e-3          5.2649599e-4   2.60768096        0.07823043                    50      0.973305           0.0034846  
     23 4521.81   -5.4892798e-3          3.3420831e-4   3.00665994         0.0901998                    50      0.951588          0.00287498  
     24 5802.32   -7.9471175e-3          2.9396095e-4   3.40587727        0.10217632                    50      0.911935          0.00310965  
     25 7242.9    -0.010727495          2.8730584e-4   3.80526057        0.11415782                    50      0.856129          0.00356166  
     26 8843.54    -0.013463878          2.6301679e-4   4.20476016         0.1261428                    50       0.78817           0.0036651  
     27 10604.2    -0.015521222          2.5310062e-4   4.60433709        0.13813011                    50      0.719608          0.00386121  
     2812525.    -0.016993983          2.5549565e-4    5.0039968         0.1501199                    50      0.653424          0.00418035  
     2914605.8    -0.017297381          2.6798795e-4   5.40369504        0.16211085                    50      0.603456          0.00472218  
     3016846.7    -0.017312523          2.8577242e-4   5.80344105        0.17410323                    50      0.558174          0.00537231  
     3119247.7    -0.016368225          3.0135741e-4   6.20322561        0.18609677                    50      0.532672          0.00617699  
     3221808.7    -0.015668849          3.2144946e-4   6.60302658         0.1980908                    50      0.505018          0.00707792  
     3324529.8    -0.015157278          3.5589713e-4   7.00285542        0.21008566                    50      0.475536          0.00829962  
     3427411.    -0.014947440          3.9623352e-4   7.40270761        0.22208123                    50      0.440819          0.00957178  
     3530452.2    -0.015147872          4.5944674e-4   7.80256676          0.234077                    50      0.397642           0.0111225  
     3633653.5    -0.015032370          5.3964070e-4   8.20244402        0.24607332                    50      0.363717           0.0132055  
     3737014.8    -0.016170897          7.0802787e-4   8.60232527        0.25806976                    50      0.302206           0.0158338  
     3840536.2    -0.014766858          8.3444978e-4   9.00222106        0.27006663                    50      0.302188           0.0204351  
     3944217.7    -0.015097771          1.1351144e-3   9.40212955        0.28206389                    50      0.263252           0.0264158  
     4048059.2    -0.015289642          1.5374507e-3   9.80203897        0.29406117                    50      0.230149           0.0339972  
     4152060.8    -0.012624691          1.7166363e-3  10.20195903        0.30605877                    50      0.268749           0.0480167  
     4256222.5   -9.5864674e-3          1.8344138e-3  10.60188851        0.31805666                    50      0.340439           0.0701945  
     4360544.2    -0.011852755          3.4220757e-3   11.0018173        0.33005452                    50      0.238197           0.0986631  
     4465026.    -0.017820748           0.013284073  11.40175425        0.34205263                    50     0.0985987            0.170273  
     4569667.9    -0.018004557           0.022631679  11.80169852        0.35405096                    50     0.0814573            0.256765  
  • example/core_shell.ses

    r5be92e8 re0721e7  
    1 name    se009051.dat                                             
    2 sample  "SiO2 pH=6 & NiO pH=10, nanosized, in D2O, Sine ++ only."                                                
    3 collimation     "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. "                                                
    4 ID      CPD                                              
    5 Date    di 27 jan 2015 18:17:50                                          
    6 Method  sine one element scan                                            
    7 Thickness [mm]  2                                                
    8 z-acceptance [nm^-1]    0.3                                              
    9 y-acceptance [nm^-1]    0.3                                              
     1FileFormatVersion       1.0 
     2DataFileTitle           se009051.dat                                             
     3Sample                  "SiO2 pH=6 & NiO pH=10, nanosized, in D2O, Sine ++ only."                                                
     4collimation             "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. "                                                
     5ID                      CPD                                              
     6Date                    di 27 jan 2015 18:17:50 
     7Method                  sine one element scan 
     8Thickness               10 
     9Thickness_unit          mm 
     10Theta_zmax              0.0168 
     11Theta_zmax_unit         radians 
     12Theta_ymax              0.0168 
     13Theta_ymax_unit         radians 
     14Orientation             Z 
     15SpinEchoLength_unit     A 
     16Depolarisation_unit     A-2 cm-1 
     17Wavelength_unit         A 
    1018                                                         
    11 spin echo length [um]   error spin-echo length  wavelength [nm] error wavelength        Sample 5        Sample 5 -error          
    12 26      5       0.161245        0.01    0.996315        0.00528445      0.161245        26 
    13 28.08   5       0.167571        0.01    0.995937        0.00522962      0.167571        28.08 
    14 30.3264 5       0.174145        0.01    0.995039        0.00371164      0.174145        30.3264 
    15 32.7525 5       0.180977        0.01    0.99448 0.00381089      0.180977        32.7525 
    16 35.3727 5       0.188076        0.01    0.992143        0.00319011      0.188076        35.3727 
    17 38.2025 5       0.195455        0.01    0.991403        0.0032496       0.195455        38.2025 
    18 41.2587 5       0.203122        0.01    0.990106        0.002978        0.203122        41.2587 
    19 44.5594 5       0.211091        0.01    0.98857 0.00297338      0.211091        44.5594 
    20 48.1242 5       0.219372        0.01    0.986236        0.00286669      0.219372        48.1242 
    21 51.9741 5       0.227978        0.01    0.982739        0.00283988      0.227978        51.9741 
    22 56.132  5       0.236922        0.01    0.979995        0.00287147      0.236922        56.132 
    23 60.6226 5       0.246217        0.01    0.976516        0.00290667      0.246217        60.6226 
    24 65.4724 5       0.255876        0.01    0.971858        0.00300969      0.255876        65.4724 
    25 70.7102 5       0.265914        0.01    0.967518        0.00310597      0.265914        70.7102 
    26 76.367  5       0.276346        0.01    0.963198        0.00310409      0.276346        76.367 
    27 82.4764 5       0.287187        0.01    0.957071        0.00306271      0.287187        82.4764 
    28 89.0745 5       0.298454        0.01    0.952784        0.00299587      0.298454        89.0745 
    29 96.2005 5       0.310162        0.01    0.947462        0.0029896       0.310162        96.2005 
    30 103.897 5       0.322331        0.01    0.941226        0.00300217      0.322331        103.897 
    31 112.208 5       0.334975        0.01    0.936151        0.00305677      0.334975        112.208 
    32 121.185 5       0.348116        0.01    0.930884        0.00307969      0.348116        121.185 
    33 130.88  5       0.361773        0.01    0.927992        0.00316999      0.361773        130.88 
    34 141.35  5       0.375965        0.01    0.926511        0.00326096      0.375965        141.35 
    35 152.658 5       0.390715        0.01    0.925098        0.00340349      0.390715        152.658 
    36 164.871 5       0.406043        0.01    0.924369        0.0035662       0.406043        164.871 
    37 178.06  5       0.421972        0.01    0.925632        0.00372489      0.421972        178.06 
    38 192.305 5       0.438526        0.01    0.922762        0.00390394      0.438526        192.305 
    39 207.69  5       0.45573 0.01    0.918753        0.00412162      0.45573 207.69 
    40 224.305 5       0.473608        0.01    0.90663 0.0043316       0.473608        224.305 
    41 242.249 5       0.492188        0.01    0.893251        0.00454417      0.492188        242.249 
    42 261.629 5       0.511497        0.01    0.878299        0.00473364      0.511497        261.629 
    43 282.559 5       0.531563        0.01    0.867266        0.00497459      0.531563        282.559 
    44 305.164 5       0.552417        0.01    0.859537        0.00521977      0.552417        305.164 
    45 329.577 5       0.574088        0.01    0.855935        0.0055573       0.574088        329.577 
    46 355.943 5       0.59661 0.01    0.851389        0.00581848      0.59661 355.943 
    47 384.419 5       0.620015        0.01    0.836311        0.00620483      0.620015        384.419 
    48 415.172 5       0.644338        0.01    0.821362        0.00668373      0.644338        415.172 
    49 448.386 5       0.669616        0.01    0.803833        0.0071725       0.669616        448.386 
    50 484.257 5       0.695886        0.01    0.792211        0.00786449      0.695886        484.257 
    51 522.998 5       0.723186        0.01    0.787407        0.00860688      0.723186        522.998 
    52 564.838 5       0.751557        0.01    0.765977        0.00938856      0.751557        564.838 
    53 610.025 5       0.781041        0.01    0.746618        0.0101567       0.781041        610.025 
    54 658.827 5       0.811682        0.01    0.734672        0.0110357       0.811682        658.827 
    55 711.533 5       0.843524        0.01    0.720351        0.0118912       0.843524        711.533 
    56 768.455 5       0.876616        0.01    0.692587        0.0130927       0.876616        768.455 
    57 829.932 5       0.911006        0.01    0.674767        0.0137763       0.911006        829.932 
    58 896.326 5       0.946745        0.01    0.664497        0.0148902       0.946745        896.326 
    59 968.032 5       0.983886        0.01    0.634464        0.0158459       0.983886        968.032 
    60 1045.47 5       1.02248 0.01    0.604792        0.0166925       1.02248 1045.47 
    61 1129.11 5       1.0626  0.01    0.591534        0.0179667       1.0626  1129.11 
    62 1219.44 5       1.10428 0.01    0.582625        0.0198223       1.10428 1219.44 
    63 1317    5       1.14761 0.01    0.549232        0.0217725       1.14761 1317 
    64 1422.36 5       1.19263 0.01    0.53863 0.024516        1.19263 1422.36 
    65 1536.15 5       1.23942 0.01    0.521353        0.0281108       1.23942 1536.15 
    66 1659.04 5       1.28804 0.01    0.460045        0.0314452       1.28804 1659.04 
    67 1791.76 5       1.33857 0.01    0.411185        0.0366899       1.33857 1791.76 
    68 1930.34 5       1.38937 0.01    0.414685        0.0478689       1.38937 1930.34 
     19BEGIN_DATA 
     20SpinEchoLength  SpinEchoLength_error  Wavelength  Wavelength_Error  Polarisation  Polarisation_error  Depolarisation  Depolarisation_error  
     21298.82023 50 2.095 0.1 0.997817556 0.00887342 -4.9779370e-5 2.0261512e-4  
     22852.57531 50 2.095 0.1 1.0026184 0.009572307 5.9579929e-5 2.1752686e-4  
     231407.9863 50 2.095 0.1 0.996013172 0.012119508 -9.1017859e-5 2.7723742e-4  
     241965.6050 50 2.095 0.1 0.991752656 0.011651473 -1.8868750e-4 2.6767598e-4  
     252525.2813 50 2.095 0.1 0.995429571 0.012193294 -1.0437182e-4 2.7908883e-4  
     263087.2927 50 2.095 0.1 0.995119819 0.009621277 -1.1146275e-4 2.2028721e-4  
     273650.3712 50 2.095 0.1 0.997497767 0.012518842 -5.7082583e-5 2.8594610e-4  
     284216.8468 50 2.095 0.1 0.994733865 0.010585926 -1.2030120e-4 2.4246770e-4  
     294785.6485 50 2.095 0.1 0.992696275 0.010677724 -1.6701950e-4 2.4507231e-4  
     305357.1552 50 2.095 0.1 0.994534294 0.012909556 -1.2487278e-4 2.9574914e-4  
     315932.4408 50 2.095 0.1 0.986412268 0.014651894 -3.1170682e-4 3.3842875e-4  
     326509.8550 50 2.095 0.1 0.988788566 0.015736468 -2.5688520e-4 3.6260666e-4  
     337092.0396 50 2.095 0.1 0.989501168 0.013611615 -2.4047103e-4 3.1341898e-4  
     347678.0994 50 2.095 0.1 0.99203933 0.016501891 -1.8210252e-4 3.7899787e-4  
     358267.6293 50 2.095 0.1 0.994062723 0.010276207 -1.3567871e-4 2.3553259e-4  
     368863.0230 50 2.095 0.1 0.979647192 0.011089421 -4.6850452e-4 2.5791174e-4  
     379462.3637 50 2.095 0.1 0.97982075 0.013168527 -4.6446835e-4 3.0621222e-4  
     3810067.841 50 2.095 0.1 0.977927712 0.014362819 -5.0853039e-4 3.3463000e-4  
     3910678.912 50 2.095 0.1 0.978351462 0.013239476 -4.9865985e-4 3.0832436e-4  
     4011295.987 50 2.095 0.1 0.981246471 0.013896898 -4.3133968e-4 3.2267975e-4  
     4111920.227 50 2.095 0.1 0.978275164 0.013747481 -5.0043677e-4 3.2017989e-4  
     4212552.713 50 2.095 0.1 0.976818842 0.013487944 -5.3437989e-4 3.1460359e-4  
     4313193.071 50 2.095 0.1 0.981546068 0.011972943 -4.2438423e-4 2.7792152e-4  
     4413841.978 50 2.095 0.1 0.96445854 0.013357572 -8.2452102e-4 3.1555561e-4  
     4514500.513 50 2.095 0.1 0.972756158 0.014756189 -6.2933879e-4 3.4562263e-4  
     4615169.786 50 2.095 0.1 0.971010184 0.014210504 -6.7027011e-4 3.3343996e-4  
     4715850.446 50 2.095 0.1 0.975009829 0.013648839 -5.7661387e-4 3.1894710e-4  
     4816543.562 50 2.095 0.1 0.969575209 0.014454336 -7.0396574e-4 3.3966327e-4  
     4917249.754 50 2.095 0.1 0.959907267 0.012806944 -9.3229353e-4 3.0398222e-4  
     5017969.627 50 2.095 0.1 0.96219203 0.011845496 -8.7812744e-4 2.8049391e-4  
     5118706.130 50 2.095 0.1 0.960157966 0.011721155 -9.2634378e-4 2.7813758e-4  
     5219459.212 50 2.095 0.1 0.95089937 0.009967287 -1.1471121e-3 2.3882201e-4  
     5320230.984 50 2.095 0.1 0.955584642 0.01104368 -1.0351259e-3 2.6331561e-4  
     5421021.899 50 2.095 0.1 0.950786938 0.011651582 -1.1498062e-3 2.7921171e-4  
     5521835.345 50 2.095 0.1 0.951403607 0.008784149 -1.1350335e-3 2.1036178e-4  
     5622671.730 50 2.095 0.1 0.94561823 0.009392383 -1.2740040e-3 2.2630383e-4  
     5723533.930 50 2.095 0.1 0.945882787 0.009468264 -1.2676305e-3 2.2806833e-4  
     5824423.504 50 2.095 0.1 0.939583158 0.009164239 -1.4198814e-3 2.2222511e-4  
     5925343.343 50 2.095 0.1 0.93731257 0.008621568 -1.4750079e-3 2.0957224e-4  
     6026294.734 50 2.095 0.1 0.936639557 0.010408821 -1.4913733e-3 2.5319842e-4  
     6127282.004 50 2.095 0.1 0.932091884 0.007112674 -1.6022666e-3 1.7386258e-4  
     6228307.499 50 2.095 0.1 0.933852763 0.009200431 -1.5592642e-3 2.2447176e-4  
     6329373.630 50 2.095 0.1 0.930283437 0.01022719 -1.6465153e-3 2.5047996e-4  
     6430485.110 50 2.095 0.1 0.930121334 0.007055747 -1.6504858e-3 1.7283645e-4  
     6531645.899 50 2.095 0.1 0.925944326 0.0078917 -1.7530356e-3 1.9418588e-4  
     6632858.438 50 2.095 0.1 0.92351759 0.005522582 -1.8128270e-3 1.3624763e-4  
     6734128.058 50 2.095 0.1 0.918228311 0.006917328 -1.9436940e-3 1.7164045e-4  
     6835460.168 50 2.095 0.1 0.908621563 0.006852313 -2.1833230e-3 1.7182490e-4  
     6936858.762 50 2.095 0.1 0.919087792 0.006995961 -1.9223776e-3 1.7342924e-4  
     7038331.748 50 2.095 0.1 0.910132214 0.006107336 -2.1454742e-3 1.5289007e-4  
     7139882.345 50 2.095 0.1 0.904491462 0.007078101 -2.2871233e-3 1.7829708e-4  
     7241519.404 50 2.095 0.1 0.90056919 0.005934428 -2.3861400e-3 1.5013907e-4  
     7343249.023 50 2.095 0.1 0.903155701 0.006249069 -2.3207959e-3 1.5764661e-4  
     7445079.637 50 2.095 0.1 0.894919827 0.006196232 -2.5295172e-3 1.5775222e-4  
     7547019.670 50 2.095 0.1 0.896358908 0.007324999 -2.4929085e-3 1.8619052e-4  
     7649077.851 50 2.095 0.1 0.88807923 0.004530465 -2.7043436e-3 1.1623128e-4  
     7751264.711 50 2.095 0.1 0.886649039 0.008924412 -2.7410654e-3 2.2932944e-4  
     7853590.598 50 2.095 0.1 0.882051947 0.006055081 -2.8595036e-3 1.5640756e-4  
     7956067.403 50 2.095 0.1 0.882573371 0.005727419 -2.8460388e-3 1.4785638e-4  
     8058709.036 50 2.095 0.1 0.876289825 0.006062306 -3.0088321e-3 1.5762389e-4  
     8161527.362 50 2.095 0.1 0.877318199 0.005307783 -2.9821094e-3 1.3784403e-4  
     8264538.932 50 2.095 0.1 0.873610284 0.006918104 -3.0786086e-3 1.8042690e-4  
     8367758.354 50 2.095 0.1 0.865954455 0.005625863 -3.2791557e-3 1.4802192e-4  
     8471205.306 50 2.095 0.1 0.868037308 0.008662141 -3.2244196e-3 2.2736248e-4  
     8574897.359 50 2.095 0.1 0.86676404 0.005761209 -3.2578647e-3 1.5144143e-4  
     8678855.560 50 2.095 0.1 0.85655072 0.005731411 -3.5279304e-3 1.5245456e-4  
     8783102.902 50 2.095 0.1 0.857936274 0.006954572 -3.4911046e-3 1.8469168e-4  
     8887663.154 50 2.095 0.1 0.85520525 0.006559622 -3.5637478e-3 1.7475934e-4  
     8992562.717 50 2.095 0.1 0.849739829 0.008907377 -3.7098230e-3 2.3883381e-4  
     9097830.316 50 2.095 0.1 0.841487161 0.006547659 -3.9321836e-3 1.7728439e-4  
     91103497.14 50 2.095 0.1 0.84801051 0.007574866 -3.7562386e-3 2.0351933e-4  
     92109596.79 50 2.095 0.1 0.843342861 0.006411993 -3.8819940e-3 1.7322909e-4  
     93116165.97 50 2.095 0.1 0.837488454 0.008141353 -4.0407107e-3 2.2148775e-4  
     94123244.47 50 2.095 0.1 0.839214922 0.009987535 -3.9937900e-3 2.7115465e-4  
     95130874.72 50 2.095 0.1 0.834336972 0.011467406 -4.1266093e-3 3.1315233e-4  
     96139102.96 50 2.095 0.1 0.828775634 0.012483365 -4.2789870e-3 3.4318369e-4  
     97147980.51 50 2.095 0.1 0.823088812 0.016035848 -4.4358638e-3 4.4389186e-4  
     98157560.98 50 2.095 0.1 0.824329178 0.014838138 -4.4015548e-3 4.1011975e-4  
     99167904.68 50 2.095 0.1 0.82023172 0.015286187 -4.5150892e-3 4.2461424e-4  
     100179076.29 50 2.095 0.1 0.827559752 0.012612741 -4.3124376e-3 3.4724985e-4  
  • example/fit.py

    rf4b36fa r1ee5ac6  
    179179 
    180180else: 
    181     print "No parameters for %s"%name 
     181    print("No parameters for %s"%name) 
    182182    sys.exit(1) 
    183183 
     
    192192else: 
    193193   problem = FitProblem(M) 
     194 
     195if __name__ == "__main__": 
     196   problem.plot() 
     197   import pylab; pylab.show() 
  • example/se008724_01.ses

    r42ade62 rb314f89  
    1 DataFileTitle   Polystyrene of Markus Strobl,  Full Sine, ++ only  
    2 Sample  Polystyrene 2 um in 53% H2O, 47% D2O  
    3 Settings        D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast  
    4 Operator        CPD  
    5 Date    ma 7 jul 2014 18:54:43  
    6 ScanType        sine one element scan  
    7 Thickness [mm]  2  
    8 Q_zmax [\AA^-1]         0.05  
    9 Q_ymax [\AA^-1]         0.05  
    10    
    11 spin echo length [nm]    error SEL       wavelength [nm]         error wavelength        polarisation    error pol  
    12 49.778  2.4889  0.211   0.01055 0.99782 0.0044367 
    13 63.041  3.152   0.211   0.01055 1.0026  0.0047862 
    14 76.487  3.8244  0.211   0.01055 0.99601 0.0060598 
    15 89.847  4.4924  0.211   0.01055 0.99175 0.0058257 
    16 103.41  5.1705  0.211   0.01055 0.99543 0.0060966 
    17 116.95  5.8475  0.211   0.01055 0.99512 0.0048106 
    18 130.61  6.5303  0.211   0.01055 0.9975  0.0062594 
    19 144.37  7.2184  0.211   0.01055 0.99473 0.005293 
    20 158.2   7.9102  0.211   0.01055 0.9927  0.0053389 
    21 172.12  8.6062  0.211   0.01055 0.99453 0.0064548 
    22 186.17  9.3087  0.211   0.01055 0.98641 0.0073259 
    23 200.28  10.014  0.211   0.01055 0.98879 0.0078682 
    24 214.46  10.723  0.211   0.01055 0.9895  0.0068058 
    25 228.73  11.436  0.211   0.01055 0.99204 0.0082509 
    26 243.12  12.156  0.211   0.01055 0.99406 0.0051381 
    27 257.55  12.878  0.211   0.01055 0.97965 0.0055447 
    28 272.22  13.611  0.211   0.01055 0.97982 0.0065843 
    29 286.9   14.345  0.211   0.01055 0.97793 0.0071814 
    30 301.73  15.087  0.211   0.01055 0.97835 0.0066197 
    31 316.75  15.838  0.211   0.01055 0.98125 0.0069484 
    32 331.82  16.591  0.211   0.01055 0.97828 0.0068737 
    33 347.16  17.358  0.211   0.01055 0.97682 0.006744 
    34 362.45  18.122  0.211   0.01055 0.98155 0.0059865 
    35 378.09  18.904  0.211   0.01055 0.96446 0.0066788 
    36 393.74  19.687  0.211   0.01055 0.97276 0.0073781 
    37 409.61  20.481  0.211   0.01055 0.97101 0.0071053 
    38 425.55  21.278  0.211   0.01055 0.97501 0.0068244 
    39 441.91  22.096  0.211   0.01055 0.96958 0.0072272 
    40 458.12  22.906  0.211   0.01055 0.95991 0.0064035 
    41 474.77  23.739  0.211   0.01055 0.96219 0.0059227 
    42 491.52  24.576  0.211   0.01055 0.96016 0.0058606 
    43 508.51  25.426  0.211   0.01055 0.9509  0.0049836 
    44 525.68  26.284  0.211   0.01055 0.95558 0.0055218 
    45 543.08  27.154  0.211   0.01055 0.95079 0.0058258 
    46 560.74  28.037  0.211   0.01055 0.9514  0.0043921 
    47 578.69  28.935  0.211   0.01055 0.94562 0.0046962 
    48 596.75  29.838  0.211   0.01055 0.94588 0.0047341 
    49 615.17  30.758  0.211   0.01055 0.93958 0.0045821 
    50 633.77  31.688  0.211   0.01055 0.93731 0.0043108 
    51 652.82  32.641  0.211   0.01055 0.93664 0.0052044 
    52 672.04  33.602  0.211   0.01055 0.93209 0.0035563 
    53 691.57  34.578  0.211   0.01055 0.93385 0.0046002 
    54 711.48  35.574  0.211   0.01055 0.93028 0.0051136 
    55 731.67  36.583  0.211   0.01055 0.93012 0.0035279 
    56 752.17  37.608  0.211   0.01055 0.92594 0.0039458 
    57 773.1   38.655  0.211   0.01055 0.92352 0.0027613 
    58 794.42  39.721  0.211   0.01055 0.91823 0.0034587 
    59 816.14  40.807  0.211   0.01055 0.90862 0.0034262 
    60 838.19  41.91   0.211   0.01055 0.91909 0.003498 
    61 860.69  43.034  0.211   0.01055 0.91013 0.0030537 
    62 883.74  44.187  0.211   0.01055 0.90449 0.0035391 
    63 907.04  45.352  0.211   0.01055 0.90057 0.0029672 
    64 930.95  46.547  0.211   0.01055 0.90316 0.0031245 
    65 955.38  47.769  0.211   0.01055 0.89492 0.0030981 
    66 980.25  49.013  0.211   0.01055 0.89636 0.0036625 
    67 1005.7  50.284  0.211   0.01055 0.88808 0.0022652 
    68 1031.7  51.586  0.211   0.01055 0.88665 0.0044622 
    69 1058.3  52.914  0.211   0.01055 0.88205 0.0030275 
    70 1085.5  54.276  0.211   0.01055 0.88257 0.0028637 
    71 1113.4  55.669  0.211   0.01055 0.87629 0.0030312 
    72 1141.9  57.093  0.211   0.01055 0.87732 0.0026539 
    73 1171    58.55   0.211   0.01055 0.87361 0.0034591 
    74 1201    60.052  0.211   0.01055 0.86595 0.0028129 
    75 1231.6  61.582  0.211   0.01055 0.86804 0.0043311 
    76 1263    63.151  0.211   0.01055 0.86676 0.0028806 
    77 1295.2  64.762  0.211   0.01055 0.85655 0.0028657 
    78 1328.4  66.419  0.211   0.01055 0.85794 0.0034773 
    79 1362.3  68.113  0.211   0.01055 0.85521 0.0032798 
    80 1397.1  69.854  0.211   0.01055 0.84974 0.0044537 
    81 1432.9  71.644  0.211   0.01055 0.84149 0.0032738 
    82 1469.7  73.487  0.211   0.01055 0.84801 0.0037874 
    83 1507.5  75.376  0.211   0.01055 0.84334 0.003206 
    84 1546.4  77.318  0.211   0.01055 0.83749 0.0040707 
    85 1586.4  79.322  0.211   0.01055 0.83921 0.0049938 
    86 1627.6  81.381  0.211   0.01055 0.83434 0.0057337 
    87 1669.9  83.497  0.211   0.01055 0.82878 0.0062417 
    88 1713.6  85.678  0.211   0.01055 0.82309 0.0080179 
    89 1758.5  87.927  0.211   0.01055 0.82433 0.0074191 
    90 1804.8  90.239  0.211   0.01055 0.82023 0.0076431 
    91 1852.5  92.627  0.211   0.01055 0.82756 0.0063064 
    92 1901.7  95.087  0.211   0.01055 0.82584 0.0052944 
    93 1952.4  97.622  0.211   0.01055 0.82799 0.0050662 
    94 2004.8  100.24  0.211   0.01055 0.82345 0.0043127 
    95 2058.8  102.94  0.211   0.01055 0.82296 0.0038 
    96 2114.5  105.72  0.211   0.01055 0.81987 0.0034072 
    97 2172    108.6   0.211   0.01055 0.82021 0.0036752 
    98 2231.4  111.57  0.211   0.01055 0.82765 0.0028851 
    99 2292.8  114.64  0.211   0.01055 0.82664 0.0038942 
    100 2356.2  117.81  0.211   0.01055 0.82702 0.0047371 
    101 2421.8  121.09  0.211   0.01055 0.81593 0.0043772 
    102 2489.6  124.48  0.211   0.01055 0.8251  0.0028026 
    103 2559.5  127.98  0.211   0.01055 0.8292  0.0024574 
    104 2631.9  131.59  0.211   0.01055 0.82626 0.0036198 
    105 2706.8  135.34  0.211   0.01055 0.8208  0.0032314 
    106 2784.3  139.22  0.211   0.01055 0.81959 0.0042731 
    107 2864.5  143.23  0.211   0.01055 0.82653 0.002699 
    108 2947.5  147.38  0.211   0.01055 0.82401 0.0036726 
    109 3033.5  151.67  0.211   0.01055 0.82361 0.0048224 
    110 3122.4  156.12  0.211   0.01055 0.82358 0.0041221 
    111 3214.5  160.73  0.211   0.01055 0.82187 0.0028807 
    112 3310.1  165.5   0.211   0.01055 0.82644 0.003516 
    113 3409    170.45  0.211   0.01055 0.82355 0.0021166 
    114 3511.4  175.57  0.211   0.01055 0.82513 0.0033911 
    115 3617.6  180.88  0.211   0.01055 0.82802 0.0015342 
    116 3727.6  186.38  0.211   0.01055 0.82663 0.0029222 
    117 3841.7  192.08  0.211   0.01055 0.82026 0.0020755 
    118 3960    198     0.211   0.01055 0.83079 0.0026394 
    119 4082.7  204.13  0.211   0.01055 0.82665 0.0027466 
    120 4209.9  210.5   0.211   0.01055 0.82774 0.0025199 
    121 4341.9  217.09  0.211   0.01055 0.83489 0.0030619 
    122 4478.7  223.94  0.211   0.01055 0.81987 0.0020988 
    123 4620.8  231.04  0.211   0.01055 0.8253  0.0023899 
    124 4768.2  238.41  0.211   0.01055 0.82653 0.0022851 
    125 4921.1  246.06  0.211   0.01055 0.82442 0.003383 
    126 5079.8  253.99  0.211   0.01055 0.82827 0.0015979 
    127 5244.7  262.23  0.211   0.01055 0.82494 0.0031129 
    128 5415.7  270.79  0.211   0.01055 0.82183 0.0030149 
    129 5593.3  279.67  0.211   0.01055 0.83217 0.0046976 
    130 5777.8  288.89  0.211   0.01055 0.82227 0.005574 
    131 5969.3  298.46  0.211   0.01055 0.82338 0.0025569 
     1FileFormatVersion       1.0 
     2DataFileTitle           Polystyrene of Markus Strobl,  Full Sine, ++ only 
     3Sample                  Polystyrene 2 um in 53% H2O, 47% D2O 
     4Settings                D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast 
     5Operator                CPD 
     6Date                    ma 7 jul 2014 18:54:43 
     7ScanType                sine one element scan 
     8Thickness               2 
     9Thickness_unit          mm 
     10Theta_zmax              0.0168 
     11Theta_zmax_unit         radians 
     12Theta_ymax              0.0168 
     13Theta_ymax_unit         radians 
     14Orientation             Z 
     15SpinEchoLength_unit     A 
     16Depolarisation_unit     A-2 cm-1 
     17Wavelength_unit         A 
     18 
     19BEGIN_DATA 
     20SpinEchoLength Depolarisation Depolarisation_error Wavelength Wavelength_error SpinEchoLength_error Polarisation  Polarisation_error 
     21497.78  -2.4509553e-4  4.9935908e-4  2.11  0.1055  24.889  0.99782  0.0044367  
     22630.41   2.9161810e-4  5.3612769e-4  2.11  0.1055   31.52   1.0026  0.0047862  
     23764.87  -4.4899949e-4  6.8328154e-4  2.11  0.1055  38.244  0.99601  0.0060598  
     24898.47  -9.3037214e-4  6.5970686e-4  2.11  0.1055  44.924  0.99175  0.0058257  
     251034.1  -5.1441728e-4  6.8783151e-4  2.11  0.1055  51.705  0.99543  0.0060966  
     261169.5  -5.4939760e-4  5.4291131e-4  2.11  0.1055  58.475  0.99512  0.0048106  
     271306.1  -2.8111792e-4  7.0473347e-4  2.11  0.1055  65.303   0.9975  0.0062594  
     281443.7  -5.9342057e-4  5.9758787e-4  2.11  0.1055  72.184  0.99473   0.005293  
     291582.  -8.2284488e-4  6.0400267e-4  2.11  0.1055  79.102   0.9927  0.0053389  
     301721.2  -6.1600315e-4  7.2890343e-4  2.11  0.1055  86.062  0.99453  0.0064548  
     311861.7  -1.5367118e-3  8.3408174e-4  2.11  0.1055  93.087  0.98641  0.0073259  
     322002.8  -1.2660661e-3  8.9366844e-4  2.11  0.1055  100.14  0.98879  0.0078682  
     332144.6  -1.1854534e-3  7.7244662e-4  2.11  0.1055  107.23   0.9895  0.0068058  
     342287.3  -8.9753711e-4  9.3406529e-4  2.11  0.1055  114.36  0.99204  0.0082509  
     352431.2  -6.6909009e-4  5.8049041e-4  2.11  0.1055  121.56  0.99406  0.0051381  
     362575.5  -2.3090130e-3  6.3564144e-4  2.11  0.1055  128.78  0.97965  0.0055447  
     372722.2  -2.2895260e-3  7.5468967e-4  2.11  0.1055  136.11  0.97982  0.0065843  
     382869.  -2.5063662e-3  8.2471984e-4  2.11  0.1055  143.45  0.97793  0.0071814  
     393017.3  -2.4581433e-3  7.5988724e-4  2.11  0.1055  150.87  0.97835  0.0066197  
     403167.5  -2.1257395e-3  7.9526201e-4  2.11  0.1055  158.38  0.98125  0.0069484  
     413318.2  -2.4661790e-3  7.8910082e-4  2.11  0.1055  165.91  0.97828  0.0068737  
     423471.6  -2.6339122e-3  7.7536843e-4  2.11  0.1055  173.58  0.97682   0.006744  
     433624.5  -2.0914090e-3  6.8496070e-4  2.11  0.1055  181.22  0.98155  0.0059865  
     443780.9  -4.0640282e-3  7.7771292e-4  2.11  0.1055  189.04  0.96446  0.0066788  
     453937.4  -3.1016697e-3  8.5181234e-4  2.11  0.1055  196.87  0.97276  0.0073781  
     464096.1  -3.3038917e-3  8.2179560e-4  2.11  0.1055  204.81  0.97101  0.0071053  
     474255.5  -2.8422039e-3  7.8606869e-4  2.11  0.1055  212.78  0.97501  0.0068244  
     484419.1  -3.4694067e-3  8.3712733e-4  2.11  0.1055  220.96  0.96958  0.0072272  
     494581.2  -4.5951067e-3  7.4919003e-4  2.11  0.1055  229.06  0.95991  0.0064035  
     504747.7  -4.3286699e-3  6.9129591e-4  2.11  0.1055  237.39  0.96219  0.0059227  
     514915.2  -4.5658612e-3  6.8549385e-4  2.11  0.1055  245.76  0.96016  0.0058606  
     525085.1  -5.6542277e-3  5.8859074e-4  2.11  0.1055  254.26   0.9509  0.0049836  
     535256.8  -5.1028496e-3  6.4896117e-4  2.11  0.1055  262.84  0.95558  0.0055218  
     545430.8  -5.6672201e-3  6.8813882e-4  2.11  0.1055  271.54  0.95079  0.0058258  
     555607.4  -5.5951905e-3  5.1845870e-4  2.11  0.1055  280.37   0.9514  0.0043921  
     565786.9  -6.2795627e-3  5.5774416e-4  2.11  0.1055  289.35  0.94562  0.0046962  
     575967.5  -6.2486880e-3  5.6209080e-4  2.11  0.1055  298.38  0.94588  0.0047341  
     586151.7  -6.9992040e-3  5.4769136e-4  2.11  0.1055  307.58  0.93958  0.0045821  
     596337.7  -7.2708619e-3  5.1651117e-4  2.11  0.1055  316.88  0.93731  0.0043108  
     606528.2  -7.3511686e-3  6.2402654e-4  2.11  0.1055  326.41  0.93664  0.0052044  
     616720.4  -7.8980596e-3  4.2849488e-4  2.11  0.1055  336.02  0.93209  0.0035563  
     626915.7  -7.6861990e-3  5.5322868e-4  2.11  0.1055  345.78  0.93385  0.0046002  
     637114.8  -8.1163566e-3  6.1733111e-4  2.11  0.1055  355.74  0.93028  0.0051136  
     647316.7  -8.1356741e-3  4.2597330e-4  2.11  0.1055  365.83  0.93012  0.0035279  
     657521.7  -8.6415221e-3  4.7858305e-4  2.11  0.1055  376.08  0.92594  0.0039458  
     667731.  -8.9354263e-3  3.3579357e-4  2.11  0.1055  386.55  0.92352  0.0027613  
     677944.2  -9.5805772e-3  4.2302546e-4  2.11  0.1055  397.21  0.91823  0.0034587  
     688161.4   -0.010762148  4.2348254e-4  2.11  0.1055  408.07  0.90862  0.0034262  
     698381.9  -9.4754418e-3  4.2743183e-4  2.11  0.1055   419.1  0.91909   0.003498  
     708606.9   -0.010575665  3.7681487e-4  2.11  0.1055  430.34  0.91013  0.0030537  
     718837.4   -0.011273784  4.3943451e-4  2.11  0.1055  441.87  0.90449  0.0035391  
     729070.4   -0.011761571  3.7002787e-4  2.11  0.1055  453.52  0.90057  0.0029672  
     739309.5   -0.011439046  3.8852675e-4  2.11  0.1055  465.47  0.90316  0.0031245  
     749553.8   -0.012468380  3.8879110e-4  2.11  0.1055  477.69  0.89492  0.0030981  
     759802.5   -0.012287815  4.5888119e-4  2.11  0.1055  490.13  0.89636  0.0036625  
     7610057.   -0.013330052  2.8645708e-4  2.11  0.1055  502.84  0.88808  0.0022652  
     7710317.   -0.013511036  5.6519968e-4  2.11  0.1055  515.86  0.88665  0.0044622  
     7810583.   -0.014095206  3.8547484e-4  2.11  0.1055  529.14  0.88205  0.0030275  
     7910855.   -0.014029017  3.6440427e-4  2.11  0.1055  542.76  0.88257  0.0028637  
     8011134.   -0.014831000  3.8848283e-4  2.11  0.1055  556.69  0.87629  0.0030312  
     8111419.   -0.014699072  3.3972822e-4  2.11  0.1055  570.93  0.87732  0.0026539  
     8211710   -0.015174999  4.4468309e-4  2.11  0.1055   585.5  0.87361  0.0034591  
     8312010   -0.016164070  3.6480986e-4  2.11  0.1055  600.52  0.86595  0.0028129  
     8412316   -0.015893340  5.6035541e-4  2.11  0.1055  615.82  0.86804  0.0043311  
     8512630   -0.016059068  3.7324087e-4  2.11  0.1055  631.51  0.86676  0.0028806  
     8612952.   -0.017389837  3.7573625e-4  2.11  0.1055  647.62  0.85655  0.0028657  
     8713284.   -0.017207735  4.5518751e-4  2.11  0.1055  664.19  0.85794  0.0034773  
     8813623.   -0.017565669  4.3070477e-4  2.11  0.1055  681.13  0.85521  0.0032798  
     8913971.   -0.018286298  5.8862675e-4  2.11  0.1055  698.54  0.84974  0.0044537  
     9014329.   -0.019381994  4.3692639e-4  2.11  0.1055  716.44  0.84149  0.0032738  
     9114697.   -0.018515178  5.0158587e-4  2.11  0.1055  734.87  0.84801  0.0037874  
     9215075.   -0.019135361  4.2693908e-4  2.11  0.1055  753.76  0.84334   0.003206  
     9315464.   -0.019917113  5.4587670e-4  2.11  0.1055  773.18  0.83749  0.0040707  
     9415864.   -0.019686699  6.6829096e-4  2.11  0.1055  793.22  0.83921  0.0049938  
     9516276.   -0.020340321  7.7178617e-4  2.11  0.1055  813.81  0.83434  0.0057337  
     9616699.   -0.021091231  8.4580203e-4  2.11  0.1055  834.97  0.82878  0.0062417  
     9717136.   -0.021864932  1.0940027e-3  2.11  0.1055  856.78  0.82309  0.0080179  
     9817585.   -0.021695868  1.0107767e-3  2.11  0.1055  879.27  0.82433  0.0074191  
     9918048.   -0.022255844  1.0464994e-3  2.11  0.1055  902.39  0.82023  0.0076431  
     10018525.   -0.021256673  8.5582923e-4  2.11  0.1055  926.27  0.82756  0.0063064  
     10119017.   -0.021490334  7.1998911e-4  2.11  0.1055  950.87  0.82584  0.0052944  
     10219524.   -0.021198334  6.8716706e-4  2.11  0.1055  976.22  0.82799  0.0050662  
     10320048.   -0.021815823  5.8818928e-4  2.11  0.1055  1002.4  0.82345  0.0043127  
     10420588.   -0.021882671  5.1857307e-4  2.11  0.1055  1029.4  0.82296     0.0038  
     10521145.   -0.022305147  4.6672141e-4  2.11  0.1055  1057.2  0.81987  0.0034072  
     10621720   -0.022258583  5.0322361e-4  2.11  0.1055   1086.  0.82021  0.0036752  
     10722314.   -0.021244460  3.9148871e-4  2.11  0.1055  1115.7  0.82765  0.0028851  
     10822928.   -0.021381594  5.2906244e-4  2.11  0.1055  1146.4  0.82664  0.0038942  
     10923562.   -0.021329979  6.4328235e-4  2.11  0.1055  1178.1  0.82702  0.0047371  
     11024218.   -0.022846153  6.0248825e-4  2.11  0.1055  1210.9  0.81593  0.0043772  
     11124896.   -0.021591012  3.8146933e-4  2.11  0.1055  1244.8   0.8251  0.0028026  
     11225595.   -0.021034332  3.3282938e-4  2.11  0.1055  1279.8   0.8292  0.0024574  
     11326319.   -0.021433232  4.9200888e-4  2.11  0.1055  1315.9  0.82626  0.0036198  
     11427068.   -0.022177827  4.4213864e-4  2.11  0.1055  1353.4   0.8208  0.0032314  
     11527843.   -0.022343508  5.8553317e-4  2.11  0.1055  1392.2  0.81959  0.0042731  
     11628645.   -0.021396539  3.6673246e-4  2.11  0.1055  1432.3  0.82653   0.002699  
     11729475.   -0.021739473  5.0054859e-4  2.11  0.1055  1473.8  0.82401  0.0036726  
     11830335.   -0.021794003  6.5757715e-4  2.11  0.1055  1516.7  0.82361  0.0048224  
     11931224.   -0.021798094  5.6210549e-4  2.11  0.1055  1561.2  0.82358  0.0041221  
     12032145.   -0.022031519  3.9364070e-4  2.11  0.1055  1607.3  0.82187  0.0028807  
     12133101.   -0.021408769  4.7779613e-4  2.11  0.1055   1655.  0.82644   0.003516  
     12234090   -0.021802185  2.8863827e-4  2.11  0.1055  1704.5  0.82355  0.0021166  
     12335114.   -0.021586929  4.6155484e-4  2.11  0.1055  1755.7  0.82513  0.0033911  
     12436176.   -0.021194265  2.0808762e-4  2.11  0.1055  1808.8  0.82802  0.0015342  
     12537276.   -0.021382952  3.9701221e-4  2.11  0.1055  1863.8  0.82663  0.0029222  
     12638417.   -0.022251737  2.8416874e-4  2.11  0.1055  1920.8  0.82026  0.0020755  
     12739600   -0.020819189  3.5679523e-4  2.11  0.1055  1980.0  0.83079  0.0026394  
     12840827.   -0.021380235  3.7314604e-4  2.11  0.1055  2041.3  0.82665  0.0027466  
     12942099.   -0.021232248  3.4189634e-4  2.11  0.1055   2105.  0.82774  0.0025199  
     13043419.   -0.020266312  4.1187633e-4  2.11  0.1055  2170.9  0.83489  0.0030619  
     13144787.   -0.022305147  2.8749557e-4  2.11  0.1055  2239.4  0.81987  0.0020988  
     13246208.   -0.021563793  3.2521680e-4  2.11  0.1055  2310.4   0.8253  0.0023899  
     13347682.   -0.021396539  3.1049291e-4  2.11  0.1055  2384.1  0.82653  0.0022851  
     13449211.   -0.021683607  4.6084892e-4  2.11  0.1055  2460.6  0.82442   0.003383  
     13550798.   -0.021160361  2.1666201e-4  2.11  0.1055  2539.9  0.82827  0.0015979  
     13652447.   -0.021612792  4.2378726e-4  2.11  0.1055  2622.3  0.82494  0.0031129  
     13754157.   -0.022036985  4.1199886e-4  2.11  0.1055  2707.9  0.82183  0.0030149  
     13855933.   -0.020632795  6.3397053e-4  2.11  0.1055  2796.7  0.83217  0.0046976  
     13957778.   -0.021976873  7.6130313e-4  2.11  0.1055  2888.9  0.82227   0.005574  
     14059693.   -0.021825370  3.4875346e-4  2.11  0.1055  2984.6  0.82338  0.0025569  
  • example/se008731_01_40pcorr.ses

    r2f8ee5f rb314f89  
    1 DataFileTitle   "Polystyrene of Markus Strobl,  Full Sine, ++ only "                                             
    2 Sample  "Polystyrene 2 um in 53% H2O, 47% D2O "                                          
    3 Settings        "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast "                                           
    4 Operator        CPD                                              
    5 Date    do 10 jul 2014 16:37:30                                                  
    6 ScanType        sine one element scan                                            
    7 Thickness [cm]  2.00E-01                                                 
    8 Q_zmax [\A^-1]  0.05                                             
    9 Q_ymax [\A^-1]  0.05                                             
    10                                                          
    11 spin echo length [A]     Depolarisation [A-2 cm-1]       error depol [A-2 cm-1]          error SEL [A]   wavelength [A]          error wavelength [A]    polarisation    error pol  
    12 391.56  0.0041929       0.0036894       19.578  2.11    0.1055  1.0037  0.0032974 
    13 1564    -0.0046571      0.0038185       78.2    2.11    0.1055  0.99586 0.003386 
    14 2735.6  -0.017007       0.0038132       136.78  2.11    0.1055  0.98497 0.0033444 
    15 3907.9  -0.033462       0.0035068       195.39  2.11    0.1055  0.97064 0.0030309 
    16 5080.2  -0.047483       0.0038208       254.01  2.11    0.1055  0.9586  0.0032613 
    17 6251.8  -0.070375       0.00376 312.59  2.11    0.1055  0.93926 0.0031446 
    18 7423.2  -0.092217       0.0037927       371.16  2.11    0.1055  0.92117 0.0031108 
    19 8595.5  -0.10238        0.004006        429.77  2.11    0.1055  0.91287 0.0032562 
    20 9767.7  -0.12672        0.0038534       488.39  2.11    0.1055  0.8933  0.0030651 
    21 10940   -0.1374 0.004243        546.98  2.11    0.1055  0.88484 0.003343 
    22 12112   -0.16072        0.0045837       605.58  2.11    0.1055  0.86666 0.0035372 
    23 13284   -0.16623        0.0045613       664.2   2.11    0.1055  0.86242 0.0035027 
    24 14456   -0.18468        0.0044918       722.79  2.11    0.1055  0.84837 0.0033931 
    25 15628   -0.19143        0.0048967       781.38  2.11    0.1055  0.84328 0.0036768 
    26 16800   -0.20029        0.0045421       840.02  2.11    0.1055  0.83666 0.0033837 
    27 17971   -0.19798        0.0046642       898.56  2.11    0.1055  0.83838 0.0034819 
    28 19143   -0.21442        0.0047052       957.17  2.11    0.1055  0.82619 0.0034614 
    29 20316   -0.20885        0.0044931       1015.8  2.11    0.1055  0.8303  0.0033218 
    30 21488   -0.21393        0.0049186       1074.4  2.11    0.1055  0.82655 0.00362 
    31 22660   -0.20685        0.004423        1133    2.11    0.1055  0.83179 0.0032758 
    32 23832   -0.20802        0.0046979       1191.6  2.11    0.1055  0.83092 0.0034758 
    33 25003   -0.19848        0.0045953       1250.2  2.11    0.1055  0.838   0.0034289 
    34 26175   -0.21117        0.0044567       1308.8  2.11    0.1055  0.82859 0.0032881 
    35 27347   -0.21283        0.004137        1367.4  2.11    0.1055  0.82736 0.0030477 
    36 28520   -0.2042 0.0044587       1426    2.11    0.1055  0.83375 0.0033101 
    37 29692   -0.2112 0.0042852       1484.6  2.11    0.1055  0.82857 0.0031615 
    38 30864   -0.20319        0.0043483       1543.2  2.11    0.1055  0.8345  0.003231 
    39 32036   -0.20752        0.0044297       1601.8  2.11    0.1055  0.83129 0.0032788 
    40 33207   -0.20654        0.0043188       1660.4  2.11    0.1055  0.83201 0.0031995 
    41 34380   -0.20126        0.0046375       1719    2.11    0.1055  0.83593 0.0034518 
    42 35551   -0.20924        0.0042871       1777.6  2.11    0.1055  0.83001 0.0031684 
    43 36724   -0.21323        0.0045471       1836.2  2.11    0.1055  0.82707 0.0033487 
    44 37895   -0.21324        0.0045354       1894.7  2.11    0.1055  0.82706 0.00334 
    45 39067   -0.19905        0.0044141       1953.4  2.11    0.1055  0.83758 0.003292 
    46 40239   -0.1991 0.0047441       2012    2.11    0.1055  0.83754 0.003538 
    47 41411   -0.20359        0.0050136       2070.5  2.11    0.1055  0.8342  0.003724 
    48 42583   -0.21032        0.0049474       2129.1  2.11    0.1055  0.82922 0.0036529 
    49 43755   -0.20689        0.0048203       2187.8  2.11    0.1055  0.83176 0.00357 
    50 44927   -0.21075        0.0052337       2246.4  2.11    0.1055  0.8289  0.0038628 
    51 46099   -0.19956        0.0047827       2304.9  2.11    0.1055  0.8372  0.0035653 
     1FileFormatVersion       1.0 
     2DataFileTitle           Polystyrene of Markus Strobl,  Full Sine, ++ only 
     3Sample                  Polystyrene 2 um in 53% H2O, 47% D2O 
     4Settings                "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast "                                           
     5Operator                CPD 
     6Date                    do 10 jul 2014 16:37:30 
     7ScanType                sine one element scan 
     8Thickness               2 
     9Thickness_unit          mm 
     10Theta_zmax              0.0168 
     11Theta_zmax_unit         radians 
     12Theta_ymax              0.0168 
     13Theta_ymax_unit         radians 
     14Orientation             Z 
     15SpinEchoLength_unit     A 
     16Depolarisation_unit     A-2 cm-1 
     17Wavelength_unit         A 
     18 
     19BEGIN_DATA 
     20SpinEchoLength Depolarisation Depolarisation_error Polarisation Polarisation_error SpinEchoLength_error Wavelength Wavelength_error  
     21391.56 0.0041929 0.0036894 1.0037 0.0032974 19.578 2.11 0.1055  
     221564 -0.0046571 0.0038185 0.99586 0.003386 78.2 2.11 0.1055  
     232735.6 -0.017007 0.0038132 0.98497 0.0033444 136.78 2.11 0.1055  
     243907.9 -0.033462 0.0035068 0.97064 0.0030309 195.39 2.11 0.1055  
     255080.2 -0.047483 0.0038208 0.9586 0.0032613 254.01 2.11 0.1055  
     266251.8 -0.070375 0.00376 0.93926 0.0031446 312.59 2.11 0.1055  
     277423.2 -0.092217 0.0037927 0.92117 0.0031108 371.16 2.11 0.1055  
     288595.5 -0.10238 0.004006 0.91287 0.0032562 429.77 2.11 0.1055  
     299767.7 -0.12672 0.0038534 0.8933 0.0030651 488.39 2.11 0.1055  
     3010940 -0.1374 0.004243 0.88484 0.003343 546.98 2.11 0.1055  
     3112112 -0.16072 0.0045837 0.86666 0.0035372 605.58 2.11 0.1055  
     3213284 -0.16623 0.0045613 0.86242 0.0035027 664.2 2.11 0.1055  
     3314456 -0.18468 0.0044918 0.84837 0.0033931 722.79 2.11 0.1055  
     3415628 -0.19143 0.0048967 0.84328 0.0036768 781.38 2.11 0.1055  
     3516800 -0.20029 0.0045421 0.83666 0.0033837 840.02 2.11 0.1055  
     3617971 -0.19798 0.0046642 0.83838 0.0034819 898.56 2.11 0.1055  
     3719143 -0.21442 0.0047052 0.82619 0.0034614 957.17 2.11 0.1055  
     3820316 -0.20885 0.0044931 0.8303 0.0033218 1015.8 2.11 0.1055  
     3921488 -0.21393 0.0049186 0.82655 0.00362 1074.4 2.11 0.1055  
     4022660 -0.20685 0.004423 0.83179 0.0032758 1133 2.11 0.1055  
     4123832 -0.20802 0.0046979 0.83092 0.0034758 1191.6 2.11 0.1055  
     4225003 -0.19848 0.0045953 0.838 0.0034289 1250.2 2.11 0.1055  
     4326175 -0.21117 0.0044567 0.82859 0.0032881 1308.8 2.11 0.1055  
     4427347 -0.21283 0.004137 0.82736 0.0030477 1367.4 2.11 0.1055  
     4528520 -0.2042 0.0044587 0.83375 0.0033101 1426 2.11 0.1055  
     4629692 -0.2112 0.0042852 0.82857 0.0031615 1484.6 2.11 0.1055  
     4730864 -0.20319 0.0043483 0.8345 0.003231 1543.2 2.11 0.1055  
     4832036 -0.20752 0.0044297 0.83129 0.0032788 1601.8 2.11 0.1055  
     4933207 -0.20654 0.0043188 0.83201 0.0031995 1660.4 2.11 0.1055  
     5034380 -0.20126 0.0046375 0.83593 0.0034518 1719 2.11 0.1055  
     5135551 -0.20924 0.0042871 0.83001 0.0031684 1777.6 2.11 0.1055  
     5236724 -0.21323 0.0045471 0.82707 0.0033487 1836.2 2.11 0.1055  
     5337895 -0.21324 0.0045354 0.82706 0.00334 1894.7 2.11 0.1055  
     5439067 -0.19905 0.0044141 0.83758 0.003292 1953.4 2.11 0.1055  
     5540239 -0.1991 0.0047441 0.83754 0.003538 2012 2.11 0.1055  
     5641411 -0.20359 0.0050136 0.8342 0.003724 2070.5 2.11 0.1055  
     5742583 -0.21032 0.0049474 0.82922 0.0036529 2129.1 2.11 0.1055  
     5843755 -0.20689 0.0048203 0.83176 0.00357 2187.8 2.11 0.1055  
     5944927 -0.21075 0.0052337 0.8289 0.0038628 2246.4 2.11 0.1055  
     6046099 -0.19956 0.0047827 0.8372 0.0035653 2304.9 2.11 0.1055  
  • example/sphere.ses

    ra98958b re0721e7  
    1 name    se009051.dat                                             
    2 sample  "SiO2 pH=6 & NiO pH=10, nanosized, in D2O, Sine ++ only."                                                
    3 collimation     "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. "                                                
    4 ID      CPD                                              
    5 Date    di 27 jan 2015 18:17:50                                          
    6 Method  sine one element scan                                            
    7 Thickness [mm]  10                                               
    8 z-acceptance [nm^-1]    0.3                                              
    9 y-acceptance [nm^-1]    0.3                                              
     1FileFormatVersion       1.0 
     2DataFileTitle           se009051.dat                                             
     3Sample                  "SiO2 pH=6 & NiO pH=10, nanosized, in D2O, Sine ++ only."                                                
     4collimation             "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. "                                                
     5ID                      CPD                                              
     6Date                    di 27 jan 2015 18:17:50 
     7Method                  sine one element scan 
     8Thickness               10 
     9Thickness_unit          mm 
     10Theta_zmax              0.0168 
     11Theta_zmax_unit         radians 
     12Theta_ymax              0.0168 
     13Theta_ymax_unit         radians 
     14Orientation             Z 
     15SpinEchoLength_unit     A 
     16Depolarisation_unit     A-2 cm-1 
     17Wavelength_unit         A 
    1018                                                         
    11 spin echo length [um]   error spin-echo length  wavelength [nm] error wavelength        Sample 5        Sample 5 -error          
    12 29.882023       5       0.2095  0.01    0.997817556     0.00887342      0.029882023     29.882023 
    13 85.257531       5       0.2095  0.01    1.0026184       0.009572307     0.085257531     85.257531 
    14 140.798634      5       0.2095  0.01    0.996013172     0.012119508     0.140798634     140.798634 
    15 196.560502      5       0.2095  0.01    0.991752656     0.011651473     0.196560502     196.560502 
    16 252.528125      5       0.2095  0.01    0.995429571     0.012193294     0.252528125     252.528125 
    17 308.729274      5       0.2095  0.01    0.995119819     0.009621277     0.308729274     308.729274 
    18 365.037117      5       0.2095  0.01    0.997497767     0.012518842     0.365037117     365.037117 
    19 421.684677      5       0.2095  0.01    0.994733865     0.010585926     0.421684677     421.684677 
    20 478.564852      5       0.2095  0.01    0.992696275     0.010677724     0.478564852     478.564852 
    21 535.715521      5       0.2095  0.01    0.994534294     0.012909556     0.535715521     535.715521 
    22 593.24408       5       0.2095  0.01    0.986412268     0.014651894     0.59324408      593.24408 
    23 650.985498      5       0.2095  0.01    0.988788566     0.015736468     0.650985498     650.985498 
    24 709.203958      5       0.2095  0.01    0.989501168     0.013611615     0.709203958     709.203958 
    25 767.809943      5       0.2095  0.01    0.99203933      0.016501891     0.767809943     767.809943 
    26 826.762925      5       0.2095  0.01    0.994062723     0.010276207     0.826762925     826.762925 
    27 886.302296      5       0.2095  0.01    0.979647192     0.011089421     0.886302296     886.302296 
    28 946.236366      5       0.2095  0.01    0.97982075      0.013168527     0.946236366     946.236366 
    29 1006.784115     5       0.2095  0.01    0.977927712     0.014362819     1.006784115     1006.784115 
    30 1067.891213     5       0.2095  0.01    0.978351462     0.013239476     1.067891213     1067.891213 
    31 1129.598722     5       0.2095  0.01    0.981246471     0.013896898     1.129598722     1129.598722 
    32 1192.022694     5       0.2095  0.01    0.978275164     0.013747481     1.192022694     1192.022694 
    33 1255.271282     5       0.2095  0.01    0.976818842     0.013487944     1.255271282     1255.271282 
    34 1319.307134     5       0.2095  0.01    0.981546068     0.011972943     1.319307134     1319.307134 
    35 1384.197821     5       0.2095  0.01    0.96445854      0.013357572     1.384197821     1384.197821 
    36 1450.051342     5       0.2095  0.01    0.972756158     0.014756189     1.450051342     1450.051342 
    37 1516.978619     5       0.2095  0.01    0.971010184     0.014210504     1.516978619     1516.978619 
    38 1585.044577     5       0.2095  0.01    0.975009829     0.013648839     1.585044577     1585.044577 
    39 1654.356194     5       0.2095  0.01    0.969575209     0.014454336     1.654356194     1654.356194 
    40 1724.97535      5       0.2095  0.01    0.959907267     0.012806944     1.72497535      1724.97535 
    41 1796.962718     5       0.2095  0.01    0.96219203      0.011845496     1.796962718     1796.962718 
    42 1870.613025     5       0.2095  0.01    0.960157966     0.011721155     1.870613025     1870.613025 
    43 1945.921163     5       0.2095  0.01    0.95089937      0.009967287     1.945921163     1945.921163 
    44 2023.09841      5       0.2095  0.01    0.955584642     0.01104368      2.02309841      2023.09841 
    45 2102.189874     5       0.2095  0.01    0.950786938     0.011651582     2.102189874     2102.189874 
    46 2183.534527     5       0.2095  0.01    0.951403607     0.008784149     2.183534527     2183.534527 
    47 2267.173021     5       0.2095  0.01    0.94561823      0.009392383     2.267173021     2267.173021 
    48 2353.393035     5       0.2095  0.01    0.945882787     0.009468264     2.353393035     2353.393035 
    49 2442.350387     5       0.2095  0.01    0.939583158     0.009164239     2.442350387     2442.350387 
    50 2534.334343     5       0.2095  0.01    0.93731257      0.008621568     2.534334343     2534.334343 
    51 2629.473374     5       0.2095  0.01    0.936639557     0.010408821     2.629473374     2629.473374 
    52 2728.20036      5       0.2095  0.01    0.932091884     0.007112674     2.72820036      2728.20036 
    53 2830.74992      5       0.2095  0.01    0.933852763     0.009200431     2.83074992      2830.74992 
    54 2937.363026     5       0.2095  0.01    0.930283437     0.01022719      2.937363026     2937.363026 
    55 3048.510973     5       0.2095  0.01    0.930121334     0.007055747     3.048510973     3048.510973 
    56 3164.589926     5       0.2095  0.01    0.925944326     0.0078917       3.164589926     3164.589926 
    57 3285.843756     5       0.2095  0.01    0.92351759      0.005522582     3.285843756     3285.843756 
    58 3412.805801     5       0.2095  0.01    0.918228311     0.006917328     3.412805801     3412.805801 
    59 3546.016792     5       0.2095  0.01    0.908621563     0.006852313     3.546016792     3546.016792 
    60 3685.876213     5       0.2095  0.01    0.919087792     0.006995961     3.685876213     3685.876213 
    61 3833.174832     5       0.2095  0.01    0.910132214     0.006107336     3.833174832     3833.174832 
    62 3988.234456     5       0.2095  0.01    0.904491462     0.007078101     3.988234456     3988.234456 
    63 4151.940408     5       0.2095  0.01    0.90056919      0.005934428     4.151940408     4151.940408 
    64 4324.902332     5       0.2095  0.01    0.903155701     0.006249069     4.324902332     4324.902332 
    65 4507.963711     5       0.2095  0.01    0.894919827     0.006196232     4.507963711     4507.963711 
    66 4701.967049     5       0.2095  0.01    0.896358908     0.007324999     4.701967049     4701.967049 
    67 4907.785127     5       0.2095  0.01    0.88807923      0.004530465     4.907785127     4907.785127 
    68 5126.471149     5       0.2095  0.01    0.886649039     0.008924412     5.126471149     5126.471149 
    69 5359.059759     5       0.2095  0.01    0.882051947     0.006055081     5.359059759     5359.059759 
    70 5606.740304     5       0.2095  0.01    0.882573371     0.005727419     5.606740304     5606.740304 
    71 5870.903592     5       0.2095  0.01    0.876289825     0.006062306     5.870903592     5870.903592 
    72 6152.736235     5       0.2095  0.01    0.877318199     0.005307783     6.152736235     6152.736235 
    73 6453.893169     5       0.2095  0.01    0.873610284     0.006918104     6.453893169     6453.893169 
    74 6775.835415     5       0.2095  0.01    0.865954455     0.005625863     6.775835415     6775.835415 
    75 7120.530581     5       0.2095  0.01    0.868037308     0.008662141     7.120530581     7120.530581 
    76 7489.735939     5       0.2095  0.01    0.86676404      0.005761209     7.489735939     7489.735939 
    77 7885.556041     5       0.2095  0.01    0.85655072      0.005731411     7.885556041     7885.556041 
    78 8310.290162     5       0.2095  0.01    0.857936274     0.006954572     8.310290162     8310.290162 
    79 8766.315359     5       0.2095  0.01    0.85520525      0.006559622     8.766315359     8766.315359 
    80 9256.271706     5       0.2095  0.01    0.849739829     0.008907377     9.256271706     9256.271706 
    81 9783.031577     5       0.2095  0.01    0.841487161     0.006547659     9.783031577     9783.031577 
    82 10349.71409     5       0.2095  0.01    0.84801051      0.007574866     10.34971409     10349.71409 
    83 10959.67923     5       0.2095  0.01    0.843342861     0.006411993     10.95967923     10959.67923 
    84 11616.5969      5       0.2095  0.01    0.837488454     0.008141353     11.6165969      11616.5969 
    85 12324.44664     5       0.2095  0.01    0.839214922     0.009987535     12.32444664     12324.44664 
    86 13087.4723      5       0.2095  0.01    0.834336972     0.011467406     13.0874723      13087.4723 
    87 13910.29604     5       0.2095  0.01    0.828775634     0.012483365     13.91029604     13910.29604 
    88 14798.05061     5       0.2095  0.01    0.823088812     0.016035848     14.79805061     14798.05061 
    89 15756.09817     5       0.2095  0.01    0.824329178     0.014838138     15.75609817     15756.09817 
    90 16790.46783     5       0.2095  0.01    0.82023172      0.015286187     16.79046783     16790.46783 
    91 17907.62866     5       0.2095  0.01    0.827559752     0.012612741     17.90762866     17907.62866 
     19BEGIN_DATA 
     20SpinEchoLength  SpinEchoLength_error  Wavelength  Wavelength_Error  Polarisation  Polarisation_error  Depolarisation  Depolarisation_error  
     21298.82023 50 2.095 0.1 0.997817556 0.00887342 -4.9779370e-5 2.0261512e-4  
     22852.57531 50 2.095 0.1 1.0026184 0.009572307 5.9579929e-5 2.1752686e-4  
     231407.9863 50 2.095 0.1 0.996013172 0.012119508 -9.1017859e-5 2.7723742e-4  
     241965.6050 50 2.095 0.1 0.991752656 0.011651473 -1.8868750e-4 2.6767598e-4  
     252525.2813 50 2.095 0.1 0.995429571 0.012193294 -1.0437182e-4 2.7908883e-4  
     263087.2927 50 2.095 0.1 0.995119819 0.009621277 -1.1146275e-4 2.2028721e-4  
     273650.3712 50 2.095 0.1 0.997497767 0.012518842 -5.7082583e-5 2.8594610e-4  
     284216.8468 50 2.095 0.1 0.994733865 0.010585926 -1.2030120e-4 2.4246770e-4  
     294785.6485 50 2.095 0.1 0.992696275 0.010677724 -1.6701950e-4 2.4507231e-4  
     305357.1552 50 2.095 0.1 0.994534294 0.012909556 -1.2487278e-4 2.9574914e-4  
     315932.4408 50 2.095 0.1 0.986412268 0.014651894 -3.1170682e-4 3.3842875e-4  
     326509.8550 50 2.095 0.1 0.988788566 0.015736468 -2.5688520e-4 3.6260666e-4  
     337092.0396 50 2.095 0.1 0.989501168 0.013611615 -2.4047103e-4 3.1341898e-4  
     347678.0994 50 2.095 0.1 0.99203933 0.016501891 -1.8210252e-4 3.7899787e-4  
     358267.6293 50 2.095 0.1 0.994062723 0.010276207 -1.3567871e-4 2.3553259e-4  
     368863.0230 50 2.095 0.1 0.979647192 0.011089421 -4.6850452e-4 2.5791174e-4  
     379462.3637 50 2.095 0.1 0.97982075 0.013168527 -4.6446835e-4 3.0621222e-4  
     3810067.841 50 2.095 0.1 0.977927712 0.014362819 -5.0853039e-4 3.3463000e-4  
     3910678.912 50 2.095 0.1 0.978351462 0.013239476 -4.9865985e-4 3.0832436e-4  
     4011295.987 50 2.095 0.1 0.981246471 0.013896898 -4.3133968e-4 3.2267975e-4  
     4111920.227 50 2.095 0.1 0.978275164 0.013747481 -5.0043677e-4 3.2017989e-4  
     4212552.713 50 2.095 0.1 0.976818842 0.013487944 -5.3437989e-4 3.1460359e-4  
     4313193.071 50 2.095 0.1 0.981546068 0.011972943 -4.2438423e-4 2.7792152e-4  
     4413841.978 50 2.095 0.1 0.96445854 0.013357572 -8.2452102e-4 3.1555561e-4  
     4514500.513 50 2.095 0.1 0.972756158 0.014756189 -6.2933879e-4 3.4562263e-4  
     4615169.786 50 2.095 0.1 0.971010184 0.014210504 -6.7027011e-4 3.3343996e-4  
     4715850.446 50 2.095 0.1 0.975009829 0.013648839 -5.7661387e-4 3.1894710e-4  
     4816543.562 50 2.095 0.1 0.969575209 0.014454336 -7.0396574e-4 3.3966327e-4  
     4917249.754 50 2.095 0.1 0.959907267 0.012806944 -9.3229353e-4 3.0398222e-4  
     5017969.627 50 2.095 0.1 0.96219203 0.011845496 -8.7812744e-4 2.8049391e-4  
     5118706.130 50 2.095 0.1 0.960157966 0.011721155 -9.2634378e-4 2.7813758e-4  
     5219459.212 50 2.095 0.1 0.95089937 0.009967287 -1.1471121e-3 2.3882201e-4  
     5320230.984 50 2.095 0.1 0.955584642 0.01104368 -1.0351259e-3 2.6331561e-4  
     5421021.899 50 2.095 0.1 0.950786938 0.011651582 -1.1498062e-3 2.7921171e-4  
     5521835.345 50 2.095 0.1 0.951403607 0.008784149 -1.1350335e-3 2.1036178e-4  
     5622671.730 50 2.095 0.1 0.94561823 0.009392383 -1.2740040e-3 2.2630383e-4  
     5723533.930 50 2.095 0.1 0.945882787 0.009468264 -1.2676305e-3 2.2806833e-4  
     5824423.504 50 2.095 0.1 0.939583158 0.009164239 -1.4198814e-3 2.2222511e-4  
     5925343.343 50 2.095 0.1 0.93731257 0.008621568 -1.4750079e-3 2.0957224e-4  
     6026294.734 50 2.095 0.1 0.936639557 0.010408821 -1.4913733e-3 2.5319842e-4  
     6127282.004 50 2.095 0.1 0.932091884 0.007112674 -1.6022666e-3 1.7386258e-4  
     6228307.499 50 2.095 0.1 0.933852763 0.009200431 -1.5592642e-3 2.2447176e-4  
     6329373.630 50 2.095 0.1 0.930283437 0.01022719 -1.6465153e-3 2.5047996e-4  
     6430485.110 50 2.095 0.1 0.930121334 0.007055747 -1.6504858e-3 1.7283645e-4  
     6531645.899 50 2.095 0.1 0.925944326 0.0078917 -1.7530356e-3 1.9418588e-4  
     6632858.438 50 2.095 0.1 0.92351759 0.005522582 -1.8128270e-3 1.3624763e-4  
     6734128.058 50 2.095 0.1 0.918228311 0.006917328 -1.9436940e-3 1.7164045e-4  
     6835460.168 50 2.095 0.1 0.908621563 0.006852313 -2.1833230e-3 1.7182490e-4  
     6936858.762 50 2.095 0.1 0.919087792 0.006995961 -1.9223776e-3 1.7342924e-4  
     7038331.748 50 2.095 0.1 0.910132214 0.006107336 -2.1454742e-3 1.5289007e-4  
     7139882.345 50 2.095 0.1 0.904491462 0.007078101 -2.2871233e-3 1.7829708e-4  
     7241519.404 50 2.095 0.1 0.90056919 0.005934428 -2.3861400e-3 1.5013907e-4  
     7343249.023 50 2.095 0.1 0.903155701 0.006249069 -2.3207959e-3 1.5764661e-4  
     7445079.637 50 2.095 0.1 0.894919827 0.006196232 -2.5295172e-3 1.5775222e-4  
     7547019.670 50 2.095 0.1 0.896358908 0.007324999 -2.4929085e-3 1.8619052e-4  
     7649077.851 50 2.095 0.1 0.88807923 0.004530465 -2.7043436e-3 1.1623128e-4  
     7751264.711 50 2.095 0.1 0.886649039 0.008924412 -2.7410654e-3 2.2932944e-4  
     7853590.598 50 2.095 0.1 0.882051947 0.006055081 -2.8595036e-3 1.5640756e-4  
     7956067.403 50 2.095 0.1 0.882573371 0.005727419 -2.8460388e-3 1.4785638e-4  
     8058709.036 50 2.095 0.1 0.876289825 0.006062306 -3.0088321e-3 1.5762389e-4  
     8161527.362 50 2.095 0.1 0.877318199 0.005307783 -2.9821094e-3 1.3784403e-4  
     8264538.932 50 2.095 0.1 0.873610284 0.006918104 -3.0786086e-3 1.8042690e-4  
     8367758.354 50 2.095 0.1 0.865954455 0.005625863 -3.2791557e-3 1.4802192e-4  
     8471205.306 50 2.095 0.1 0.868037308 0.008662141 -3.2244196e-3 2.2736248e-4  
     8574897.359 50 2.095 0.1 0.86676404 0.005761209 -3.2578647e-3 1.5144143e-4  
     8678855.560 50 2.095 0.1 0.85655072 0.005731411 -3.5279304e-3 1.5245456e-4  
     8783102.902 50 2.095 0.1 0.857936274 0.006954572 -3.4911046e-3 1.8469168e-4  
     8887663.154 50 2.095 0.1 0.85520525 0.006559622 -3.5637478e-3 1.7475934e-4  
     8992562.717 50 2.095 0.1 0.849739829 0.008907377 -3.7098230e-3 2.3883381e-4  
     9097830.316 50 2.095 0.1 0.841487161 0.006547659 -3.9321836e-3 1.7728439e-4  
     91103497.14 50 2.095 0.1 0.84801051 0.007574866 -3.7562386e-3 2.0351933e-4  
     92109596.79 50 2.095 0.1 0.843342861 0.006411993 -3.8819940e-3 1.7322909e-4  
     93116165.97 50 2.095 0.1 0.837488454 0.008141353 -4.0407107e-3 2.2148775e-4  
     94123244.47 50 2.095 0.1 0.839214922 0.009987535 -3.9937900e-3 2.7115465e-4  
     95130874.72 50 2.095 0.1 0.834336972 0.011467406 -4.1266093e-3 3.1315233e-4  
     96139102.96 50 2.095 0.1 0.828775634 0.012483365 -4.2789870e-3 3.4318369e-4  
     97147980.51 50 2.095 0.1 0.823088812 0.016035848 -4.4358638e-3 4.4389186e-4  
     98157560.98 50 2.095 0.1 0.824329178 0.014838138 -4.4015548e-3 4.1011975e-4  
     99167904.68 50 2.095 0.1 0.82023172 0.015286187 -4.5150892e-3 4.2461424e-4  
     100179076.29 50 2.095 0.1 0.827559752 0.012612741 -4.3124376e-3 3.4724985e-4  
  • explore/asymint.py

    ra1c32c2 rc462169  
    3030import numpy as np 
    3131import mpmath as mp 
    32 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 
     32from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos 
    3333from numpy.polynomial.legendre import leggauss 
    3434from scipy.integrate import dblquad, simps, romb, romberg 
     
    4141shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
    4242Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
     43 
     44DTYPE = 'd' 
    4345 
    4446class MPenv: 
     
    7274    sas_sinx_x = staticmethod(sp.sas_sinx_x) 
    7375    pi = np.pi 
    74     mpf = staticmethod(float) 
     76    #mpf = staticmethod(float) 
     77    mpf = staticmethod(lambda x: np.array(x, DTYPE)) 
    7578 
    7679SLD = 3 
     
    220223    #A, B, C = 4450, 14000, 47 
    221224    #A, B, C = 445, 140, 47  # integer for the sake of mpf 
    222     A, B, C = 6800, 114, 1380 
    223     DA, DB, DC = 2300, 21, 58 
     225    A, B, C = 114, 1380, 6800 
     226    DA, DB, DC = 21, 58, 2300 
    224227    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 
    225232    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
    226233    #SLD_SOLVENT,CONTRAST = 0, 4 
     
    233240        DA, DC = DC, DA 
    234241        SLDA, SLDC = SLDC, SLDA 
     242    #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    235243    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    236244    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
     
    348356    return n**2, Iq 
    349357 
     358def gauss_quad_usub(q, n=150, dtype=DTYPE): 
     359    """ 
     360    Compute the integral using gaussian quadrature for n = 20, 76 or 150. 
     361 
     362    Use *u = sin theta* substitution, and restrict integration over a single 
     363    quadrant for shapes that are mirror symmetric about AB, AC and BC planes. 
     364 
     365    Note that this doesn't work for fcc/bcc paracrystals, which instead step 
     366    over the entire 4 pi surface uniformly in theta-phi. 
     367    """ 
     368    z, w = leggauss(n) 
     369    cos_theta = 0.5 * (z + 1) 
     370    theta = arccos(cos_theta) 
     371    phi = pi/2*(0.5 * (z + 1)) 
     372    Atheta, Aphi = np.meshgrid(theta, phi) 
     373    Aw = w[None, :] * w[:, None] 
     374    q, Atheta, Aphi, Aw = [np.asarray(v, dtype=dtype) for v in (q, Atheta, Aphi, Aw)] 
     375    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi) 
     376    Iq = np.sum(Zq*Aw)*0.25 
     377    return n**2, Iq 
     378 
    350379def gridded_2d(q, n=300): 
    351380    """ 
     
    395424    print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
    396425    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)) 
    397429    #gridded_2d(Q, n=2**8+1) 
    398430    gridded_2d(Q, n=2**10+1) 
  • explore/realspace.py

    r8fb2a94 r362a66f  
    33import time 
    44from copy import copy 
     5import os 
     6import argparse 
     7from collections import OrderedDict 
    58 
    69import numpy as np 
    710from numpy import pi, radians, sin, cos, sqrt 
    8 from numpy.random import poisson, uniform 
     11from numpy.random import poisson, uniform, randn, rand 
    912from numpy.polynomial.legendre import leggauss 
    1013from scipy.integrate import simps 
    1114from scipy.special import j1 as J1 
     15 
     16try: 
     17    import numba 
     18    USE_NUMBA = True 
     19except ImportError: 
     20    USE_NUMBA = False 
    1221 
    1322# Definition of rotation matrices comes from wikipedia: 
     
    8291        raise NotImplementedError() 
    8392 
     93    def dims(self): 
     94        # type: () -> float, float, float 
     95        raise NotImplementedError() 
     96 
    8497    def rotate(self, theta, phi, psi): 
    8598        self.rotation = rotation(theta, phi, psi) * self.rotation 
     
    116129                     for s2 in shapes] 
    117130        self.r_max = max(distances + [s.r_max for s in shapes]) 
    118  
    119     def volume(self): 
    120         return sum(shape.volume() for shape in self.shapes) 
     131        self.volume = sum(shape.volume for shape in self.shapes) 
    121132 
    122133    def sample(self, density): 
     
    133144        self._scale = np.array([a/2, b/2, c/2])[None, :] 
    134145        self.r_max = sqrt(a**2 + b**2 + c**2) 
    135  
    136     def volume(self): 
    137         return self.a*self.b*self.c 
     146        self.dims = a, b, c 
     147        self.volume = a*b*c 
    138148 
    139149    def sample(self, density): 
    140         num_points = poisson(density*self.a*self.b*self.c) 
     150        num_points = poisson(density*self.volume) 
    141151        points = self._scale*uniform(-1, 1, size=(num_points, 3)) 
    142152        values = self.value.repeat(points.shape[0]) 
     
    152162        self._scale = np.array([ra, rb, length/2])[None, :] 
    153163        self.r_max = sqrt(4*max(ra, rb)**2 + length**2) 
    154  
    155     def volume(self): 
    156         return pi*self.ra*self.rb*self.length 
     164        self.dims = 2*ra, 2*rb, length 
     165        self.volume = pi*ra*rb*length 
    157166 
    158167    def sample(self, density): 
    159         # density of the bounding box 
     168        # randomly sample from a box of side length 2*r, excluding anything 
     169        # not in the cylinder 
    160170        num_points = poisson(density*4*self.ra*self.rb*self.length) 
    161171        points = uniform(-1, 1, size=(num_points, 3)) 
    162172        radius = points[:, 0]**2 + points[:, 1]**2 
    163         points = self._scale*points[radius <= 1] 
     173        points = points[radius <= 1] 
    164174        values = self.value.repeat(points.shape[0]) 
    165         return values, self._adjust(points) 
     175        return values, self._adjust(self._scale*points) 
     176 
     177class EllipticalBicelle(Shape): 
     178    def __init__(self, ra, rb, length, 
     179                 thick_rim, thick_face, 
     180                 value_core, value_rim, value_face, 
     181                 center=(0, 0, 0), orientation=(0, 0, 0)): 
     182        self.rotate(*orientation) 
     183        self.shift(*center) 
     184        self.value = value_core 
     185        self.ra, self.rb, self.length = ra, rb, length 
     186        self.thick_rim, self.thick_face = thick_rim, thick_face 
     187        self.value_rim, self.value_face = value_rim, value_face 
     188 
     189        # reset cylinder to outer dimensions for calculating scale, etc. 
     190        ra = self.ra + self.thick_rim 
     191        rb = self.rb + self.thick_rim 
     192        length = self.length + 2*self.thick_face 
     193        self._scale = np.array([ra, rb, length/2])[None, :] 
     194        self.r_max = sqrt(4*max(ra, rb)**2 + length**2) 
     195        self.dims = 2*ra, 2*rb, length 
     196        self.volume = pi*ra*rb*length 
     197 
     198    def sample(self, density): 
     199        # randomly sample from a box of side length 2*r, excluding anything 
     200        # not in the cylinder 
     201        ra = self.ra + self.thick_rim 
     202        rb = self.rb + self.thick_rim 
     203        length = self.length + 2*self.thick_face 
     204        num_points = poisson(density*4*ra*rb*length) 
     205        points = uniform(-1, 1, size=(num_points, 3)) 
     206        radius = points[:, 0]**2 + points[:, 1]**2 
     207        points = points[radius <= 1] 
     208        # set all to core value first 
     209        values = np.ones_like(points[:, 0])*self.value 
     210        # then set value to face value if |z| > face/(length/2)) 
     211        values[abs(points[:, 2]) > self.length/(self.length + 2*self.thick_face)] = self.value_face 
     212        # finally set value to rim value if outside the core ellipse 
     213        radius = (points[:, 0]**2*(1 + self.thick_rim/self.ra)**2 
     214                  + points[:, 1]**2*(1 + self.thick_rim/self.rb)**2) 
     215        values[radius>1] = self.value_rim 
     216        return values, self._adjust(self._scale*points) 
    166217 
    167218class TriaxialEllipsoid(Shape): 
     
    174225        self._scale = np.array([ra, rb, rc])[None, :] 
    175226        self.r_max = 2*max(ra, rb, rc) 
    176  
    177     def volume(self): 
    178         return 4*pi/3 * self.ra * self.rb * self.rc 
     227        self.dims = 2*ra, 2*rb, 2*rc 
     228        self.volume = 4*pi/3 * ra * rb * rc 
    179229 
    180230    def sample(self, density): 
     
    194244        self.rotate(*orientation) 
    195245        self.shift(*center) 
     246        helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 
     247        total_radius = self.helix_radius + self.tube_radius 
    196248        self.helix_radius, self.helix_pitch = helix_radius, helix_pitch 
    197249        self.tube_radius, self.tube_length = tube_radius, tube_length 
    198         helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 
    199         self.r_max = sqrt((2*helix_radius + 2*tube_radius)*2 
    200                           + (helix_length + 2*tube_radius)**2) 
    201  
    202     def volume(self): 
     250        self.r_max = sqrt(4*total_radius + (helix_length + 2*tube_radius)**2) 
     251        self.dims = 2*total_radius, 2*total_radius, helix_length 
    203252        # small tube radius approximation; for larger tubes need to account 
    204253        # for the fact that the inner length is much shorter than the outer 
    205254        # length 
    206         return pi*self.tube_radius**2*self.tube_length 
     255        self.volume = pi*self.tube_radius**2*self.tube_length 
    207256 
    208257    def points(self, density): 
     
    244293        return values, self._adjust(points) 
    245294 
    246 NUMBA = False 
    247 if NUMBA: 
     295def csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
     296    core = Box(a, b, c, sld_core) 
     297    side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0)) 
     298    side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0)) 
     299    side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2)) 
     300    side_a2 = copy(side_a).shift(-a-da, 0, 0) 
     301    side_b2 = copy(side_b).shift(0, -b-db, 0) 
     302    side_c2 = copy(side_c).shift(0, 0, -c-dc) 
     303    shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
     304    shape.dims = 2*da+a, 2*db+b, 2*dc+c 
     305    return shape 
     306 
     307def _Iqxy(values, x, y, z, qa, qb, qc): 
     308    """I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r)""" 
     309    Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 
     310            for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 
     311    return Iq 
     312 
     313if USE_NUMBA: 
     314    # Override simple numpy solution with numba if available 
    248315    from numba import njit 
    249316    @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])") 
     
    252319        for j in range(len(Iq)): 
    253320            total = 0. + 0j 
    254             for k in range(len(Iq)): 
     321            for k in range(len(values)): 
    255322                total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k])) 
    256323            Iq[j] = abs(total)**2 
     
    263330    values = rho*volume 
    264331    x, y, z = points.T 
     332    values, x, y, z, qa, qb, qc = [np.asarray(v, 'd') 
     333                                   for v in (values, x, y, z, qa, qb, qc)] 
    265334 
    266335    # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r) 
    267     if NUMBA: 
    268         Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 
    269     else: 
    270         Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 
    271               for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 
     336    Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 
    272337    return np.asarray(Iq).reshape(qx.shape) / np.sum(volume) 
    273338 
     
    338403""" 
    339404 
    340 if NUMBA: 
     405if USE_NUMBA: 
     406    # Override simple numpy solution with numba if available 
    341407    @njit("f8[:](f8[:], f8[:], f8[:,:])") 
    342     def _calc_Pr_uniform_jit(r, rho, points): 
     408    def _calc_Pr_uniform(r, rho, points): 
    343409        dr = r[0] 
    344410        n_max = len(r) 
     
    362428    # P(r) with uniform steps in r is 3x faster; check if we are uniform 
    363429    # before continuing 
     430    r, rho, points = [np.asarray(v, 'd') for v in (r, rho, points)] 
    364431    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 
    365432        Pr = _calc_Pr_nonuniform(r, rho, points) 
    366433    else: 
    367         if NUMBA: 
    368             Pr = _calc_Pr_uniform_jit(r, rho, points) 
    369         else: 
    370             Pr = _calc_Pr_uniform(r, rho, points) 
     434        Pr = _calc_Pr_uniform(r, rho, points) 
    371435    return Pr / Pr.max() 
    372436 
     
    399463    return edges 
    400464 
    401 def plot_calc(r, Pr, q, Iq, theory=None): 
     465# -------------- plotters ---------------- 
     466def plot_calc(r, Pr, q, Iq, theory=None, title=None): 
    402467    import matplotlib.pyplot as plt 
    403468    plt.subplot(211) 
     
    405470    plt.xlabel('r (A)') 
    406471    plt.ylabel('Pr (1/A^2)') 
     472    if title is not None: 
     473        plt.title(title) 
    407474    plt.subplot(212) 
    408475    plt.loglog(q, Iq, '-', label='from Pr') 
     
    410477    plt.ylabel('Iq') 
    411478    if theory is not None: 
    412         plt.loglog(theory[0], theory[1], '-', label='analytic') 
     479        plt.loglog(theory[0], theory[1]/theory[1][0], '-', label='analytic') 
    413480        plt.legend() 
    414481 
    415 def plot_calc_2d(qx, qy, Iqxy, theory=None): 
     482def plot_calc_2d(qx, qy, Iqxy, theory=None, title=None): 
    416483    import matplotlib.pyplot as plt 
    417484    qx, qy = bin_edges(qx), bin_edges(qy) 
     
    419486    if theory is not None: 
    420487        plt.subplot(121) 
    421     plt.pcolormesh(qx, qy, np.log10(Iqxy)) 
     488    #plt.pcolor(qx, qy, np.log10(Iqxy)) 
     489    extent = [qx[0], qx[-1], qy[0], qy[-1]] 
     490    plt.imshow(np.log10(Iqxy), extent=extent, interpolation="nearest", 
     491               origin='lower') 
    422492    plt.xlabel('qx (1/A)') 
    423493    plt.ylabel('qy (1/A)') 
     494    plt.axis('equal') 
     495    plt.axis(extent) 
     496    #plt.grid(True) 
     497    if title is not None: 
     498        plt.title(title) 
    424499    if theory is not None: 
    425500        plt.subplot(122) 
    426         plt.pcolormesh(qx, qy, np.log10(theory)) 
     501        plt.imshow(np.log10(theory), extent=extent, interpolation="nearest", 
     502                   origin='lower') 
     503        plt.axis('equal') 
     504        plt.axis(extent) 
    427505        plt.xlabel('qx (1/A)') 
    428506 
     
    440518    index = np.random.choice(n, size=500) if n > 500 else slice(None, None) 
    441519    ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index]) 
     520    # make square axes 
     521    minmax = np.array([points.min(), points.max()]) 
     522    ax.scatter(minmax, minmax, minmax, c='w') 
    442523    #low, high = points.min(axis=0), points.max(axis=0) 
    443524    #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]]) 
     525    ax.set_xlabel("x") 
     526    ax.set_ylabel("y") 
     527    ax.set_zlabel("z") 
    444528    ax.autoscale(True) 
    445529 
    446 def check_shape(shape, fn=None): 
    447     rho_solvent = 0 
    448     q = np.logspace(-3, 0, 200) 
    449     r = shape.r_bins(q, r_step=0.01) 
    450     sampling_density = 6*5000 / shape.volume() 
    451     rho, points = shape.sample(sampling_density) 
    452     t0 = time.time() 
    453     Pr = calc_Pr(r, rho-rho_solvent, points) 
    454     print("calc Pr time", time.time() - t0) 
    455     Iq = calc_Iq(q, r, Pr) 
    456     theory = (q, fn(q)) if fn is not None else None 
    457  
    458     import pylab 
    459     #plot_points(rho, points); pylab.figure() 
    460     plot_calc(r, Pr, q, Iq, theory=theory) 
    461     pylab.show() 
    462  
    463 def check_shape_2d(shape, fn=None, view=(0, 0, 0)): 
    464     rho_solvent = 0 
    465     nq, qmax = 100, 1.0 
    466     qx = np.linspace(0.0, qmax, nq) 
    467     qy = np.linspace(0.0, qmax, nq) 
    468     Qx, Qy = np.meshgrid(qx, qy) 
    469     sampling_density = 50000 / shape.volume() 
    470     #t0 = time.time() 
    471     rho, points = shape.sample(sampling_density) 
    472     #print("sample time", time.time() - t0) 
    473     t0 = time.time() 
    474     Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
    475     print("calc time", time.time() - t0) 
    476     theory = fn(Qx, Qy) if fn is not None else None 
    477     Iqxy += 0.001 * Iqxy.max() 
    478     if theory is not None: 
    479         theory += 0.001 * theory.max() 
    480  
    481     import pylab 
    482     #plot_points(rho, points); pylab.figure() 
    483     plot_calc_2d(qx, qy, Iqxy, theory=theory) 
    484     pylab.show() 
    485  
     530# ----------- Analytic models -------------- 
    486531def sas_sinx_x(x): 
    487532    with np.errstate(all='ignore'): 
     
    510555    for k, qk in enumerate(q): 
    511556        qab, qc = qk*sin_alpha, qk*cos_alpha 
    512         Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2) 
     557        Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 
    513558        Iq[k] = np.sum(w*Fq**2) 
    514     Iq = Iq/Iq[0] 
     559    Iq = Iq 
    515560    return Iq 
    516561 
    517562def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)): 
    518563    qa, qb, qc = invert_view(qx, qy, view) 
    519     qab = np.sqrt(qa**2 + qb**2) 
    520     Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2) 
     564    qab = sqrt(qa**2 + qb**2) 
     565    Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 
    521566    Iq = Fq**2 
    522567    return Iq.reshape(qx.shape) 
     
    524569def sphere_Iq(q, radius): 
    525570    Iq = sas_3j1x_x(q*radius)**2 
    526     return Iq/Iq[0] 
     571    return Iq 
     572 
     573def box_Iq(q, a, b, c): 
     574    z, w = leggauss(76) 
     575    outer_sum = np.zeros_like(q) 
     576    for cos_alpha, outer_w in zip((z+1)/2, w): 
     577        sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 
     578        qc = q*cos_alpha 
     579        siC = c*sas_sinx_x(c*qc/2) 
     580        inner_sum = np.zeros_like(q) 
     581        for beta, inner_w in zip((z + 1)*pi/4, w): 
     582            qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta) 
     583            siA = a*sas_sinx_x(a*qa/2) 
     584            siB = b*sas_sinx_x(b*qb/2) 
     585            Fq = siA*siB*siC 
     586            inner_sum += inner_w * Fq**2 
     587        outer_sum += outer_w * inner_sum 
     588    Iq = outer_sum / 4  # = outer*um*zm*8.0/(4.0*M_PI) 
     589    return Iq 
     590 
     591def box_Iqxy(qx, qy, a, b, c, view=(0, 0, 0)): 
     592    qa, qb, qc = invert_view(qx, qy, view) 
     593    sia = sas_sinx_x(qa*a/2) 
     594    sib = sas_sinx_x(qb*b/2) 
     595    sic = sas_sinx_x(qc*c/2) 
     596    Fq = sia*sib*sic 
     597    Iq = Fq**2 
     598    return Iq.reshape(qx.shape) 
    527599 
    528600def csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core): 
     
    539611        sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 
    540612        qc = q*cos_alpha 
    541         siC = c*j0(c*qc/2) 
    542         siCt = tC*j0(tC*qc/2) 
     613        siC = c*sas_sinx_x(c*qc/2) 
     614        siCt = tC*sas_sinx_x(tC*qc/2) 
    543615        inner_sum = np.zeros_like(q) 
    544616        for beta, inner_w in zip((z + 1)*pi/4, w): 
    545617            qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta) 
    546             siA = a*j0(a*qa/2) 
    547             siB = b*j0(b*qb/2) 
    548             siAt = tA*j0(tA*qa/2) 
    549             siBt = tB*j0(tB*qb/2) 
     618            siA = a*sas_sinx_x(a*qa/2) 
     619            siB = b*sas_sinx_x(b*qb/2) 
     620            siAt = tA*sas_sinx_x(tA*qa/2) 
     621            siBt = tB*sas_sinx_x(tB*qb/2) 
    550622            if overlapping: 
    551623                Fq = (dr0*siA*siB*siC 
     
    584656    return Iq.reshape(qx.shape) 
    585657 
    586 def check_cylinder(radius=25, length=125, rho=2.): 
     658def sasmodels_Iq(kernel, q, pars): 
     659    from sasmodels.data import empty_data1D 
     660    from sasmodels.direct_model import DirectModel 
     661    data = empty_data1D(q) 
     662    calculator = DirectModel(data, kernel) 
     663    Iq = calculator(**pars) 
     664    return Iq 
     665 
     666def sasmodels_Iqxy(kernel, qx, qy, pars, view): 
     667    from sasmodels.data import Data2D 
     668    from sasmodels.direct_model import DirectModel 
     669    Iq = 100 * np.ones_like(qx) 
     670    data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq)) 
     671    data.x_bins = qx[0,:] 
     672    data.y_bins = qy[:,0] 
     673    data.filename = "fake data" 
     674 
     675    calculator = DirectModel(data, kernel) 
     676    pars_plus_view = pars.copy() 
     677    pars_plus_view.update(theta=view[0], phi=view[1], psi=view[2]) 
     678    Iqxy = calculator(**pars_plus_view) 
     679    return Iqxy.reshape(qx.shape) 
     680 
     681def wrap_sasmodel(name, **pars): 
     682    from sasmodels.core import load_model 
     683    kernel = load_model(name) 
     684    fn = lambda q: sasmodels_Iq(kernel, q, pars) 
     685    fn_xy = lambda qx, qy, view: sasmodels_Iqxy(kernel, qx, qy, pars, view) 
     686    return fn, fn_xy 
     687 
     688 
     689# --------- Test cases ----------- 
     690 
     691def build_cylinder(radius=25, length=125, rho=2.): 
    587692    shape = EllipticalCylinder(radius, radius, length, rho) 
    588     fn = lambda q: cylinder_Iq(q, radius, length) 
    589     check_shape(shape, fn) 
    590  
    591 def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)): 
    592     shape = EllipticalCylinder(radius, radius, length, rho) 
    593     fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    594     check_shape_2d(shape, fn, view=view) 
    595  
    596 def check_cylinder_2d_lattice(radius=25, length=125, rho=2., 
    597                               view=(0, 0, 0)): 
    598     nx, dx = 1, 2*radius 
    599     ny, dy = 30, 2*radius 
    600     nz, dz = 30, length 
    601     dx, dy, dz = 2*dx, 2*dy, 2*dz 
    602     def center(*args): 
    603         sigma = 0.333 
    604         space = 2 
    605         return [(space*n+np.random.randn()*sigma)*x for n, x in args] 
    606     shapes = [EllipticalCylinder(radius, radius, length, rho, 
    607                                  #center=(ix*dx, iy*dy, iz*dz) 
    608                                  orientation=np.random.randn(3)*0, 
    609                                  center=center((ix, dx), (iy, dy), (iz, dz)) 
    610                                 ) 
     693    fn = lambda q: cylinder_Iq(q, radius, length)*rho**2 
     694    fn_xy = lambda qx, qy, view: cylinder_Iqxy(qx, qy, radius, length, view=view)*rho**2 
     695    return shape, fn, fn_xy 
     696 
     697def build_sphere(radius=125, rho=2): 
     698    shape = TriaxialEllipsoid(radius, radius, radius, rho) 
     699    fn = lambda q: sphere_Iq(q, radius)*rho**2 
     700    fn_xy = lambda qx, qy, view: sphere_Iq(np.sqrt(qx**2+qy**2), radius)*rho**2 
     701    return shape, fn, fn_xy 
     702 
     703def build_box(a=10, b=20, c=30, rho=2.): 
     704    shape = Box(a, b, c, rho) 
     705    fn = lambda q: box_Iq(q, a, b, c)*rho**2 
     706    fn_xy = lambda qx, qy, view: box_Iqxy(qx, qy, a, b, c, view=view)*rho**2 
     707    return shape, fn, fn_xy 
     708 
     709def build_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
     710    shape = csbox(a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
     711    fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
     712    fn_xy = lambda qx, qy, view: csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 
     713                                            slda, sldb, sldc, sld_core, view=view) 
     714    return shape, fn, fn_xy 
     715 
     716def build_ellcyl(ra=25, rb=50, length=125, rho=2.): 
     717    shape = EllipticalCylinder(ra, rb, length, rho) 
     718    fn, fn_xy = wrap_sasmodel( 
     719        'elliptical_cylinder', 
     720        scale=1, 
     721        background=0, 
     722        radius_minor=ra, 
     723        axis_ratio=rb/ra, 
     724        length=length, 
     725        sld=rho, 
     726        sld_solvent=0, 
     727    ) 
     728    return shape, fn, fn_xy 
     729 
     730def build_cscyl(ra=30, rb=90, length=30, thick_rim=8, thick_face=14, 
     731                sld_core=4, sld_rim=1, sld_face=7): 
     732    shape = EllipticalBicelle( 
     733        ra=ra, rb=rb, length=length, 
     734        thick_rim=thick_rim, thick_face=thick_face, 
     735        value_core=sld_core, value_rim=sld_rim, value_face=sld_face, 
     736        ) 
     737    fn, fn_xy = wrap_sasmodel( 
     738        'core_shell_bicelle_elliptical', 
     739        scale=1, 
     740        background=0, 
     741        radius=ra, 
     742        x_core=rb/ra, 
     743        length=length, 
     744        thick_rim=thick_rim, 
     745        thick_face=thick_face, 
     746        sld_core=sld_core, 
     747        sld_face=sld_face, 
     748        sld_rim=sld_rim, 
     749        sld_solvent=0, 
     750    ) 
     751    return shape, fn, fn_xy 
     752 
     753def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     754                  shuffle=0, rotate=0): 
     755    a, b, c = shape.dims 
     756    shapes = [copy(shape) 
     757              .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     758                     (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     759                     (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     760              .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
    611761              for ix in range(nx) 
    612762              for iy in range(ny) 
    613763              for iz in range(nz)] 
    614     shape = Composite(shapes) 
    615     fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    616     check_shape_2d(shape, fn, view=view) 
    617  
    618 def check_sphere(radius=125, rho=2): 
    619     shape = TriaxialEllipsoid(radius, radius, radius, rho) 
    620     fn = lambda q: sphere_Iq(q, radius) 
    621     check_shape(shape, fn) 
    622  
    623 def check_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
    624     core = Box(a, b, c, sld_core) 
    625     side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0)) 
    626     side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0)) 
    627     side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2)) 
    628     side_a2 = copy(side_a).shift(-a-da, 0, 0) 
    629     side_b2 = copy(side_b).shift(0, -b-db, 0) 
    630     side_c2 = copy(side_c).shift(0, 0, -c-dc) 
    631     shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
    632     def fn(q): 
    633         return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
    634     #check_shape(shape, fn) 
    635  
    636     view = (20, 30, 40) 
    637     def fn_xy(qx, qy): 
    638         return csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 
    639                           slda, sldb, sldc, sld_core, view=view) 
    640     check_shape_2d(shape, fn_xy, view=view) 
     764    lattice = Composite(shapes) 
     765    return lattice 
     766 
     767 
     768SHAPE_FUNCTIONS = OrderedDict([ 
     769    ("cyl", build_cylinder), 
     770    ("ellcyl", build_ellcyl), 
     771    ("sphere", build_sphere), 
     772    ("box", build_box), 
     773    ("csbox", build_csbox), 
     774    ("cscyl", build_cscyl), 
     775]) 
     776SHAPES = list(SHAPE_FUNCTIONS.keys()) 
     777 
     778def check_shape(title, shape, fn=None, show_points=False, 
     779                mesh=100, qmax=1.0, r_step=0.01, samples=5000): 
     780    rho_solvent = 0 
     781    qmin = qmax/100. 
     782    q = np.logspace(np.log10(qmin), np.log10(qmax), mesh) 
     783    r = shape.r_bins(q, r_step=r_step) 
     784    sampling_density = samples / shape.volume 
     785    rho, points = shape.sample(sampling_density) 
     786    t0 = time.time() 
     787    Pr = calc_Pr(r, rho-rho_solvent, points) 
     788    print("calc Pr time", time.time() - t0) 
     789    Iq = calc_Iq(q, r, Pr) 
     790    theory = (q, fn(q)) if fn is not None else None 
     791 
     792    import pylab 
     793    if show_points: 
     794         plot_points(rho, points); pylab.figure() 
     795    plot_calc(r, Pr, q, Iq, theory=theory, title=title) 
     796    pylab.gcf().canvas.set_window_title(title) 
     797    pylab.show() 
     798 
     799def check_shape_2d(title, shape, fn=None, view=(0, 0, 0), show_points=False, 
     800                   mesh=100, qmax=1.0, samples=5000): 
     801    rho_solvent = 0 
     802    #qx = np.linspace(0.0, qmax, mesh) 
     803    #qy = np.linspace(0.0, qmax, mesh) 
     804    qx = np.linspace(-qmax, qmax, mesh) 
     805    qy = np.linspace(-qmax, qmax, mesh) 
     806    Qx, Qy = np.meshgrid(qx, qy) 
     807    sampling_density = samples / shape.volume 
     808    t0 = time.time() 
     809    rho, points = shape.sample(sampling_density) 
     810    print("point generation time", time.time() - t0) 
     811    t0 = time.time() 
     812    Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
     813    print("calc Iqxy time", time.time() - t0) 
     814    t0 = time.time() 
     815    theory = fn(Qx, Qy, view) if fn is not None else None 
     816    print("calc theory time", time.time() - t0) 
     817    Iqxy += 0.001 * Iqxy.max() 
     818    if theory is not None: 
     819        theory += 0.001 * theory.max() 
     820 
     821    import pylab 
     822    if show_points: 
     823        plot_points(rho, points); pylab.figure() 
     824    plot_calc_2d(qx, qy, Iqxy, theory=theory, title=title) 
     825    pylab.gcf().canvas.set_window_title(title) 
     826    pylab.show() 
     827 
     828def main(): 
     829    parser = argparse.ArgumentParser( 
     830        description="Compute scattering from realspace sampling", 
     831        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     832        ) 
     833    parser.add_argument('-d', '--dim', type=int, default=1, 
     834                        help='dimension 1 or 2') 
     835    parser.add_argument('-m', '--mesh', type=int, default=100, 
     836                        help='number of mesh points') 
     837    parser.add_argument('-s', '--samples', type=int, default=5000, 
     838                        help="number of sample points") 
     839    parser.add_argument('-q', '--qmax', type=float, default=0.5, 
     840                        help='max q') 
     841    parser.add_argument('-v', '--view', type=str, default='0,0,0', 
     842                        help='theta,phi,psi angles') 
     843    parser.add_argument('-n', '--lattice', type=str, default='1,1,1', 
     844                        help='lattice size') 
     845    parser.add_argument('-z', '--spacing', type=str, default='2,2,2', 
     846                        help='lattice spacing') 
     847    parser.add_argument('-r', '--rotate', type=float, default=0., 
     848                        help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") 
     849    parser.add_argument('-w', '--shuffle', type=float, default=0., 
     850                        help="position relative to lattice, gaussian < 0.3, uniform otherwise") 
     851    parser.add_argument('-p', '--plot', action='store_true', 
     852                        help='plot points') 
     853    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], 
     854                        help='oriented shape') 
     855    parser.add_argument('pars', type=str, nargs='*', help='shape parameters') 
     856    opts = parser.parse_args() 
     857    pars = {key: float(value) for p in opts.pars for key, value in [p.split('=')]} 
     858    nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 
     859    dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 
     860    shuffle, rotate = opts.shuffle, opts.rotate 
     861    shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 
     862    if nx > 1 or ny > 1 or nz > 1: 
     863        shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 
     864    title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 
     865    if opts.dim == 1: 
     866        check_shape(title, shape, fn, show_points=opts.plot, 
     867                    mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
     868    else: 
     869        view = tuple(float(v) for v in opts.view.split(',')) 
     870        check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 
     871                       mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
     872 
    641873 
    642874if __name__ == "__main__": 
    643     check_cylinder(radius=10, length=40) 
    644     #check_cylinder_2d(radius=10, length=40, view=(90,30,0)) 
    645     #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0)) 
    646     #check_sphere() 
    647     #check_csbox() 
    648     #check_csbox(da=100, db=200, dc=300) 
     875    main() 
  • sasmodels/bumps_model.py

    r2d81cfe r49d1f8b8  
    137137    """ 
    138138    _cache = None # type: Dict[str, np.ndarray] 
    139     def __init__(self, data, model, cutoff=1e-5, name=None): 
     139    def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): 
    140140        # type: (Data, Model, float) -> None 
    141141        # remember inputs so we can inspect from outside 
     
    145145        self._interpret_data(data, model.sasmodel) 
    146146        self._cache = {} 
     147        self.extra_pars = extra_pars 
    147148 
    148149    def update(self): 
     
    166167        Return a dictionary of parameters 
    167168        """ 
    168         return self.model.parameters() 
     169        pars = self.model.parameters() 
     170        if self.extra_pars: 
     171            pars.update(self.extra_pars) 
     172        return pars 
    169173 
    170174    def theory(self): 
  • sasmodels/compare.py

    r2a7e20e r1fbadb2  
    4040from . import core 
    4141from . import kerneldll 
     42from . import kernelcl 
    4243from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4344from .direct_model import DirectModel, get_mesh 
     
    623624 
    624625 
     626def time_calculation(calculator, pars, evals=1): 
     627    # type: (Calculator, ParameterSet, int) -> Tuple[np.ndarray, float] 
     628    """ 
     629    Compute the average calculation time over N evaluations. 
     630 
     631    An additional call is generated without polydispersity in order to 
     632    initialize the calculation engine, and make the average more stable. 
     633    """ 
     634    # initialize the code so time is more accurate 
     635    if evals > 1: 
     636        calculator(**suppress_pd(pars)) 
     637    toc = tic() 
     638    # make sure there is at least one eval 
     639    value = calculator(**pars) 
     640    for _ in range(evals-1): 
     641        value = calculator(**pars) 
     642    average_time = toc()*1000. / evals 
     643    #print("I(q)",value) 
     644    return value, average_time 
     645 
     646def make_data(opts): 
     647    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
     648    """ 
     649    Generate an empty dataset, used with the model to set Q points 
     650    and resolution. 
     651 
     652    *opts* contains the options, with 'qmax', 'nq', 'res', 
     653    'accuracy', 'is2d' and 'view' parsed from the command line. 
     654    """ 
     655    qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res'] 
     656    if opts['is2d']: 
     657        q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray 
     658        data = empty_data2D(q, resolution=res) 
     659        data.accuracy = opts['accuracy'] 
     660        set_beam_stop(data, qmin) 
     661        index = ~data.mask 
     662    else: 
     663        if opts['view'] == 'log' and not opts['zero']: 
     664            q = np.logspace(math.log10(qmin), math.log10(qmax), nq) 
     665        else: 
     666            q = np.linspace(qmin, qmax, nq) 
     667        if opts['zero']: 
     668            q = np.hstack((0, q)) 
     669        data = empty_data1D(q, resolution=res) 
     670        index = slice(None, None) 
     671    return data, index 
     672 
    625673DTYPE_MAP = { 
    626674    'half': '16', 
     
    643691    Return a model calculator using the OpenCL calculation engine. 
    644692    """ 
    645     if not core.HAVE_OPENCL: 
    646         raise RuntimeError("OpenCL not available") 
    647     model = core.build_model(model_info, dtype=dtype, platform="ocl") 
    648     calculator = DirectModel(data, model, cutoff=cutoff) 
    649     calculator.engine = "OCL%s"%DTYPE_MAP[str(model.dtype)] 
    650     return calculator 
    651693 
    652694def eval_ctypes(model_info, data, dtype='double', cutoff=0.): 
     
    660702    return calculator 
    661703 
    662 def time_calculation(calculator, pars, evals=1): 
    663     # type: (Calculator, ParameterSet, int) -> Tuple[np.ndarray, float] 
    664     """ 
    665     Compute the average calculation time over N evaluations. 
    666  
    667     An additional call is generated without polydispersity in order to 
    668     initialize the calculation engine, and make the average more stable. 
    669     """ 
    670     # initialize the code so time is more accurate 
    671     if evals > 1: 
    672         calculator(**suppress_pd(pars)) 
    673     toc = tic() 
    674     # make sure there is at least one eval 
    675     value = calculator(**pars) 
    676     for _ in range(evals-1): 
    677         value = calculator(**pars) 
    678     average_time = toc()*1000. / evals 
    679     #print("I(q)",value) 
    680     return value, average_time 
    681  
    682 def make_data(opts): 
    683     # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
    684     """ 
    685     Generate an empty dataset, used with the model to set Q points 
    686     and resolution. 
    687  
    688     *opts* contains the options, with 'qmax', 'nq', 'res', 
    689     'accuracy', 'is2d' and 'view' parsed from the command line. 
    690     """ 
    691     qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res'] 
    692     if opts['is2d']: 
    693         q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray 
    694         data = empty_data2D(q, resolution=res) 
    695         data.accuracy = opts['accuracy'] 
    696         set_beam_stop(data, qmin) 
    697         index = ~data.mask 
    698     else: 
    699         if opts['view'] == 'log' and not opts['zero']: 
    700             q = np.logspace(math.log10(qmin), math.log10(qmax), nq) 
    701         else: 
    702             q = np.linspace(qmin, qmax, nq) 
    703         if opts['zero']: 
    704             q = np.hstack((0, q)) 
    705         data = empty_data1D(q, resolution=res) 
    706         index = slice(None, None) 
    707     return data, index 
    708  
    709704def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
    710705    # type: (ModelInfo, Data, str, float) -> Calculator 
     
    718713        set_integration_size(model_info, ngauss) 
    719714 
    720     if dtype is None or not dtype.endswith('!'): 
    721         return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    722     else: 
    723         return eval_ctypes(model_info, data, dtype=dtype[:-1], cutoff=cutoff) 
     715    if dtype != "default" and not dtype.endswith('!') and not kernelcl.use_opencl(): 
     716        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) 
     717 
     718    model = core.build_model(model_info, dtype=dtype, platform="ocl") 
     719    calculator = DirectModel(data, model, cutoff=cutoff) 
     720    engine_type = calculator._model.__class__.__name__.replace('Model', '').upper() 
     721    bits = calculator._model.dtype.itemsize*8 
     722    precision = "fast" if getattr(calculator._model, 'fast', False) else str(bits) 
     723    calculator.engine = "%s[%s]" % (engine_type, precision) 
     724    return calculator 
    724725 
    725726def _show_invalid(data, theory): 
     
    751752    """ 
    752753    for k in range(opts['sets']): 
    753         if k > 1: 
     754        if k > 0: 
    754755            # print a separate seed for each dataset for better reproducibility 
    755756            new_seed = np.random.randint(1000000) 
    756             print("Set %d uses -random=%i"%(k+1, new_seed)) 
     757            print("=== Set %d uses -random=%i ==="%(k+1, new_seed)) 
    757758            np.random.seed(new_seed) 
    758759        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
     
    761762        result = run_models(opts, verbose=True) 
    762763        if opts['plot']: 
     764            if opts['is2d'] and k > 0: 
     765                import matplotlib.pyplot as plt 
     766                plt.figure() 
    763767            limits = plot_models(opts, result, limits=limits, setnum=k) 
    764768        if opts['show_weights']: 
     
    13281332 
    13291333    # Evaluate preset parameter expressions 
     1334    # Note: need to replace ':' with '_' in parameter names and expressions 
     1335    # in order to support math on magnetic parameters. 
    13301336    context = MATH.copy() 
    13311337    context['np'] = np 
    1332     context.update(pars) 
     1338    context.update((k.replace(':', '_'), v) for k, v in pars.items()) 
    13331339    context.update((k, v) for k, v in presets.items() if isinstance(v, float)) 
     1340    #for k,v in sorted(context.items()): print(k, v) 
    13341341    for k, v in presets.items(): 
    13351342        if not isinstance(v, float) and not k.endswith('_type'): 
    1336             presets[k] = eval(v, context) 
     1343            presets[k] = eval(v.replace(':', '_'), context) 
    13371344    context.update(presets) 
    1338     context.update((k, v) for k, v in presets2.items() if isinstance(v, float)) 
     1345    context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float)) 
    13391346    for k, v in presets2.items(): 
    13401347        if not isinstance(v, float) and not k.endswith('_type'): 
    1341             presets2[k] = eval(v, context) 
     1348            presets2[k] = eval(v.replace(':', '_'), context) 
    13421349 
    13431350    # update parameters with presets 
     
    13691376    path = os.path.dirname(info.filename) 
    13701377    url = "file://" + path.replace("\\", "/")[2:] + "/" 
    1371     rst2html.view_html_qtapp(html, url) 
     1378    rst2html.view_html_wxapp(html, url) 
    13721379 
    13731380def explore(opts): 
     
    14901497            vmin, vmax = limits 
    14911498            self.limits = vmax*1e-7, 1.3*vmax 
    1492             import pylab; pylab.clf() 
     1499            import pylab 
     1500            pylab.clf() 
    14931501            plot_models(self.opts, result, limits=self.limits) 
    14941502 
  • sasmodels/convert.py

    r2d81cfe ra69d8cd  
    633633def test_backward_forward(): 
    634634    from .core import list_models 
     635    L = lambda name: _check_one(name, seed=1) 
    635636    for name in list_models('all'): 
    636         L = lambda: _check_one(name, seed=1) 
    637         L.description = name 
    638         yield L 
     637        yield L, name 
  • sasmodels/core.py

    r2d81cfe r3221de0  
    2121from . import mixture 
    2222from . import kernelpy 
     23from . import kernelcl 
    2324from . import kerneldll 
    2425from . import custom 
    25  
    26 if os.environ.get("SAS_OPENCL", "").lower() == "none": 
    27     HAVE_OPENCL = False 
    28 else: 
    29     try: 
    30         from . import kernelcl 
    31         HAVE_OPENCL = True 
    32     except Exception: 
    33         HAVE_OPENCL = False 
    34  
    35 CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "") 
    36 if CUSTOM_MODEL_PATH == "": 
    37     CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
    38     if not os.path.isdir(CUSTOM_MODEL_PATH): 
    39         os.makedirs(CUSTOM_MODEL_PATH) 
    4026 
    4127# pylint: disable=unused-import 
     
    4733    pass 
    4834# pylint: enable=unused-import 
     35 
     36CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "") 
     37if CUSTOM_MODEL_PATH == "": 
     38    CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
     39    if not os.path.isdir(CUSTOM_MODEL_PATH): 
     40        os.makedirs(CUSTOM_MODEL_PATH) 
    4941 
    5042# TODO: refactor composite model support 
     
    148140    used with functions in generate to build the docs or extract model info. 
    149141    """ 
    150     if '@' in model_string: 
    151         parts = model_string.split('@') 
    152         if len(parts) != 2: 
    153             raise ValueError("Use P@S to apply a structure factor S to model P") 
    154         P_info, Q_info = [load_model_info(part) for part in parts] 
    155         return product.make_product_info(P_info, Q_info) 
    156  
    157     product_parts = [] 
    158     addition_parts = [] 
    159  
    160     addition_parts_names = model_string.split('+') 
    161     if len(addition_parts_names) >= 2: 
    162         addition_parts = [load_model_info(part) for part in addition_parts_names] 
    163     elif len(addition_parts_names) == 1: 
    164         product_parts_names = model_string.split('*') 
    165         if len(product_parts_names) >= 2: 
    166             product_parts = [load_model_info(part) for part in product_parts_names] 
    167         elif len(product_parts_names) == 1: 
    168             if "custom." in product_parts_names[0]: 
    169                 # Extract ModelName from "custom.ModelName" 
    170                 pattern = "custom.([A-Za-z0-9_-]+)" 
    171                 result = re.match(pattern, product_parts_names[0]) 
    172                 if result is None: 
    173                     raise ValueError("Model name in invalid format: " + product_parts_names[0]) 
    174                 model_name = result.group(1) 
    175                 # Use ModelName to find the path to the custom model file 
    176                 model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py") 
    177                 if not os.path.isfile(model_path): 
    178                     raise ValueError("The model file {} doesn't exist".format(model_path)) 
    179                 kernel_module = custom.load_custom_kernel_module(model_path) 
    180                 return modelinfo.make_model_info(kernel_module) 
    181             # Model is a core model 
    182             kernel_module = generate.load_kernel_module(product_parts_names[0]) 
    183             return modelinfo.make_model_info(kernel_module) 
    184  
    185     model = None 
    186     if len(product_parts) > 1: 
    187         model = mixture.make_mixture_info(product_parts, operation='*') 
    188     if len(addition_parts) > 1: 
    189         if model is not None: 
    190             addition_parts.append(model) 
    191         model = mixture.make_mixture_info(addition_parts, operation='+') 
    192     return model 
     142    if "+" in model_string: 
     143        parts = [load_model_info(part) 
     144                 for part in model_string.split("+")] 
     145        return mixture.make_mixture_info(parts, operation='+') 
     146    elif "*" in model_string: 
     147        parts = [load_model_info(part) 
     148                 for part in model_string.split("*")] 
     149        return mixture.make_mixture_info(parts, operation='*') 
     150    elif "@" in model_string: 
     151        p_info, q_info = [load_model_info(part) 
     152                          for part in model_string.split("@")] 
     153        return product.make_product_info(p_info, q_info) 
     154    # We are now dealing with a pure model 
     155    elif "custom." in model_string: 
     156        pattern = "custom.([A-Za-z0-9_-]+)" 
     157        result = re.match(pattern, model_string) 
     158        if result is None: 
     159            raise ValueError("Model name in invalid format: " + model_string) 
     160        model_name = result.group(1) 
     161        # Use ModelName to find the path to the custom model file 
     162        model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py") 
     163        if not os.path.isfile(model_path): 
     164            raise ValueError("The model file {} doesn't exist".format(model_path)) 
     165        kernel_module = custom.load_custom_kernel_module(model_path) 
     166        return modelinfo.make_model_info(kernel_module) 
     167    kernel_module = generate.load_kernel_module(model_string) 
     168    return modelinfo.make_model_info(kernel_module) 
    193169 
    194170 
     
    292268    if platform is None: 
    293269        platform = "ocl" 
    294     if platform == "ocl" and not HAVE_OPENCL or not model_info.opencl: 
     270    if not kernelcl.use_opencl() or not model_info.opencl: 
    295271        platform = "dll" 
    296272 
     
    336312    print("\n".join(list_models(kind))) 
    337313 
     314def test_composite_order(): 
     315    def test_models(fst, snd): 
     316        """Confirm that two models produce the same parameters""" 
     317        fst = load_model(fst) 
     318        snd = load_model(snd) 
     319        # Un-disambiguate parameter names so that we can check if the same 
     320        # parameters are in a pair of composite models. Since each parameter in 
     321        # the mixture model is tagged as e.g., A_sld, we ought to use a 
     322        # regex subsitution s/^[A-Z]+_/_/, but removing all uppercase letters 
     323        # is good enough. 
     324        fst = [[x for x in p.name if x == x.lower()] for p in fst.info.parameters.kernel_parameters] 
     325        snd = [[x for x in p.name if x == x.lower()] for p in snd.info.parameters.kernel_parameters] 
     326        assert sorted(fst) == sorted(snd), "{} != {}".format(fst, snd) 
     327 
     328    def build_test(first, second): 
     329        test = lambda description: test_models(first, second) 
     330        description = first + " vs. " + second 
     331        return test, description 
     332 
     333    yield build_test( 
     334        "cylinder+sphere", 
     335        "sphere+cylinder") 
     336    yield build_test( 
     337        "cylinder*sphere", 
     338        "sphere*cylinder") 
     339    yield build_test( 
     340        "cylinder@hardsphere*sphere", 
     341        "sphere*cylinder@hardsphere") 
     342    yield build_test( 
     343        "barbell+sphere*cylinder@hardsphere", 
     344        "sphere*cylinder@hardsphere+barbell") 
     345    yield build_test( 
     346        "barbell+cylinder@hardsphere*sphere", 
     347        "cylinder@hardsphere*sphere+barbell") 
     348    yield build_test( 
     349        "barbell+sphere*cylinder@hardsphere", 
     350        "barbell+cylinder@hardsphere*sphere") 
     351    yield build_test( 
     352        "sphere*cylinder@hardsphere+barbell", 
     353        "cylinder@hardsphere*sphere+barbell") 
     354    yield build_test( 
     355        "barbell+sphere*cylinder@hardsphere", 
     356        "cylinder@hardsphere*sphere+barbell") 
     357    yield build_test( 
     358        "barbell+cylinder@hardsphere*sphere", 
     359        "sphere*cylinder@hardsphere+barbell") 
     360 
     361def test_composite(): 
     362    # type: () -> None 
     363    """Check that model load works""" 
     364    #Test the the model produces the parameters that we would expect 
     365    model = load_model("cylinder@hardsphere*sphere") 
     366    actual = [p.name for p in model.info.parameters.kernel_parameters] 
     367    target = ("sld sld_solvent radius length theta phi volfraction" 
     368              " A_sld A_sld_solvent A_radius").split() 
     369    assert target == actual, "%s != %s"%(target, actual) 
     370 
    338371if __name__ == "__main__": 
    339372    list_models_main() 
  • sasmodels/data.py

    r2d81cfe rd86f0fc  
    706706    else: 
    707707        vmin_scaled, vmax_scaled = vmin, vmax 
    708     plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)), 
     708    #nx, ny = len(data.x_bins), len(data.y_bins) 
     709    x_bins, y_bins, image = _build_matrix(data, plottable) 
     710    plt.imshow(image, 
    709711               interpolation='nearest', aspect=1, origin='lower', 
    710712               extent=[xmin, xmax, ymin, ymax], 
     
    714716    return vmin, vmax 
    715717 
     718 
     719# === The following is modified from sas.sasgui.plottools.PlotPanel 
     720def _build_matrix(self, plottable): 
     721    """ 
     722    Build a matrix for 2d plot from a vector 
     723    Returns a matrix (image) with ~ square binning 
     724    Requirement: need 1d array formats of 
     725    self.data, self.qx_data, and self.qy_data 
     726    where each one corresponds to z, x, or y axis values 
     727 
     728    """ 
     729    # No qx or qy given in a vector format 
     730    if self.qx_data is None or self.qy_data is None \ 
     731            or self.qx_data.ndim != 1 or self.qy_data.ndim != 1: 
     732        return self.x_bins, self.y_bins, plottable 
     733 
     734    # maximum # of loops to fillup_pixels 
     735    # otherwise, loop could never stop depending on data 
     736    max_loop = 1 
     737    # get the x and y_bin arrays. 
     738    x_bins, y_bins = _get_bins(self) 
     739    # set zero to None 
     740 
     741    #Note: Can not use scipy.interpolate.Rbf: 
     742    # 'cause too many data points (>10000)<=JHC. 
     743    # 1d array to use for weighting the data point averaging 
     744    #when they fall into a same bin. 
     745    weights_data = np.ones([self.data.size]) 
     746    # get histogram of ones w/len(data); this will provide 
     747    #the weights of data on each bins 
     748    weights, xedges, yedges = np.histogram2d(x=self.qy_data, 
     749                                             y=self.qx_data, 
     750                                             bins=[y_bins, x_bins], 
     751                                             weights=weights_data) 
     752    # get histogram of data, all points into a bin in a way of summing 
     753    image, xedges, yedges = np.histogram2d(x=self.qy_data, 
     754                                           y=self.qx_data, 
     755                                           bins=[y_bins, x_bins], 
     756                                           weights=plottable) 
     757    # Now, normalize the image by weights only for weights>1: 
     758    # If weight == 1, there is only one data point in the bin so 
     759    # that no normalization is required. 
     760    image[weights > 1] = image[weights > 1] / weights[weights > 1] 
     761    # Set image bins w/o a data point (weight==0) as None (was set to zero 
     762    # by histogram2d.) 
     763    image[weights == 0] = None 
     764 
     765    # Fill empty bins with 8 nearest neighbors only when at least 
     766    #one None point exists 
     767    loop = 0 
     768 
     769    # do while loop until all vacant bins are filled up up 
     770    #to loop = max_loop 
     771    while (weights == 0).any(): 
     772        if loop >= max_loop:  # this protects never-ending loop 
     773            break 
     774        image = _fillup_pixels(image=image, weights=weights) 
     775        loop += 1 
     776 
     777    return x_bins, y_bins, image 
     778 
     779def _get_bins(self): 
     780    """ 
     781    get bins 
     782    set x_bins and y_bins into self, 1d arrays of the index with 
     783    ~ square binning 
     784    Requirement: need 1d array formats of 
     785    self.qx_data, and self.qy_data 
     786    where each one corresponds to  x, or y axis values 
     787    """ 
     788    # find max and min values of qx and qy 
     789    xmax = self.qx_data.max() 
     790    xmin = self.qx_data.min() 
     791    ymax = self.qy_data.max() 
     792    ymin = self.qy_data.min() 
     793 
     794    # calculate the range of qx and qy: this way, it is a little 
     795    # more independent 
     796    x_size = xmax - xmin 
     797    y_size = ymax - ymin 
     798 
     799    # estimate the # of pixels on each axes 
     800    npix_y = int(np.floor(np.sqrt(len(self.qy_data)))) 
     801    npix_x = int(np.floor(len(self.qy_data) / npix_y)) 
     802 
     803    # bin size: x- & y-directions 
     804    xstep = x_size / (npix_x - 1) 
     805    ystep = y_size / (npix_y - 1) 
     806 
     807    # max and min taking account of the bin sizes 
     808    xmax = xmax + xstep / 2.0 
     809    xmin = xmin - xstep / 2.0 
     810    ymax = ymax + ystep / 2.0 
     811    ymin = ymin - ystep / 2.0 
     812 
     813    # store x and y bin centers in q space 
     814    x_bins = np.linspace(xmin, xmax, npix_x) 
     815    y_bins = np.linspace(ymin, ymax, npix_y) 
     816 
     817    return x_bins, y_bins 
     818 
     819def _fillup_pixels(image=None, weights=None): 
     820    """ 
     821    Fill z values of the empty cells of 2d image matrix 
     822    with the average over up-to next nearest neighbor points 
     823 
     824    :param image: (2d matrix with some zi = None) 
     825 
     826    :return: image (2d array ) 
     827 
     828    :TODO: Find better way to do for-loop below 
     829 
     830    """ 
     831    # No image matrix given 
     832    if image is None or np.ndim(image) != 2 \ 
     833            or np.isfinite(image).all() \ 
     834            or weights is None: 
     835        return image 
     836    # Get bin size in y and x directions 
     837    len_y = len(image) 
     838    len_x = len(image[1]) 
     839    temp_image = np.zeros([len_y, len_x]) 
     840    weit = np.zeros([len_y, len_x]) 
     841    # do for-loop for all pixels 
     842    for n_y in range(len(image)): 
     843        for n_x in range(len(image[1])): 
     844            # find only null pixels 
     845            if weights[n_y][n_x] > 0 or np.isfinite(image[n_y][n_x]): 
     846                continue 
     847            else: 
     848                # find 4 nearest neighbors 
     849                # check where or not it is at the corner 
     850                if n_y != 0 and np.isfinite(image[n_y - 1][n_x]): 
     851                    temp_image[n_y][n_x] += image[n_y - 1][n_x] 
     852                    weit[n_y][n_x] += 1 
     853                if n_x != 0 and np.isfinite(image[n_y][n_x - 1]): 
     854                    temp_image[n_y][n_x] += image[n_y][n_x - 1] 
     855                    weit[n_y][n_x] += 1 
     856                if n_y != len_y - 1 and np.isfinite(image[n_y + 1][n_x]): 
     857                    temp_image[n_y][n_x] += image[n_y + 1][n_x] 
     858                    weit[n_y][n_x] += 1 
     859                if n_x != len_x - 1 and np.isfinite(image[n_y][n_x + 1]): 
     860                    temp_image[n_y][n_x] += image[n_y][n_x + 1] 
     861                    weit[n_y][n_x] += 1 
     862                # go 4 next nearest neighbors when no non-zero 
     863                # neighbor exists 
     864                if n_y != 0 and n_x != 0 and \ 
     865                        np.isfinite(image[n_y - 1][n_x - 1]): 
     866                    temp_image[n_y][n_x] += image[n_y - 1][n_x - 1] 
     867                    weit[n_y][n_x] += 1 
     868                if n_y != len_y - 1 and n_x != 0 and \ 
     869                        np.isfinite(image[n_y + 1][n_x - 1]): 
     870                    temp_image[n_y][n_x] += image[n_y + 1][n_x - 1] 
     871                    weit[n_y][n_x] += 1 
     872                if n_y != len_y and n_x != len_x - 1 and \ 
     873                        np.isfinite(image[n_y - 1][n_x + 1]): 
     874                    temp_image[n_y][n_x] += image[n_y - 1][n_x + 1] 
     875                    weit[n_y][n_x] += 1 
     876                if n_y != len_y - 1 and n_x != len_x - 1 and \ 
     877                        np.isfinite(image[n_y + 1][n_x + 1]): 
     878                    temp_image[n_y][n_x] += image[n_y + 1][n_x + 1] 
     879                    weit[n_y][n_x] += 1 
     880 
     881    # get it normalized 
     882    ind = (weit > 0) 
     883    image[ind] = temp_image[ind] / weit[ind] 
     884 
     885    return image 
     886 
     887 
    716888def demo(): 
    717889    # type: () -> None 
  • sasmodels/details.py

    r108e70e r885753a  
    243243    offset = np.cumsum(np.hstack((0, length))) 
    244244    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    245     # Pad value array to a 32 value boundaryd 
     245    # Pad value array to a 32 value boundary 
    246246    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    247247    extra = (32 - data_len%32)%32 
     
    250250    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
    251251    #call_details.show() 
     252    #print("data", data) 
    252253    return call_details, data, is_magnetic 
    253254 
     
    296297    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 
    297298    mag = mag.reshape(-1, 3) 
    298     scale = mag[:, 0] 
    299     if np.any(scale): 
     299    if np.any(mag[:, 0] != 0.0): 
     300        M0 = mag[:, 0].copy() 
    300301        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    301         cos_theta = cos(theta) 
    302         mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
    303         mag[:, 1] = scale*sin(theta)  # my 
    304         mag[:, 2] = -scale*cos_theta*sin(phi)  # mz 
     302        mag[:, 0] = +M0*cos(theta)*cos(phi)  # mx 
     303        mag[:, 1] = +M0*sin(theta) # my 
     304        mag[:, 2] = -M0*cos(theta)*sin(phi)  # mz 
    305305        return True 
    306306    else: 
  • sasmodels/direct_model.py

    r2d81cfe rd18d6dd  
    345345            self._kernel = self._model.make_kernel(self._kernel_inputs) 
    346346 
     347        # Need to pull background out of resolution for multiple scattering 
     348        background = pars.get('background', 0.) 
     349        pars = pars.copy() 
     350        pars['background'] = 0. 
     351 
    347352        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
    348353        # Storing the calculated Iq values so that they can be plotted. 
     
    357362                np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 
    358363            ) 
    359         return result 
     364        return result + background 
    360365 
    361366 
  • sasmodels/generate.py

    r108e70e rd86f0fc  
    169169 
    170170import sys 
    171 from os.path import abspath, dirname, join as joinpath, exists, getmtime 
     171from os import environ 
     172from os.path import abspath, dirname, join as joinpath, exists, getmtime, sep 
    172173import re 
    173174import string 
     
    214215 
    215216EXTERNAL_DIR = 'sasmodels-data' 
    216 DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_template.c') 
     217DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_iq.c') 
    217218MODEL_PATH = joinpath(DATA_PATH, 'models') 
    218219 
     
    289290    import loops. 
    290291    """ 
    291     if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)): 
    292         import os.path 
     292    if info.source and any(lib.startswith('lib/gauss') for lib in info.source): 
    293293        from .gengauss import gengauss 
    294         path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n) 
    295         if not os.path.exists(path): 
     294        path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n) 
     295        if not exists(path): 
    296296            gengauss(n, path) 
    297297        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 
    298                         else lib for lib in info.source] 
     298                       else lib for lib in info.source] 
    299299 
    300300def format_units(units): 
     
    320320                     for w, h in zip(column_widths, PARTABLE_HEADERS)] 
    321321 
    322     sep = " ".join("="*w for w in column_widths) 
     322    underbar = " ".join("="*w for w in column_widths) 
    323323    lines = [ 
    324         sep, 
     324        underbar, 
    325325        " ".join("%-*s" % (w, h) 
    326326                 for w, h in zip(column_widths, PARTABLE_HEADERS)), 
    327         sep, 
     327        underbar, 
    328328        ] 
    329329    for p in pars: 
     
    334334            "%*g" % (column_widths[3], p.default), 
    335335            ])) 
    336     lines.append(sep) 
     336    lines.append(underbar) 
    337337    return "\n".join(lines) 
    338338 
     
    389389    # TODO: fails DRY; templates appear two places. 
    390390    model_templates = [joinpath(DATA_PATH, filename) 
    391                        for filename in ('kernel_header.c', 'kernel_iq.cl')] 
     391                       for filename in ('kernel_header.c', 'kernel_iq.c')] 
    392392    source_files = (model_sources(model_info) 
    393393                    + model_templates 
     
    612612    """ 
    613613    spaces = " "*depth 
    614     sep = "\n" + spaces 
    615     return spaces + sep.join(s.split("\n")) 
     614    interline_separator = "\n" + spaces 
     615    return spaces + interline_separator.join(s.split("\n")) 
    616616 
    617617 
     
    619619def load_template(filename): 
    620620    # type: (str) -> str 
     621    """ 
     622    Load template file from sasmodels resource directory. 
     623    """ 
    621624    path = joinpath(DATA_PATH, filename) 
    622625    mtime = getmtime(path) 
     
    900903        kernel_module = load_custom_kernel_module(model_name) 
    901904    else: 
    902         from sasmodels import models 
    903         __import__('sasmodels.models.'+model_name) 
    904         kernel_module = getattr(models, model_name, None) 
     905        try: 
     906            from sasmodels import models 
     907            __import__('sasmodels.models.'+model_name) 
     908            kernel_module = getattr(models, model_name, None) 
     909        except ImportError: 
     910            # If the model isn't a built in model, try the plugin directory 
     911            plugin_path = environ.get('SAS_MODELPATH', None) 
     912            if plugin_path is not None: 
     913                file_name = model_name.split(sep)[-1] 
     914                model_name = plugin_path + sep + file_name + ".py" 
     915                kernel_module = load_custom_kernel_module(model_name) 
     916            else: 
     917                raise 
    905918    return kernel_module 
    906919 
  • sasmodels/guyou.py

    r9d9fcbd rb3703f5  
    3131# 
    3232# 2017-11-01 Paul Kienzle 
    33 # * converted to python, using degrees rather than radians 
     33# * converted to python, with degrees rather than radians 
     34""" 
     35Convert between latitude-longitude and Guyou map coordinates. 
     36""" 
     37 
    3438from __future__ import division, print_function 
    3539 
    3640import numpy as np 
    37 from numpy import sqrt, pi, tan, cos, sin, log, exp, arctan2 as atan2, sign, radians, degrees 
     41from numpy import sqrt, pi, tan, cos, sin, sign, radians, degrees 
     42from numpy import sinh, arctan as atan 
     43 
     44# scipy version of special functions 
     45from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF 
    3846 
    3947_ = """ 
     
    5967""" 
    6068 
    61 # scipy version of special functions 
    62 from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF 
    63 from numpy import sinh, sign, arctan as atan 
    64  
    6569def ellipticJi(u, v, m): 
    6670    scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m) 
    6771    u, v, m = np.broadcast_arrays(u, v, m) 
    68     result = np.empty_like([u,u,u], 'D') 
    69     real = v==0 
    70     imag = u==0 
     72    result = np.empty_like([u, u, u], 'D') 
     73    real = (v == 0) 
     74    imag = (u == 0) 
    7175    mixed = ~(real|imag) 
    7276    result[:, real] = _ellipticJi_real(u[real], m[real]) 
    7377    result[:, imag] = _ellipticJi_imag(v[imag], m[imag]) 
    7478    result[:, mixed] = _ellipticJi(u[mixed], v[mixed], m[mixed]) 
    75     return result[0,:] if scalar else result 
     79    return result[0, :] if scalar else result 
    7680 
    7781def _ellipticJi_real(u, m): 
     
    104108        result = np.empty_like(phi, 'D') 
    105109        index = (phi == 0) 
    106         result[index] = ellipticF(atan(sinh(abs(phi[index]))), 1-m[index]) * sign(psi[index]) 
     110        result[index] = ellipticF(atan(sinh(abs(phi[index]))), 
     111                                  1-m[index]) * sign(psi[index]) 
    107112        result[~index] = ellipticFi(phi[~index], psi[~index], m[~index]) 
    108113        return result.reshape(1)[0] if scalar else result 
     
    117122    cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2 
    118123    re = ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi) 
    119     im = ellipticF(atan(sqrt(np.maximum(0,(cotlambda2 / cotphi2 - 1) / m))), 1 - m) * sign(psi) 
     124    im = ellipticF(atan(sqrt(np.maximum(0, (cotlambda2 / cotphi2 - 1) / m))), 
     125                   1 - m) * sign(psi) 
    120126    return re + 1j*im 
    121127 
    122 sqrt2 = sqrt(2) 
     128SQRT2 = sqrt(2) 
    123129 
    124130# [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion 
     
    127133# K = 3.165103454447431823666270142140819753058976299237578486994... 
    128134def guyou(lam, phi): 
     135    """Transform from (latitude, longitude) to point (x, y)""" 
    129136    # [PAK] wrap into [-pi/2, pi/2] radians 
    130137    x, y = np.asarray(lam), np.asarray(phi) 
     
    135142 
    136143    # Compute constant K 
    137     cos_u = (sqrt2 - 1) / (sqrt2 + 1) 
     144    cos_u = (SQRT2 - 1) / (SQRT2 + 1) 
    138145    sinsq_u = 1 - cos_u**2 
    139146    K = ellipticF(pi/2, sinsq_u) 
     
    144151    at = atan(r * (cos(lam) - 1j*sin(lam))) 
    145152    t = ellipticFi(at.real, at.imag, sinsq_u) 
    146     x, y = (-t.imag, sign(phi + (phi==0))*(0.5 * K - t.real)) 
     153    x, y = (-t.imag, sign(phi + (phi == 0))*(0.5 * K - t.real)) 
    147154 
    148155    # [PAK] convert to degrees, and return to original tile 
     
    150157 
    151158def guyou_invert(x, y): 
     159    """Transform from point (x, y) on plot to (latitude, longitude)""" 
    152160    # [PAK] wrap into [-pi/2, pi/2] radians 
    153161    x, y = np.asarray(x), np.asarray(y) 
     
    158166 
    159167    # compute constant K 
    160     cos_u = (sqrt2 - 1) / (sqrt2 + 1) 
     168    cos_u = (SQRT2 - 1) / (SQRT2 + 1) 
    161169    sinsq_u = 1 - cos_u**2 
    162170    K = ellipticF(pi/2, sinsq_u) 
     
    174182 
    175183def plot_grid(): 
     184    """Plot the latitude-longitude grid for Guyou transform""" 
    176185    import matplotlib.pyplot as plt 
    177186    from numpy import linspace 
     
    186195        plt.plot(x, y, 'g') 
    187196 
    188     for long in range(-limit, limit+1, step): 
    189         x, y = guyou(scale*long, scale*long_line) 
     197    for longitude in range(-limit, limit+1, step): 
     198        x, y = guyou(scale*longitude, scale*long_line) 
    190199        plt.plot(x, y, 'b') 
    191200    #plt.xlabel('longitude') 
    192201    plt.ylabel('latitude') 
     202    plt.title('forward transform') 
    193203 
    194204    plt.subplot(212) 
     
    202212    plt.xlabel('longitude') 
    203213    plt.ylabel('latitude') 
     214    plt.title('inverse transform') 
     215 
     216def main(): 
     217    """Show the Guyou transformation""" 
     218    plot_grid() 
     219    import matplotlib.pyplot as plt 
     220    plt.show() 
    204221 
    205222if __name__ == "__main__": 
    206     plot_grid() 
    207     import matplotlib.pyplot as plt; plt.show() 
    208  
     223    main() 
    209224 
    210225_ = """ 
  • sasmodels/jitter.py

    r8cfb486 rb3703f5  
    11#!/usr/bin/env python 
    22""" 
    3 Application to explore the difference between sasview 3.x orientation 
    4 dispersity and possible replacement algorithms. 
     3Jitter Explorer 
     4=============== 
     5 
     6Application to explore orientation angle and angular dispersity. 
    57""" 
    68from __future__ import division, print_function 
    79 
    8 import sys, os 
    9 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
    10  
    1110import argparse 
    1211 
    13 import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
     12try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d 
     13    import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot 
     14except ImportError: 
     15    pass 
     16 
    1417import matplotlib.pyplot as plt 
    15 from matplotlib.widgets import Slider, CheckButtons 
    16 from matplotlib import cm 
     18from matplotlib.widgets import Slider 
    1719import numpy as np 
    1820from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1921 
    20 def draw_beam(ax, view=(0, 0)): 
     22def draw_beam(axes, view=(0, 0)): 
    2123    """ 
    2224    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
    2325    """ 
    24     #ax.plot([0,0],[0,0],[1,-1]) 
    25     #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
     26    #axes.plot([0,0],[0,0],[1,-1]) 
     27    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    2628 
    2729    steps = 25 
     
    4042    x, y, z = [v.reshape(shape) for v in points] 
    4143 
    42     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    43  
    44 def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)): 
     44    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
     45 
     46def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1): 
     47    """Draw an ellipsoid.""" 
     48    a, b, c = size 
     49    u = np.linspace(0, 2 * np.pi, steps) 
     50    v = np.linspace(0, np.pi, steps) 
     51    x = a*np.outer(np.cos(u), np.sin(v)) 
     52    y = b*np.outer(np.sin(u), np.sin(v)) 
     53    z = c*np.outer(np.ones_like(u), np.cos(v)) 
     54    x, y, z = transform_xyz(view, jitter, x, y, z) 
     55 
     56    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
     57 
     58    draw_labels(axes, view, jitter, [ 
     59        ('c+', [+0, +0, +c], [+1, +0, +0]), 
     60        ('c-', [+0, +0, -c], [+0, +0, -1]), 
     61        ('a+', [+a, +0, +0], [+0, +0, +1]), 
     62        ('a-', [-a, +0, +0], [+0, +0, -1]), 
     63        ('b+', [+0, +b, +0], [-1, +0, +0]), 
     64        ('b-', [+0, -b, +0], [-1, +0, +0]), 
     65    ]) 
     66 
     67def draw_sc(axes, size, view, jitter, steps=None, alpha=1): 
     68    """Draw points for simple cubic paracrystal""" 
     69    atoms = _build_sc() 
     70    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     71 
     72def draw_fcc(axes, size, view, jitter, steps=None, alpha=1): 
     73    """Draw points for face-centered cubic paracrystal""" 
     74    # Build the simple cubic crystal 
     75    atoms = _build_sc() 
     76    # Define the centers for each face 
     77    # x planes at -1, 0, 1 have four centers per plane, at +/- 0.5 in y and z 
     78    x, y, z = ( 
     79        [-1]*4 + [0]*4 + [1]*4, 
     80        ([-0.5]*2 + [0.5]*2)*3, 
     81        [-0.5, 0.5]*12, 
     82    ) 
     83    # y and z planes can be generated by substituting x for y and z respectively 
     84    atoms.extend(zip(x+y+z, y+z+x, z+x+y)) 
     85    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     86 
     87def draw_bcc(axes, size, view, jitter, steps=None, alpha=1): 
     88    """Draw points for body-centered cubic paracrystal""" 
     89    # Build the simple cubic crystal 
     90    atoms = _build_sc() 
     91    # Define the centers for each octant 
     92    # x plane at +/- 0.5 have four centers per plane at +/- 0.5 in y and z 
     93    x, y, z = ( 
     94        [-0.5]*4 + [0.5]*4, 
     95        ([-0.5]*2 + [0.5]*2)*2, 
     96        [-0.5, 0.5]*8, 
     97    ) 
     98    atoms.extend(zip(x, y, z)) 
     99    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     100 
     101def _draw_crystal(axes, size, view, jitter, atoms=None): 
     102    atoms, size = np.asarray(atoms, 'd').T, np.asarray(size, 'd') 
     103    x, y, z = atoms*size[:, None] 
     104    x, y, z = transform_xyz(view, jitter, x, y, z) 
     105    axes.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o') 
     106    axes.scatter(x[1:], y[1:], z[1:], c='r', marker='o') 
     107 
     108def _build_sc(): 
     109    # three planes of 9 dots for x at -1, 0 and 1 
     110    x, y, z = ( 
     111        [-1]*9 + [0]*9 + [1]*9, 
     112        ([-1]*3 + [0]*3 + [1]*3)*3, 
     113        [-1, 0, 1]*9, 
     114    ) 
     115    atoms = list(zip(x, y, z)) 
     116    #print(list(enumerate(atoms))) 
     117    # Pull the dot at (0, 0, 1) to the front of the list 
     118    # It will be highlighted in the view 
     119    index = 14 
     120    highlight = atoms[index] 
     121    del atoms[index] 
     122    atoms.insert(0, highlight) 
     123    return atoms 
     124 
     125def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1): 
     126    """Draw a parallelepiped.""" 
     127    a, b, c = size 
     128    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     129    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     130    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
     131    tri = np.array([ 
     132        # counter clockwise triangles 
     133        # z: up/down, x: right/left, y: front/back 
     134        [0, 1, 2], [3, 2, 1], # top face 
     135        [6, 5, 4], [5, 6, 7], # bottom face 
     136        [0, 2, 6], [6, 4, 0], # right face 
     137        [1, 5, 7], [7, 3, 1], # left face 
     138        [2, 3, 6], [7, 6, 3], # front face 
     139        [4, 1, 0], [5, 1, 4], # back face 
     140    ]) 
     141 
     142    x, y, z = transform_xyz(view, jitter, x, y, z) 
     143    axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
     144 
     145    # Draw pink face on box. 
     146    # Since I can't control face color, instead draw a thin box situated just 
     147    # in front of the "c+" face.  Use the c face so that rotations about psi 
     148    # rotate that face. 
     149    if 1: 
     150        x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     151        y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     152        z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
     153        x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) 
     154        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) 
     155 
     156    draw_labels(axes, view, jitter, [ 
     157        ('c+', [+0, +0, +c], [+1, +0, +0]), 
     158        ('c-', [+0, +0, -c], [+0, +0, -1]), 
     159        ('a+', [+a, +0, +0], [+0, +0, +1]), 
     160        ('a-', [-a, +0, +0], [+0, +0, -1]), 
     161        ('b+', [+0, +b, +0], [-1, +0, +0]), 
     162        ('b-', [+0, -b, +0], [-1, +0, +0]), 
     163    ]) 
     164 
     165def draw_sphere(axes, radius=10., steps=100): 
     166    """Draw a sphere""" 
     167    u = np.linspace(0, 2 * np.pi, steps) 
     168    v = np.linspace(0, np.pi, steps) 
     169 
     170    x = radius * np.outer(np.cos(u), np.sin(v)) 
     171    y = radius * np.outer(np.sin(u), np.sin(v)) 
     172    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
     173    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
     174 
     175def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), 
     176                draw_shape=draw_parallelepiped): 
    45177    """ 
    46178    Represent jitter as a set of shapes at different orientations. 
     
    49181    scale = 0.95/sqrt(sum(v**2 for v in size)) 
    50182    size = tuple(scale*v for v in size) 
    51     draw_shape = draw_parallelepiped 
    52     #draw_shape = draw_ellipsoid 
    53183 
    54184    #np.random.seed(10) 
     
    56186    cloud = [ 
    57187        [-1, -1, -1], 
    58         [-1, -1,  0], 
    59         [-1, -1,  1], 
    60         [-1,  0, -1], 
    61         [-1,  0,  0], 
    62         [-1,  0,  1], 
    63         [-1,  1, -1], 
    64         [-1,  1,  0], 
    65         [-1,  1,  1], 
    66         [ 0, -1, -1], 
    67         [ 0, -1,  0], 
    68         [ 0, -1,  1], 
    69         [ 0,  0, -1], 
    70         [ 0,  0,  0], 
    71         [ 0,  0,  1], 
    72         [ 0,  1, -1], 
    73         [ 0,  1,  0], 
    74         [ 0,  1,  1], 
    75         [ 1, -1, -1], 
    76         [ 1, -1,  0], 
    77         [ 1, -1,  1], 
    78         [ 1,  0, -1], 
    79         [ 1,  0,  0], 
    80         [ 1,  0,  1], 
    81         [ 1,  1, -1], 
    82         [ 1,  1,  0], 
    83         [ 1,  1,  1], 
     188        [-1, -1, +0], 
     189        [-1, -1, +1], 
     190        [-1, +0, -1], 
     191        [-1, +0, +0], 
     192        [-1, +0, +1], 
     193        [-1, +1, -1], 
     194        [-1, +1, +0], 
     195        [-1, +1, +1], 
     196        [+0, -1, -1], 
     197        [+0, -1, +0], 
     198        [+0, -1, +1], 
     199        [+0, +0, -1], 
     200        [+0, +0, +0], 
     201        [+0, +0, +1], 
     202        [+0, +1, -1], 
     203        [+0, +1, +0], 
     204        [+0, +1, +1], 
     205        [+1, -1, -1], 
     206        [+1, -1, +0], 
     207        [+1, -1, +1], 
     208        [+1, +0, -1], 
     209        [+1, +0, +0], 
     210        [+1, +0, +1], 
     211        [+1, +1, -1], 
     212        [+1, +1, +0], 
     213        [+1, +1, +1], 
    84214    ] 
    85215    dtheta, dphi, dpsi = jitter 
     
    90220    if dpsi == 0: 
    91221        cloud = [v for v in cloud if v[2] == 0] 
    92     draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
     222    draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    93223    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    94224    for point in cloud: 
    95225        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
    96         draw_shape(ax, size, view, delta, alpha=0.8) 
     226        draw_shape(axes, size, view, delta, alpha=0.8) 
    97227    for v in 'xyz': 
    98228        a, b, c = size 
    99         lim = np.sqrt(a**2+b**2+c**2) 
    100         getattr(ax, 'set_'+v+'lim')([-lim, lim]) 
    101         getattr(ax, v+'axis').label.set_text(v) 
    102  
    103 def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
    104     """Draw an ellipsoid.""" 
    105     a,b,c = size 
    106     u = np.linspace(0, 2 * np.pi, steps) 
    107     v = np.linspace(0, np.pi, steps) 
    108     x = a*np.outer(np.cos(u), np.sin(v)) 
    109     y = b*np.outer(np.sin(u), np.sin(v)) 
    110     z = c*np.outer(np.ones_like(u), np.cos(v)) 
    111     x, y, z = transform_xyz(view, jitter, x, y, z) 
    112  
    113     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    114  
    115     draw_labels(ax, view, jitter, [ 
    116          ('c+', [ 0, 0, c], [ 1, 0, 0]), 
    117          ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
    118          ('a+', [ a, 0, 0], [ 0, 0, 1]), 
    119          ('a-', [-a, 0, 0], [ 0, 0,-1]), 
    120          ('b+', [ 0, b, 0], [-1, 0, 0]), 
    121          ('b-', [ 0,-b, 0], [-1, 0, 0]), 
    122     ]) 
    123  
    124 def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
    125     """Draw a parallelepiped.""" 
    126     a,b,c = size 
    127     x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    128     y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    129     z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
    130     tri = np.array([ 
    131         # counter clockwise triangles 
    132         # z: up/down, x: right/left, y: front/back 
    133         [0,1,2], [3,2,1], # top face 
    134         [6,5,4], [5,6,7], # bottom face 
    135         [0,2,6], [6,4,0], # right face 
    136         [1,5,7], [7,3,1], # left face 
    137         [2,3,6], [7,6,3], # front face 
    138         [4,1,0], [5,1,4], # back face 
    139     ]) 
    140  
    141     x, y, z = transform_xyz(view, jitter, x, y, z) 
    142     ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    143  
    144     draw_labels(ax, view, jitter, [ 
    145          ('c+', [ 0, 0, c], [ 1, 0, 0]), 
    146          ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
    147          ('a+', [ a, 0, 0], [ 0, 0, 1]), 
    148          ('a-', [-a, 0, 0], [ 0, 0,-1]), 
    149          ('b+', [ 0, b, 0], [-1, 0, 0]), 
    150          ('b-', [ 0,-b, 0], [-1, 0, 0]), 
    151     ]) 
    152  
    153 def draw_sphere(ax, radius=10., steps=100): 
    154     """Draw a sphere""" 
    155     u = np.linspace(0, 2 * np.pi, steps) 
    156     v = np.linspace(0, np.pi, steps) 
    157  
    158     x = radius * np.outer(np.cos(u), np.sin(v)) 
    159     y = radius * np.outer(np.sin(u), np.sin(v)) 
    160     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    161     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
     229        lim = np.sqrt(a**2 + b**2 + c**2) 
     230        getattr(axes, 'set_'+v+'lim')([-lim, lim]) 
     231        getattr(axes, v+'axis').label.set_text(v) 
    162232 
    163233PROJECTIONS = [ 
     
    167237    'azimuthal_equal_area', 
    168238] 
    169 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 
     239def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
    170240              projection='equirectangular'): 
    171241    """ 
     
    226296        Should allow free movement in theta, but phi is distorted. 
    227297    """ 
    228     theta, phi, psi = view 
    229     dtheta, dphi, dpsi = jitter 
    230  
    231     t = np.linspace(-1, 1, n) 
    232     weights = np.ones_like(t) 
     298    # TODO: try Kent distribution instead of a gaussian warped by projection 
     299 
     300    dist_x = np.linspace(-1, 1, n) 
     301    weights = np.ones_like(dist_x) 
    233302    if dist == 'gaussian': 
    234         t *= 3 
    235         weights = exp(-0.5*t**2) 
     303        dist_x *= 3 
     304        weights = exp(-0.5*dist_x**2) 
    236305    elif dist == 'rectangle': 
    237306        # Note: uses sasmodels ridiculous definition of rectangle width 
    238         t *= sqrt(3) 
     307        dist_x *= sqrt(3) 
    239308    elif dist == 'uniform': 
    240309        pass 
     
    243312 
    244313    if projection == 'equirectangular':  #define PROJECTION 1 
    245         def rotate(theta_i, phi_j): 
     314        def _rotate(theta_i, phi_j): 
    246315            return Rx(phi_j)*Ry(theta_i) 
    247         def weight(theta_i, phi_j, wi, wj): 
    248             return wi*wj*abs(cos(radians(theta_i))) 
     316        def _weight(theta_i, phi_j, w_i, w_j): 
     317            return w_i*w_j*abs(cos(radians(theta_i))) 
    249318    elif projection == 'sinusoidal':  #define PROJECTION 2 
    250         def rotate(theta_i, phi_j): 
     319        def _rotate(theta_i, phi_j): 
    251320            latitude = theta_i 
    252321            scale = cos(radians(latitude)) 
     
    254323            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    255324            return Rx(longitude)*Ry(latitude) 
    256         def weight(theta_i, phi_j, wi, wj): 
     325        def _weight(theta_i, phi_j, w_i, w_j): 
    257326            latitude = theta_i 
    258327            scale = cos(radians(latitude)) 
    259             w = 1 if abs(phi_j) < abs(scale)*180 else 0 
    260             return w*wi*wj 
     328            active = 1 if abs(phi_j) < abs(scale)*180 else 0 
     329            return active*w_i*w_j 
    261330    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
    262         def rotate(theta_i, phi_j): 
    263             from guyou import guyou_invert 
     331        def _rotate(theta_i, phi_j): 
     332            from .guyou import guyou_invert 
    264333            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
    265334            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
    266335            return Rx(longitude[0])*Ry(latitude[0]) 
    267         def weight(theta_i, phi_j, wi, wj): 
    268             return wi*wj 
     336        def _weight(theta_i, phi_j, w_i, w_j): 
     337            return w_i*w_j 
    269338    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
    270         def rotate(theta_i, phi_j): 
     339        def _rotate(theta_i, phi_j): 
    271340            latitude = sqrt(theta_i**2 + phi_j**2) 
    272341            longitude = degrees(np.arctan2(phi_j, theta_i)) 
    273342            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    274343            return Rz(longitude)*Ry(latitude) 
    275         def weight(theta_i, phi_j, wi, wj): 
     344        def _weight(theta_i, phi_j, w_i, w_j): 
    276345            # Weighting for each point comes from the integral: 
    277346            #     \int\int I(q, lat, log) sin(lat) dlat dlog 
     
    303372            # the entire sphere, and treats theta and phi identically. 
    304373            latitude = sqrt(theta_i**2 + phi_j**2) 
    305             w = sin(radians(latitude))/latitude if latitude != 0 else 1 
    306             return w*wi*wj if latitude < 180 else 0 
     374            weight = sin(radians(latitude))/latitude if latitude != 0 else 1 
     375            return weight*w_i*w_j if latitude < 180 else 0 
    307376    elif projection == 'azimuthal_equal_area': 
    308         def rotate(theta_i, phi_j): 
    309             R = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
    310             latitude = 180-degrees(2*np.arccos(R)) 
     377        def _rotate(theta_i, phi_j): 
     378            radius = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
     379            latitude = 180-degrees(2*np.arccos(radius)) 
    311380            longitude = degrees(np.arctan2(phi_j, theta_i)) 
    312381            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    313382            return Rz(longitude)*Ry(latitude) 
    314         def weight(theta_i, phi_j, wi, wj): 
     383        def _weight(theta_i, phi_j, w_i, w_j): 
    315384            latitude = sqrt(theta_i**2 + phi_j**2) 
    316             w = sin(radians(latitude))/latitude if latitude != 0 else 1 
    317             return w*wi*wj if latitude < 180 else 0 
     385            weight = sin(radians(latitude))/latitude if latitude != 0 else 1 
     386            return weight*w_i*w_j if latitude < 180 else 0 
    318387    else: 
    319388        raise ValueError("unknown projection %r"%projection) 
    320389 
    321390    # mesh in theta, phi formed by rotating z 
     391    dtheta, dphi, dpsi = jitter 
    322392    z = np.matrix([[0], [0], [radius]]) 
    323     points = np.hstack([rotate(theta_i, phi_j)*z 
    324                         for theta_i in dtheta*t 
    325                         for phi_j in dphi*t]) 
    326     # select just the active points (i.e., those with phi < 180 
    327     w = np.array([weight(theta_i, phi_j, wi, wj) 
    328                   for wi, theta_i in zip(weights, dtheta*t) 
    329                   for wj, phi_j in zip(weights, dphi*t)]) 
    330     #print(max(w), min(w), min(w[w>0])) 
    331     points = points[:, w>0] 
    332     w = w[w>0] 
    333     w /= max(w) 
    334  
    335     if 0: # Kent distribution 
    336         points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t]) 
    337         xp, yp, zp = [np.array(v).flatten() for v in points] 
    338         kappa = max(1e6, radians(dtheta)/(2*pi)) 
    339         beta = 1/max(1e-6, radians(dphi)/(2*pi))/kappa 
    340         w = exp(kappa*zp) #+ beta*(xp**2 + yp**2) 
    341         print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w)) 
    342         #w /= abs(cos(radians( 
    343         #w /= sum(w) 
     393    points = np.hstack([_rotate(theta_i, phi_j)*z 
     394                        for theta_i in dtheta*dist_x 
     395                        for phi_j in dphi*dist_x]) 
     396    dist_w = np.array([_weight(theta_i, phi_j, w_i, w_j) 
     397                       for w_i, theta_i in zip(weights, dtheta*dist_x) 
     398                       for w_j, phi_j in zip(weights, dphi*dist_x)]) 
     399    #print(max(dist_w), min(dist_w), min(dist_w[dist_w > 0])) 
     400    points = points[:, dist_w > 0] 
     401    dist_w = dist_w[dist_w > 0] 
     402    dist_w /= max(dist_w) 
    344403 
    345404    # rotate relative to beam 
     
    348407    x, y, z = [np.array(v).flatten() for v in points] 
    349408    #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51)) 
    350     ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    351  
    352 def draw_labels(ax, view, jitter, text): 
     409    axes.scatter(x, y, z, c=dist_w, marker='o', vmin=0., vmax=1.) 
     410 
     411def draw_labels(axes, view, jitter, text): 
    353412    """ 
    354413    Draw text at a particular location. 
     
    364423    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
    365424        zdir = np.asarray(zdir).flatten() 
    366         ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     425        axes.text(p[0], p[1], p[2], label, zdir=zdir) 
    367426 
    368427# Definition of rotation matrices comes from wikipedia: 
     
    370429def Rx(angle): 
    371430    """Construct a matrix to rotate points about *x* by *angle* degrees.""" 
    372     a = radians(angle) 
    373     R = [[1, 0, 0], 
    374          [0, +cos(a), -sin(a)], 
    375          [0, +sin(a), +cos(a)]] 
    376     return np.matrix(R) 
     431    angle = radians(angle) 
     432    rot = [[1, 0, 0], 
     433           [0, +cos(angle), -sin(angle)], 
     434           [0, +sin(angle), +cos(angle)]] 
     435    return np.matrix(rot) 
    377436 
    378437def Ry(angle): 
    379438    """Construct a matrix to rotate points about *y* by *angle* degrees.""" 
    380     a = radians(angle) 
    381     R = [[+cos(a), 0, +sin(a)], 
    382          [0, 1, 0], 
    383          [-sin(a), 0, +cos(a)]] 
    384     return np.matrix(R) 
     439    angle = radians(angle) 
     440    rot = [[+cos(angle), 0, +sin(angle)], 
     441           [0, 1, 0], 
     442           [-sin(angle), 0, +cos(angle)]] 
     443    return np.matrix(rot) 
    385444 
    386445def Rz(angle): 
    387446    """Construct a matrix to rotate points about *z* by *angle* degrees.""" 
    388     a = radians(angle) 
    389     R = [[+cos(a), -sin(a), 0], 
    390          [+sin(a), +cos(a), 0], 
    391          [0, 0, 1]] 
    392     return np.matrix(R) 
     447    angle = radians(angle) 
     448    rot = [[+cos(angle), -sin(angle), 0], 
     449           [+sin(angle), +cos(angle), 0], 
     450           [0, 0, 1]] 
     451    return np.matrix(rot) 
    393452 
    394453def transform_xyz(view, jitter, x, y, z): 
     
    398457    x, y, z = [np.asarray(v) for v in (x, y, z)] 
    399458    shape = x.shape 
    400     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     459    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
    401460    points = apply_jitter(jitter, points) 
    402461    points = orient_relative_to_beam(view, points) 
     
    455514        return data[offset], data[-1] 
    456515 
    457 def draw_scattering(calculator, ax, view, jitter, dist='gaussian'): 
     516def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): 
    458517    """ 
    459518    Plot the scattering for the particular view. 
    460519 
    461     *calculator* is returned from :func:`build_model`.  *ax* are the 3D axes 
     520    *calculator* is returned from :func:`build_model`.  *axes* are the 3D axes 
    462521    on which the data will be plotted.  *view* and *jitter* are the current 
    463522    orientation and orientation dispersity.  *dist* is one of the sasmodels 
     
    472531    theta, phi, psi = view 
    473532    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter] 
    474     theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)] 
     533    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd > 0, phi_pd > 0, psi_pd > 0)] 
    475534    ## increase pd_n for testing jitter integration rather than simple viz 
    476535    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] 
     
    500559    if 0: 
    501560        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
    502         level[level<0] = 0 
     561        level[level < 0] = 0 
    503562        colors = plt.get_cmap()(level) 
    504         ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
     563        axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
    505564    elif 1: 
    506         ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
    507                     levels=np.linspace(vmin, vmax, 24)) 
     565        axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
     566                      levels=np.linspace(vmin, vmax, 24)) 
    508567    else: 
    509         ax.pcolormesh(qx, qy, Iqxy) 
     568        axes.pcolormesh(qx, qy, Iqxy) 
    510569 
    511570def build_model(model_name, n=150, qmax=0.5, **pars): 
     
    524583    for details. 
    525584    """ 
    526     from sasmodels.core import load_model_info, build_model 
     585    from sasmodels.core import load_model_info, build_model as build_sasmodel 
    527586    from sasmodels.data import empty_data2D 
    528587    from sasmodels.direct_model import DirectModel 
    529588 
    530589    model_info = load_model_info(model_name) 
    531     model = build_model(model_info) #, dtype='double!') 
     590    model = build_sasmodel(model_info) #, dtype='double!') 
    532591    q = np.linspace(-qmax, qmax, n) 
    533592    data = empty_data2D(q, q) 
     
    547606    return calculator 
    548607 
    549 def select_calculator(model_name, n=150, size=(10,40,100)): 
     608def select_calculator(model_name, n=150, size=(10, 40, 100)): 
    550609    """ 
    551610    Create a model calculator for the given shape. 
     
    561620    """ 
    562621    a, b, c = size 
     622    d_factor = 0.06  # for paracrystal models 
    563623    if model_name == 'sphere': 
    564624        calculator = build_model('sphere', n=n, radius=c) 
    565625        a = b = c 
     626    elif model_name == 'sc_paracrystal': 
     627        a = b = c 
     628        dnn = c 
     629        radius = 0.5*c 
     630        calculator = build_model('sc_paracrystal', n=n, dnn=dnn, 
     631                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     632                                 background=0) 
     633    elif model_name == 'fcc_paracrystal': 
     634        a = b = c 
     635        # nearest neigbour distance dnn should be 2 radius, but I think the 
     636        # model uses lattice spacing rather than dnn in its calculations 
     637        dnn = 0.5*c 
     638        radius = sqrt(2)/4 * c 
     639        calculator = build_model('fcc_paracrystal', n=n, dnn=dnn, 
     640                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     641                                 background=0) 
    566642    elif model_name == 'bcc_paracrystal': 
    567         calculator = build_model('bcc_paracrystal', n=n, dnn=c, 
    568                                   d_factor=0.06, radius=40) 
    569643        a = b = c 
     644        # nearest neigbour distance dnn should be 2 radius, but I think the 
     645        # model uses lattice spacing rather than dnn in its calculations 
     646        dnn = 0.5*c 
     647        radius = sqrt(3)/2 * c 
     648        calculator = build_model('bcc_paracrystal', n=n, dnn=dnn, 
     649                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     650                                 background=0) 
    570651    elif model_name == 'cylinder': 
    571652        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) 
     
    588669 
    589670SHAPES = [ 
    590     'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal', 
    591     'cylinder', 'ellipsoid', 
    592     'sphere', 
    593  ] 
     671    'parallelepiped', 
     672    'sphere', 'ellipsoid', 'triaxial_ellipsoid', 
     673    'cylinder', 
     674    'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal', 
     675] 
     676 
     677DRAW_SHAPES = { 
     678    'fcc_paracrystal': draw_fcc, 
     679    'bcc_paracrystal': draw_bcc, 
     680    'sc_paracrystal': draw_sc, 
     681    'parallelepiped': draw_parallelepiped, 
     682} 
    594683 
    595684DISTRIBUTIONS = [ 
     
    608697    Show an interactive orientation and jitter demo. 
    609698 
    610     *model_name* is one of the models available in :func:`select_model`. 
     699    *model_name* is one of: sphere, ellipsoid, triaxial_ellipsoid, 
     700    parallelepiped, cylinder, or sc/fcc/bcc_paracrystal 
     701 
     702    *size* gives the dimensions (a, b, c) of the shape. 
     703 
     704    *dist* is the type of dispersition: gaussian, rectangle, or uniform. 
     705 
     706    *mesh* is the number of points in the dispersion mesh. 
     707 
     708    *projection* is the map projection to use for the mesh: equirectangular, 
     709    sinusoidal, guyou, azimuthal_equidistance, or azimuthal_equal_area. 
    611710    """ 
    612711    # projection number according to 1-order position in list, but 
     
    620719    # set up calculator 
    621720    calculator, size = select_calculator(model_name, n=150, size=size) 
     721    draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped) 
    622722 
    623723    ## uncomment to set an independent the colour range for every view 
     
    639739    plt.gcf().canvas.set_window_title(projection) 
    640740    #gs = gridspec.GridSpec(2,1,height_ratios=[4,1]) 
    641     #ax = plt.subplot(gs[0], projection='3d') 
    642     ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') 
     741    #axes = plt.subplot(gs[0], projection='3d') 
     742    axes = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') 
    643743    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection 
    644         ax.axis('square') 
     744        axes.axis('square') 
    645745    except Exception: 
    646746        pass 
     
    649749 
    650750    ## add control widgets to plot 
    651     axtheta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    652     axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    653     axpsi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
    654     stheta = Slider(axtheta, 'Theta', -90, 90, valinit=theta) 
    655     sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    656     spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
    657  
    658     axdtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    659     axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    660     axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
     751    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
     752    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
     753    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
     754    stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
     755    sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
     756    spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
     757 
     758    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
     759    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
     760    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
    661761    # Note: using ridiculous definition of rectangle distribution, whose width 
    662762    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
    663763    # the maximum width to 90. 
    664764    dlimit = DIST_LIMITS[dist] 
    665     sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
    666     sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
    667     sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
     765    sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
     766    sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
     767    sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
     768 
    668769 
    669770    ## callback to draw the new view 
     
    675776        limit = [0, 0, 0.5, 5][dims] 
    676777        jitter = [0 if v < limit else v for v in jitter] 
    677         ax.cla() 
    678         draw_beam(ax, (0, 0)) 
    679         draw_jitter(ax, view, jitter, dist=dist, size=size) 
    680         #draw_jitter(ax, view, (0,0,0)) 
    681         draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection) 
    682         draw_scattering(calculator, ax, view, jitter, dist=dist) 
     778        axes.cla() 
     779        draw_beam(axes, (0, 0)) 
     780        draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape) 
     781        #draw_jitter(axes, view, (0,0,0)) 
     782        draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection) 
     783        draw_scattering(calculator, axes, view, jitter, dist=dist) 
    683784        plt.gcf().canvas.draw() 
    684785 
    685786    ## bind control widgets to view updater 
    686     stheta.on_changed(lambda v: update(v,'theta')) 
     787    stheta.on_changed(lambda v: update(v, 'theta')) 
    687788    sphi.on_changed(lambda v: update(v, 'phi')) 
    688789    spsi.on_changed(lambda v: update(v, 'psi')) 
     
    702803        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
    703804        ) 
    704     parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection') 
    705     parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths') 
    706     parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution') 
    707     parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh') 
    708     parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape') 
     805    parser.add_argument('-p', '--projection', choices=PROJECTIONS, 
     806                        default=PROJECTIONS[0], 
     807                        help='coordinate projection') 
     808    parser.add_argument('-s', '--size', type=str, default='10,40,100', 
     809                        help='a,b,c lengths') 
     810    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, 
     811                        default=DISTRIBUTIONS[0], 
     812                        help='jitter distribution') 
     813    parser.add_argument('-m', '--mesh', type=int, default=30, 
     814                        help='#points in theta-phi mesh') 
     815    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], 
     816                        help='oriented shape') 
    709817    opts = parser.parse_args() 
    710818    size = tuple(int(v) for v in opts.size.split(',')) 
  • sasmodels/kernel_iq.c

    r924a119 rdc6f601  
    7979//     du * (m_sigma_y + 1j*m_sigma_z); 
    8080// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag 
    81 static void set_spin_weights(double in_spin, double out_spin, double spins[4]) 
     81static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 
    8282{ 
    8383  in_spin = clip(in_spin, 0.0, 1.0); 
    8484  out_spin = clip(out_spin, 0.0, 1.0); 
    85   spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd 
    86   spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real 
    87   spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real 
    88   spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu 
    89   spins[4] = spins[1]; // du imag 
    90   spins[5] = spins[2]; // ud imag 
     85  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that 
     86  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
     87  // which is likely to be the case for simple models. 
     88  weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd 
     89  weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real 
     90  weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real 
     91  weight[3] = sqrt(in_spin * out_spin);             // uu 
     92  weight[4] = weight[1]; // du.imag 
     93  weight[5] = weight[2]; // ud.imag 
    9194} 
    9295 
    9396// Compute the magnetic sld 
    9497static double mag_sld( 
    95   const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
     98  const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag 
    9699  const double qx, const double qy, 
    97100  const double px, const double py, 
     
    103106    const double perp = qy*mx - qx*my; 
    104107    switch (xs) { 
     108      default: // keep compiler happy; condition ensures xs in [0,1,2,3] 
    105109      case 0: // uu => sld - D M_perpx 
    106110          return sld - px*perp; 
    107       case 1: // ud real => -D M_perpy 
     111      case 1: // ud.real => -D M_perpy 
    108112          return py*perp; 
    109       case 2: // du real => -D M_perpy 
     113      case 2: // du.real => -D M_perpy 
    110114          return py*perp; 
    111       case 3: // dd real => sld + D M_perpx 
     115      case 3: // dd => sld + D M_perpx 
    112116          return sld + px*perp; 
    113117    } 
    114118  } else { 
    115119    if (xs== 4) { 
    116       return -mz;  // ud imag => -D M_perpz 
     120      return -mz;  // du.imag => +D M_perpz 
    117121    } else { // index == 5 
    118       return mz;   // du imag => D M_perpz 
     122      return +mz;  // ud.imag => -D M_perpz 
    119123    } 
    120124  } 
     
    172176// Apply the rotation matrix returned from qac_rotation to the point (qx,qy), 
    173177// returning R*[qx,qy]' = [qa,qc]' 
    174 static double 
     178static void 
    175179qac_apply( 
    176180    QACRotation *rotation, 
     
    245249// Apply the rotation matrix returned from qabc_rotation to the point (qx,qy), 
    246250// returning R*[qx,qy]' = [qa,qb,qc]' 
    247 static double 
     251static void 
    248252qabc_apply( 
    249253    QABCRotation *rotation, 
     
    311315  //     up_angle = values[NUM_PARS+4]; 
    312316  // TODO: could precompute more magnetism parameters before calling the kernel. 
    313   double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
     317  double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
    314318  double cos_mspin, sin_mspin; 
    315   set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
     319  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights); 
    316320  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    317321#endif // MAGNETIC 
     
    659663 
    660664            // loop over uu, ud real, du real, dd, ud imag, du imag 
    661             for (int xs=0; xs<6; xs++) { 
    662               const double xs_weight = spins[xs]; 
     665            for (unsigned int xs=0; xs<6; xs++) { 
     666              const double xs_weight = xs_weights[xs]; 
    663667              if (xs_weight > 1.e-8) { 
    664668                // Since the cross section weight is significant, set the slds 
     
    673677                  local_values.vector[sld_index] = 
    674678                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); 
     679//if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n", 
     680//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); 
    675681                } 
    676682                scattering += xs_weight * CALL_KERNEL(); 
  • sasmodels/kernelcl.py

    r2d81cfe rd86f0fc  
    5858import numpy as np  # type: ignore 
    5959 
     60 
     61# Attempt to setup opencl. This may fail if the opencl package is not 
     62# installed or if it is installed but there are no devices available. 
    6063try: 
    61     #raise NotImplementedError("OpenCL not yet implemented for new kernel template") 
    6264    import pyopencl as cl  # type: ignore 
     65    from pyopencl import mem_flags as mf 
     66    from pyopencl.characterize import get_fast_inaccurate_build_options 
    6367    # Ask OpenCL for the default context so that we know that one exists 
    6468    cl.create_some_context(interactive=False) 
     69    HAVE_OPENCL = True 
     70    OPENCL_ERROR = "" 
    6571except Exception as exc: 
    66     warnings.warn("OpenCL startup failed with ***" 
    67                   + str(exc) + "***; using C compiler instead") 
    68     raise RuntimeError("OpenCL not available") 
    69  
    70 from pyopencl import mem_flags as mf 
    71 from pyopencl.characterize import get_fast_inaccurate_build_options 
     72    HAVE_OPENCL = False 
     73    OPENCL_ERROR = str(exc) 
    7274 
    7375from . import generate 
     
    102104        cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS] 
    103105 
    104 fix_pyopencl_include() 
    105  
     106if HAVE_OPENCL: 
     107    fix_pyopencl_include() 
    106108 
    107109# The max loops number is limited by the amount of local memory available 
     
    128130""" 
    129131 
     132def use_opencl(): 
     133    return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none" 
    130134 
    131135ENV = None 
     136def reset_environment(): 
     137    """ 
     138    Call to create a new OpenCL context, such as after a change to SAS_OPENCL. 
     139    """ 
     140    global ENV 
     141    ENV = GpuEnvironment() if use_opencl() else None 
     142 
    132143def environment(): 
    133144    # type: () -> "GpuEnvironment" 
     
    137148    This provides an OpenCL context and one queue per device. 
    138149    """ 
    139     global ENV 
    140150    if ENV is None: 
    141         ENV = GpuEnvironment() 
     151        if not HAVE_OPENCL: 
     152            raise RuntimeError("OpenCL startup failed with ***" 
     153                               + OPENCL_ERROR + "***; using C compiler instead") 
     154        reset_environment() 
     155        if ENV is None: 
     156            raise RuntimeError("SAS_OPENCL=None in environment") 
    142157    return ENV 
    143158 
     
    193208    Build a model to run on the gpu. 
    194209 
    195     Returns the compiled program and its type.  The returned type will 
    196     be float32 even if the desired type is float64 if any of the 
    197     devices in the context do not support the cl_khr_fp64 extension. 
     210    Returns the compiled program and its type. 
     211 
     212    Raises an error if the desired precision is not available. 
    198213    """ 
    199214    dtype = np.dtype(dtype) 
  • sasmodels/kerneldll.py

    r1ddb794 r1662ebe  
    123123    # add openmp support if not running on a mac 
    124124    if sys.platform != "darwin": 
    125         CC.append("-fopenmp") 
     125        # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9) 
     126        # Shut it off for all unix until we can investigate. 
     127        #CC.append("-fopenmp") 
     128        pass 
    126129    def compile_command(source, output): 
    127130        """unix compiler command""" 
     
    158161        return CC + [source, "-o", output, "-lm"] 
    159162 
    160 # Windows-specific solution 
    161 if os.name == 'nt': 
    162     # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
    163     DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    164     if not os.path.exists(DLL_PATH): 
    165         os.makedirs(DLL_PATH) 
    166 else: 
    167     # Set up the default path for compiled modules. 
    168     DLL_PATH = tempfile.gettempdir() 
     163# Assume the default location of module DLLs is in .sasmodels/compiled_models. 
     164DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    169165 
    170166ALLOW_SINGLE_PRECISION_DLLS = True 
     
    234230 
    235231    Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
    236     The default is the system temporary directory. 
     232    The default is in ~/.sasmodels/compiled_models. 
    237233    """ 
    238234    if dtype == F16: 
     
    251247        need_recompile = dll_time < newest_source 
    252248    if need_recompile: 
     249        # Make sure the DLL path exists 
     250        if not os.path.exists(DLL_PATH): 
     251            os.makedirs(DLL_PATH) 
    253252        basename = splitext(os.path.basename(dll))[0] + "_" 
    254253        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
  • sasmodels/mixture.py

    r2d81cfe recf895e  
    9494            # If model is a sum model, each constituent model gets its own scale parameter 
    9595            scale_prefix = prefix 
    96             if prefix == '' and part.operation == '*': 
     96            if prefix == '' and getattr(part, "operation", '') == '*': 
    9797                # `part` is a composition product model. Find the prefixes of 
    98                 # it's parameters to form a new prefix for the scale, eg: 
    99                 # a model with A*B*C will have ABC_scale 
     98                # it's parameters to form a new prefix for the scale. 
     99                # For example, a model with A*B*C will have ABC_scale. 
    100100                sub_prefixes = [] 
    101101                for param in part.parameters.kernel_parameters: 
     
    233233        return self 
    234234 
    235     def next(self): 
     235    def __next__(self): 
    236236        # type: () -> Tuple[List[Callable], CallDetails, np.ndarray] 
    237237        if self.part_num >= len(self.parts): 
     
    251251 
    252252        return kernel, call_details, values 
     253 
     254    # CRUFT: py2 support 
     255    next = __next__ 
    253256 
    254257    def _part_details(self, info, par_index): 
  • sasmodels/model_test.py

    r2d81cfe r3221de0  
    6262from .exception import annotate_exception 
    6363from .modelinfo import expand_pars 
     64from .kernelcl import use_opencl 
    6465 
    6566# pylint: disable=unused-import 
     
    123124 
    124125        if is_py:  # kernel implemented in python 
    125             test_name = "Model: %s, Kernel: python"%model_name 
     126            test_name = "%s-python"%model_name 
    126127            test_method_name = "test_%s_python" % model_info.id 
    127128            test = ModelTestCase(test_name, model_info, 
     
    134135 
    135136            # test using dll if desired 
    136             if 'dll' in loaders or not core.HAVE_OPENCL: 
    137                 test_name = "Model: %s, Kernel: dll"%model_name 
     137            if 'dll' in loaders or not use_opencl(): 
     138                test_name = "%s-dll"%model_name 
    138139                test_method_name = "test_%s_dll" % model_info.id 
    139140                test = ModelTestCase(test_name, model_info, 
     
    145146 
    146147            # test using opencl if desired and available 
    147             if 'opencl' in loaders and core.HAVE_OPENCL: 
    148                 test_name = "Model: %s, Kernel: OpenCL"%model_name 
     148            if 'opencl' in loaders and use_opencl(): 
     149                test_name = "%s-opencl"%model_name 
    149150                test_method_name = "test_%s_opencl" % model_info.id 
    150151                # Using dtype=None so that the models that are only 
     
    160161 
    161162    return suite 
    162  
    163163 
    164164def _hide_model_case_from_nose(): 
     
    368368 
    369369    # Build a test suite containing just the model 
    370     loaders = ['opencl'] if core.HAVE_OPENCL else ['dll'] 
     370    loaders = ['opencl'] if use_opencl() else ['dll'] 
    371371    models = [model] 
    372372    try: 
     
    425425        verbosity = 1 
    426426    if models and models[0] == 'opencl': 
    427         if not core.HAVE_OPENCL: 
     427        if not use_opencl(): 
    428428            print("opencl is not available") 
    429429            return 1 
     
    435435        models = models[1:] 
    436436    elif models and models[0] == 'opencl_and_dll': 
    437         loaders = ['opencl', 'dll'] if core.HAVE_OPENCL else ['dll'] 
     437        loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    438438        models = models[1:] 
    439439    else: 
    440         loaders = ['opencl', 'dll'] if core.HAVE_OPENCL else ['dll'] 
     440        loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    441441    if not models: 
    442442        print("""\ 
     
    467467    Run "nosetests sasmodels" on the command line to invoke it. 
    468468    """ 
    469     loaders = ['opencl', 'dll'] if core.HAVE_OPENCL else ['dll'] 
     469    loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    470470    tests = make_suite(loaders, ['all']) 
    471     for test_i in tests: 
    472         # In order for nosetest to see the correct test name, need to set 
    473         # the description attribute of the returned function.  Since we 
    474         # can't do this for the returned instance, wrap it in a lambda and 
    475         # set the description on the lambda.  Otherwise we could just do: 
    476         #    yield test_i.run_all 
    477         L = lambda: test_i.run_all() 
    478         L.description = test_i.test_name 
    479         yield L 
     471    def build_test(test): 
     472        # In order for nosetest to show the test name, wrap the test.run_all 
     473        # instance in function that takes the test name as a parameter which 
     474        # will be displayed when the test is run.  Do this as a function so 
     475        # that it properly captures the context for tests that captured and 
     476        # run later.  If done directly in the for loop, then the looping 
     477        # variable test will be shared amongst all the tests, and we will be 
     478        # repeatedly testing vesicle. 
     479 
     480        # Note: in sasview sas.sasgui.perspectives.fitting.gpu_options 
     481        # requires that the test.description field be set. 
     482        wrap = lambda: test.run_all() 
     483        wrap.description = test.test_name 
     484        return wrap 
     485        # The following would work with nosetests and pytest: 
     486        #     return lambda name: test.run_all(), test.test_name 
     487 
     488    for test in tests: 
     489        yield build_test(test) 
    480490 
    481491 
  • sasmodels/modelinfo.py

    r5ab99b7 r1662ebe  
    7474        processed.append(parse_parameter(*p)) 
    7575    partable = ParameterTable(processed) 
     76    partable.check_angles() 
    7677    return partable 
    7778 
     
    426427        # type: (List[Parameter]) -> None 
    427428        self.kernel_parameters = parameters 
    428         self._check_angles() 
    429429        self._set_vector_lengths() 
    430430 
     
    476476        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
    477477 
    478     def _check_angles(self): 
     478    def check_angles(self): 
     479        """ 
     480        Check that orientation angles are theta, phi and possibly psi. 
     481        """ 
    479482        theta = phi = psi = -1 
    480483        for k, p in enumerate(self.kernel_parameters): 
     
    499502            if psi >= 0 and psi != phi+1: 
    500503                raise TypeError("psi must follow phi") 
    501             #if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): 
    502             #    raise TypeError("orientation parameters must appear at the " 
    503             #                    "end of the parameter table") 
     504            if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): 
     505                raise TypeError("orientation parameters must appear at the " 
     506                                "end of the parameter table") 
    504507        elif theta >= 0 or phi >= 0 or psi >= 0: 
    505508            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
  • sasmodels/models/core_shell_bicelle_elliptical.py

    r2d81cfe rfc3ae1b  
    99scattering length densities. The form factor is normalized by the total 
    1010particle volume. 
    11  
    1211 
    1312.. figure:: img/core_shell_bicelle_geometry.png 
     
    9190for further information. 
    9291 
    93  
    9492.. figure:: img/elliptical_cylinder_angle_definition.png 
    9593 
    9694    Definition of the angles for the oriented core_shell_bicelle_elliptical particles. 
    9795 
    98  
     96Model verified using Monte Carlo simulation for 1D and 2D scattering. 
    9997 
    10098References 
     
    108106* **Author:** Richard Heenan **Date:** December 14, 2016 
    109107* **Last Modified by:**  Richard Heenan **Date:** December 14, 2016 
    110 * **Last Reviewed by:**  Richard Heenan BEWARE 2d data yet to be checked **Date:** December 14, 2016 
     108* **Last Reviewed by:**  Paul Kienzle **Date:** Feb 28, 2018 
    111109""" 
    112110 
     
    130128#             ["name", "units", default, [lower, upper], "type", "description"], 
    131129parameters = [ 
    132     ["radius",         "Ang",       30, [0, inf],    "volume",      "Cylinder core radius"], 
    133     ["x_core",        "None",       3,  [0, inf],    "volume",      "axial ratio of core, X = r_polar/r_equatorial"], 
     130    ["radius",         "Ang",       30, [0, inf],    "volume",      "Cylinder core radius r_minor"], 
     131    ["x_core",        "None",       3,  [0, inf],    "volume",      "Axial ratio of core, X = r_major/r_minor"], 
    134132    ["thick_rim",  "Ang",            8, [0, inf],    "volume",      "Rim shell thickness"], 
    135133    ["thick_face", "Ang",           14, [0, inf],    "volume",      "Cylinder face thickness"], 
     
    139137    ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    140138    ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"], 
    141     ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    142     ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
    143     ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"] 
     139    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "Cylinder axis to beam angle"], 
     140    ["phi",         "degrees",    0,    [-360, 360], "orientation", "Rotation about beam"], 
     141    ["psi",         "degrees",    0,    [-360, 360], "orientation", "Rotation about cylinder axis"] 
    144142    ] 
    145143 
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r108e70e rfc3ae1b  
    140140#             ["name", "units", default, [lower, upper], "type", "description"], 
    141141parameters = [ 
    142     ["radius",         "Ang",       30, [0, inf],    "volume",      "Cylinder core radius"], 
    143     ["x_core",        "None",       3,  [0, inf],    "volume",      "axial ratio of core, X = r_polar/r_equatorial"], 
     142    ["radius",         "Ang",       30, [0, inf],    "volume",      "Cylinder core radius r_minor"], 
     143    ["x_core",        "None",       3,  [0, inf],    "volume",      "Axial ratio of core, X = r_major/r_minor"], 
    144144    ["thick_rim",  "Ang",            8, [0, inf],    "volume",      "Rim or belt shell thickness"], 
    145145    ["thick_face", "Ang",           14, [0, inf],    "volume",      "Cylinder face thickness"], 
     
    149149    ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    150150    ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"], 
    151     ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"], 
    152     ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    153     ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
    154     ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"], 
     151    ["sigma",       "Ang",        0,    [0, inf],    "",            "Interfacial roughness"], 
     152    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "Cylinder axis to beam angle"], 
     153    ["phi",         "degrees",    0,    [-360, 360], "orientation", "Rotation about beam"], 
     154    ["psi",         "degrees",    0,    [-360, 360], "orientation", "Rotation about cylinder axis"], 
    155155    ] 
    156156 
  • sasmodels/models/core_shell_cylinder.c

    r108e70e rd86f0fc  
    4848 
    4949 
    50 double Iqac(double qab, double qc, 
     50static double 
     51Iqac(double qab, double qc, 
    5152    double core_sld, 
    5253    double shell_sld, 
  • sasmodels/models/hollow_rectangular_prism.c

    r108e70e rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio, double thickness); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    63{ 
    74    double length_b = length_a * b2a_ratio; 
     
    1613} 
    1714 
    18 double Iq(double q, 
     15static double 
     16Iq(double q, 
    1917    double sld, 
    2018    double solvent_sld, 
     
    8583} 
    8684 
    87 double Iqabc(double qa, double qb, double qc, 
     85static double 
     86Iqabc(double qa, double qb, double qc, 
    8887    double sld, 
    8988    double solvent_sld, 
  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r74768cb rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    63{ 
    74    double length_b = length_a * b2a_ratio; 
     
    118} 
    129 
    13 double Iq(double q, 
     10static double 
     11Iq(double q, 
    1412    double sld, 
    1513    double solvent_sld, 
  • sasmodels/models/polymer_excl_volume.py

    r2d81cfe raa90015  
    1717$a$ is the statistical segment length of the polymer chain, 
    1818and $n$ is the degree of polymerization. 
    19 This integral was later put into an almost analytical form as follows 
     19 
     20This integral was put into an almost analytical form as follows 
    2021(Hammouda, 1993) 
     22 
     23.. math:: 
     24 
     25    P(Q)=\frac{1}{\nu U^{1/2\nu}} 
     26    \left\{ 
     27        \gamma\left(\frac{1}{2\nu},U\right) - 
     28        \frac{1}{U^{1/2\nu}}\gamma\left(\frac{1}{\nu},U\right) 
     29    \right\} 
     30 
     31and later recast as (for example, Hore, 2013; Hammouda & Kim, 2017) 
    2132 
    2233.. math:: 
     
    2940.. math:: 
    3041 
    31     \gamma(x,U)=\int_0^{U}dt\ exp(-t)t^{x-1} 
     42    \gamma(x,U)=\int_0^{U}dt\ \exp(-t)t^{x-1} 
    3243 
    3344and the variable $U$ is given in terms of the scattering vector $Q$ as 
     
    3748    U=\frac{Q^2a^2n^{2\nu}}{6} = \frac{Q^2R_{g}^2(2\nu+1)(2\nu+2)}{6} 
    3849 
     50The two analytic forms are equivalent. In the 1993 paper 
     51 
     52.. math:: 
     53 
     54    \frac{1}{\nu U^{1/2\nu}} 
     55 
     56has been factored out. 
     57 
     58**SasView implements the 1993 expression**. 
     59 
    3960The square of the radius-of-gyration is defined as 
    4061 
     
    4364    R_{g}^2 = \frac{a^2n^{2\nu}}{(2\nu+1)(2\nu+2)} 
    4465 
    45 Note that this model applies only in the mass fractal range (ie, $5/3<=m<=3$ ) 
    46 and **does not apply** to surface fractals ( $3<m<=4$ ). 
    47 It also does not reproduce the rigid rod limit (m=1) because it assumes chain 
    48 flexibility from the outset. It may cover a portion of the semi-flexible chain 
    49 range ( $1<m<5/3$ ). 
     66.. note:: 
     67    This model applies only in the mass fractal range (ie, $5/3<=m<=3$ ) 
     68    and **does not apply** to surface fractals ( $3<m<=4$ ). 
     69    It also does not reproduce the rigid rod limit (m=1) because it assumes chain 
     70    flexibility from the outset. It may cover a portion of the semi-flexible chain 
     71    range ( $1<m<5/3$ ). 
    5072 
    5173A low-Q expansion yields the Guinier form and a high-Q expansion yields the 
     
    7395.. math:: 
    7496 
    75     P(Q) = \frac{2}{Q^4R_{g}^4} \left[exp(-Q^2R_{g}^2) - 1 + Q^2R_{g}^2 \right] 
     97    P(Q) = \frac{2}{Q^4R_{g}^4} \left[\exp(-Q^2R_{g}^2) - 1 + Q^2R_{g}^2 \right] 
    7698 
    7799For 2D data: The 2D scattering intensity is calculated in the same way as 1D, 
     
    89111 
    90112B Hammouda, *SANS from Homogeneous Polymer Mixtures - A Unified Overview, 
    91 Advances in Polym. Sci.* 106(1993) 87-133 
     113Advances in Polym. Sci.* 106 (1993) 87-133 
     114 
     115M Hore et al, *Co-Nonsolvency of Poly(n-isopropylacrylamide) in Deuterated 
     116Water/Ethanol Mixtures* 46 (2013) 7894-7901 
     117 
     118B Hammouda & M-H Kim, *The empirical core-chain model* 247 (2017) 434-440 
    92119""" 
    93120 
     
    124151    with errstate(divide='ignore', invalid='ignore'): 
    125152        upow = power(usub, -0.5*porod_exp) 
     153        # Note: scipy gammainc is "regularized", being gamma(s,x)/Gamma(s), 
     154        # so need to scale by Gamma(s) to recover gamma(s, x). 
    126155        result = (porod_exp*upow * 
    127156                  (gamma(0.5*porod_exp)*gammainc(0.5*porod_exp, usub) - 
  • sasmodels/models/rectangular_prism.c

    r108e70e rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    63{ 
    74    return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio); 
    85} 
    96 
    10 double Iq(double q, 
     7static double 
     8Iq(double q, 
    119    double sld, 
    1210    double solvent_sld, 
     
    7169 
    7270 
    73 double Iqabc(double qa, double qb, double qc, 
     71static double 
     72Iqabc(double qa, double qb, double qc, 
    7473    double sld, 
    7574    double solvent_sld, 
  • sasmodels/rst2html.py

    r2d81cfe r1fbadb2  
    5656    # others don't work properly with math_output! 
    5757    if math_output == "mathjax": 
    58         settings = {"math_output": math_output} 
     58        # TODO: this is copied from docs/conf.py; there should be only one 
     59        mathjax_path = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML" 
     60        settings = {"math_output": math_output + " " + mathjax_path} 
    5961    else: 
    6062        settings = {"math-output": math_output} 
  • sasmodels/sasview_model.py

    r2d81cfe rd533590  
    593593            # Check whether we have a list of ndarrays [qx,qy] 
    594594            qx, qy = qdist 
    595             if not self._model_info.parameters.has_2d: 
    596                 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    597             else: 
    598                 return self.calculate_Iq(qx, qy) 
     595            return self.calculate_Iq(qx, qy) 
    599596 
    600597        elif isinstance(qdist, np.ndarray): 
     
    665662 
    666663    def _calculate_Iq(self, qx, qy=None): 
    667         #core.HAVE_OPENCL = False 
    668664        if self._model is None: 
    669665            self._model = core.build_model(self._model_info) 
     
    678674        call_details, values, is_magnetic = make_kernel_args(calculator, pairs) 
    679675        #call_details.show() 
    680         #print("pairs", pairs) 
     676        #print("================ parameters ==================") 
     677        #for p, v in zip(parameters.call_parameters, pairs): print(p.name, v[0]) 
    681678        #for k, p in enumerate(self._model_info.parameters.call_parameters): 
    682679        #    print(k, p.name, *pairs[k]) 
     
    689686        self._intermediate_results = getattr(calculator, 'results', None) 
    690687        calculator.release() 
    691         self._model.release() 
     688        #self._model.release() 
    692689        return result 
    693690 
     
    722719 
    723720    def set_dispersion(self, parameter, dispersion): 
    724         # type: (str, weights.Dispersion) -> Dict[str, Any] 
     721        # type: (str, weights.Dispersion) -> None 
    725722        """ 
    726723        Set the dispersion object for a model parameter 
     
    745742            self.dispersion[parameter] = dispersion.get_pars() 
    746743        else: 
    747             raise ValueError("%r is not a dispersity or orientation parameter") 
     744            raise ValueError("%r is not a dispersity or orientation parameter" 
     745                             % parameter) 
    748746 
    749747    def _dispersion_mesh(self): 
     
    859857    # type: () -> None 
    860858    """ 
    861     Load and run cylinder model from sas.models.CylinderModel 
     859    Load and run cylinder model as sas-models-CylinderModel 
    862860    """ 
    863861    if not SUPPORT_OLD_STYLE_PLUGINS: 
     
    872870    CylinderModel().evalDistribution([0.1, 0.1]) 
    873871 
     872def magnetic_demo(): 
     873    Model = _make_standard_model('sphere') 
     874    model = Model() 
     875    model.setParam('M0:sld', 8) 
     876    q = np.linspace(-0.35, 0.35, 500) 
     877    qx, qy = np.meshgrid(q, q) 
     878    result = model.calculate_Iq(qx.flatten(), qy.flatten()) 
     879    result = result.reshape(qx.shape) 
     880 
     881    import pylab 
     882    pylab.imshow(np.log(result + 0.001)) 
     883    pylab.show() 
     884 
    874885if __name__ == "__main__": 
    875886    print("cylinder(0.1,0.1)=%g"%test_cylinder()) 
     887    #magnetic_demo() 
    876888    #test_product() 
    877889    #test_structure_factor() 
  • sasmodels/weights.py

    r2d81cfe r3d58247  
    9797        return x, px 
    9898 
     99class UniformDispersion(Dispersion): 
     100    r""" 
     101    Uniform dispersion, with width $\sigma$. 
     102 
     103    .. math:: 
     104 
     105        w = 1 
     106    """ 
     107    type = "uniform" 
     108    default = dict(npts=35, width=0, nsigmas=None) 
     109    def _weights(self, center, sigma, lb, ub): 
     110        x = np.linspace(center-sigma, center+sigma, self.npts) 
     111        x = x[(x >= lb) & (x <= ub)] 
     112        return x, np.ones_like(x) 
    99113 
    100114class RectangleDispersion(Dispersion): 
     
    107121    """ 
    108122    type = "rectangle" 
    109     default = dict(npts=35, width=0, nsigmas=1.70325) 
    110     def _weights(self, center, sigma, lb, ub): 
    111         x = self._linspace(center, sigma, lb, ub) 
    112         x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] 
    113         return x, np.ones_like(x) 
    114  
     123    default = dict(npts=35, width=0, nsigmas=1.73205) 
     124    def _weights(self, center, sigma, lb, ub): 
     125         x = self._linspace(center, sigma, lb, ub) 
     126         x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] 
     127         return x, np.ones_like(x) 
    115128 
    116129class LogNormalDispersion(Dispersion): 
     
    190203        return x, px 
    191204 
     205class BoltzmannDispersion(Dispersion): 
     206    r""" 
     207    Boltzmann dispersion, with $\sigma=k T/E$. 
     208 
     209    .. math:: 
     210 
     211        w = \exp\left( -|x - c|/\sigma\right) 
     212    """ 
     213    type = "boltzmann" 
     214    default = dict(npts=35, width=0, nsigmas=3) 
     215    def _weights(self, center, sigma, lb, ub): 
     216        x = self._linspace(center, sigma, lb, ub) 
     217        px = np.exp(-np.abs(x-center) / np.abs(sigma)) 
     218        return x, px 
    192219 
    193220# dispersion name -> disperser lookup table. 
     
    196223MODELS = OrderedDict((d.type, d) for d in ( 
    197224    RectangleDispersion, 
     225    UniformDispersion, 
    198226    ArrayDispersion, 
    199227    LogNormalDispersion, 
    200228    GaussianDispersion, 
    201229    SchulzDispersion, 
     230    BoltzmannDispersion 
    202231)) 
    203232 
  • setup.py

    r32e3c9b r1f991d6  
    1 try: 
    2     from setuptools import setup 
    3 except ImportError: 
    4     from distutils.core import setup 
     1import sys 
     2from setuptools import setup 
     3from setuptools.command.test import test as TestCommand 
     4 
     5class PyTest(TestCommand): 
     6    user_options = [('pytest-args=', 'a', "Arguments to pass to pytest")] 
     7 
     8    def initialize_options(self): 
     9        TestCommand.initialize_options(self) 
     10        self.pytest_args = '' 
     11 
     12    def run_tests(self): 
     13        import shlex 
     14        #import here, cause outside the eggs aren't loaded 
     15        import pytest 
     16        errno = pytest.main(shlex.split(self.pytest_args)) 
     17        sys.exit(errno) 
    518 
    619def find_version(package): 
     
    4861        'sasmodels': ['*.c', '*.cl'], 
    4962    }, 
    50     install_requires = [ 
     63    install_requires=[ 
    5164    ], 
    52     extras_require = { 
     65    extras_require={ 
    5366        'OpenCL': ["pyopencl"], 
    5467        'Bumps': ["bumps"], 
    5568        'TinyCC': ["tinycc"], 
    5669    }, 
     70    build_requires=['setuptools'], 
     71    test_requires=['pytest'], 
     72    cmdclass = {'test': PyTest}, 
    5773) 
  • .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_*/ 
Note: See TracChangeset for help on using the changeset viewer.