Changeset 1662ebe in sasmodels
- Timestamp:
- May 2, 2018 5:02:29 PM (7 years ago)
- 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. - Files:
-
- 9 added
- 40 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
.travis.yml
r2d09df1 r335271e 37 37 - conda update --yes conda 38 38 - conda info -a 39 - conda install --yes python=$PY numpy scipy cython mako cffi39 - conda install --yes python=$PY numpy scipy matplotlib docutils setuptools pytest 40 40 install: 41 41 - pip install bumps 42 - pip install unittest-xml-reporting 42 - pip install unittest-xml-reporting tinycc 43 43 script: 44 44 - 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 46 48 before_deploy: 47 49 - echo -e "Host danse.chem.utk.edu\n\tStrictHostKeyChecking no\n" >> ~/.ssh/config -
appveyor.yml
r1258e32 r335271e 4 4 # /E:ON and /V:ON options are not enabled in the batch script interpreter 5 5 # 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" 7 8 8 9 # Workaround for https://github.com/conda/conda-build/issues/636 … … 22 23 - TARGET_ARCH: "x86" 23 24 CONDA_PY: "35" 24 PY_CONDITION: "python >=3.5 "25 PY_CONDITION: "python >=3.5,<3.6" 25 26 CONDA_INSTALL_LOCN: "C:\\Miniconda35" 26 27 - TARGET_ARCH: "x64" … … 34 35 - TARGET_ARCH: "x64" 35 36 CONDA_PY: "35" 36 PY_CONDITION: "python >=3.5 "37 PY_CONDITION: "python >=3.5,<3.6" 37 38 CONDA_INSTALL_LOCN: "C:\\Miniconda35-x64" 38 39 … … 43 44 44 45 install: 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. 46 48 - cmd: set CONDA_NPY=19 47 49 48 50 # 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% 50 56 51 57 # Use the pre-installed Miniconda for the desired arch … … 56 62 # so that we can update it. Then we remove it so that 57 63 # 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%"60 64 - cmd: conda update --yes --quiet conda python 61 - cmd: set "PATH=%OLDPATH%"62 65 - cmd: call %CONDA_INSTALL_LOCN%\Scripts\activate.bat 63 64 - cmd: conda config --add channels conda-forge65 - cmd: conda config --set show_channel_urls true66 - cmd: conda update --yes --quiet conda67 - cmd: conda install --yes --quiet obvious-ci68 - cmd: conda install --yes --quiet numpy toolchain scipy cython cffi66 #- 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 69 72 #- 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 70 74 - cmd: pip install bumps unittest-xml-reporting tinycc 71 75 72 76 build_script: 73 77 # 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 75 82 76 83 test_script: 77 84 # 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 190 190 *sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for 191 191 equirectangular and *PROJECTION=2* for sinusoidal. The more complicated 192 Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh192 Guyou and Postel projections are not implemented. See jitter.draw_mesh 193 193 for details. 194 194 … … 274 274 code. 275 275 276 * explore/jitter.py* is for exploring different options for handling277 orientation and orientation dispersity. It uses * explore/guyou.py* to276 *sasmodels/jitter.py* is for exploring different options for handling 277 orientation and orientation dispersity. It uses *sasmodels/guyou.py* to 278 278 generate the Guyou projection. 279 279 -
doc/genapi.py
r706f466 r0d5a655 69 69 ('exception', 'Annotate exceptions'), 70 70 ('generate', 'Model parser'), 71 ('jitter', 'Orientation explorer'), 72 ('guyou', 'Guyou map projection'), 71 73 ('kernel', 'Evaluator type definitions'), 72 74 ('kernelcl', 'OpenCL model evaluator'), -
doc/guide/orientation/orientation.rst
r5fb0634 r0d5a655 52 52 yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the 53 53 $a$-axis. 54 55 You can explore the view and jitter angles interactively using 56 :func:`sasmodels.jitter.run`. Enter the following into the python 57 interpreter:: 58 59 from sasmodels import jitter 60 jitter.run() 54 61 55 62 More formally, starting with axes $a$-$b$-$c$ of the particle aligned -
doc/guide/pd/polydispersity.rst
reda8b30 r29afc50 23 23 average over the size distribution. 24 24 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. 25 Each distribution is characterized by a center value $\bar x$ or 26 $x_\text{med}$, a width parameter $\sigma$ (note this is *not necessarily* 27 the standard deviation, so read the description carefully), the number of 28 sigmas $N_\sigma$ to include from the tails of the distribution, and the 29 number of points used to compute the average. The center of the distribution 30 is set by the value of the model parameter. The meaning of a polydispersity 31 parameter *PD* (not to be confused with a molecular weight distributions 32 in polymer science) in a model depends on the type of parameter it is being 33 applied too. 34 35 The distribution width applied to *volume* (ie, shape-describing) parameters 36 is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$. 37 However, the distribution width applied to *orientation* (ie, angle-describing) 38 parameters is just $\sigma = \mathrm{PD}$. 39 40 $N_\sigma$ determines how far into the tails to evaluate the distribution, 41 with larger values of $N_\sigma$ required for heavier tailed distributions. 34 42 The scattering in general falls rapidly with $qr$ so the usual assumption 35 43 that $G(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)G(r - 3\sigma_r)$ … … 42 50 calculations are generally more robust with more data points or more angles. 43 51 44 The following five distribution functions are provided: 45 52 The following distribution functions are provided: 53 54 * *Uniform Distribution* 46 55 * *Rectangular Distribution* 47 56 * *Gaussian Distribution* 57 * *Boltzmann Distribution* 48 58 * *Lognormal Distribution* 49 59 * *Schulz Distribution* … … 52 62 These are all implemented as *number-average* distributions. 53 63 64 Additional distributions are under consideration. 65 66 Suggested Applications 67 ^^^^^^^^^^^^^^^^^^^^^^ 68 69 If applying polydispersion to parameters describing particle sizes, use 70 the Lognormal or Schulz distributions. 71 72 If applying polydispersion to parameters describing interfacial thicknesses 73 or angular orientations, use the Gaussian or Boltzmann distributions. 74 75 If applying polydispersion to parameters describing angles, use the Uniform 76 distribution. Beware of using distributions that are always positive (eg, the 77 Lognormal) because angles can be negative! 78 79 The array distribution allows a user-defined distribution to be applied. 80 81 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 82 83 Uniform Distribution 84 ^^^^^^^^^^^^^^^^^^^^ 85 86 The 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 96 where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 97 distribution, $\sigma$ is the half-width, and *Norm* is a normalization 98 factor which is determined during the numerical calculation. 99 100 The 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 108 The value $N_\sigma$ is ignored for this distribution. 109 54 110 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 55 111 … … 63 119 f(x) = \frac{1}{\text{Norm}} 64 120 \begin{cases} 65 1 & \text{for } |x - \bar x| \leq w \\66 0 & \text{for } |x - \bar x| > w121 1 & \text{for } |x - \bar x| \leq w \\ 122 0 & \text{for } |x - \bar x| > w 67 123 \end{cases} 68 124 69 where $\bar x$ is the mean of the distribution, $w$ is the half-width, and70 *Norm* is a normalization factor which is determined during the numerical 71 calculation.125 where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 126 distribution, $w$ is the half-width, and *Norm* is a normalization 127 factor which is determined during the numerical calculation. 72 128 73 129 Note that the standard deviation and the half width $w$ are different! … … 77 133 .. math:: \sigma = w / \sqrt{3} 78 134 79 whilst the polydispersity i s135 whilst the polydispersity in sasmodels is given by 80 136 81 137 .. math:: \text{PD} = \sigma / \bar x … … 84 140 85 141 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. 86 146 87 147 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ … … 95 155 96 156 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 159 where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 160 distribution and *Norm* is a normalization factor which is determined 161 during the numerical calculation. 162 163 The polydispersity in sasmodels is given by 103 164 104 165 .. math:: \text{PD} = \sigma / \bar x … … 107 168 108 169 Normal distribution. 170 171 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 172 173 Boltzmann Distribution 174 ^^^^^^^^^^^^^^^^^^^^^^ 175 176 The 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 183 where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 184 distribution and *Norm* is a normalization factor which is determined 185 during the numerical calculation. 186 187 The width is defined as 188 189 .. math:: \sigma=\frac{k T}{E} 190 191 which 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. 109 197 110 198 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ … … 113 201 ^^^^^^^^^^^^^^^^^^^^^^ 114 202 203 The Lognormal Distribution describes a function of $x$ where $\ln (x)$ has 204 a normal distribution. The result is a distribution that is skewed towards 205 larger values of $x$. 206 115 207 The Lognormal Distribution is defined as 116 208 117 209 .. math:: 118 210 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 215 where *Norm* is a normalization factor which will be determined during 216 the numerical calculation, $\mu=\ln(x_\text{med})$ and $x_\text{med}$ 217 is the *median* value of the *lognormal* distribution, but $\sigma$ is 218 a parameter describing the width of the underlying *normal* distribution. 219 220 $x_\text{med}$ will be the value given for the respective size parameter 221 in sasmodels, for example, *radius=60*. 222 223 The polydispersity in sasmodels is given by 224 225 .. math:: \text{PD} = \sigma = p / x_\text{med} 226 227 The mean value of the distribution is given by $\bar x = \exp(\mu+ \sigma^2/2)$ 228 and the peak value by $\max x = \exp(\mu - \sigma^2)$. 229 230 The variance (the square of the standard deviation) of the *lognormal* 231 distribution is given by 232 233 .. math:: 234 235 \nu = [\exp({\sigma}^2) - 1] \exp({2\mu + \sigma^2}) 236 237 Note that larger values of PD might need a larger number of points 238 and $N_\sigma$. 139 239 140 240 .. figure:: pd_lognormal.jpg 141 241 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 244 For further information on the Lognormal distribution see: 245 http://en.wikipedia.org/wiki/Log-normal_distribution and 246 http://mathworld.wolfram.com/LogNormalDistribution.html 147 247 148 248 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ … … 151 251 ^^^^^^^^^^^^^^^^^^^ 152 252 253 The Schulz (sometimes written Schultz) distribution is similar to the 254 Lognormal distribution, in that it is also skewed towards larger values of 255 $x$, but which has computational advantages over the Lognormal distribution. 256 153 257 The Schulz distribution is defined as 154 258 155 259 .. math:: 156 260 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 264 where $\bar x$ ($x_\text{mean}$ in the figure) is the mean of the 265 distribution, *Norm* is a normalization factor which is determined 266 during the numerical calculation, and $z$ is a measure of the width 267 of the distribution such that 163 268 164 269 .. math:: z = (1-p^2) / p^2 165 270 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. 271 where $p$ is the polydispersity in sasmodels given by 272 273 .. math:: PD = p = \sigma / \bar x 274 275 and $\sigma$ is the RMS deviation from $\bar x$. 276 277 Note that larger values of PD might need a larger number of points 278 and $N_\sigma$. For example, for PD=0.7 with radius=60 |Ang|, at least 279 Npts>=160 and Nsigmas>=15 are required. 172 280 173 281 .. figure:: pd_schulz.jpg … … 176 284 177 285 For further information on the Schulz distribution see: 178 M Kotlarchyk & S-H Chen, *J Chem Phys*, (1983), 79, 2461. 286 M Kotlarchyk & S-H Chen, *J Chem Phys*, (1983), 79, 2461 and 287 M Kotlarchyk, RB Stephens, and JS Huang, *J Phys Chem*, (1988), 92, 1533 179 288 180 289 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ … … 183 292 ^^^^^^^^^^^^^^^^^^ 184 293 185 This user-definable distribution should be given as a s asimple ASCII text294 This user-definable distribution should be given as a simple ASCII text 186 295 file where the array is defined by two columns of numbers: $x$ and $f(x)$. 187 296 The $f(x)$ will be normalized to 1 during the computation. … … 202 311 given for the model will have no affect, and will be ignored when computing 203 312 the average. This means that any parameter with an array distribution will 204 not be fit table.313 not be fitable. 205 314 206 315 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ … … 216 325 related) except when the DLS polydispersity parameter is <0.13. 217 326 327 .. math:: 328 329 p_{DLS} = \sqrt(\nu / \bar x^2) 330 331 where $\nu$ is the variance of the distribution and $\bar x$ is the mean 332 value of $x$. 333 218 334 For more information see: 219 335 S King, C Washington & R Heenan, *Phys Chem Chem Phys*, (2005), 7, 143 … … 225 341 | 2015-05-01 Steve King 226 342 | 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 1 FileFormatVersion 1.0 2 DataFileTitle SiO2 stober 3 Sample SiO2 100pcH2O 4 collimation unknown 5 ID CPD 6 Date di 27 jan 2015 18:17:50 7 Method Offspec 8 Thickness 1 9 Thickness_unit mm 10 Theta_zmax 0.0100 11 Theta_zmax_unit radians 12 Theta_ymax 0.0100 13 Theta_ymax_unit radians 14 Orientation Z 15 SpinEchoLength_unit A 16 Depolarisation_unit A-2 cm-1 17 Wavelength_unit A 10 18 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 19 BEGIN_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 28 12525. -0.016993983 2.5549565e-4 5.0039968 0.1501199 50 0.653424 0.00418035 29 14605.8 -0.017297381 2.6798795e-4 5.40369504 0.16211085 50 0.603456 0.00472218 30 16846.7 -0.017312523 2.8577242e-4 5.80344105 0.17410323 50 0.558174 0.00537231 31 19247.7 -0.016368225 3.0135741e-4 6.20322561 0.18609677 50 0.532672 0.00617699 32 21808.7 -0.015668849 3.2144946e-4 6.60302658 0.1980908 50 0.505018 0.00707792 33 24529.8 -0.015157278 3.5589713e-4 7.00285542 0.21008566 50 0.475536 0.00829962 34 27411. -0.014947440 3.9623352e-4 7.40270761 0.22208123 50 0.440819 0.00957178 35 30452.2 -0.015147872 4.5944674e-4 7.80256676 0.234077 50 0.397642 0.0111225 36 33653.5 -0.015032370 5.3964070e-4 8.20244402 0.24607332 50 0.363717 0.0132055 37 37014.8 -0.016170897 7.0802787e-4 8.60232527 0.25806976 50 0.302206 0.0158338 38 40536.2 -0.014766858 8.3444978e-4 9.00222106 0.27006663 50 0.302188 0.0204351 39 44217.7 -0.015097771 1.1351144e-3 9.40212955 0.28206389 50 0.263252 0.0264158 40 48059.2 -0.015289642 1.5374507e-3 9.80203897 0.29406117 50 0.230149 0.0339972 41 52060.8 -0.012624691 1.7166363e-3 10.20195903 0.30605877 50 0.268749 0.0480167 42 56222.5 -9.5864674e-3 1.8344138e-3 10.60188851 0.31805666 50 0.340439 0.0701945 43 60544.2 -0.011852755 3.4220757e-3 11.0018173 0.33005452 50 0.238197 0.0986631 44 65026. -0.017820748 0.013284073 11.40175425 0.34205263 50 0.0985987 0.170273 45 69667.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 1 FileFormatVersion 1.0 2 DataFileTitle se009051.dat 3 Sample "SiO2 pH=6 & NiO pH=10, nanosized, in D2O, Sine ++ only." 4 collimation "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. " 5 ID CPD 6 Date di 27 jan 2015 18:17:50 7 Method sine one element scan 8 Thickness 10 9 Thickness_unit mm 10 Theta_zmax 0.0168 11 Theta_zmax_unit radians 12 Theta_ymax 0.0168 13 Theta_ymax_unit radians 14 Orientation Z 15 SpinEchoLength_unit A 16 Depolarisation_unit A-2 cm-1 17 Wavelength_unit A 10 18 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 19 BEGIN_DATA 20 SpinEchoLength SpinEchoLength_error Wavelength Wavelength_Error Polarisation Polarisation_error Depolarisation Depolarisation_error 21 298.82023 50 2.095 0.1 0.997817556 0.00887342 -4.9779370e-5 2.0261512e-4 22 852.57531 50 2.095 0.1 1.0026184 0.009572307 5.9579929e-5 2.1752686e-4 23 1407.9863 50 2.095 0.1 0.996013172 0.012119508 -9.1017859e-5 2.7723742e-4 24 1965.6050 50 2.095 0.1 0.991752656 0.011651473 -1.8868750e-4 2.6767598e-4 25 2525.2813 50 2.095 0.1 0.995429571 0.012193294 -1.0437182e-4 2.7908883e-4 26 3087.2927 50 2.095 0.1 0.995119819 0.009621277 -1.1146275e-4 2.2028721e-4 27 3650.3712 50 2.095 0.1 0.997497767 0.012518842 -5.7082583e-5 2.8594610e-4 28 4216.8468 50 2.095 0.1 0.994733865 0.010585926 -1.2030120e-4 2.4246770e-4 29 4785.6485 50 2.095 0.1 0.992696275 0.010677724 -1.6701950e-4 2.4507231e-4 30 5357.1552 50 2.095 0.1 0.994534294 0.012909556 -1.2487278e-4 2.9574914e-4 31 5932.4408 50 2.095 0.1 0.986412268 0.014651894 -3.1170682e-4 3.3842875e-4 32 6509.8550 50 2.095 0.1 0.988788566 0.015736468 -2.5688520e-4 3.6260666e-4 33 7092.0396 50 2.095 0.1 0.989501168 0.013611615 -2.4047103e-4 3.1341898e-4 34 7678.0994 50 2.095 0.1 0.99203933 0.016501891 -1.8210252e-4 3.7899787e-4 35 8267.6293 50 2.095 0.1 0.994062723 0.010276207 -1.3567871e-4 2.3553259e-4 36 8863.0230 50 2.095 0.1 0.979647192 0.011089421 -4.6850452e-4 2.5791174e-4 37 9462.3637 50 2.095 0.1 0.97982075 0.013168527 -4.6446835e-4 3.0621222e-4 38 10067.841 50 2.095 0.1 0.977927712 0.014362819 -5.0853039e-4 3.3463000e-4 39 10678.912 50 2.095 0.1 0.978351462 0.013239476 -4.9865985e-4 3.0832436e-4 40 11295.987 50 2.095 0.1 0.981246471 0.013896898 -4.3133968e-4 3.2267975e-4 41 11920.227 50 2.095 0.1 0.978275164 0.013747481 -5.0043677e-4 3.2017989e-4 42 12552.713 50 2.095 0.1 0.976818842 0.013487944 -5.3437989e-4 3.1460359e-4 43 13193.071 50 2.095 0.1 0.981546068 0.011972943 -4.2438423e-4 2.7792152e-4 44 13841.978 50 2.095 0.1 0.96445854 0.013357572 -8.2452102e-4 3.1555561e-4 45 14500.513 50 2.095 0.1 0.972756158 0.014756189 -6.2933879e-4 3.4562263e-4 46 15169.786 50 2.095 0.1 0.971010184 0.014210504 -6.7027011e-4 3.3343996e-4 47 15850.446 50 2.095 0.1 0.975009829 0.013648839 -5.7661387e-4 3.1894710e-4 48 16543.562 50 2.095 0.1 0.969575209 0.014454336 -7.0396574e-4 3.3966327e-4 49 17249.754 50 2.095 0.1 0.959907267 0.012806944 -9.3229353e-4 3.0398222e-4 50 17969.627 50 2.095 0.1 0.96219203 0.011845496 -8.7812744e-4 2.8049391e-4 51 18706.130 50 2.095 0.1 0.960157966 0.011721155 -9.2634378e-4 2.7813758e-4 52 19459.212 50 2.095 0.1 0.95089937 0.009967287 -1.1471121e-3 2.3882201e-4 53 20230.984 50 2.095 0.1 0.955584642 0.01104368 -1.0351259e-3 2.6331561e-4 54 21021.899 50 2.095 0.1 0.950786938 0.011651582 -1.1498062e-3 2.7921171e-4 55 21835.345 50 2.095 0.1 0.951403607 0.008784149 -1.1350335e-3 2.1036178e-4 56 22671.730 50 2.095 0.1 0.94561823 0.009392383 -1.2740040e-3 2.2630383e-4 57 23533.930 50 2.095 0.1 0.945882787 0.009468264 -1.2676305e-3 2.2806833e-4 58 24423.504 50 2.095 0.1 0.939583158 0.009164239 -1.4198814e-3 2.2222511e-4 59 25343.343 50 2.095 0.1 0.93731257 0.008621568 -1.4750079e-3 2.0957224e-4 60 26294.734 50 2.095 0.1 0.936639557 0.010408821 -1.4913733e-3 2.5319842e-4 61 27282.004 50 2.095 0.1 0.932091884 0.007112674 -1.6022666e-3 1.7386258e-4 62 28307.499 50 2.095 0.1 0.933852763 0.009200431 -1.5592642e-3 2.2447176e-4 63 29373.630 50 2.095 0.1 0.930283437 0.01022719 -1.6465153e-3 2.5047996e-4 64 30485.110 50 2.095 0.1 0.930121334 0.007055747 -1.6504858e-3 1.7283645e-4 65 31645.899 50 2.095 0.1 0.925944326 0.0078917 -1.7530356e-3 1.9418588e-4 66 32858.438 50 2.095 0.1 0.92351759 0.005522582 -1.8128270e-3 1.3624763e-4 67 34128.058 50 2.095 0.1 0.918228311 0.006917328 -1.9436940e-3 1.7164045e-4 68 35460.168 50 2.095 0.1 0.908621563 0.006852313 -2.1833230e-3 1.7182490e-4 69 36858.762 50 2.095 0.1 0.919087792 0.006995961 -1.9223776e-3 1.7342924e-4 70 38331.748 50 2.095 0.1 0.910132214 0.006107336 -2.1454742e-3 1.5289007e-4 71 39882.345 50 2.095 0.1 0.904491462 0.007078101 -2.2871233e-3 1.7829708e-4 72 41519.404 50 2.095 0.1 0.90056919 0.005934428 -2.3861400e-3 1.5013907e-4 73 43249.023 50 2.095 0.1 0.903155701 0.006249069 -2.3207959e-3 1.5764661e-4 74 45079.637 50 2.095 0.1 0.894919827 0.006196232 -2.5295172e-3 1.5775222e-4 75 47019.670 50 2.095 0.1 0.896358908 0.007324999 -2.4929085e-3 1.8619052e-4 76 49077.851 50 2.095 0.1 0.88807923 0.004530465 -2.7043436e-3 1.1623128e-4 77 51264.711 50 2.095 0.1 0.886649039 0.008924412 -2.7410654e-3 2.2932944e-4 78 53590.598 50 2.095 0.1 0.882051947 0.006055081 -2.8595036e-3 1.5640756e-4 79 56067.403 50 2.095 0.1 0.882573371 0.005727419 -2.8460388e-3 1.4785638e-4 80 58709.036 50 2.095 0.1 0.876289825 0.006062306 -3.0088321e-3 1.5762389e-4 81 61527.362 50 2.095 0.1 0.877318199 0.005307783 -2.9821094e-3 1.3784403e-4 82 64538.932 50 2.095 0.1 0.873610284 0.006918104 -3.0786086e-3 1.8042690e-4 83 67758.354 50 2.095 0.1 0.865954455 0.005625863 -3.2791557e-3 1.4802192e-4 84 71205.306 50 2.095 0.1 0.868037308 0.008662141 -3.2244196e-3 2.2736248e-4 85 74897.359 50 2.095 0.1 0.86676404 0.005761209 -3.2578647e-3 1.5144143e-4 86 78855.560 50 2.095 0.1 0.85655072 0.005731411 -3.5279304e-3 1.5245456e-4 87 83102.902 50 2.095 0.1 0.857936274 0.006954572 -3.4911046e-3 1.8469168e-4 88 87663.154 50 2.095 0.1 0.85520525 0.006559622 -3.5637478e-3 1.7475934e-4 89 92562.717 50 2.095 0.1 0.849739829 0.008907377 -3.7098230e-3 2.3883381e-4 90 97830.316 50 2.095 0.1 0.841487161 0.006547659 -3.9321836e-3 1.7728439e-4 91 103497.14 50 2.095 0.1 0.84801051 0.007574866 -3.7562386e-3 2.0351933e-4 92 109596.79 50 2.095 0.1 0.843342861 0.006411993 -3.8819940e-3 1.7322909e-4 93 116165.97 50 2.095 0.1 0.837488454 0.008141353 -4.0407107e-3 2.2148775e-4 94 123244.47 50 2.095 0.1 0.839214922 0.009987535 -3.9937900e-3 2.7115465e-4 95 130874.72 50 2.095 0.1 0.834336972 0.011467406 -4.1266093e-3 3.1315233e-4 96 139102.96 50 2.095 0.1 0.828775634 0.012483365 -4.2789870e-3 3.4318369e-4 97 147980.51 50 2.095 0.1 0.823088812 0.016035848 -4.4358638e-3 4.4389186e-4 98 157560.98 50 2.095 0.1 0.824329178 0.014838138 -4.4015548e-3 4.1011975e-4 99 167904.68 50 2.095 0.1 0.82023172 0.015286187 -4.5150892e-3 4.2461424e-4 100 179076.29 50 2.095 0.1 0.827559752 0.012612741 -4.3124376e-3 3.4724985e-4 -
example/fit.py
rf4b36fa r1ee5ac6 179 179 180 180 else: 181 print "No parameters for %s"%name181 print("No parameters for %s"%name) 182 182 sys.exit(1) 183 183 … … 192 192 else: 193 193 problem = FitProblem(M) 194 195 if __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 1 FileFormatVersion 1.0 2 DataFileTitle Polystyrene of Markus Strobl, Full Sine, ++ only 3 Sample Polystyrene 2 um in 53% H2O, 47% D2O 4 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 5 Operator CPD 6 Date ma 7 jul 2014 18:54:43 7 ScanType sine one element scan 8 Thickness 2 9 Thickness_unit mm 10 Theta_zmax 0.0168 11 Theta_zmax_unit radians 12 Theta_ymax 0.0168 13 Theta_ymax_unit radians 14 Orientation Z 15 SpinEchoLength_unit A 16 Depolarisation_unit A-2 cm-1 17 Wavelength_unit A 18 19 BEGIN_DATA 20 SpinEchoLength Depolarisation Depolarisation_error Wavelength Wavelength_error SpinEchoLength_error Polarisation Polarisation_error 21 497.78 -2.4509553e-4 4.9935908e-4 2.11 0.1055 24.889 0.99782 0.0044367 22 630.41 2.9161810e-4 5.3612769e-4 2.11 0.1055 31.52 1.0026 0.0047862 23 764.87 -4.4899949e-4 6.8328154e-4 2.11 0.1055 38.244 0.99601 0.0060598 24 898.47 -9.3037214e-4 6.5970686e-4 2.11 0.1055 44.924 0.99175 0.0058257 25 1034.1 -5.1441728e-4 6.8783151e-4 2.11 0.1055 51.705 0.99543 0.0060966 26 1169.5 -5.4939760e-4 5.4291131e-4 2.11 0.1055 58.475 0.99512 0.0048106 27 1306.1 -2.8111792e-4 7.0473347e-4 2.11 0.1055 65.303 0.9975 0.0062594 28 1443.7 -5.9342057e-4 5.9758787e-4 2.11 0.1055 72.184 0.99473 0.005293 29 1582. -8.2284488e-4 6.0400267e-4 2.11 0.1055 79.102 0.9927 0.0053389 30 1721.2 -6.1600315e-4 7.2890343e-4 2.11 0.1055 86.062 0.99453 0.0064548 31 1861.7 -1.5367118e-3 8.3408174e-4 2.11 0.1055 93.087 0.98641 0.0073259 32 2002.8 -1.2660661e-3 8.9366844e-4 2.11 0.1055 100.14 0.98879 0.0078682 33 2144.6 -1.1854534e-3 7.7244662e-4 2.11 0.1055 107.23 0.9895 0.0068058 34 2287.3 -8.9753711e-4 9.3406529e-4 2.11 0.1055 114.36 0.99204 0.0082509 35 2431.2 -6.6909009e-4 5.8049041e-4 2.11 0.1055 121.56 0.99406 0.0051381 36 2575.5 -2.3090130e-3 6.3564144e-4 2.11 0.1055 128.78 0.97965 0.0055447 37 2722.2 -2.2895260e-3 7.5468967e-4 2.11 0.1055 136.11 0.97982 0.0065843 38 2869. -2.5063662e-3 8.2471984e-4 2.11 0.1055 143.45 0.97793 0.0071814 39 3017.3 -2.4581433e-3 7.5988724e-4 2.11 0.1055 150.87 0.97835 0.0066197 40 3167.5 -2.1257395e-3 7.9526201e-4 2.11 0.1055 158.38 0.98125 0.0069484 41 3318.2 -2.4661790e-3 7.8910082e-4 2.11 0.1055 165.91 0.97828 0.0068737 42 3471.6 -2.6339122e-3 7.7536843e-4 2.11 0.1055 173.58 0.97682 0.006744 43 3624.5 -2.0914090e-3 6.8496070e-4 2.11 0.1055 181.22 0.98155 0.0059865 44 3780.9 -4.0640282e-3 7.7771292e-4 2.11 0.1055 189.04 0.96446 0.0066788 45 3937.4 -3.1016697e-3 8.5181234e-4 2.11 0.1055 196.87 0.97276 0.0073781 46 4096.1 -3.3038917e-3 8.2179560e-4 2.11 0.1055 204.81 0.97101 0.0071053 47 4255.5 -2.8422039e-3 7.8606869e-4 2.11 0.1055 212.78 0.97501 0.0068244 48 4419.1 -3.4694067e-3 8.3712733e-4 2.11 0.1055 220.96 0.96958 0.0072272 49 4581.2 -4.5951067e-3 7.4919003e-4 2.11 0.1055 229.06 0.95991 0.0064035 50 4747.7 -4.3286699e-3 6.9129591e-4 2.11 0.1055 237.39 0.96219 0.0059227 51 4915.2 -4.5658612e-3 6.8549385e-4 2.11 0.1055 245.76 0.96016 0.0058606 52 5085.1 -5.6542277e-3 5.8859074e-4 2.11 0.1055 254.26 0.9509 0.0049836 53 5256.8 -5.1028496e-3 6.4896117e-4 2.11 0.1055 262.84 0.95558 0.0055218 54 5430.8 -5.6672201e-3 6.8813882e-4 2.11 0.1055 271.54 0.95079 0.0058258 55 5607.4 -5.5951905e-3 5.1845870e-4 2.11 0.1055 280.37 0.9514 0.0043921 56 5786.9 -6.2795627e-3 5.5774416e-4 2.11 0.1055 289.35 0.94562 0.0046962 57 5967.5 -6.2486880e-3 5.6209080e-4 2.11 0.1055 298.38 0.94588 0.0047341 58 6151.7 -6.9992040e-3 5.4769136e-4 2.11 0.1055 307.58 0.93958 0.0045821 59 6337.7 -7.2708619e-3 5.1651117e-4 2.11 0.1055 316.88 0.93731 0.0043108 60 6528.2 -7.3511686e-3 6.2402654e-4 2.11 0.1055 326.41 0.93664 0.0052044 61 6720.4 -7.8980596e-3 4.2849488e-4 2.11 0.1055 336.02 0.93209 0.0035563 62 6915.7 -7.6861990e-3 5.5322868e-4 2.11 0.1055 345.78 0.93385 0.0046002 63 7114.8 -8.1163566e-3 6.1733111e-4 2.11 0.1055 355.74 0.93028 0.0051136 64 7316.7 -8.1356741e-3 4.2597330e-4 2.11 0.1055 365.83 0.93012 0.0035279 65 7521.7 -8.6415221e-3 4.7858305e-4 2.11 0.1055 376.08 0.92594 0.0039458 66 7731. -8.9354263e-3 3.3579357e-4 2.11 0.1055 386.55 0.92352 0.0027613 67 7944.2 -9.5805772e-3 4.2302546e-4 2.11 0.1055 397.21 0.91823 0.0034587 68 8161.4 -0.010762148 4.2348254e-4 2.11 0.1055 408.07 0.90862 0.0034262 69 8381.9 -9.4754418e-3 4.2743183e-4 2.11 0.1055 419.1 0.91909 0.003498 70 8606.9 -0.010575665 3.7681487e-4 2.11 0.1055 430.34 0.91013 0.0030537 71 8837.4 -0.011273784 4.3943451e-4 2.11 0.1055 441.87 0.90449 0.0035391 72 9070.4 -0.011761571 3.7002787e-4 2.11 0.1055 453.52 0.90057 0.0029672 73 9309.5 -0.011439046 3.8852675e-4 2.11 0.1055 465.47 0.90316 0.0031245 74 9553.8 -0.012468380 3.8879110e-4 2.11 0.1055 477.69 0.89492 0.0030981 75 9802.5 -0.012287815 4.5888119e-4 2.11 0.1055 490.13 0.89636 0.0036625 76 10057. -0.013330052 2.8645708e-4 2.11 0.1055 502.84 0.88808 0.0022652 77 10317. -0.013511036 5.6519968e-4 2.11 0.1055 515.86 0.88665 0.0044622 78 10583. -0.014095206 3.8547484e-4 2.11 0.1055 529.14 0.88205 0.0030275 79 10855. -0.014029017 3.6440427e-4 2.11 0.1055 542.76 0.88257 0.0028637 80 11134. -0.014831000 3.8848283e-4 2.11 0.1055 556.69 0.87629 0.0030312 81 11419. -0.014699072 3.3972822e-4 2.11 0.1055 570.93 0.87732 0.0026539 82 11710 -0.015174999 4.4468309e-4 2.11 0.1055 585.5 0.87361 0.0034591 83 12010 -0.016164070 3.6480986e-4 2.11 0.1055 600.52 0.86595 0.0028129 84 12316 -0.015893340 5.6035541e-4 2.11 0.1055 615.82 0.86804 0.0043311 85 12630 -0.016059068 3.7324087e-4 2.11 0.1055 631.51 0.86676 0.0028806 86 12952. -0.017389837 3.7573625e-4 2.11 0.1055 647.62 0.85655 0.0028657 87 13284. -0.017207735 4.5518751e-4 2.11 0.1055 664.19 0.85794 0.0034773 88 13623. -0.017565669 4.3070477e-4 2.11 0.1055 681.13 0.85521 0.0032798 89 13971. -0.018286298 5.8862675e-4 2.11 0.1055 698.54 0.84974 0.0044537 90 14329. -0.019381994 4.3692639e-4 2.11 0.1055 716.44 0.84149 0.0032738 91 14697. -0.018515178 5.0158587e-4 2.11 0.1055 734.87 0.84801 0.0037874 92 15075. -0.019135361 4.2693908e-4 2.11 0.1055 753.76 0.84334 0.003206 93 15464. -0.019917113 5.4587670e-4 2.11 0.1055 773.18 0.83749 0.0040707 94 15864. -0.019686699 6.6829096e-4 2.11 0.1055 793.22 0.83921 0.0049938 95 16276. -0.020340321 7.7178617e-4 2.11 0.1055 813.81 0.83434 0.0057337 96 16699. -0.021091231 8.4580203e-4 2.11 0.1055 834.97 0.82878 0.0062417 97 17136. -0.021864932 1.0940027e-3 2.11 0.1055 856.78 0.82309 0.0080179 98 17585. -0.021695868 1.0107767e-3 2.11 0.1055 879.27 0.82433 0.0074191 99 18048. -0.022255844 1.0464994e-3 2.11 0.1055 902.39 0.82023 0.0076431 100 18525. -0.021256673 8.5582923e-4 2.11 0.1055 926.27 0.82756 0.0063064 101 19017. -0.021490334 7.1998911e-4 2.11 0.1055 950.87 0.82584 0.0052944 102 19524. -0.021198334 6.8716706e-4 2.11 0.1055 976.22 0.82799 0.0050662 103 20048. -0.021815823 5.8818928e-4 2.11 0.1055 1002.4 0.82345 0.0043127 104 20588. -0.021882671 5.1857307e-4 2.11 0.1055 1029.4 0.82296 0.0038 105 21145. -0.022305147 4.6672141e-4 2.11 0.1055 1057.2 0.81987 0.0034072 106 21720 -0.022258583 5.0322361e-4 2.11 0.1055 1086. 0.82021 0.0036752 107 22314. -0.021244460 3.9148871e-4 2.11 0.1055 1115.7 0.82765 0.0028851 108 22928. -0.021381594 5.2906244e-4 2.11 0.1055 1146.4 0.82664 0.0038942 109 23562. -0.021329979 6.4328235e-4 2.11 0.1055 1178.1 0.82702 0.0047371 110 24218. -0.022846153 6.0248825e-4 2.11 0.1055 1210.9 0.81593 0.0043772 111 24896. -0.021591012 3.8146933e-4 2.11 0.1055 1244.8 0.8251 0.0028026 112 25595. -0.021034332 3.3282938e-4 2.11 0.1055 1279.8 0.8292 0.0024574 113 26319. -0.021433232 4.9200888e-4 2.11 0.1055 1315.9 0.82626 0.0036198 114 27068. -0.022177827 4.4213864e-4 2.11 0.1055 1353.4 0.8208 0.0032314 115 27843. -0.022343508 5.8553317e-4 2.11 0.1055 1392.2 0.81959 0.0042731 116 28645. -0.021396539 3.6673246e-4 2.11 0.1055 1432.3 0.82653 0.002699 117 29475. -0.021739473 5.0054859e-4 2.11 0.1055 1473.8 0.82401 0.0036726 118 30335. -0.021794003 6.5757715e-4 2.11 0.1055 1516.7 0.82361 0.0048224 119 31224. -0.021798094 5.6210549e-4 2.11 0.1055 1561.2 0.82358 0.0041221 120 32145. -0.022031519 3.9364070e-4 2.11 0.1055 1607.3 0.82187 0.0028807 121 33101. -0.021408769 4.7779613e-4 2.11 0.1055 1655. 0.82644 0.003516 122 34090 -0.021802185 2.8863827e-4 2.11 0.1055 1704.5 0.82355 0.0021166 123 35114. -0.021586929 4.6155484e-4 2.11 0.1055 1755.7 0.82513 0.0033911 124 36176. -0.021194265 2.0808762e-4 2.11 0.1055 1808.8 0.82802 0.0015342 125 37276. -0.021382952 3.9701221e-4 2.11 0.1055 1863.8 0.82663 0.0029222 126 38417. -0.022251737 2.8416874e-4 2.11 0.1055 1920.8 0.82026 0.0020755 127 39600 -0.020819189 3.5679523e-4 2.11 0.1055 1980.0 0.83079 0.0026394 128 40827. -0.021380235 3.7314604e-4 2.11 0.1055 2041.3 0.82665 0.0027466 129 42099. -0.021232248 3.4189634e-4 2.11 0.1055 2105. 0.82774 0.0025199 130 43419. -0.020266312 4.1187633e-4 2.11 0.1055 2170.9 0.83489 0.0030619 131 44787. -0.022305147 2.8749557e-4 2.11 0.1055 2239.4 0.81987 0.0020988 132 46208. -0.021563793 3.2521680e-4 2.11 0.1055 2310.4 0.8253 0.0023899 133 47682. -0.021396539 3.1049291e-4 2.11 0.1055 2384.1 0.82653 0.0022851 134 49211. -0.021683607 4.6084892e-4 2.11 0.1055 2460.6 0.82442 0.003383 135 50798. -0.021160361 2.1666201e-4 2.11 0.1055 2539.9 0.82827 0.0015979 136 52447. -0.021612792 4.2378726e-4 2.11 0.1055 2622.3 0.82494 0.0031129 137 54157. -0.022036985 4.1199886e-4 2.11 0.1055 2707.9 0.82183 0.0030149 138 55933. -0.020632795 6.3397053e-4 2.11 0.1055 2796.7 0.83217 0.0046976 139 57778. -0.021976873 7.6130313e-4 2.11 0.1055 2888.9 0.82227 0.005574 140 59693. -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 1 FileFormatVersion 1.0 2 DataFileTitle Polystyrene of Markus Strobl, Full Sine, ++ only 3 Sample Polystyrene 2 um in 53% H2O, 47% D2O 4 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 " 5 Operator CPD 6 Date do 10 jul 2014 16:37:30 7 ScanType sine one element scan 8 Thickness 2 9 Thickness_unit mm 10 Theta_zmax 0.0168 11 Theta_zmax_unit radians 12 Theta_ymax 0.0168 13 Theta_ymax_unit radians 14 Orientation Z 15 SpinEchoLength_unit A 16 Depolarisation_unit A-2 cm-1 17 Wavelength_unit A 18 19 BEGIN_DATA 20 SpinEchoLength Depolarisation Depolarisation_error Polarisation Polarisation_error SpinEchoLength_error Wavelength Wavelength_error 21 391.56 0.0041929 0.0036894 1.0037 0.0032974 19.578 2.11 0.1055 22 1564 -0.0046571 0.0038185 0.99586 0.003386 78.2 2.11 0.1055 23 2735.6 -0.017007 0.0038132 0.98497 0.0033444 136.78 2.11 0.1055 24 3907.9 -0.033462 0.0035068 0.97064 0.0030309 195.39 2.11 0.1055 25 5080.2 -0.047483 0.0038208 0.9586 0.0032613 254.01 2.11 0.1055 26 6251.8 -0.070375 0.00376 0.93926 0.0031446 312.59 2.11 0.1055 27 7423.2 -0.092217 0.0037927 0.92117 0.0031108 371.16 2.11 0.1055 28 8595.5 -0.10238 0.004006 0.91287 0.0032562 429.77 2.11 0.1055 29 9767.7 -0.12672 0.0038534 0.8933 0.0030651 488.39 2.11 0.1055 30 10940 -0.1374 0.004243 0.88484 0.003343 546.98 2.11 0.1055 31 12112 -0.16072 0.0045837 0.86666 0.0035372 605.58 2.11 0.1055 32 13284 -0.16623 0.0045613 0.86242 0.0035027 664.2 2.11 0.1055 33 14456 -0.18468 0.0044918 0.84837 0.0033931 722.79 2.11 0.1055 34 15628 -0.19143 0.0048967 0.84328 0.0036768 781.38 2.11 0.1055 35 16800 -0.20029 0.0045421 0.83666 0.0033837 840.02 2.11 0.1055 36 17971 -0.19798 0.0046642 0.83838 0.0034819 898.56 2.11 0.1055 37 19143 -0.21442 0.0047052 0.82619 0.0034614 957.17 2.11 0.1055 38 20316 -0.20885 0.0044931 0.8303 0.0033218 1015.8 2.11 0.1055 39 21488 -0.21393 0.0049186 0.82655 0.00362 1074.4 2.11 0.1055 40 22660 -0.20685 0.004423 0.83179 0.0032758 1133 2.11 0.1055 41 23832 -0.20802 0.0046979 0.83092 0.0034758 1191.6 2.11 0.1055 42 25003 -0.19848 0.0045953 0.838 0.0034289 1250.2 2.11 0.1055 43 26175 -0.21117 0.0044567 0.82859 0.0032881 1308.8 2.11 0.1055 44 27347 -0.21283 0.004137 0.82736 0.0030477 1367.4 2.11 0.1055 45 28520 -0.2042 0.0044587 0.83375 0.0033101 1426 2.11 0.1055 46 29692 -0.2112 0.0042852 0.82857 0.0031615 1484.6 2.11 0.1055 47 30864 -0.20319 0.0043483 0.8345 0.003231 1543.2 2.11 0.1055 48 32036 -0.20752 0.0044297 0.83129 0.0032788 1601.8 2.11 0.1055 49 33207 -0.20654 0.0043188 0.83201 0.0031995 1660.4 2.11 0.1055 50 34380 -0.20126 0.0046375 0.83593 0.0034518 1719 2.11 0.1055 51 35551 -0.20924 0.0042871 0.83001 0.0031684 1777.6 2.11 0.1055 52 36724 -0.21323 0.0045471 0.82707 0.0033487 1836.2 2.11 0.1055 53 37895 -0.21324 0.0045354 0.82706 0.00334 1894.7 2.11 0.1055 54 39067 -0.19905 0.0044141 0.83758 0.003292 1953.4 2.11 0.1055 55 40239 -0.1991 0.0047441 0.83754 0.003538 2012 2.11 0.1055 56 41411 -0.20359 0.0050136 0.8342 0.003724 2070.5 2.11 0.1055 57 42583 -0.21032 0.0049474 0.82922 0.0036529 2129.1 2.11 0.1055 58 43755 -0.20689 0.0048203 0.83176 0.00357 2187.8 2.11 0.1055 59 44927 -0.21075 0.0052337 0.8289 0.0038628 2246.4 2.11 0.1055 60 46099 -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 1 FileFormatVersion 1.0 2 DataFileTitle se009051.dat 3 Sample "SiO2 pH=6 & NiO pH=10, nanosized, in D2O, Sine ++ only." 4 collimation "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. " 5 ID CPD 6 Date di 27 jan 2015 18:17:50 7 Method sine one element scan 8 Thickness 10 9 Thickness_unit mm 10 Theta_zmax 0.0168 11 Theta_zmax_unit radians 12 Theta_ymax 0.0168 13 Theta_ymax_unit radians 14 Orientation Z 15 SpinEchoLength_unit A 16 Depolarisation_unit A-2 cm-1 17 Wavelength_unit A 10 18 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 19 BEGIN_DATA 20 SpinEchoLength SpinEchoLength_error Wavelength Wavelength_Error Polarisation Polarisation_error Depolarisation Depolarisation_error 21 298.82023 50 2.095 0.1 0.997817556 0.00887342 -4.9779370e-5 2.0261512e-4 22 852.57531 50 2.095 0.1 1.0026184 0.009572307 5.9579929e-5 2.1752686e-4 23 1407.9863 50 2.095 0.1 0.996013172 0.012119508 -9.1017859e-5 2.7723742e-4 24 1965.6050 50 2.095 0.1 0.991752656 0.011651473 -1.8868750e-4 2.6767598e-4 25 2525.2813 50 2.095 0.1 0.995429571 0.012193294 -1.0437182e-4 2.7908883e-4 26 3087.2927 50 2.095 0.1 0.995119819 0.009621277 -1.1146275e-4 2.2028721e-4 27 3650.3712 50 2.095 0.1 0.997497767 0.012518842 -5.7082583e-5 2.8594610e-4 28 4216.8468 50 2.095 0.1 0.994733865 0.010585926 -1.2030120e-4 2.4246770e-4 29 4785.6485 50 2.095 0.1 0.992696275 0.010677724 -1.6701950e-4 2.4507231e-4 30 5357.1552 50 2.095 0.1 0.994534294 0.012909556 -1.2487278e-4 2.9574914e-4 31 5932.4408 50 2.095 0.1 0.986412268 0.014651894 -3.1170682e-4 3.3842875e-4 32 6509.8550 50 2.095 0.1 0.988788566 0.015736468 -2.5688520e-4 3.6260666e-4 33 7092.0396 50 2.095 0.1 0.989501168 0.013611615 -2.4047103e-4 3.1341898e-4 34 7678.0994 50 2.095 0.1 0.99203933 0.016501891 -1.8210252e-4 3.7899787e-4 35 8267.6293 50 2.095 0.1 0.994062723 0.010276207 -1.3567871e-4 2.3553259e-4 36 8863.0230 50 2.095 0.1 0.979647192 0.011089421 -4.6850452e-4 2.5791174e-4 37 9462.3637 50 2.095 0.1 0.97982075 0.013168527 -4.6446835e-4 3.0621222e-4 38 10067.841 50 2.095 0.1 0.977927712 0.014362819 -5.0853039e-4 3.3463000e-4 39 10678.912 50 2.095 0.1 0.978351462 0.013239476 -4.9865985e-4 3.0832436e-4 40 11295.987 50 2.095 0.1 0.981246471 0.013896898 -4.3133968e-4 3.2267975e-4 41 11920.227 50 2.095 0.1 0.978275164 0.013747481 -5.0043677e-4 3.2017989e-4 42 12552.713 50 2.095 0.1 0.976818842 0.013487944 -5.3437989e-4 3.1460359e-4 43 13193.071 50 2.095 0.1 0.981546068 0.011972943 -4.2438423e-4 2.7792152e-4 44 13841.978 50 2.095 0.1 0.96445854 0.013357572 -8.2452102e-4 3.1555561e-4 45 14500.513 50 2.095 0.1 0.972756158 0.014756189 -6.2933879e-4 3.4562263e-4 46 15169.786 50 2.095 0.1 0.971010184 0.014210504 -6.7027011e-4 3.3343996e-4 47 15850.446 50 2.095 0.1 0.975009829 0.013648839 -5.7661387e-4 3.1894710e-4 48 16543.562 50 2.095 0.1 0.969575209 0.014454336 -7.0396574e-4 3.3966327e-4 49 17249.754 50 2.095 0.1 0.959907267 0.012806944 -9.3229353e-4 3.0398222e-4 50 17969.627 50 2.095 0.1 0.96219203 0.011845496 -8.7812744e-4 2.8049391e-4 51 18706.130 50 2.095 0.1 0.960157966 0.011721155 -9.2634378e-4 2.7813758e-4 52 19459.212 50 2.095 0.1 0.95089937 0.009967287 -1.1471121e-3 2.3882201e-4 53 20230.984 50 2.095 0.1 0.955584642 0.01104368 -1.0351259e-3 2.6331561e-4 54 21021.899 50 2.095 0.1 0.950786938 0.011651582 -1.1498062e-3 2.7921171e-4 55 21835.345 50 2.095 0.1 0.951403607 0.008784149 -1.1350335e-3 2.1036178e-4 56 22671.730 50 2.095 0.1 0.94561823 0.009392383 -1.2740040e-3 2.2630383e-4 57 23533.930 50 2.095 0.1 0.945882787 0.009468264 -1.2676305e-3 2.2806833e-4 58 24423.504 50 2.095 0.1 0.939583158 0.009164239 -1.4198814e-3 2.2222511e-4 59 25343.343 50 2.095 0.1 0.93731257 0.008621568 -1.4750079e-3 2.0957224e-4 60 26294.734 50 2.095 0.1 0.936639557 0.010408821 -1.4913733e-3 2.5319842e-4 61 27282.004 50 2.095 0.1 0.932091884 0.007112674 -1.6022666e-3 1.7386258e-4 62 28307.499 50 2.095 0.1 0.933852763 0.009200431 -1.5592642e-3 2.2447176e-4 63 29373.630 50 2.095 0.1 0.930283437 0.01022719 -1.6465153e-3 2.5047996e-4 64 30485.110 50 2.095 0.1 0.930121334 0.007055747 -1.6504858e-3 1.7283645e-4 65 31645.899 50 2.095 0.1 0.925944326 0.0078917 -1.7530356e-3 1.9418588e-4 66 32858.438 50 2.095 0.1 0.92351759 0.005522582 -1.8128270e-3 1.3624763e-4 67 34128.058 50 2.095 0.1 0.918228311 0.006917328 -1.9436940e-3 1.7164045e-4 68 35460.168 50 2.095 0.1 0.908621563 0.006852313 -2.1833230e-3 1.7182490e-4 69 36858.762 50 2.095 0.1 0.919087792 0.006995961 -1.9223776e-3 1.7342924e-4 70 38331.748 50 2.095 0.1 0.910132214 0.006107336 -2.1454742e-3 1.5289007e-4 71 39882.345 50 2.095 0.1 0.904491462 0.007078101 -2.2871233e-3 1.7829708e-4 72 41519.404 50 2.095 0.1 0.90056919 0.005934428 -2.3861400e-3 1.5013907e-4 73 43249.023 50 2.095 0.1 0.903155701 0.006249069 -2.3207959e-3 1.5764661e-4 74 45079.637 50 2.095 0.1 0.894919827 0.006196232 -2.5295172e-3 1.5775222e-4 75 47019.670 50 2.095 0.1 0.896358908 0.007324999 -2.4929085e-3 1.8619052e-4 76 49077.851 50 2.095 0.1 0.88807923 0.004530465 -2.7043436e-3 1.1623128e-4 77 51264.711 50 2.095 0.1 0.886649039 0.008924412 -2.7410654e-3 2.2932944e-4 78 53590.598 50 2.095 0.1 0.882051947 0.006055081 -2.8595036e-3 1.5640756e-4 79 56067.403 50 2.095 0.1 0.882573371 0.005727419 -2.8460388e-3 1.4785638e-4 80 58709.036 50 2.095 0.1 0.876289825 0.006062306 -3.0088321e-3 1.5762389e-4 81 61527.362 50 2.095 0.1 0.877318199 0.005307783 -2.9821094e-3 1.3784403e-4 82 64538.932 50 2.095 0.1 0.873610284 0.006918104 -3.0786086e-3 1.8042690e-4 83 67758.354 50 2.095 0.1 0.865954455 0.005625863 -3.2791557e-3 1.4802192e-4 84 71205.306 50 2.095 0.1 0.868037308 0.008662141 -3.2244196e-3 2.2736248e-4 85 74897.359 50 2.095 0.1 0.86676404 0.005761209 -3.2578647e-3 1.5144143e-4 86 78855.560 50 2.095 0.1 0.85655072 0.005731411 -3.5279304e-3 1.5245456e-4 87 83102.902 50 2.095 0.1 0.857936274 0.006954572 -3.4911046e-3 1.8469168e-4 88 87663.154 50 2.095 0.1 0.85520525 0.006559622 -3.5637478e-3 1.7475934e-4 89 92562.717 50 2.095 0.1 0.849739829 0.008907377 -3.7098230e-3 2.3883381e-4 90 97830.316 50 2.095 0.1 0.841487161 0.006547659 -3.9321836e-3 1.7728439e-4 91 103497.14 50 2.095 0.1 0.84801051 0.007574866 -3.7562386e-3 2.0351933e-4 92 109596.79 50 2.095 0.1 0.843342861 0.006411993 -3.8819940e-3 1.7322909e-4 93 116165.97 50 2.095 0.1 0.837488454 0.008141353 -4.0407107e-3 2.2148775e-4 94 123244.47 50 2.095 0.1 0.839214922 0.009987535 -3.9937900e-3 2.7115465e-4 95 130874.72 50 2.095 0.1 0.834336972 0.011467406 -4.1266093e-3 3.1315233e-4 96 139102.96 50 2.095 0.1 0.828775634 0.012483365 -4.2789870e-3 3.4318369e-4 97 147980.51 50 2.095 0.1 0.823088812 0.016035848 -4.4358638e-3 4.4389186e-4 98 157560.98 50 2.095 0.1 0.824329178 0.014838138 -4.4015548e-3 4.1011975e-4 99 167904.68 50 2.095 0.1 0.82023172 0.015286187 -4.5150892e-3 4.2461424e-4 100 179076.29 50 2.095 0.1 0.827559752 0.012612741 -4.3124376e-3 3.4724985e-4 -
explore/asymint.py
ra1c32c2 rc462169 30 30 import numpy as np 31 31 import mpmath as mp 32 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 32 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos 33 33 from numpy.polynomial.legendre import leggauss 34 34 from scipy.integrate import dblquad, simps, romb, romberg … … 41 41 shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 42 42 Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 43 44 DTYPE = 'd' 43 45 44 46 class MPenv: … … 72 74 sas_sinx_x = staticmethod(sp.sas_sinx_x) 73 75 pi = np.pi 74 mpf = staticmethod(float) 76 #mpf = staticmethod(float) 77 mpf = staticmethod(lambda x: np.array(x, DTYPE)) 75 78 76 79 SLD = 3 … … 220 223 #A, B, C = 4450, 14000, 47 221 224 #A, B, C = 445, 140, 47 # integer for the sake of mpf 222 A, B, C = 6800, 114, 1380223 DA, DB, DC = 2 300, 21, 58225 A, B, C = 114, 1380, 6800 226 DA, DB, DC = 21, 58, 2300 224 227 SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 228 ## default parameters from sasmodels 229 #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 230 ## swap A-B-C to C-B-A 231 #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 225 232 #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 226 233 #SLD_SOLVENT,CONTRAST = 0, 4 … … 233 240 DA, DC = DC, DA 234 241 SLDA, SLDC = SLDC, SLDA 242 #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 235 243 NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 236 244 NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) … … 348 356 return n**2, Iq 349 357 358 def 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 350 379 def gridded_2d(q, n=300): 351 380 """ … … 395 424 print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 396 425 print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 426 print("gauss-20 usub", *gauss_quad_usub(Q, n=20)) 427 print("gauss-76 usub", *gauss_quad_usub(Q, n=76)) 428 print("gauss-150 usub", *gauss_quad_usub(Q, n=150)) 397 429 #gridded_2d(Q, n=2**8+1) 398 430 gridded_2d(Q, n=2**10+1) -
explore/realspace.py
r8fb2a94 r362a66f 3 3 import time 4 4 from copy import copy 5 import os 6 import argparse 7 from collections import OrderedDict 5 8 6 9 import numpy as np 7 10 from numpy import pi, radians, sin, cos, sqrt 8 from numpy.random import poisson, uniform 11 from numpy.random import poisson, uniform, randn, rand 9 12 from numpy.polynomial.legendre import leggauss 10 13 from scipy.integrate import simps 11 14 from scipy.special import j1 as J1 15 16 try: 17 import numba 18 USE_NUMBA = True 19 except ImportError: 20 USE_NUMBA = False 12 21 13 22 # Definition of rotation matrices comes from wikipedia: … … 82 91 raise NotImplementedError() 83 92 93 def dims(self): 94 # type: () -> float, float, float 95 raise NotImplementedError() 96 84 97 def rotate(self, theta, phi, psi): 85 98 self.rotation = rotation(theta, phi, psi) * self.rotation … … 116 129 for s2 in shapes] 117 130 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) 121 132 122 133 def sample(self, density): … … 133 144 self._scale = np.array([a/2, b/2, c/2])[None, :] 134 145 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 138 148 139 149 def sample(self, density): 140 num_points = poisson(density*self. a*self.b*self.c)150 num_points = poisson(density*self.volume) 141 151 points = self._scale*uniform(-1, 1, size=(num_points, 3)) 142 152 values = self.value.repeat(points.shape[0]) … … 152 162 self._scale = np.array([ra, rb, length/2])[None, :] 153 163 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 157 166 158 167 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 160 170 num_points = poisson(density*4*self.ra*self.rb*self.length) 161 171 points = uniform(-1, 1, size=(num_points, 3)) 162 172 radius = points[:, 0]**2 + points[:, 1]**2 163 points = self._scale*points[radius <= 1]173 points = points[radius <= 1] 164 174 values = self.value.repeat(points.shape[0]) 165 return values, self._adjust(points) 175 return values, self._adjust(self._scale*points) 176 177 class 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) 166 217 167 218 class TriaxialEllipsoid(Shape): … … 174 225 self._scale = np.array([ra, rb, rc])[None, :] 175 226 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 179 229 180 230 def sample(self, density): … … 194 244 self.rotate(*orientation) 195 245 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 196 248 self.helix_radius, self.helix_pitch = helix_radius, helix_pitch 197 249 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 203 252 # small tube radius approximation; for larger tubes need to account 204 253 # for the fact that the inner length is much shorter than the outer 205 254 # length 206 returnpi*self.tube_radius**2*self.tube_length255 self.volume = pi*self.tube_radius**2*self.tube_length 207 256 208 257 def points(self, density): … … 244 293 return values, self._adjust(points) 245 294 246 NUMBA = False 247 if NUMBA: 295 def 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 307 def _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 313 if USE_NUMBA: 314 # Override simple numpy solution with numba if available 248 315 from numba import njit 249 316 @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])") … … 252 319 for j in range(len(Iq)): 253 320 total = 0. + 0j 254 for k in range(len( Iq)):321 for k in range(len(values)): 255 322 total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k])) 256 323 Iq[j] = abs(total)**2 … … 263 330 values = rho*volume 264 331 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)] 265 334 266 335 # 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()) 272 337 return np.asarray(Iq).reshape(qx.shape) / np.sum(volume) 273 338 … … 338 403 """ 339 404 340 if NUMBA: 405 if USE_NUMBA: 406 # Override simple numpy solution with numba if available 341 407 @njit("f8[:](f8[:], f8[:], f8[:,:])") 342 def _calc_Pr_uniform _jit(r, rho, points):408 def _calc_Pr_uniform(r, rho, points): 343 409 dr = r[0] 344 410 n_max = len(r) … … 362 428 # P(r) with uniform steps in r is 3x faster; check if we are uniform 363 429 # before continuing 430 r, rho, points = [np.asarray(v, 'd') for v in (r, rho, points)] 364 431 if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 365 432 Pr = _calc_Pr_nonuniform(r, rho, points) 366 433 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) 371 435 return Pr / Pr.max() 372 436 … … 399 463 return edges 400 464 401 def plot_calc(r, Pr, q, Iq, theory=None): 465 # -------------- plotters ---------------- 466 def plot_calc(r, Pr, q, Iq, theory=None, title=None): 402 467 import matplotlib.pyplot as plt 403 468 plt.subplot(211) … … 405 470 plt.xlabel('r (A)') 406 471 plt.ylabel('Pr (1/A^2)') 472 if title is not None: 473 plt.title(title) 407 474 plt.subplot(212) 408 475 plt.loglog(q, Iq, '-', label='from Pr') … … 410 477 plt.ylabel('Iq') 411 478 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') 413 480 plt.legend() 414 481 415 def plot_calc_2d(qx, qy, Iqxy, theory=None ):482 def plot_calc_2d(qx, qy, Iqxy, theory=None, title=None): 416 483 import matplotlib.pyplot as plt 417 484 qx, qy = bin_edges(qx), bin_edges(qy) … … 419 486 if theory is not None: 420 487 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') 422 492 plt.xlabel('qx (1/A)') 423 493 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) 424 499 if theory is not None: 425 500 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) 427 505 plt.xlabel('qx (1/A)') 428 506 … … 440 518 index = np.random.choice(n, size=500) if n > 500 else slice(None, None) 441 519 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') 442 523 #low, high = points.min(axis=0), points.max(axis=0) 443 524 #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") 444 528 ax.autoscale(True) 445 529 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 -------------- 486 531 def sas_sinx_x(x): 487 532 with np.errstate(all='ignore'): … … 510 555 for k, qk in enumerate(q): 511 556 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) 513 558 Iq[k] = np.sum(w*Fq**2) 514 Iq = Iq /Iq[0]559 Iq = Iq 515 560 return Iq 516 561 517 562 def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)): 518 563 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) 521 566 Iq = Fq**2 522 567 return Iq.reshape(qx.shape) … … 524 569 def sphere_Iq(q, radius): 525 570 Iq = sas_3j1x_x(q*radius)**2 526 return Iq/Iq[0] 571 return Iq 572 573 def 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 591 def 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) 527 599 528 600 def csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core): … … 539 611 sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 540 612 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) 543 615 inner_sum = np.zeros_like(q) 544 616 for beta, inner_w in zip((z + 1)*pi/4, w): 545 617 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) 550 622 if overlapping: 551 623 Fq = (dr0*siA*siB*siC … … 584 656 return Iq.reshape(qx.shape) 585 657 586 def check_cylinder(radius=25, length=125, rho=2.): 658 def 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 666 def 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 681 def 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 691 def build_cylinder(radius=25, length=125, rho=2.): 587 692 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 697 def 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 703 def 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 709 def 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 716 def 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 730 def 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 753 def 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)) 611 761 for ix in range(nx) 612 762 for iy in range(ny) 613 763 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 768 SHAPE_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 ]) 776 SHAPES = list(SHAPE_FUNCTIONS.keys()) 777 778 def 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 799 def 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 828 def 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 641 873 642 874 if __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 137 137 """ 138 138 _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): 140 140 # type: (Data, Model, float) -> None 141 141 # remember inputs so we can inspect from outside … … 145 145 self._interpret_data(data, model.sasmodel) 146 146 self._cache = {} 147 self.extra_pars = extra_pars 147 148 148 149 def update(self): … … 166 167 Return a dictionary of parameters 167 168 """ 168 return self.model.parameters() 169 pars = self.model.parameters() 170 if self.extra_pars: 171 pars.update(self.extra_pars) 172 return pars 169 173 170 174 def theory(self): -
sasmodels/compare.py
r2a7e20e r1fbadb2 40 40 from . import core 41 41 from . import kerneldll 42 from . import kernelcl 42 43 from .data import plot_theory, empty_data1D, empty_data2D, load_data 43 44 from .direct_model import DirectModel, get_mesh … … 623 624 624 625 626 def 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 646 def 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 625 673 DTYPE_MAP = { 626 674 'half': '16', … … 643 691 Return a model calculator using the OpenCL calculation engine. 644 692 """ 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 calculator651 693 652 694 def eval_ctypes(model_info, data, dtype='double', cutoff=0.): … … 660 702 return calculator 661 703 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 to668 initialize the calculation engine, and make the average more stable.669 """670 # initialize the code so time is more accurate671 if evals > 1:672 calculator(**suppress_pd(pars))673 toc = tic()674 # make sure there is at least one eval675 value = calculator(**pars)676 for _ in range(evals-1):677 value = calculator(**pars)678 average_time = toc()*1000. / evals679 #print("I(q)",value)680 return value, average_time681 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 points686 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.ndarray694 data = empty_data2D(q, resolution=res)695 data.accuracy = opts['accuracy']696 set_beam_stop(data, qmin)697 index = ~data.mask698 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, index708 709 704 def make_engine(model_info, data, dtype, cutoff, ngauss=0): 710 705 # type: (ModelInfo, Data, str, float) -> Calculator … … 718 713 set_integration_size(model_info, ngauss) 719 714 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 724 725 725 726 def _show_invalid(data, theory): … … 751 752 """ 752 753 for k in range(opts['sets']): 753 if k > 1:754 if k > 0: 754 755 # print a separate seed for each dataset for better reproducibility 755 756 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)) 757 758 np.random.seed(new_seed) 758 759 opts['pars'] = parse_pars(opts, maxdim=maxdim) … … 761 762 result = run_models(opts, verbose=True) 762 763 if opts['plot']: 764 if opts['is2d'] and k > 0: 765 import matplotlib.pyplot as plt 766 plt.figure() 763 767 limits = plot_models(opts, result, limits=limits, setnum=k) 764 768 if opts['show_weights']: … … 1328 1332 1329 1333 # Evaluate preset parameter expressions 1334 # Note: need to replace ':' with '_' in parameter names and expressions 1335 # in order to support math on magnetic parameters. 1330 1336 context = MATH.copy() 1331 1337 context['np'] = np 1332 context.update( pars)1338 context.update((k.replace(':', '_'), v) for k, v in pars.items()) 1333 1339 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) 1334 1341 for k, v in presets.items(): 1335 1342 if not isinstance(v, float) and not k.endswith('_type'): 1336 presets[k] = eval(v , context)1343 presets[k] = eval(v.replace(':', '_'), context) 1337 1344 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)) 1339 1346 for k, v in presets2.items(): 1340 1347 if not isinstance(v, float) and not k.endswith('_type'): 1341 presets2[k] = eval(v , context)1348 presets2[k] = eval(v.replace(':', '_'), context) 1342 1349 1343 1350 # update parameters with presets … … 1369 1376 path = os.path.dirname(info.filename) 1370 1377 url = "file://" + path.replace("\\", "/")[2:] + "/" 1371 rst2html.view_html_ qtapp(html, url)1378 rst2html.view_html_wxapp(html, url) 1372 1379 1373 1380 def explore(opts): … … 1490 1497 vmin, vmax = limits 1491 1498 self.limits = vmax*1e-7, 1.3*vmax 1492 import pylab; pylab.clf() 1499 import pylab 1500 pylab.clf() 1493 1501 plot_models(self.opts, result, limits=self.limits) 1494 1502 -
sasmodels/convert.py
r2d81cfe ra69d8cd 633 633 def test_backward_forward(): 634 634 from .core import list_models 635 L = lambda name: _check_one(name, seed=1) 635 636 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 21 21 from . import mixture 22 22 from . import kernelpy 23 from . import kernelcl 23 24 from . import kerneldll 24 25 from . import custom 25 26 if os.environ.get("SAS_OPENCL", "").lower() == "none":27 HAVE_OPENCL = False28 else:29 try:30 from . import kernelcl31 HAVE_OPENCL = True32 except Exception:33 HAVE_OPENCL = False34 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)40 26 41 27 # pylint: disable=unused-import … … 47 33 pass 48 34 # pylint: enable=unused-import 35 36 CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "") 37 if 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) 49 41 50 42 # TODO: refactor composite model support … … 148 140 used with functions in generate to build the docs or extract model info. 149 141 """ 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) 193 169 194 170 … … 292 268 if platform is None: 293 269 platform = "ocl" 294 if platform == "ocl" and not HAVE_OPENCLor not model_info.opencl:270 if not kernelcl.use_opencl() or not model_info.opencl: 295 271 platform = "dll" 296 272 … … 336 312 print("\n".join(list_models(kind))) 337 313 314 def 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 361 def 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 338 371 if __name__ == "__main__": 339 372 list_models_main() -
sasmodels/data.py
r2d81cfe rd86f0fc 706 706 else: 707 707 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, 709 711 interpolation='nearest', aspect=1, origin='lower', 710 712 extent=[xmin, xmax, ymin, ymax], … … 714 716 return vmin, vmax 715 717 718 719 # === The following is modified from sas.sasgui.plottools.PlotPanel 720 def _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 779 def _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 819 def _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 716 888 def demo(): 717 889 # type: () -> None -
sasmodels/details.py
r108e70e r885753a 243 243 offset = np.cumsum(np.hstack((0, length))) 244 244 call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 245 # Pad value array to a 32 value boundary d245 # Pad value array to a 32 value boundary 246 246 data_len = nvalues + 2*sum(len(v) for v in dispersity) 247 247 extra = (32 - data_len%32)%32 … … 250 250 is_magnetic = convert_magnetism(kernel.info.parameters, data) 251 251 #call_details.show() 252 #print("data", data) 252 253 return call_details, data, is_magnetic 253 254 … … 296 297 mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 297 298 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() 300 301 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 305 305 return True 306 306 else: -
sasmodels/direct_model.py
r2d81cfe rd18d6dd 345 345 self._kernel = self._model.make_kernel(self._kernel_inputs) 346 346 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 347 352 Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 348 353 # Storing the calculated Iq values so that they can be plotted. … … 357 362 np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 358 363 ) 359 return result 364 return result + background 360 365 361 366 -
sasmodels/generate.py
r108e70e rd86f0fc 169 169 170 170 import sys 171 from os.path import abspath, dirname, join as joinpath, exists, getmtime 171 from os import environ 172 from os.path import abspath, dirname, join as joinpath, exists, getmtime, sep 172 173 import re 173 174 import string … … 214 215 215 216 EXTERNAL_DIR = 'sasmodels-data' 216 DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_ template.c')217 DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_iq.c') 217 218 MODEL_PATH = joinpath(DATA_PATH, 'models') 218 219 … … 289 290 import loops. 290 291 """ 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): 293 293 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): 296 296 gengauss(n, path) 297 297 info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 298 298 else lib for lib in info.source] 299 299 300 300 def format_units(units): … … 320 320 for w, h in zip(column_widths, PARTABLE_HEADERS)] 321 321 322 sep= " ".join("="*w for w in column_widths)322 underbar = " ".join("="*w for w in column_widths) 323 323 lines = [ 324 sep,324 underbar, 325 325 " ".join("%-*s" % (w, h) 326 326 for w, h in zip(column_widths, PARTABLE_HEADERS)), 327 sep,327 underbar, 328 328 ] 329 329 for p in pars: … … 334 334 "%*g" % (column_widths[3], p.default), 335 335 ])) 336 lines.append( sep)336 lines.append(underbar) 337 337 return "\n".join(lines) 338 338 … … 389 389 # TODO: fails DRY; templates appear two places. 390 390 model_templates = [joinpath(DATA_PATH, filename) 391 for filename in ('kernel_header.c', 'kernel_iq.c l')]391 for filename in ('kernel_header.c', 'kernel_iq.c')] 392 392 source_files = (model_sources(model_info) 393 393 + model_templates … … 612 612 """ 613 613 spaces = " "*depth 614 sep= "\n" + spaces615 return spaces + sep.join(s.split("\n"))614 interline_separator = "\n" + spaces 615 return spaces + interline_separator.join(s.split("\n")) 616 616 617 617 … … 619 619 def load_template(filename): 620 620 # type: (str) -> str 621 """ 622 Load template file from sasmodels resource directory. 623 """ 621 624 path = joinpath(DATA_PATH, filename) 622 625 mtime = getmtime(path) … … 900 903 kernel_module = load_custom_kernel_module(model_name) 901 904 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 905 918 return kernel_module 906 919 -
sasmodels/guyou.py
r9d9fcbd rb3703f5 31 31 # 32 32 # 2017-11-01 Paul Kienzle 33 # * converted to python, using degrees rather than radians 33 # * converted to python, with degrees rather than radians 34 """ 35 Convert between latitude-longitude and Guyou map coordinates. 36 """ 37 34 38 from __future__ import division, print_function 35 39 36 40 import numpy as np 37 from numpy import sqrt, pi, tan, cos, sin, log, exp, arctan2 as atan2, sign, radians, degrees 41 from numpy import sqrt, pi, tan, cos, sin, sign, radians, degrees 42 from numpy import sinh, arctan as atan 43 44 # scipy version of special functions 45 from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF 38 46 39 47 _ = """ … … 59 67 """ 60 68 61 # scipy version of special functions62 from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF63 from numpy import sinh, sign, arctan as atan64 65 69 def ellipticJi(u, v, m): 66 70 scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m) 67 71 u, v, m = np.broadcast_arrays(u, v, m) 68 result = np.empty_like([u, u,u], 'D')69 real = v==070 imag = u==072 result = np.empty_like([u, u, u], 'D') 73 real = (v == 0) 74 imag = (u == 0) 71 75 mixed = ~(real|imag) 72 76 result[:, real] = _ellipticJi_real(u[real], m[real]) 73 77 result[:, imag] = _ellipticJi_imag(v[imag], m[imag]) 74 78 result[:, mixed] = _ellipticJi(u[mixed], v[mixed], m[mixed]) 75 return result[0, :] if scalar else result79 return result[0, :] if scalar else result 76 80 77 81 def _ellipticJi_real(u, m): … … 104 108 result = np.empty_like(phi, 'D') 105 109 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]) 107 112 result[~index] = ellipticFi(phi[~index], psi[~index], m[~index]) 108 113 return result.reshape(1)[0] if scalar else result … … 117 122 cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2 118 123 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) 120 126 return re + 1j*im 121 127 122 sqrt2 = sqrt(2)128 SQRT2 = sqrt(2) 123 129 124 130 # [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion … … 127 133 # K = 3.165103454447431823666270142140819753058976299237578486994... 128 134 def guyou(lam, phi): 135 """Transform from (latitude, longitude) to point (x, y)""" 129 136 # [PAK] wrap into [-pi/2, pi/2] radians 130 137 x, y = np.asarray(lam), np.asarray(phi) … … 135 142 136 143 # Compute constant K 137 cos_u = ( sqrt2 - 1) / (sqrt2 + 1)144 cos_u = (SQRT2 - 1) / (SQRT2 + 1) 138 145 sinsq_u = 1 - cos_u**2 139 146 K = ellipticF(pi/2, sinsq_u) … … 144 151 at = atan(r * (cos(lam) - 1j*sin(lam))) 145 152 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)) 147 154 148 155 # [PAK] convert to degrees, and return to original tile … … 150 157 151 158 def guyou_invert(x, y): 159 """Transform from point (x, y) on plot to (latitude, longitude)""" 152 160 # [PAK] wrap into [-pi/2, pi/2] radians 153 161 x, y = np.asarray(x), np.asarray(y) … … 158 166 159 167 # compute constant K 160 cos_u = ( sqrt2 - 1) / (sqrt2 + 1)168 cos_u = (SQRT2 - 1) / (SQRT2 + 1) 161 169 sinsq_u = 1 - cos_u**2 162 170 K = ellipticF(pi/2, sinsq_u) … … 174 182 175 183 def plot_grid(): 184 """Plot the latitude-longitude grid for Guyou transform""" 176 185 import matplotlib.pyplot as plt 177 186 from numpy import linspace … … 186 195 plt.plot(x, y, 'g') 187 196 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) 190 199 plt.plot(x, y, 'b') 191 200 #plt.xlabel('longitude') 192 201 plt.ylabel('latitude') 202 plt.title('forward transform') 193 203 194 204 plt.subplot(212) … … 202 212 plt.xlabel('longitude') 203 213 plt.ylabel('latitude') 214 plt.title('inverse transform') 215 216 def main(): 217 """Show the Guyou transformation""" 218 plot_grid() 219 import matplotlib.pyplot as plt 220 plt.show() 204 221 205 222 if __name__ == "__main__": 206 plot_grid() 207 import matplotlib.pyplot as plt; plt.show() 208 223 main() 209 224 210 225 _ = """ -
sasmodels/jitter.py
r8cfb486 rb3703f5 1 1 #!/usr/bin/env python 2 2 """ 3 Application to explore the difference between sasview 3.x orientation 4 dispersity and possible replacement algorithms. 3 Jitter Explorer 4 =============== 5 6 Application to explore orientation angle and angular dispersity. 5 7 """ 6 8 from __future__ import division, print_function 7 9 8 import sys, os9 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))10 11 10 import argparse 12 11 13 import mpl_toolkits.mplot3d # Adds projection='3d' option to subplot 12 try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d 13 import mpl_toolkits.mplot3d # Adds projection='3d' option to subplot 14 except ImportError: 15 pass 16 14 17 import matplotlib.pyplot as plt 15 from matplotlib.widgets import Slider, CheckButtons 16 from matplotlib import cm 18 from matplotlib.widgets import Slider 17 19 import numpy as np 18 20 from numpy import pi, cos, sin, sqrt, exp, degrees, radians 19 21 20 def draw_beam(ax , view=(0, 0)):22 def draw_beam(axes, view=(0, 0)): 21 23 """ 22 24 Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 23 25 """ 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) 26 28 27 29 steps = 25 … … 40 42 x, y, z = [v.reshape(shape) for v in points] 41 43 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 46 def 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 67 def 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 72 def 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 87 def 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 101 def _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 108 def _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 125 def 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 165 def 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 175 def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), 176 draw_shape=draw_parallelepiped): 45 177 """ 46 178 Represent jitter as a set of shapes at different orientations. … … 49 181 scale = 0.95/sqrt(sum(v**2 for v in size)) 50 182 size = tuple(scale*v for v in size) 51 draw_shape = draw_parallelepiped52 #draw_shape = draw_ellipsoid53 183 54 184 #np.random.seed(10) … … 56 186 cloud = [ 57 187 [-1, -1, -1], 58 [-1, -1, 59 [-1, -1, 60 [-1, 61 [-1, 0,0],62 [-1, 0,1],63 [-1, 64 [-1, 1,0],65 [-1, 1,1],66 [ 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 [ 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], 84 214 ] 85 215 dtheta, dphi, dpsi = jitter … … 90 220 if dpsi == 0: 91 221 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) 93 223 scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 94 224 for point in cloud: 95 225 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) 97 227 for v in 'xyz': 98 228 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) 162 232 163 233 PROJECTIONS = [ … … 167 237 'azimuthal_equal_area', 168 238 ] 169 def draw_mesh(ax , view, jitter, radius=1.2, n=11, dist='gaussian',239 def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 170 240 projection='equirectangular'): 171 241 """ … … 226 296 Should allow free movement in theta, but phi is distorted. 227 297 """ 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) 233 302 if dist == 'gaussian': 234 t*= 3235 weights = exp(-0.5* t**2)303 dist_x *= 3 304 weights = exp(-0.5*dist_x**2) 236 305 elif dist == 'rectangle': 237 306 # Note: uses sasmodels ridiculous definition of rectangle width 238 t*= sqrt(3)307 dist_x *= sqrt(3) 239 308 elif dist == 'uniform': 240 309 pass … … 243 312 244 313 if projection == 'equirectangular': #define PROJECTION 1 245 def rotate(theta_i, phi_j):314 def _rotate(theta_i, phi_j): 246 315 return Rx(phi_j)*Ry(theta_i) 247 def weight(theta_i, phi_j, wi, wj):248 return w i*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))) 249 318 elif projection == 'sinusoidal': #define PROJECTION 2 250 def rotate(theta_i, phi_j):319 def _rotate(theta_i, phi_j): 251 320 latitude = theta_i 252 321 scale = cos(radians(latitude)) … … 254 323 #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 255 324 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): 257 326 latitude = theta_i 258 327 scale = cos(radians(latitude)) 259 w= 1 if abs(phi_j) < abs(scale)*180 else 0260 return w*wi*wj328 active = 1 if abs(phi_j) < abs(scale)*180 else 0 329 return active*w_i*w_j 261 330 elif projection == 'guyou': #define PROJECTION 3 (eventually?) 262 def rotate(theta_i, phi_j):263 from guyou import guyou_invert331 def _rotate(theta_i, phi_j): 332 from .guyou import guyou_invert 264 333 #latitude, longitude = guyou_invert([theta_i], [phi_j]) 265 334 longitude, latitude = guyou_invert([phi_j], [theta_i]) 266 335 return Rx(longitude[0])*Ry(latitude[0]) 267 def weight(theta_i, phi_j, wi, wj):268 return w i*wj336 def _weight(theta_i, phi_j, w_i, w_j): 337 return w_i*w_j 269 338 elif projection == 'azimuthal_equidistance': # Note: Rz Ry, not Rx Ry 270 def rotate(theta_i, phi_j):339 def _rotate(theta_i, phi_j): 271 340 latitude = sqrt(theta_i**2 + phi_j**2) 272 341 longitude = degrees(np.arctan2(phi_j, theta_i)) 273 342 #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 274 343 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): 276 345 # Weighting for each point comes from the integral: 277 346 # \int\int I(q, lat, log) sin(lat) dlat dlog … … 303 372 # the entire sphere, and treats theta and phi identically. 304 373 latitude = sqrt(theta_i**2 + phi_j**2) 305 w = sin(radians(latitude))/latitude if latitude != 0 else 1306 return w *wi*wj if latitude < 180 else 0374 weight = sin(radians(latitude))/latitude if latitude != 0 else 1 375 return weight*w_i*w_j if latitude < 180 else 0 307 376 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)) 311 380 longitude = degrees(np.arctan2(phi_j, theta_i)) 312 381 #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 313 382 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): 315 384 latitude = sqrt(theta_i**2 + phi_j**2) 316 w = sin(radians(latitude))/latitude if latitude != 0 else 1317 return w *wi*wj if latitude < 180 else 0385 weight = sin(radians(latitude))/latitude if latitude != 0 else 1 386 return weight*w_i*w_j if latitude < 180 else 0 318 387 else: 319 388 raise ValueError("unknown projection %r"%projection) 320 389 321 390 # mesh in theta, phi formed by rotating z 391 dtheta, dphi, dpsi = jitter 322 392 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) 344 403 345 404 # rotate relative to beam … … 348 407 x, y, z = [np.array(v).flatten() for v in points] 349 408 #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 411 def draw_labels(axes, view, jitter, text): 353 412 """ 354 413 Draw text at a particular location. … … 364 423 for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 365 424 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) 367 426 368 427 # Definition of rotation matrices comes from wikipedia: … … 370 429 def Rx(angle): 371 430 """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) 377 436 378 437 def Ry(angle): 379 438 """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) 385 444 386 445 def Rz(angle): 387 446 """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) 393 452 394 453 def transform_xyz(view, jitter, x, y, z): … … 398 457 x, y, z = [np.asarray(v) for v in (x, y, z)] 399 458 shape = x.shape 400 points = np.matrix([x.flatten(), y.flatten(),z.flatten()])459 points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 401 460 points = apply_jitter(jitter, points) 402 461 points = orient_relative_to_beam(view, points) … … 455 514 return data[offset], data[-1] 456 515 457 def draw_scattering(calculator, ax , view, jitter, dist='gaussian'):516 def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): 458 517 """ 459 518 Plot the scattering for the particular view. 460 519 461 *calculator* is returned from :func:`build_model`. *ax * are the 3D axes520 *calculator* is returned from :func:`build_model`. *axes* are the 3D axes 462 521 on which the data will be plotted. *view* and *jitter* are the current 463 522 orientation and orientation dispersity. *dist* is one of the sasmodels … … 472 531 theta, phi, psi = view 473 532 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)] 475 534 ## increase pd_n for testing jitter integration rather than simple viz 476 535 #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] … … 500 559 if 0: 501 560 level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 502 level[level <0] = 0561 level[level < 0] = 0 503 562 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) 505 564 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)) 508 567 else: 509 ax .pcolormesh(qx, qy, Iqxy)568 axes.pcolormesh(qx, qy, Iqxy) 510 569 511 570 def build_model(model_name, n=150, qmax=0.5, **pars): … … 524 583 for details. 525 584 """ 526 from sasmodels.core import load_model_info, build_model 585 from sasmodels.core import load_model_info, build_model as build_sasmodel 527 586 from sasmodels.data import empty_data2D 528 587 from sasmodels.direct_model import DirectModel 529 588 530 589 model_info = load_model_info(model_name) 531 model = build_ model(model_info) #, dtype='double!')590 model = build_sasmodel(model_info) #, dtype='double!') 532 591 q = np.linspace(-qmax, qmax, n) 533 592 data = empty_data2D(q, q) … … 547 606 return calculator 548 607 549 def select_calculator(model_name, n=150, size=(10, 40,100)):608 def select_calculator(model_name, n=150, size=(10, 40, 100)): 550 609 """ 551 610 Create a model calculator for the given shape. … … 561 620 """ 562 621 a, b, c = size 622 d_factor = 0.06 # for paracrystal models 563 623 if model_name == 'sphere': 564 624 calculator = build_model('sphere', n=n, radius=c) 565 625 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) 566 642 elif model_name == 'bcc_paracrystal': 567 calculator = build_model('bcc_paracrystal', n=n, dnn=c,568 d_factor=0.06, radius=40)569 643 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) 570 651 elif model_name == 'cylinder': 571 652 calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) … … 588 669 589 670 SHAPES = [ 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 677 DRAW_SHAPES = { 678 'fcc_paracrystal': draw_fcc, 679 'bcc_paracrystal': draw_bcc, 680 'sc_paracrystal': draw_sc, 681 'parallelepiped': draw_parallelepiped, 682 } 594 683 595 684 DISTRIBUTIONS = [ … … 608 697 Show an interactive orientation and jitter demo. 609 698 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. 611 710 """ 612 711 # projection number according to 1-order position in list, but … … 620 719 # set up calculator 621 720 calculator, size = select_calculator(model_name, n=150, size=size) 721 draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped) 622 722 623 723 ## uncomment to set an independent the colour range for every view … … 639 739 plt.gcf().canvas.set_window_title(projection) 640 740 #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') 643 743 try: # CRUFT: not all versions of matplotlib accept 'square' 3d projection 644 ax .axis('square')744 axes.axis('square') 645 745 except Exception: 646 746 pass … … 649 749 650 750 ## add control widgets to plot 651 ax theta= plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor)652 ax phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor)653 ax psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor)654 stheta = Slider(ax theta, 'Theta', -90, 90, valinit=theta)655 sphi = Slider(ax phi, 'Phi', -180, 180, valinit=phi)656 spsi = Slider(ax psi, 'Psi', -180, 180, valinit=psi)657 658 ax dtheta= plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor)659 ax dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor)660 ax dpsi= 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) 661 761 # Note: using ridiculous definition of rectangle distribution, whose width 662 762 # in sasmodels is sqrt(3) times the given width. Divide by sqrt(3) to keep 663 763 # the maximum width to 90. 664 764 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 668 769 669 770 ## callback to draw the new view … … 675 776 limit = [0, 0, 0.5, 5][dims] 676 777 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) 683 784 plt.gcf().canvas.draw() 684 785 685 786 ## bind control widgets to view updater 686 stheta.on_changed(lambda v: update(v, 'theta'))787 stheta.on_changed(lambda v: update(v, 'theta')) 687 788 sphi.on_changed(lambda v: update(v, 'phi')) 688 789 spsi.on_changed(lambda v: update(v, 'psi')) … … 702 803 formatter_class=argparse.ArgumentDefaultsHelpFormatter, 703 804 ) 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') 709 817 opts = parser.parse_args() 710 818 size = tuple(int(v) for v in opts.size.split(',')) -
sasmodels/kernel_iq.c
r924a119 rdc6f601 79 79 // du * (m_sigma_y + 1j*m_sigma_z); 80 80 // 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])81 static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 82 82 { 83 83 in_spin = clip(in_spin, 0.0, 1.0); 84 84 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 91 94 } 92 95 93 96 // Compute the magnetic sld 94 97 static double mag_sld( 95 const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=upimag98 const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag 96 99 const double qx, const double qy, 97 100 const double px, const double py, … … 103 106 const double perp = qy*mx - qx*my; 104 107 switch (xs) { 108 default: // keep compiler happy; condition ensures xs in [0,1,2,3] 105 109 case 0: // uu => sld - D M_perpx 106 110 return sld - px*perp; 107 case 1: // ud 111 case 1: // ud.real => -D M_perpy 108 112 return py*perp; 109 case 2: // du 113 case 2: // du.real => -D M_perpy 110 114 return py*perp; 111 case 3: // dd real=> sld + D M_perpx115 case 3: // dd => sld + D M_perpx 112 116 return sld + px*perp; 113 117 } 114 118 } else { 115 119 if (xs== 4) { 116 return -mz; // ud imag => -D M_perpz120 return -mz; // du.imag => +D M_perpz 117 121 } else { // index == 5 118 return mz; // du imag =>D M_perpz122 return +mz; // ud.imag => -D M_perpz 119 123 } 120 124 } … … 172 176 // Apply the rotation matrix returned from qac_rotation to the point (qx,qy), 173 177 // returning R*[qx,qy]' = [qa,qc]' 174 static double178 static void 175 179 qac_apply( 176 180 QACRotation *rotation, … … 245 249 // Apply the rotation matrix returned from qabc_rotation to the point (qx,qy), 246 250 // returning R*[qx,qy]' = [qa,qb,qc]' 247 static double251 static void 248 252 qabc_apply( 249 253 QABCRotation *rotation, … … 311 315 // up_angle = values[NUM_PARS+4]; 312 316 // 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, fill317 double xs_weights[8]; // uu, ud real, du real, dd, ud imag, du imag, fill, fill 314 318 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); 316 320 SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 317 321 #endif // MAGNETIC … … 659 663 660 664 // 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]; 663 667 if (xs_weight > 1.e-8) { 664 668 // Since the cross section weight is significant, set the slds … … 673 677 local_values.vector[sld_index] = 674 678 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); 675 681 } 676 682 scattering += xs_weight * CALL_KERNEL(); -
sasmodels/kernelcl.py
r2d81cfe rd86f0fc 58 58 import numpy as np # type: ignore 59 59 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. 60 63 try: 61 #raise NotImplementedError("OpenCL not yet implemented for new kernel template")62 64 import pyopencl as cl # type: ignore 65 from pyopencl import mem_flags as mf 66 from pyopencl.characterize import get_fast_inaccurate_build_options 63 67 # Ask OpenCL for the default context so that we know that one exists 64 68 cl.create_some_context(interactive=False) 69 HAVE_OPENCL = True 70 OPENCL_ERROR = "" 65 71 except 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) 72 74 73 75 from . import generate … … 102 104 cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS] 103 105 104 fix_pyopencl_include() 105 106 if HAVE_OPENCL: 107 fix_pyopencl_include() 106 108 107 109 # The max loops number is limited by the amount of local memory available … … 128 130 """ 129 131 132 def use_opencl(): 133 return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none" 130 134 131 135 ENV = None 136 def 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 132 143 def environment(): 133 144 # type: () -> "GpuEnvironment" … … 137 148 This provides an OpenCL context and one queue per device. 138 149 """ 139 global ENV140 150 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") 142 157 return ENV 143 158 … … 193 208 Build a model to run on the gpu. 194 209 195 Returns the compiled program and its type. The returned type will196 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. 198 213 """ 199 214 dtype = np.dtype(dtype) -
sasmodels/kerneldll.py
r1ddb794 r1662ebe 123 123 # add openmp support if not running on a mac 124 124 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 126 129 def compile_command(source, output): 127 130 """unix compiler command""" … … 158 161 return CC + [source, "-o", output, "-lm"] 159 162 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. 164 DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 169 165 170 166 ALLOW_SINGLE_PRECISION_DLLS = True … … 234 230 235 231 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. 237 233 """ 238 234 if dtype == F16: … … 251 247 need_recompile = dll_time < newest_source 252 248 if need_recompile: 249 # Make sure the DLL path exists 250 if not os.path.exists(DLL_PATH): 251 os.makedirs(DLL_PATH) 253 252 basename = splitext(os.path.basename(dll))[0] + "_" 254 253 system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) -
sasmodels/mixture.py
r2d81cfe recf895e 94 94 # If model is a sum model, each constituent model gets its own scale parameter 95 95 scale_prefix = prefix 96 if prefix == '' and part.operation== '*':96 if prefix == '' and getattr(part, "operation", '') == '*': 97 97 # `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_scale98 # it's parameters to form a new prefix for the scale. 99 # For example, a model with A*B*C will have ABC_scale. 100 100 sub_prefixes = [] 101 101 for param in part.parameters.kernel_parameters: … … 233 233 return self 234 234 235 def next(self):235 def __next__(self): 236 236 # type: () -> Tuple[List[Callable], CallDetails, np.ndarray] 237 237 if self.part_num >= len(self.parts): … … 251 251 252 252 return kernel, call_details, values 253 254 # CRUFT: py2 support 255 next = __next__ 253 256 254 257 def _part_details(self, info, par_index): -
sasmodels/model_test.py
r2d81cfe r3221de0 62 62 from .exception import annotate_exception 63 63 from .modelinfo import expand_pars 64 from .kernelcl import use_opencl 64 65 65 66 # pylint: disable=unused-import … … 123 124 124 125 if is_py: # kernel implemented in python 125 test_name = " Model: %s, Kernel:python"%model_name126 test_name = "%s-python"%model_name 126 127 test_method_name = "test_%s_python" % model_info.id 127 128 test = ModelTestCase(test_name, model_info, … … 134 135 135 136 # test using dll if desired 136 if 'dll' in loaders or not core.HAVE_OPENCL:137 test_name = " Model: %s, Kernel:dll"%model_name137 if 'dll' in loaders or not use_opencl(): 138 test_name = "%s-dll"%model_name 138 139 test_method_name = "test_%s_dll" % model_info.id 139 140 test = ModelTestCase(test_name, model_info, … … 145 146 146 147 # test using opencl if desired and available 147 if 'opencl' in loaders and core.HAVE_OPENCL:148 test_name = " Model: %s, Kernel: OpenCL"%model_name148 if 'opencl' in loaders and use_opencl(): 149 test_name = "%s-opencl"%model_name 149 150 test_method_name = "test_%s_opencl" % model_info.id 150 151 # Using dtype=None so that the models that are only … … 160 161 161 162 return suite 162 163 163 164 164 def _hide_model_case_from_nose(): … … 368 368 369 369 # Build a test suite containing just the model 370 loaders = ['opencl'] if core.HAVE_OPENCLelse ['dll']370 loaders = ['opencl'] if use_opencl() else ['dll'] 371 371 models = [model] 372 372 try: … … 425 425 verbosity = 1 426 426 if models and models[0] == 'opencl': 427 if not core.HAVE_OPENCL:427 if not use_opencl(): 428 428 print("opencl is not available") 429 429 return 1 … … 435 435 models = models[1:] 436 436 elif models and models[0] == 'opencl_and_dll': 437 loaders = ['opencl', 'dll'] if core.HAVE_OPENCLelse ['dll']437 loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 438 438 models = models[1:] 439 439 else: 440 loaders = ['opencl', 'dll'] if core.HAVE_OPENCLelse ['dll']440 loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 441 441 if not models: 442 442 print("""\ … … 467 467 Run "nosetests sasmodels" on the command line to invoke it. 468 468 """ 469 loaders = ['opencl', 'dll'] if core.HAVE_OPENCLelse ['dll']469 loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 470 470 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) 480 490 481 491 -
sasmodels/modelinfo.py
r5ab99b7 r1662ebe 74 74 processed.append(parse_parameter(*p)) 75 75 partable = ParameterTable(processed) 76 partable.check_angles() 76 77 return partable 77 78 … … 426 427 # type: (List[Parameter]) -> None 427 428 self.kernel_parameters = parameters 428 self._check_angles()429 429 self._set_vector_lengths() 430 430 … … 476 476 self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 477 477 478 def _check_angles(self): 478 def check_angles(self): 479 """ 480 Check that orientation angles are theta, phi and possibly psi. 481 """ 479 482 theta = phi = psi = -1 480 483 for k, p in enumerate(self.kernel_parameters): … … 499 502 if psi >= 0 and psi != phi+1: 500 503 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") 504 507 elif theta >= 0 or phi >= 0 or psi >= 0: 505 508 raise TypeError("oriented shapes must have both theta and phi and maybe psi") -
sasmodels/models/core_shell_bicelle_elliptical.py
r2d81cfe rfc3ae1b 9 9 scattering length densities. The form factor is normalized by the total 10 10 particle volume. 11 12 11 13 12 .. figure:: img/core_shell_bicelle_geometry.png … … 91 90 for further information. 92 91 93 94 92 .. figure:: img/elliptical_cylinder_angle_definition.png 95 93 96 94 Definition of the angles for the oriented core_shell_bicelle_elliptical particles. 97 95 98 96 Model verified using Monte Carlo simulation for 1D and 2D scattering. 99 97 100 98 References … … 108 106 * **Author:** Richard Heenan **Date:** December 14, 2016 109 107 * **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, 2016108 * **Last Reviewed by:** Paul Kienzle **Date:** Feb 28, 2018 111 109 """ 112 110 … … 130 128 # ["name", "units", default, [lower, upper], "type", "description"], 131 129 parameters = [ 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"], 134 132 ["thick_rim", "Ang", 8, [0, inf], "volume", "Rim shell thickness"], 135 133 ["thick_face", "Ang", 14, [0, inf], "volume", "Cylinder face thickness"], … … 139 137 ["sld_rim", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Cylinder rim scattering length density"], 140 138 ["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"] 144 142 ] 145 143 -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py
r108e70e rfc3ae1b 140 140 # ["name", "units", default, [lower, upper], "type", "description"], 141 141 parameters = [ 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"], 144 144 ["thick_rim", "Ang", 8, [0, inf], "volume", "Rim or belt shell thickness"], 145 145 ["thick_face", "Ang", 14, [0, inf], "volume", "Cylinder face thickness"], … … 149 149 ["sld_rim", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Cylinder rim scattering length density"], 150 150 ["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"], 155 155 ] 156 156 -
sasmodels/models/core_shell_cylinder.c
r108e70e rd86f0fc 48 48 49 49 50 double Iqac(double qab, double qc, 50 static double 51 Iqac(double qab, double qc, 51 52 double core_sld, 52 53 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) 1 static double 2 form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 6 3 { 7 4 double length_b = length_a * b2a_ratio; … … 16 13 } 17 14 18 double Iq(double q, 15 static double 16 Iq(double q, 19 17 double sld, 20 18 double solvent_sld, … … 85 83 } 86 84 87 double Iqabc(double qa, double qb, double qc, 85 static double 86 Iqabc(double qa, double qb, double qc, 88 87 double sld, 89 88 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) 1 static double 2 form_volume(double length_a, double b2a_ratio, double c2a_ratio) 6 3 { 7 4 double length_b = length_a * b2a_ratio; … … 11 8 } 12 9 13 double Iq(double q, 10 static double 11 Iq(double q, 14 12 double sld, 15 13 double solvent_sld, -
sasmodels/models/polymer_excl_volume.py
r2d81cfe raa90015 17 17 $a$ is the statistical segment length of the polymer chain, 18 18 and $n$ is the degree of polymerization. 19 This integral was later put into an almost analytical form as follows 19 20 This integral was put into an almost analytical form as follows 20 21 (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 31 and later recast as (for example, Hore, 2013; Hammouda & Kim, 2017) 21 32 22 33 .. math:: … … 29 40 .. math:: 30 41 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} 32 43 33 44 and the variable $U$ is given in terms of the scattering vector $Q$ as … … 37 48 U=\frac{Q^2a^2n^{2\nu}}{6} = \frac{Q^2R_{g}^2(2\nu+1)(2\nu+2)}{6} 38 49 50 The two analytic forms are equivalent. In the 1993 paper 51 52 .. math:: 53 54 \frac{1}{\nu U^{1/2\nu}} 55 56 has been factored out. 57 58 **SasView implements the 1993 expression**. 59 39 60 The square of the radius-of-gyration is defined as 40 61 … … 43 64 R_{g}^2 = \frac{a^2n^{2\nu}}{(2\nu+1)(2\nu+2)} 44 65 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$ ). 50 72 51 73 A low-Q expansion yields the Guinier form and a high-Q expansion yields the … … 73 95 .. math:: 74 96 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] 76 98 77 99 For 2D data: The 2D scattering intensity is calculated in the same way as 1D, … … 89 111 90 112 B Hammouda, *SANS from Homogeneous Polymer Mixtures - A Unified Overview, 91 Advances in Polym. Sci.* 106(1993) 87-133 113 Advances in Polym. Sci.* 106 (1993) 87-133 114 115 M Hore et al, *Co-Nonsolvency of Poly(n-isopropylacrylamide) in Deuterated 116 Water/Ethanol Mixtures* 46 (2013) 7894-7901 117 118 B Hammouda & M-H Kim, *The empirical core-chain model* 247 (2017) 434-440 92 119 """ 93 120 … … 124 151 with errstate(divide='ignore', invalid='ignore'): 125 152 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). 126 155 result = (porod_exp*upow * 127 156 (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) 1 static double 2 form_volume(double length_a, double b2a_ratio, double c2a_ratio) 6 3 { 7 4 return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio); 8 5 } 9 6 10 double Iq(double q, 7 static double 8 Iq(double q, 11 9 double sld, 12 10 double solvent_sld, … … 71 69 72 70 73 double Iqabc(double qa, double qb, double qc, 71 static double 72 Iqabc(double qa, double qb, double qc, 74 73 double sld, 75 74 double solvent_sld, -
sasmodels/rst2html.py
r2d81cfe r1fbadb2 56 56 # others don't work properly with math_output! 57 57 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} 59 61 else: 60 62 settings = {"math-output": math_output} -
sasmodels/sasview_model.py
r2d81cfe rd533590 593 593 # Check whether we have a list of ndarrays [qx,qy] 594 594 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) 599 596 600 597 elif isinstance(qdist, np.ndarray): … … 665 662 666 663 def _calculate_Iq(self, qx, qy=None): 667 #core.HAVE_OPENCL = False668 664 if self._model is None: 669 665 self._model = core.build_model(self._model_info) … … 678 674 call_details, values, is_magnetic = make_kernel_args(calculator, pairs) 679 675 #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]) 681 678 #for k, p in enumerate(self._model_info.parameters.call_parameters): 682 679 # print(k, p.name, *pairs[k]) … … 689 686 self._intermediate_results = getattr(calculator, 'results', None) 690 687 calculator.release() 691 self._model.release()688 #self._model.release() 692 689 return result 693 690 … … 722 719 723 720 def set_dispersion(self, parameter, dispersion): 724 # type: (str, weights.Dispersion) -> Dict[str, Any]721 # type: (str, weights.Dispersion) -> None 725 722 """ 726 723 Set the dispersion object for a model parameter … … 745 742 self.dispersion[parameter] = dispersion.get_pars() 746 743 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) 748 746 749 747 def _dispersion_mesh(self): … … 859 857 # type: () -> None 860 858 """ 861 Load and run cylinder model from sas.models.CylinderModel859 Load and run cylinder model as sas-models-CylinderModel 862 860 """ 863 861 if not SUPPORT_OLD_STYLE_PLUGINS: … … 872 870 CylinderModel().evalDistribution([0.1, 0.1]) 873 871 872 def 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 874 885 if __name__ == "__main__": 875 886 print("cylinder(0.1,0.1)=%g"%test_cylinder()) 887 #magnetic_demo() 876 888 #test_product() 877 889 #test_structure_factor() -
sasmodels/weights.py
r2d81cfe r3d58247 97 97 return x, px 98 98 99 class 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) 99 113 100 114 class RectangleDispersion(Dispersion): … … 107 121 """ 108 122 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) 115 128 116 129 class LogNormalDispersion(Dispersion): … … 190 203 return x, px 191 204 205 class 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 192 219 193 220 # dispersion name -> disperser lookup table. … … 196 223 MODELS = OrderedDict((d.type, d) for d in ( 197 224 RectangleDispersion, 225 UniformDispersion, 198 226 ArrayDispersion, 199 227 LogNormalDispersion, 200 228 GaussianDispersion, 201 229 SchulzDispersion, 230 BoltzmannDispersion 202 231 )) 203 232 -
setup.py
r32e3c9b r1f991d6 1 try: 2 from setuptools import setup 3 except ImportError: 4 from distutils.core import setup 1 import sys 2 from setuptools import setup 3 from setuptools.command.test import test as TestCommand 4 5 class 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) 5 18 6 19 def find_version(package): … … 48 61 'sasmodels': ['*.c', '*.cl'], 49 62 }, 50 install_requires =[63 install_requires=[ 51 64 ], 52 extras_require ={65 extras_require={ 53 66 'OpenCL': ["pyopencl"], 54 67 'Bumps': ["bumps"], 55 68 'TinyCC': ["tinycc"], 56 69 }, 70 build_requires=['setuptools'], 71 test_requires=['pytest'], 72 cmdclass = {'test': PyTest}, 57 73 ) -
.gitignore
re9ed2de r6ceca44 8 8 *.so 9 9 *.obj 10 *.o 10 11 /doc/_build/ 11 12 /doc/api/ … … 19 20 /.pydevproject 20 21 /.idea 22 .vscode 21 23 /sasmodels.egg-info/ 22 24 /example/Fit_*/
Note: See TracChangeset
for help on using the changeset viewer.