Changes in / [7a516d0:49284e1] in sasmodels


Ignore:
Files:
4 added
2 deleted
50 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    r8a5f021 re9ed2de  
    2222/example/Fit_*/ 
    2323/example/batch_fit.csv 
     24# Note: changes to gauss20, gauss76 and gauss150 are still tracked since they 
     25# are part of the repo and required by some models. 
     26/sasmodels/models/lib/gauss*.c 

  • .travis.yml

    rce8c388 r2d09df1  
    66    env: 
    77    - PY=2.7 
     8    - DEPLOY=True 
    89  - os: linux 
    910    env: 
    … …  
    5152  on: 
    5253    branch: master 
     54    condition: $DEPLOY = True 
    5355notifications: 
    5456  slack: 

  • appveyor.yml

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

  • doc/developer/overview.rst

    r3d40839 r0d5a655  
    185185jitter applied before particle orientation. 
    186186 
     187When computing the orientation dispersity integral, the weights for 
     188the individual points depends on the map projection used to translate jitter 
     189angles into latitude/longitude.  The choice of projection is set by 
     190*sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for 
     191equirectangular and *PROJECTION=2* for sinusoidal.  The more complicated 
     192Guyou and Postel projections are not implemented. See jitter.draw_mesh 
     193for details. 
     194 
    187195For numerical integration within form factors etc. sasmodels is mostly using 
    188196Gaussian quadrature with 20, 76 or 150 points depending on the model. It also 
    … …  
    199207Useful testing routines include: 
    200208 
    201 :mod:`asymint` a direct implementation of the surface integral for certain 
    202 models to get a more trusted value for the 1D integral using a 
    203 reimplementation of the 2D kernel in python and mpmath (which computes math 
    204 functions to arbitrary precision). It uses $\theta$ ranging from 0 to $\pi$ 
    205 and $\phi$ ranging from 0 to $2\pi$. It perhaps would benefit from including 
    206 the U-substitution for $\theta$. 
    207  
    208 :mod:`check1d` uses sasmodels 1D integration and compares that with a 
     209The *sascomp* utility is used to view and compare models with different 
     210parameters and calculation engines. The usual case is to simply plot a 
     211model that you are developing:: 
     212 
     213    python sascomp path/to/model.py 
     214 
     215Once the obvious problems are addressed, check the numerical precision 
     216across a variety of randomly generated inputs:: 
     217 
     218    python sascomp -engine=single,double path/to/model.py -sets=10 
     219 
     220You can compare different parameter values for the same or different models. 
     221For example when looking along the long axis of a cylinder ($\theta=0$), 
     222dispersity in $\theta$ should be equivalent to dispersity in $\phi$ 
     223when $\phi=90$:: 
     224 
     225    python sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle \\ 
     226    phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10 
     227 
     228It turns out that they are not because the equirectangular map projection 
     229weights the points by $\cos(\theta)$ so $\Delta\theta$ is not identical 
     230to $\Delta\phi$.  Setting *PROJECTION=2* in :mod:`sasmodels.generate` helps 
     231somewhat.  Postel would help even more in this case, though leading 
     232to distortions elsewhere.  See :mod:`sasmodels.compare` for many more details. 
     233 
     234*sascomp -ngauss=n* allows you to set the number of quadrature points used 
     235for the 1D integration for any model.  For example, a carbon nanotube with 
     236length 10 $\mu$\ m and radius 1 nm is not computed correctly at high $q$:: 
     237 
     238    python sascomp cylinder length=100000 radius=10 -ngauss=76,500 -double -highq 
     239 
     240Note: ticket 702 gives some forms for long rods and thin disks that may 
     241be used in place of the integration, depending on $q$, radius and length; if 
     242the cylinder model is updated with these corrections then above call will show 
     243no difference. 
     244 
     245*explore/check1d.py* uses sasmodels 1D integration and compares that with a 
    209246rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to 
    210247$\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle 
    … …  
    214251similar reasoning.] This should rotate the sample through the entire 
    215252$\theta$-$\phi$ surface according to the pattern that you see in jitter.py when 
    216 you modify it to use 'rectangle' rather than 'gaussian' for its distribution 
    217 without changing the viewing angle. In order to match the 1-D pattern for 
    218 an arbitrary viewing angle on triaxial shapes, we need to integrate 
     253you use 'rectangle' rather than 'gaussian' for its distribution without 
     254changing the viewing angle. In order to match the 1-D pattern for an arbitrary 
     255viewing angle on triaxial shapes, we need to integrate 
    219256over $\theta$, $\phi$ and $\Psi$. 
    220257 
    221 When computing the dispersity integral, weights are scaled by 
    222 $|\cos(\delta \theta)|$ to account for the points in $\phi$ getting closer 
    223 together as $\delta \theta$ increases. 
    224 [This will probably change so that instead of adjusting the weights, we will 
    225 adjust $\delta\theta$-$\delta\phi$ mesh so that the point density in 
    226 $\delta\phi$ is lower at larger $\delta\theta$. The flag USE_SCALED_PHI in 
    227 *kernel_iq.c* selects an alternative algorithm.] 
    228  
    229 The integrated dispersion is computed at a set of $(qx, qy)$ points $(q 
    230 \cos(\alpha), q \sin(\alpha))$ at some angle $\alpha$ (currently angle=0) for 
    231 each $q$ used in the 1-D integration. The individual $q$ points should be 
    232 equivalent to asymint rect-n when the viewing angle is set to 
    233 $(\theta,\phi,\Psi) = (90, 0, 0)$. Such tests can help to validate that 2d 
    234 models are consistent with 1d models. 
    235  
    236 :mod:`sascomp -sphere=n` uses the same rectangular distribution as check1d to 
    237 compute the pattern of the $q_x$-$q_y$ grid. 
    238  
    239 The :mod:`sascomp` utility can be used for 2d as well as 1d calculations to 
    240 compare results for two sets of parameters or processor types, for example 
    241 these two oriented cylinders here should be equivalent. 
    242  
    243 :mod:`\./sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10` 
    244  
     258*sascomp -sphere=n* uses the same rectangular distribution as check1d to 
     259compute the pattern of the $q_x$-$q_y$ grid.  This ought to produce a 
     260spherically symmetric pattern. 
     261 
     262*explore/precision.py* investigates the accuracy of individual functions 
     263on the different execution platforms.  This includes the basic special 
     264functions as well as more complex expressions used in scattering.  In many 
     265cases the OpenCL function in sasmodels will use a polynomial approximation 
     266over part of the range to improve accuracy, as shown in:: 
     267 
     268    python explore/precision.py 3j1/x 
     269 
     270*explore/realspace.py* allows you to set up a Monte Carlo simulation of your 
     271model by sampling random points within and computing the 1D and 2D scattering 
     272patterns.  This was used to check the core shell parallelepiped example.  This 
     273is similar to the general sas calculator in sasview, though it uses different 
     274code. 
     275 
     276*sasmodels/jitter.py* is for exploring different options for handling 
     277orientation and orientation dispersity.  It uses *sasmodels/guyou.py* to 
     278generate the Guyou projection. 
     279 
     280*explore/asymint.py* is a direct implementation of the 1D integration for 
     281a number of different models.  It calculates the integral for a particular 
     282$q$ using several different integration schemes, including mpmath with 100 
     283bits of precision (33 digits), so you can use it to check the target values 
     284for the 1D model tests.  This is not a general purpose tool; you will need to 
     285edit the file to change the model parameters. It does not currently 
     286apply the $u=cos(\theta)$ substitution that many models are using 
     287internally so the 76-point gaussian quadrature may not match the value 
     288produced by the eqivalent model from sasmodels. 
     289 
     290*explore/symint.py* is for rotationally symmetric models (just cylinder for 
     291now), so it can compute an entire curve rather than a single $q$ point.  It 
     292includes code to compare the long cylinder approximation to cylinder. 
     293 
     294*explore/rpa.py* is for checking different implementations of the RPA model 
     295against calculations for specific blends.  This is a work in (suspended) 
     296progress. 
     297 
     298There are a few modules left over in *explore* that can probably be removed. 
     299*angular_pd.py* was an early version of *jitter.py*.  *sc.py* and *sc.c* 
     300was used to explore different calculations for paracrystal models; it has 
     301been absorbed into *asymint.py*. *transform_angles.py* translates orientation 
     302parameters from the SasView 3.x definition to sasmodels. 
     303 
     304*explore/angles.py* generates code for the view and jitter transformations. 
     305Keep this around since it may be needed if we add new projections. 
    245306 
    246307Testing 

  • doc/genapi.py

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

  • doc/guide/orientation/orientation.rst

    r82592da r0d5a655  
    44================== 
    55 
    6 With two dimensional small angle diffraction data SasView will calculate 
     6With two dimensional small angle diffraction data sasmodels will calculate 
    77scattering from oriented particles, applicable for example to shear flow 
    88or orientation in a magnetic field. 
    99 
    1010In general we first need to define the reference orientation 
    11 of the particles with respect to the incoming neutron or X-ray beam. This 
    12 is done using three angles: $\theta$ and $\phi$ define the orientation of 
    13 the axis of the particle, angle $\Psi$ is defined as the orientation of 
    14 the major axis of the particle cross section with respect to its starting 
    15 position along the beam direction. The figures below are for an elliptical 
    16 cross section cylinder, but may be applied analogously to other shapes of 
    17 particle. 
     11of the particle's $a$-$b$-$c$ axes with respect to the incoming 
     12neutron or X-ray beam. This is done using three angles: $\theta$ and $\phi$ 
     13define the orientation of the $c$-axis of the particle, and angle $\Psi$ is 
     14defined as the orientation of the major axis of the particle cross section 
     15with respect to its starting position along the beam direction (or 
     16equivalently, as rotation about the $c$ axis). There is an unavoidable 
     17ambiguity when $c$ is aligned with $z$ in that $\phi$ and $\Psi$ both 
     18serve to rotate the particle about $c$, but this symmetry is destroyed 
     19when $\theta$ is not a multiple of 180. 
     20 
     21The figures below are for an elliptical cross section cylinder, but may 
     22be applied analogously to other shapes of particle. 
    1823 
    1924.. note:: 
    … …  
    2934 
    3035    Definition of angles for oriented elliptical cylinder, where axis_ratio 
    31     b/a is shown >1, Note that rotation $\theta$, initially in the $x$-$z$ 
     36    b/a is shown >1. Note that rotation $\theta$, initially in the $x$-$z$ 
    3237    plane, is carried out first, then rotation $\phi$ about the $z$-axis, 
    3338    finally rotation $\Psi$ is around the axis of the cylinder. The neutron 
    34     or X-ray beam is along the $z$ axis. 
     39    or X-ray beam is along the $-z$ axis. 
    3540 
    3641.. figure:: 
    … …  
    4045    with $\Psi$ = 0. 
    4146 
    42 Having established the mean direction of the particle we can then apply 
    43 angular orientation distributions. This is done by a numerical integration 
    44 over a range of angles in a similar way to particle size dispersity. 
    45 In the current version of sasview the orientational dispersity is defined 
    46 with respect to the axes of the particle. 
     47Having established the mean direction of the particle (the view) we can then 
     48apply angular orientation distributions (jitter). This is done by a numerical 
     49integration over a range of angles in a similar way to particle size 
     50dispersity. The orientation dispersity is defined with respect to the 
     51$a$-$b$-$c$ axes of the particle, with roll angle $\Psi$ about the $c$-axis, 
     52yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the 
     53$a$-axis. 
     54 
     55You can explore the view and jitter angles interactively using 
     56:func:`sasmodels.jitter.run`.  Enter the following into the python 
     57interpreter:: 
     58 
     59    from sasmodels import jitter 
     60    jitter.run() 
     61 
     62More formally, starting with axes $a$-$b$-$c$ of the particle aligned 
     63with axes $x$-$y$-$z$ of the laboratory frame, the orientation dispersity 
     64is applied first, using the 
     65`Tait-Bryan <https://en.wikipedia.org/wiki/Euler_angles#Conventions_2>`_ 
     66$x$-$y'$-$z''$ convention with angles $\Delta\phi$-$\Delta\theta$-$\Delta\Psi$. 
     67The reference orientation then follows, using the 
     68`Euler angles <https://en.wikipedia.org/wiki/Euler_angles#Conventions>`_ 
     69$z$-$y'$-$z''$ with angles $\phi$-$\theta$-$\Psi$.  This is implemented 
     70using rotation matrices as 
     71 
     72.. math:: 
     73 
     74    R = R_z(\phi)\, R_y(\theta)\, R_z(\Psi)\, 
     75        R_x(\Delta\phi)\, R_y(\Delta\theta)\, R_z(\Delta\Psi) 
     76 
     77To transform detector $(q_x, q_y)$ values into $(q_a, q_b, q_c)$ for the 
     78shape in its canonical orientation, use 
     79 
     80.. math:: 
     81 
     82    [q_a, q_b, q_c]^T = R^{-1} \, [q_x, q_y, 0]^T 
     83 
     84 
     85The inverse rotation is easily calculated by rotating the opposite directions 
     86in the reverse order, so 
     87 
     88.. math:: 
     89 
     90    R^{-1} = R_z(-\Delta\Psi)\, R_y(-\Delta\theta)\, R_x(-\Delta\phi)\, 
     91             R_z(-\Psi)\, R_y(-\theta)\, R_z(-\phi) 
     92 
    4793 
    4894The $\theta$ and $\phi$ orientation parameters for the cylinder only appear 
    49 when fitting 2d data. On introducing "Orientational Distribution" in 
    50 the angles, "distribution of theta" and "distribution of phi" parameters will 
     95when fitting 2d data. On introducing "Orientational Distribution" in the 
     96angles, "distribution of theta" and "distribution of phi" parameters will 
    5197appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ 
    52 of the cylinder, the $b$ and $a$ axes of the cylinder cross section. (When 
    53 $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the 
    54 instrument.) The third orientation distribution, in $\Psi$, is about the $c$ 
    55 axis of the particle. Some experimentation may be required to understand the 
    56 2d patterns fully. A number of different shapes of distribution are 
    57 available, as described for polydispersity, see :ref:`polydispersityhelp` . 
    58  
    59 Earlier versions of SasView had numerical integration issues in some 
    60 circumstances when distributions passed through 90 degrees. The distributions 
    61 in particle coordinates are more robust, but should still be approached with 
    62 care for large ranges of angle. 
     98of the cylinder, which correspond to the $b$ and $a$ axes of the cylinder 
     99cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and 
     100$X$ axes of the instrument.) The third orientation distribution, in $\Psi$, 
     101is about the $c$ axis of the particle. Some experimentation may be required 
     102to understand the 2d patterns fully. A number of different shapes of 
     103distribution are available, as described for size dispersity, see 
     104:ref:`polydispersityhelp`. 
     105 
     106Given that the angular dispersion distribution is defined in cartesian space, 
     107over a cube defined by 
     108 
     109.. math:: 
     110 
     111    [-\Delta \theta, \Delta \theta] \times 
     112    [-\Delta \phi, \Delta \phi] \times 
     113    [-\Delta \Psi, \Delta \Psi] 
     114 
     115but the orientation is defined over a sphere, we are left with a 
     116`map projection <https://en.wikipedia.org/wiki/List_of_map_projections>`_ 
     117problem, with different tradeoffs depending on how values in $\Delta\theta$ 
     118and $\Delta\phi$ are translated into latitude/longitude on the sphere. 
     119 
     120Sasmodels is using the 
     121`equirectangular projection <https://en.wikipedia.org/wiki/Equirectangular_projection>`_. 
     122In this projection, square patches in angular dispersity become wedge-shaped 
     123patches on the sphere. To correct for the changing point density, there is a 
     124scale factor of $\sin(\Delta\theta)$ that applies to each point in the 
     125integral. This is not enough, though. Consider a shape which is tumbling 
     126freely around the $b$ axis, with $\Delta\theta$ uniform in $[-180, 180]$. At 
     127$\pm 90$, all points in $\Delta\phi$ map to the pole, so the jitter will have 
     128a distinct angular preference. If the spin axis is along the beam (which 
     129will be the case for $\theta=90$ and $\Psi=90$) the scattering pattern 
     130should be circularly symmetric, but it will go to zero at $q_x = 0$ due to the 
     131$\sin(\Delta\theta)$ correction. This problem does not appear for a shape 
     132that is tumbling freely around the $a$ axis, with $\Delta\phi$ uniform in 
     133$[-180, 180]$, so swap the $a$ and $b$ axes so $\Delta\theta < \Delta\phi$ 
     134and adjust $\Psi$ by 90. This works with the current sasmodels shapes due to 
     135symmetry. 
     136 
     137Alternative projections were considered. 
     138The `sinusoidal projection <https://en.wikipedia.org/wiki/Sinusoidal_projection>`_ 
     139works by scaling $\Delta\phi$ as $\Delta\theta$ increases, and dropping those 
     140points outside $[-180, 180]$. The distortions are a little less for middle 
     141ranges of $\Delta\theta$, but they are still severe for large $\Delta\theta$ 
     142and the model is much harder to explain. 
     143The `azimuthal equidistance projection <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>`_ 
     144also improves on the equirectangular projection by extending the range of 
     145reasonable values for the $\Delta\theta$ range, with $\Delta\phi$ forming a 
     146wedge that cuts to the opposite side of the sphere rather than cutting to the 
     147pole. This projection has the nice property that distance from the center are 
     148preserved, and that $\Delta\theta$ and $\Delta\phi$ act the same. 
     149The `azimuthal equal area projection <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_ 
     150is like the azimuthal equidistance projection, but it preserves area instead 
     151of distance. It also has the same behaviour for $\Delta\theta$ and $\Delta\phi$. 
     152The `Guyou projection <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>`_ 
     153has an excellent balance with reasonable distortion in both $\Delta\theta$ 
     154and $\Delta\phi$, as well as preserving small patches. However, it requires 
     155considerably more computational overhead, and we have not yet derived the 
     156formula for the distortion correction, measuring the degree of stretch at 
     157the point $(\Delta\theta, \Delta\phi)$ on the map. 
    63158 
    64159.. note:: 
    65     Note that the form factors for oriented particles are also performing 
    66     numerical integrations over one or more variables, so care should be taken, 
    67     especially with very large particles or more extreme aspect ratios. In such  
    68     cases results may not be accurate, particularly at very high Q, unless the model 
    69     has been specifically coded to use limiting forms of the scattering equations. 
    70      
    71     For best numerical results keep the $\theta$ distribution narrower than the $\phi$  
    72     distribution. Thus for asymmetric particles, such as elliptical_cylinder, you may  
    73     need to reorder the sizes of the three axes to acheive the desired result.  
    74     This is due to the issues of mapping a rectangular distribution onto the  
    75     surface of a sphere. 
    76  
    77 Users can experiment with the values of *Npts* and *Nsigs*, the number of steps  
    78 used in the integration and the range spanned in number of standard deviations.  
    79 The standard deviation is entered in units of degrees. For a "rectangular"  
    80 distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations.  
    81 The new "uniform" distribution avoids this by letting you directly specify the  
     160    Note that the form factors for oriented particles are performing 
     161    numerical integrations over one or more variables, so care should be 
     162    taken, especially with very large particles or more extreme aspect 
     163    ratios. In such cases results may not be accurate, particularly at very 
     164    high Q, unless the model has been specifically coded to use limiting 
     165    forms of the scattering equations. 
     166 
     167    For best numerical results keep the $\theta$ distribution narrower than 
     168    the $\phi$ distribution. Thus for asymmetric particles, such as 
     169    elliptical_cylinder, you may need to reorder the sizes of the three axes 
     170    to acheive the desired result. This is due to the issues of mapping a 
     171    rectanglar distribution onto the surface of a sphere. 
     172 
     173Users can experiment with the values of *Npts* and *Nsigs*, the number of steps 
     174used in the integration and the range spanned in number of standard deviations. 
     175The standard deviation is entered in units of degrees. For a "rectangular" 
     176distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations. 
     177The new "uniform" distribution avoids this by letting you directly specify the 
    82178half width. 
    83179 
    84 The angular distributions will be truncated outside of the range -180 to +180  
    85 degrees, so beware of using saying a broad Gaussian distribution with large value 
    86 of *Nsigs*, as the array of *Npts* may be truncated to many fewer points than would  
    87 give a good integration,as well as becoming rather meaningless. (At some point  
    88 in the future the actual polydispersity arrays may be made available to the user  
    89 for inspection.) 
     180The angular distributions may be truncated outside of the range -180 to +180 
     181degrees, so beware of using saying a broad Gaussian distribution with large 
     182value of *Nsigs*, as the array of *Npts* may be truncated to many fewer 
     183points than would give a good integration,as well as becoming rather 
     184meaningless. (At some point in the future the actual dispersion arrays may be 
     185made available to the user for inspection.) 
    90186 
    91187Some more detailed technical notes are provided in the developer section of 
    92188this manual :ref:`orientation_developer` . 
    93189 
     190This definition of orientation is new to SasView 4.2.  In earlier versions, 
     191the orientation distribution appeared as a distribution of view angles. 
     192This led to strange effects when $c$ was aligned with $z$, where changes 
     193to the $\phi$ angle served only to rotate the shape about $c$, rather than 
     194having a consistent interpretation as the pitch of the shape relative to 
     195the flow field defining the reference orientation.  Prior to SasView 4.1, 
     196the reference orientation was defined using a Tait-Bryan convention, making 
     197it difficult to control.  Now, rotation in $\theta$ modifies the spacings 
     198in the refraction pattern, and rotation in $\phi$ rotates it in the detector 
     199plane. 
     200 
     201 
    94202*Document History* 
    95203 
    96204| 2017-11-06 Richard Heenan 
     205| 2017-12-20 Paul Kienzle 

  • doc/guide/plugin.rst

    r3048ec6 r7e6bc45e  
    292292**Note: The order of the parameters in the definition will be the order of the 
    293293parameters in the user interface and the order of the parameters in Iq(), 
    294 Iqxy() and form_volume(). And** *scale* **and** *background* **parameters are 
    295 implicit to all models, so they do not need to be included in the parameter table.** 
     294Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background* 
     295**parameters are implicit to all models, so they do not need to be included 
     296in the parameter table.** 
    296297 
    297298- **"name"** is the name of the parameter shown on the FitPage. 
    … …  
    362363    scattered intensity. 
    363364 
    364   - "volume" parameters are passed to Iq(), Iqxy(), and form_volume(), and 
    365     have polydispersity loops generated automatically. 
    366  
    367   - "orientation" parameters are only passed to Iqxy(), and have angular 
    368     dispersion. 
     365  - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(), 
     366    and have polydispersity loops generated automatically. 
     367 
     368  - "orientation" parameters are not passed, but instead are combined with 
     369    orientation dispersity to translate *qx* and *qy* to *qa*, *qb* and *qc*. 
     370    These parameters should appear at the end of the table with the specific 
     371    names *theta*, *phi* and for asymmetric shapes *psi*, in that order. 
    369372 
    370373Some models will have integer parameters, such as number of pearls in the 
    … …  
    419422That is, the individual models do not need to include polydispersity 
    420423calculations, but instead rely on numerical integration to compute the 
    421 appropriately smeared pattern.   Angular dispersion values over polar angle 
    422 $\theta$ requires an additional $\cos \theta$ weighting due to decreased 
    423 arc length for the equatorial angle $\phi$ with increasing latitude. 
     424appropriately smeared pattern. 
    424425 
    425426Python Models 
    … …  
    468469barbell).  If I(q; pars) is NaN for any $q$, then those parameters will be 
    469470ignored, and not included in the calculation of the weighted polydispersity. 
    470  
    471 Similar to *Iq*, you can define *Iqxy(qx, qy, par1, par2, ...)* where the 
    472 parameter list includes any orientation parameters.  If *Iqxy* is not defined, 
    473 then it will default to *Iqxy = Iq(sqrt(qx**2+qy**2), par1, par2, ...)*. 
    474471 
    475472Models should define *form_volume(par1, par2, ...)* where the parameter 
    … …  
    497494    } 
    498495 
    499 *Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*, 
    500 and it includes orientation parameters. 
    501  
    502496*form_volume* defines the volume of the shape. As in python models, it 
    503497includes only the volume parameters. 
    504498 
    505 *Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and 
    506 *form_volume* will default to 1.0. 
    507  
    508499**source=['fn.c', ...]** includes the listed C source files in the 
    509 program before *Iq* and *Iqxy* are defined. This allows you to extend the 
    510 library of C functions available to your model. 
     500program before *Iq* and *form_volume* are defined. This allows you to 
     501extend the library of C functions available to your model. 
     502 
     503*c_code* includes arbitrary C code into your kernel, which can be 
     504handy for defining helper functions for *Iq* and *form_volume*. Note that 
     505you can put the full function definition for *Iq* and *form_volume* 
     506(include function declaration) into *c_code* as well, or put them into an 
     507external C file and add that file to the list of sources. 
    511508 
    512509Models are defined using double precision declarations for the 
    … …  
    533530    #define INVALID(v) (v.bell_radius < v.radius) 
    534531 
     532The INVALID define can go into *Iq*, or *c_code*, or an external C file 
     533listed in *source*. 
     534 
     535Oriented Shapes 
     536............... 
     537 
     538If the scattering is dependent on the orientation of the shape, then you 
     539will need to include *orientation* parameters *theta*, *phi* and *psi* 
     540at the end of the parameter table.  As described in the section 
     541:ref:`orientation`, the individual $(q_x, q_y)$ points on the detector will 
     542be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its 
     543canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the 
     544laboratory frame and beam travelling along $-z$. 
     545 
     546The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
     547*par1*, etc. are the parameters to the model.  If the shape is rotationally 
     548symmetric about *c* then *psi* is not needed, and the model is called 
     549as *Iqac(qab, qc, par1, par2, ...)*.  In either case, the orientation 
     550parameters are not included in the function call. 
     551 
     552For 1D oriented shapes, an integral over all angles is usually needed for 
     553the *Iq* function. Given symmetry and the substitution $u = \cos(\alpha)$, 
     554$du = -\sin(\alpha)\,d\alpha$ this becomes 
     555 
     556.. math:: 
     557 
     558    I(q) &= \frac{1}{4\pi} \int_{-\pi/2}^{pi/2} \int_{-pi}^{pi} 
     559            F(q_a, q_b, q_c)^2 \sin(\alpha)\,d\beta\,d\alpha \\ 
     560        &= \frac{8}{4\pi} \int_{0}^{pi/2} \int_{0}^{\pi/2} 
     561            F^2 \sin(\alpha)\,d\beta\,d\alpha \\ 
     562        &= \frac{8}{4\pi} \int_1^0 \int_{0}^{\pi/2} - F^2 \,d\beta\,du \\ 
     563        &= \frac{8}{4\pi} \int_0^1 \int_{0}^{\pi/2} F^2 \,d\beta\,du 
     564 
     565for 
     566 
     567.. math:: 
     568 
     569    q_a &= q \sin(\alpha)\sin(\beta) = q \sqrt{1-u^2} \sin(\beta) \\ 
     570    q_b &= q \sin(\alpha)\cos(\beta) = q \sqrt{1-u^2} \cos(\beta) \\ 
     571    q_c &= q \cos(\alpha) = q u 
     572 
     573Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the 
     574numerical integration is then:: 
     575 
     576    double outer_sum = 0.0; 
     577    for (int i = 0; i < GAUSS_N; i++) { 
     578        const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5; 
     579        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
     580        const double qc = cos_alpha * q; 
     581        double inner_sum = 0.0; 
     582        for (int j = 0; j < GAUSS_N; j++) { 
     583            const double beta = M_PI_4 * GAUSS_Z[j] + M_PI_4; 
     584            double sin_beta, cos_beta; 
     585            SINCOS(beta, sin_beta, cos_beta); 
     586            const double qa = sin_alpha * sin_beta * q; 
     587            const double qb = sin_alpha * cos_beta * q; 
     588            const double form = Fq(qa, qb, qc, ...); 
     589            inner_sum += GAUSS_W[j] * form * form; 
     590        } 
     591        outer_sum += GAUSS_W[i] * inner_sum; 
     592    } 
     593    outer_sum *= 0.25; // = 8/(4 pi) * outer_sum * (pi/2) / 4 
     594 
     595The *z* values for the Gauss-Legendre integration extends from -1 to 1, so 
     596the double sum of *w[i]w[j]* explains the factor of 4.  Correcting for the 
     597average *dz[i]dz[j]* gives $(1-0) \cdot (\pi/2-0) = \pi/2$.  The $8/(4 \pi)$ 
     598factor comes from the integral over the quadrant.  With less symmetry (eg., 
     599in the bcc and fcc paracrystal models), then an integral over the entire 
     600sphere may be necessary. 
     601 
     602For simpler models which are rotationally symmetric a single integral 
     603suffices: 
     604 
     605.. math:: 
     606 
     607    I(q) &= \frac{1}{\pi}\int_{-\pi/2}^{\pi/2} 
     608            F(q_{ab}, q_c)^2 \sin(\alpha)\,d\alpha/\pi \\ 
     609        &= \frac{2}{\pi} \int_0^1 F^2\,du 
     610 
     611for 
     612 
     613.. math:: 
     614 
     615    q_{ab} &= q \sin(\alpha) = q \sqrt{1 - u^2} \\ 
     616    q_c &= q \cos(\alpha) = q u 
     617 
     618 
     619with integration loop:: 
     620 
     621    double sum = 0.0; 
     622    for (int i = 0; i < GAUSS_N; i++) { 
     623        const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5; 
     624        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
     625        const double qab = sin_alpha * q; 
     626        const double qc = cos_alpha * q; 
     627        const double form = Fq(qab, qc, ...); 
     628        sum += GAUSS_W[j] * form * form; 
     629    } 
     630    sum *= 0.5; // = 2/pi * sum * (pi/2) / 2 
     631 
     632Magnetism 
     633......... 
     634 
     635Magnetism is supported automatically for all shapes by modifying the 
     636effective SLD of particle according to the Halpern-Johnson vector 
     637describing the interaction between neutron spin and magnetic field.  All 
     638parameters marked as type *sld* in the parameter table are treated as 
     639possibly magnetic particles with magnitude *M0* and direction 
     640*mtheta* and *mphi*.  Polarization parameters are also provided 
     641automatically for magnetic models to set the spin state of the measurement. 
     642 
     643For more complicated systems where magnetism is not uniform throughout 
     644the individual particles, you will need to write your own models. 
     645You should not mark the nuclear sld as type *sld*, but instead leave 
     646them unmarked and provide your own magnetism and polarization parameters. 
     647For 2D measurements you will need $(q_x, q_y)$ values for the measurement 
     648to compute the proper magnetism and orientation, which you can implement 
     649using *Iqxy(qx, qy, par1, par2, ...)*. 
     650 
    535651Special Functions 
    536652................. 
    … …  
    543659    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: 
    544660        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ 
    545     exp, log, pow(x,y), expm1, sqrt: 
    546         Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$. 
    547         The function expm1(x) is accurate across all $x$, including $x$ 
    548         very close to zero. 
     661    exp, log, pow(x,y), expm1, log1p, sqrt, cbrt: 
     662        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$, 
     663        $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x) 
     664        are accurate across all $x$, including $x$ very close to zero. 
    549665    sin, cos, tan, asin, acos, atan: 
    550666        Trigonometry functions and inverses, operating on radians. 
    … …  
    557673        atan(y/x) would return a value in quadrant I. Similarly for 
    558674        quadrants II and IV when $x$ and $y$ have opposite sign. 
    559     fmin(x,y), fmax(x,y), trunc, rint: 
     675    fabs(x), fmin(x,y), fmax(x,y), trunc, rint: 
    560676        Floating point functions.  rint(x) returns the nearest integer. 
    561677    NAN: 
    562678        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that 
    563679        you cannot use :code:`x == NAN` to test for NaN values since that 
    564         will always return false.  NAN does not equal NAN! 
     680        will always return false.  NAN does not equal NAN!  The alternative, 
     681        :code:`x != x` may fail if the compiler optimizes the test away. 
    565682    INFINITY: 
    566683        $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x) 
    … …  
    734851        Similar arrays are available in :code:`gauss20.c` for 20-point 
    735852        quadrature and in :code:`gauss150.c` for 150-point quadrature. 
     853        The macros :code:`GAUSS_N`, :code:`GAUSS_Z` and :code:`GAUSS_W` are 
     854        defined so that you can change the order of the integration by 
     855        selecting an different source without touching the C code. 
    736856 
    737857        :code:`source = ["lib/gauss76.c", ...]` 
    … …  
    791911can be a model that runs 1000x faster on a good card.  Even your laptop may 
    792912show a 50x improvement or more over the equivalent pure python model. 
    793  
    794 External C Models 
    795 ................. 
    796  
    797 External C models are very much like embedded C models, except that 
    798 *Iq*, *Iqxy* and *form_volume* are defined in an external source file 
    799 loaded using the *source=[...]* statement. You need to supply the function 
    800 declarations for each of these that you need instead of building them 
    801 automatically from the parameter table. 
    802913 
    803914 
    … …  
    10021113variable name *Rg* for example because $R_g$ is the right name for the model 
    10031114parameter then ignore the lint errors.  Also, ignore *missing-docstring* 
    1004 for standard model functions *Iq*, *Iqxy*, etc. 
     1115for standard model functions *Iq*, *Iqac*, etc. 
    10051116 
    10061117We will have delinting sessions at the SasView Code Camps, where we can 

  • explore/asymint.py

    r1820208 ra1c32c2  
    8686    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 
    8787    def Fq(qa, qb, qc): 
    88         siA = env.sas_sinx_x(0.5*a*qa/2) 
    89         siB = env.sas_sinx_x(0.5*b*qb/2) 
    90         siC = env.sas_sinx_x(0.5*c*qc/2) 
     88        siA = env.sas_sinx_x(a*qa/2) 
     89        siB = env.sas_sinx_x(b*qb/2) 
     90        siC = env.sas_sinx_x(c*qc/2) 
    9191        return siA * siB * siC 
    9292    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c) 
    9393    volume = a*b*c 
    9494    norm = CONTRAST**2*volume/10000 
     95    return norm, Fq 
     96 
     97def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv): 
     98    overlapping = False 
     99    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 
     100    da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc) 
     101    slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc) 
     102    dr0 = CONTRAST 
     103    drA, drB, drC = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT 
     104    tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc 
     105    def Fq(qa, qb, qc): 
     106        siA = a*env.sas_sinx_x(a*qa/2) 
     107        siB = b*env.sas_sinx_x(b*qb/2) 
     108        siC = c*env.sas_sinx_x(c*qc/2) 
     109        siAt = tA*env.sas_sinx_x(tA*qa/2) 
     110        siBt = tB*env.sas_sinx_x(tB*qb/2) 
     111        siCt = tC*env.sas_sinx_x(tC*qc/2) 
     112        if overlapping: 
     113            return (dr0*siA*siB*siC 
     114                    + drA*(siAt-siA)*siB*siC 
     115                    + drB*siAt*(siBt-siB)*siC 
     116                    + drC*siAt*siBt*(siCt-siC)) 
     117        else: 
     118            return (dr0*siA*siB*siC 
     119                    + drA*(siAt-siA)*siB*siC 
     120                    + drB*siA*(siBt-siB)*siC 
     121                    + drC*siA*siB*(siCt-siC)) 
     122    Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c) 
     123    if overlapping: 
     124        volume = a*b*c + 2*da*b*c + 2*tA*db*c + 2*tA*tB*dc 
     125    else: 
     126        volume = a*b*c + 2*da*b*c + 2*a*db*c + 2*a*b*dc 
     127    norm = 1/(volume*10000) 
    95128    return norm, Fq 
    96129 
    … …  
    184217    NORM, KERNEL = make_parallelepiped(A, B, C) 
    185218    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
     219elif shape == 'core_shell_parallelepiped': 
     220    #A, B, C = 4450, 14000, 47 
     221    #A, B, C = 445, 140, 47  # integer for the sake of mpf 
     222    A, B, C = 6800, 114, 1380 
     223    DA, DB, DC = 2300, 21, 58 
     224    SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 
     225    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
     226    #SLD_SOLVENT,CONTRAST = 0, 4 
     227    if 1: # C shortest 
     228        B, C = C, B 
     229        DB, DC = DC, DB 
     230        SLDB, SLDC = SLDC, SLDB 
     231    elif 0: # C longest 
     232        A, C = C, A 
     233        DA, DC = DC, DA 
     234        SLDA, SLDC = SLDC, SLDA 
     235    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
     236    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
    186237elif shape == 'paracrystal': 
    187238    LATTICE = 'bcc' 
    … …  
    342393    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    343394    print("gauss-500", *gauss_quad_2d(Q, n=500)) 
     395    print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
     396    print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 
    344397    #gridded_2d(Q, n=2**8+1) 
    345398    gridded_2d(Q, n=2**10+1) 
    346     #gridded_2d(Q, n=2**13+1) 
     399    #gridded_2d(Q, n=2**12+1) 
    347400    #gridded_2d(Q, n=2**15+1) 
    348     if shape != 'paracrystal':  # adaptive forms are too slow! 
     401    if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
     402        # adaptive forms on models for which the calculations are fast enough 
    349403        print("dblquad", *scipy_dblquad_2d(Q)) 
    350404        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 

  • explore/precision.py

    r2a602c7 r2a7e20e  
    345345) 
    346346add_function( 
    347     name="(1/2+(1-cos(x))/x^2-sin(x)/x)/x", 
     347    name="(1/2-sin(x)/x+(1-cos(x))/x^2)/x", 
    348348    mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x, 
    349349    np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x, 
    … …  
    609609    names = ", ".join(sorted(ALL_FUNCTIONS)) 
    610610    print("""\ 
    611 usage: precision.py [-f/a/r] [-x<range>] name... 
     611usage: precision.py [-f/a/r] [-x<range>] "name" ... 
    612612where 
    613613    -f indicates that the function value should be plotted, 
    … …  
    620620      zoom indicates linear stepping in [1000, 1010] 
    621621      neg indicates linear stepping in [-100.1, 100.1] 
    622 and name is "all [first]" or one of: 
     622and name is "all" or one of: 
    623623    """+names) 
    624624    sys.exit(1) 

  • sasmodels/compare.py

    r2d81cfe r2a7e20e  
    4242from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4343from .direct_model import DirectModel, get_mesh 
    44 from .generate import FLOAT_RE 
     44from .generate import FLOAT_RE, set_integration_size 
    4545from .weights import plot_weights 
    4646 
    … …  
    9292    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh 
    9393    -neval=1 sets the number of evals for more accurate timing 
     94    -ngauss=0 overrides the number of points in the 1-D gaussian quadrature 
    9495 
    9596    === precision options === 
    … …  
    693694        data = empty_data2D(q, resolution=res) 
    694695        data.accuracy = opts['accuracy'] 
    695         set_beam_stop(data, 0.0004) 
     696        set_beam_stop(data, qmin) 
    696697        index = ~data.mask 
    697698    else: 
    … …  
    706707    return data, index 
    707708 
    708 def make_engine(model_info, data, dtype, cutoff): 
     709def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
    709710    # type: (ModelInfo, Data, str, float) -> Calculator 
    710711    """ 
    … …  
    714715    than OpenCL. 
    715716    """ 
     717    if ngauss: 
     718        set_integration_size(model_info, ngauss) 
     719 
    716720    if dtype is None or not dtype.endswith('!'): 
    717721        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    … …  
    954958    'poly', 'mono', 'cutoff=', 
    955959    'magnetic', 'nonmagnetic', 
    956     'accuracy=', 
     960    'accuracy=', 'ngauss=', 
    957961    'neval=',  # for timing... 
    958962 
    … …  
    10891093        'show_weights' : False, 
    10901094        'sphere'    : 0, 
     1095        'ngauss'    : '0', 
    10911096    } 
    10921097    for arg in flags: 
    … …  
    11151120        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
    11161121        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
     1122        elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:] 
    11171123        elif arg.startswith('-random='): 
    11181124            opts['seed'] = int(arg[8:]) 
    … …  
    11691175 
    11701176    comparison = any(PAR_SPLIT in v for v in values) 
     1177 
    11711178    if PAR_SPLIT in name: 
    11721179        names = name.split(PAR_SPLIT, 2) 
    … …  
    11811188        return None 
    11821189 
     1190    if PAR_SPLIT in opts['ngauss']: 
     1191        opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)] 
     1192        comparison = True 
     1193    else: 
     1194        opts['ngauss'] = [int(opts['ngauss'])]*2 
     1195 
    11831196    if PAR_SPLIT in opts['engine']: 
    11841197        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    … …  
    11991212        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12001213 
    1201     base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
     1214    base = make_engine(model_info[0], data, opts['engine'][0], 
     1215                       opts['cutoff'][0], opts['ngauss'][0]) 
    12021216    if comparison: 
    1203         comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
     1217        comp = make_engine(model_info[1], data, opts['engine'][1], 
     1218                           opts['cutoff'][1], opts['ngauss'][1]) 
    12041219    else: 
    12051220        comp = None 
    … …  
    12741289        if model_info != model_info2: 
    12751290            pars2 = randomize_pars(model_info2, pars2) 
    1276             limit_dimensions(model_info, pars2, maxdim) 
     1291            limit_dimensions(model_info2, pars2, maxdim) 
    12771292            # Share values for parameters with the same name 
    12781293            for k, v in pars.items(): 

  • sasmodels/details.py

    r2d81cfe r108e70e  
    258258    # type: (...) -> Sequence[np.ndarray] 
    259259    """ 
     260    **Deprecated** Theta weights will be computed in the kernel wrapper if 
     261    they are needed. 
     262 
    260263    If there is a theta parameter, update the weights of that parameter so that 
    261264    the cosine weighting required for polar integration is preserved. 
    … …  
    272275    Returns updated weights vectors 
    273276    """ 
    274     # TODO: explain in a comment why scale and background are missing 
    275277    # Apparently the parameters.theta_offset similarly skips scale and 
    276278    # and background, so the indexing works out, but they are still shipped 
    … …  
    279281        index = parameters.theta_offset 
    280282        theta = dispersity[index] 
    281         # TODO: modify the dispersity vector to avoid the theta=-90,90,270,... 
    282283        theta_weight = abs(cos(radians(theta))) 
    283284        weights = tuple(theta_weight*w if k == index else w 

  • sasmodels/generate.py

    r2d81cfe r108e70e  
    77    particular dimensions averaged over all orientations. 
    88 
    9     *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form 
    10     with particular dimensions for a single orientation. 
    11  
    12     *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the 
    13     polarized neutron spin states (up-up, up-down, down-up, down-down) for 
    14     a form with particular dimensions for a single orientation. 
     9    *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc 
     10    for a rotationally symmetric form with particular dimensions. 
     11    qab, qc are determined from shape orientation and scattering angles. 
     12    This call is used if the shape has orientation parameters theta and phi. 
     13 
     14    *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc 
     15    for a form with particular dimensions.  qa, qb, qc are determined from 
     16    shape orientation and scattering angles. This call is used if the shape 
     17    has orientation parameters theta, phi and psi. 
     18 
     19    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy.  Use this 
     20    to create an arbitrary 2D theory function, needed for q-dependent 
     21    background functions and for models with non-uniform magnetism. 
    1522 
    1623    *form_volume(p1, p2, ...)* returns the volume of the form with particular 
    … …  
    3138scale and background parameters for each model. 
    3239 
    33 *Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 
    34 functions written for OpenCL.  All functions need prototype declarations 
    35 even if the are defined before they are used.  OpenCL does not support 
    36 *#include* preprocessor directives, so instead the list of includes needs 
    37 to be given as part of the metadata in the kernel module definition. 
    38 The included files should be listed using a path relative to the kernel 
    39 module, or if using "lib/file.c" if it is one of the standard includes 
    40 provided with the sasmodels source.  The includes need to be listed in 
    41 order so that functions are defined before they are used. 
     40C code should be stylized C-99 functions written for OpenCL. All functions 
     41need prototype declarations even if the are defined before they are used. 
     42Although OpenCL supports *#include* preprocessor directives, the list of 
     43includes should be given as part of the metadata in the kernel module 
     44definition. The included files should be listed using a path relative to the 
     45kernel module, or if using "lib/file.c" if it is one of the standard includes 
     46provided with the sasmodels source. The includes need to be listed in order 
     47so that functions are defined before they are used. 
    4248 
    4349Floating point values should be declared as *double*.  For single precision 
    … …  
    107113    present, the volume ratio is 1. 
    108114 
    109     *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the 
    110     C source code for the body of the volume, Iq, and Iqxy functions 
     115    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
     116    the C source code for the body of the volume, Iq, and Iqac functions 
    111117    respectively.  These can also be defined in the last source file. 
    112118 
    113     *Iq* and *Iqxy* also be instead be python functions defining the 
     119    *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the 
    114120    kernel.  If they are marked as *Iq.vectorized = True* then the 
    115121    kernel is passed the entire *q* vector at once, otherwise it is 
    … …  
    168174from zlib import crc32 
    169175from inspect import currentframe, getframeinfo 
     176import logging 
    170177 
    171178import numpy as np  # type: ignore 
    … …  
    181188    pass 
    182189# pylint: enable=unused-import 
     190 
     191logger = logging.getLogger(__name__) 
    183192 
    184193# jitter projection to use in the kernel code.  See explore/jitter.py 
    … …  
    270279""" 
    271280 
     281 
     282def set_integration_size(info, n): 
     283    # type: (ModelInfo, int) -> None 
     284    """ 
     285    Update the model definition, replacing the gaussian integration with 
     286    a gaussian integration of a different size. 
     287 
     288    Note: this really ought to be a method in modelinfo, but that leads to 
     289    import loops. 
     290    """ 
     291    if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)): 
     292        import os.path 
     293        from .gengauss import gengauss 
     294        path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n) 
     295        if not os.path.exists(path): 
     296            gengauss(n, path) 
     297        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 
     298                        else lib for lib in info.source] 
    272299 
    273300def format_units(units): 
    … …  
    608635 
    609636""" 
    610 def _gen_fn(name, pars, body, filename, line): 
    611     # type: (str, List[Parameter], str, str, int) -> str 
     637def _gen_fn(model_info, name, pars): 
     638    # type: (ModelInfo, str, List[Parameter]) -> str 
    612639    """ 
    613640    Generate a function given pars and body. 
    … …  
    621648    """ 
    622649    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 
     650    body = getattr(model_info, name) 
     651    filename = model_info.filename 
     652    # Note: if symbol is defined strangely in the module then default it to 1 
     653    lineno = model_info.lineno.get(name, 1) 
    623654    return _FN_TEMPLATE % { 
    624655        'name': name, 'pars': par_decl, 'body': body, 
    625         'filename': filename.replace('\\', '\\\\'), 'line': line, 
     656        'filename': filename.replace('\\', '\\\\'), 'line': lineno, 
    626657    } 
    627658 
    … …  
    638669 
    639670# type in IQXY pattern could be single, float, double, long double, ... 
    640 _IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)", 
     671_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
    641672                           flags=re.MULTILINE) 
    642 def _have_Iqxy(sources): 
     673def find_xy_mode(source): 
    643674    # type: (List[str]) -> bool 
    644675    """ 
    645     Return true if any file defines Iqxy. 
     676    Return the xy mode as qa, qac, qabc or qxy. 
    646677 
    647678    Note this is not a C parser, and so can be easily confused by 
    648679    non-standard syntax.  Also, it will incorrectly identify the following 
    649     as having Iqxy:: 
     680    as having 2D models:: 
    650681 
    651682        /* 
    652         double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     683        double Iqac(qab, qc, ...) { ... fill this in later ... } 
    653684        */ 
    654685 
    655     If you want to comment out an Iqxy function, use // on the front of the 
    656     line instead. 
    657     """ 
    658     for _path, code in sources: 
    659         if _IQXY_PATTERN.search(code): 
    660             return True 
    661     return False 
    662  
    663  
    664 def _add_source(source, code, path): 
     686    If you want to comment out the function, use // on the front of the 
     687    line:: 
     688 
     689        /* 
     690        // double Iqac(qab, qc, ...) { ... fill this in later ... } 
     691        */ 
     692 
     693    """ 
     694    for code in source: 
     695        m = _IQXY_PATTERN.search(code) 
     696        if m is not None: 
     697            return m.group('mode') 
     698    return 'qa' 
     699 
     700 
     701def _add_source(source, code, path, lineno=1): 
    665702    """ 
    666703    Add a file to the list of source code chunks, tagged with path and line. 
    667704    """ 
    668705    path = path.replace('\\', '\\\\') 
    669     source.append('#line 1 "%s"' % path) 
     706    source.append('#line %d "%s"' % (lineno, path)) 
    670707    source.append(code) 
    671708 
    … …  
    698735    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    699736 
    700     # What kind of 2D model do we need? 
    701     xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str) 
    702                else 'qac' if not partable.is_asymmetric 
    703                else 'qabc') 
    704  
    705737    # Build initial sources 
    706738    source = [] 
    … …  
    709741        _add_source(source, code, path) 
    710742 
     743    if model_info.c_code: 
     744        _add_source(source, model_info.c_code, model_info.filename, 
     745                    lineno=model_info.lineno.get('c_code', 1)) 
     746 
    711747    # Make parameters for q, qx, qy so that we can use them in declarations 
    712     q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     748    q, qx, qy, qab, qa, qb, qc \ 
     749        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
    713750    # Generate form_volume function, etc. from body only 
    714751    if isinstance(model_info.form_volume, str): 
    715752        pars = partable.form_volume_parameters 
    716         source.append(_gen_fn('form_volume', pars, model_info.form_volume, 
    717                               model_info.filename, model_info._form_volume_line)) 
     753        source.append(_gen_fn(model_info, 'form_volume', pars)) 
    718754    if isinstance(model_info.Iq, str): 
    719755        pars = [q] + partable.iq_parameters 
    720         source.append(_gen_fn('Iq', pars, model_info.Iq, 
    721                               model_info.filename, model_info._Iq_line)) 
     756        source.append(_gen_fn(model_info, 'Iq', pars)) 
    722757    if isinstance(model_info.Iqxy, str): 
    723         pars = [qx, qy] + partable.iqxy_parameters 
    724         source.append(_gen_fn('Iqxy', pars, model_info.Iqxy, 
    725                               model_info.filename, model_info._Iqxy_line)) 
     758        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters 
     759        source.append(_gen_fn(model_info, 'Iqxy', pars)) 
     760    if isinstance(model_info.Iqac, str): 
     761        pars = [qab, qc] + partable.iq_parameters 
     762        source.append(_gen_fn(model_info, 'Iqac', pars)) 
     763    if isinstance(model_info.Iqabc, str): 
     764        pars = [qa, qb, qc] + partable.iq_parameters 
     765        source.append(_gen_fn(model_info, 'Iqabc', pars)) 
     766 
     767    # What kind of 2D model do we need?  Is it consistent with the parameters? 
     768    xy_mode = find_xy_mode(source) 
     769    if xy_mode == 'qabc' and not partable.is_asymmetric: 
     770        raise ValueError("asymmetric oriented models need to define Iqabc") 
     771    elif xy_mode == 'qac' and partable.is_asymmetric: 
     772        raise ValueError("symmetric oriented models need to define Iqac") 
     773    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'): 
     774        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode) 
     775    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'): 
     776        if xy_mode == 'qxy': 
     777            logger.warn("oriented shapes should define Iqac or Iqabc") 
     778        else: 
     779            raise ValueError("Expected function Iqac or Iqabc for oriented shape") 
    726780 
    727781    # Define the parameter table 
    … …  
    749803    if xy_mode == 'qabc': 
    750804        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    751         call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
     805        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars 
    752806        clear_iqxy = "#undef CALL_IQ_ABC" 
    753807    elif xy_mode == 'qac': 
    754808        pars = ",".join(["_qa", "_qc"] + model_refs) 
    755         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
     809        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    756810        clear_iqxy = "#undef CALL_IQ_AC" 
    757     else:  # xy_mode == 'qa' 
     811    elif xy_mode == 'qa': 
    758812        pars = ",".join(["_qa"] + model_refs) 
    759813        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    760814        clear_iqxy = "#undef CALL_IQ_A" 
     815    elif xy_mode == 'qxy': 
     816        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
     817        pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs) 
     818        call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars 
     819        clear_iqxy = "#undef CALL_IQ_XY" 
     820        if partable.orientation_parameters: 
     821            call_iqxy += "\n#define HAVE_THETA" 
     822            clear_iqxy += "\n#undef HAVE_THETA" 
     823        if partable.is_asymmetric: 
     824            call_iqxy += "\n#define HAVE_PSI" 
     825            clear_iqxy += "\n#undef HAVE_PSI" 
     826 
    761827 
    762828    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 

  • sasmodels/kernel_header.c

    r8698a0d r108e70e  
    150150inline double cube(double x) { return x*x*x; } 
    151151inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 
     152 
     153// CRUFT: support old style models with orientation received qx, qy and angles 
     154 
     155// To rotate from the canonical position to theta, phi, psi, first rotate by 
     156// psi about the major axis, oriented along z, which is a rotation in the 
     157// detector plane xy. Next rotate by theta about the y axis, aligning the major 
     158// axis in the xz plane. Finally, rotate by phi in the detector plane xy. 
     159// To compute the scattering, undo these rotations in reverse order: 
     160//     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi 
     161// The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit 
     162// vector in the q direction. 
     163// To change between counterclockwise and clockwise rotation, change the 
     164// sign of phi and psi. 
     165 
     166#if 1 
     167//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
     168#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     169    SINCOS(phi*M_PI_180, sn, cn); \ 
     170    q = sqrt(qx*qx + qy*qy); \ 
     171    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
     172    sn = sqrt(1 - cn*cn); \ 
     173    } while (0) 
     174#else 
     175// SasView 3.x definition of orientation 
     176#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     177    SINCOS(theta*M_PI_180, sn, cn); \ 
     178    q = sqrt(qx*qx + qy*qy);\ 
     179    cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ 
     180    sn = sqrt(1 - cn*cn); \ 
     181    } while (0) 
     182#endif 
     183 
     184#if 1 
     185#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \ 
     186    q = sqrt(qx*qx + qy*qy); \ 
     187    const double qxhat = qx/q; \ 
     188    const double qyhat = qy/q; \ 
     189    double sin_theta, cos_theta; \ 
     190    double sin_phi, cos_phi; \ 
     191    double sin_psi, cos_psi; \ 
     192    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     193    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     194    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     195    xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \ 
     196         + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \ 
     197    yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \ 
     198         + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \ 
     199    zhat = qxhat*(-sin_theta*cos_phi) \ 
     200         + qyhat*(-sin_theta*sin_phi); \ 
     201    } while (0) 
     202#else 
     203// SasView 3.x definition of orientation 
     204#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
     205    q = sqrt(qx*qx + qy*qy); \ 
     206    const double qxhat = qx/q; \ 
     207    const double qyhat = qy/q; \ 
     208    double sin_theta, cos_theta; \ 
     209    double sin_phi, cos_phi; \ 
     210    double sin_psi, cos_psi; \ 
     211    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     212    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     213    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     214    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 
     215    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 
     216    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
     217    } while (0) 
     218#endif 

  • sasmodels/kernel_iq.c

    r6aee3ab rec8d4ac  
    3131//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
    3232//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
     33//  CALL_IQ_XY(qx, qy, table) : call the Iqxy function for arbitrary models 
    3334//  INVALID(table) : test if the current point is feesible to calculate.  This 
    3435//      will be defined in the kernel definition file. 
    … …  
    173174static double 
    174175qac_apply( 
    175     QACRotation rotation, 
     176    QACRotation *rotation, 
    176177    double qx, double qy, 
    177178    double *qa_out, double *qc_out) 
    178179{ 
    179     const double dqc = rotation.R31*qx + rotation.R32*qy; 
     180    const double dqc = rotation->R31*qx + rotation->R32*qy; 
    180181    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    181182    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
    … …  
    246247static double 
    247248qabc_apply( 
    248     QABCRotation rotation, 
     249    QABCRotation *rotation, 
    249250    double qx, double qy, 
    250251    double *qa_out, double *qb_out, double *qc_out) 
    251252{ 
    252     *qa_out = rotation.R11*qx + rotation.R12*qy; 
    253     *qb_out = rotation.R21*qx + rotation.R22*qy; 
    254     *qc_out = rotation.R31*qx + rotation.R32*qy; 
     253    *qa_out = rotation->R11*qx + rotation->R12*qy; 
     254    *qb_out = rotation->R21*qx + rotation->R22*qy; 
     255    *qc_out = rotation->R31*qx + rotation->R32*qy; 
    255256} 
    256257 
    … …  
    453454  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code. 
    454455  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi); 
    455   #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 
     456  #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc) 
    456457  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) 
    457458 
    … …  
    467468  local_values.table.psi = 0.; 
    468469  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi) 
    469   #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
     470  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    470471  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
    471 #endif 
    472  
    473 // Doing jitter projection code outside the previous if block so that we don't 
    474 // need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This 
    475 // will become more important if we implement more projections, or more 
    476 // complicated projections. 
    477 #if defined(CALL_IQ) || defined(CALL_IQ_A) 
     472#elif defined(CALL_IQ_XY) 
     473  // direct call to qx,qy calculator 
     474  double qx, qy; 
     475  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     476  #define BUILD_ROTATION() do {} while(0) 
     477  #define APPLY_ROTATION() do {} while(0) 
     478  #define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table) 
     479#endif 
     480 
     481// Define APPLY_PROJECTION depending on model symmetries. We do this outside 
     482// the previous if block so that we don't need to repeat the identical 
     483// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
     484// if we implement more projections, or more complicated projections. 
     485#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation 
    478486  #define APPLY_PROJECTION() const double weight=weight0 
    479 #else // !spherosymmetric projection 
     487#elif defined(CALL_IQ_XY) // pass orientation to the model 
     488  // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc 
     489  // Need to plug the values for the orientation angles back into parameter 
     490  // table in case they were overridden by the orientation offset.  This 
     491  // means that orientation dispersity will not work for these models, but 
     492  // it was broken anyway, so no matter.  Still want to provide Iqxy in case 
     493  // the user model wants full control of orientation/magnetism. 
     494  #if defined(HAVE_PSI) 
     495    const double theta = values[details->theta_par+2]; 
     496    const double phi = values[details->theta_par+3]; 
     497    const double psi = values[details->theta_par+4]; 
     498    double weight; 
     499    #define APPLY_PROJECTION() do { \ 
     500      local_values.table.theta = theta; \ 
     501      local_values.table.phi = phi; \ 
     502      local_values.table.psi = psi; \ 
     503      weight=weight0; \ 
     504    } while (0) 
     505  #elif defined(HAVE_THETA) 
     506    const double theta = values[details->theta_par+2]; 
     507    const double phi = values[details->theta_par+3]; 
     508    double weight; 
     509    #define APPLY_PROJECTION() do { \ 
     510      local_values.table.theta = theta; \ 
     511      local_values.table.phi = phi; \ 
     512      weight=weight0; \ 
     513    } while (0) 
     514  #else 
     515    #define APPLY_PROJECTION() const double weight=weight0 
     516  #endif 
     517#else // apply jitter and view before calling the model 
    480518  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 
    481519  const double theta = values[details->theta_par+2]; 
    … …  
    488526  // we go through the mesh. 
    489527  double dtheta, dphi, weight; 
    490   #if PROJECTION == 1 
     528  #if PROJECTION == 1 // equirectangular 
    491529    #define APPLY_PROJECTION() do { \ 
    492530      dtheta = local_values.table.theta; \ 
    … …  
    494532      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \ 
    495533    } while (0) 
    496   #elif PROJECTION == 2 
     534  #elif PROJECTION == 2 // sinusoidal 
    497535    #define APPLY_PROJECTION() do { \ 
    498536      dtheta = local_values.table.theta; \ 
    … …  
    504542    } while (0) 
    505543  #endif 
    506 #endif // !spherosymmetric projection 
     544#endif // done defining APPLY_PROJECTION 
    507545 
    508546// ** define looping macros ** 

  • sasmodels/kernelpy.py

    r2d81cfe r108e70e  
    2626# pylint: enable=unused-import 
    2727 
     28logger = logging.getLogger(__name__) 
     29 
    2830class PyModel(KernelModel): 
    2931    """ 
    … …  
    3133    """ 
    3234    def __init__(self, model_info): 
    33         # Make sure Iq and Iqxy are available and vectorized 
     35        # Make sure Iq is available and vectorized 
    3436        _create_default_functions(model_info) 
    3537        self.info = model_info 
    3638        self.dtype = np.dtype('d') 
    37         logging.info("load python model " + self.info.name) 
     39        logger.info("load python model " + self.info.name) 
    3840 
    3941    def make_kernel(self, q_vectors): 
    4042        q_input = PyInput(q_vectors, dtype=F64) 
    41         kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
    42         return PyKernel(kernel, self.info, q_input) 
     43        return PyKernel(self.info, q_input) 
    4344 
    4445    def release(self): 
    … …  
    8990    Callable SAS kernel. 
    9091 
    91     *kernel* is the DllKernel object to call. 
     92    *kernel* is the kernel object to call. 
    9293 
    9394    *model_info* is the module information 
    … …  
    104105    Call :meth:`release` when done with the kernel instance. 
    105106    """ 
    106     def __init__(self, kernel, model_info, q_input): 
     107    def __init__(self, model_info, q_input): 
    107108        # type: (callable, ModelInfo, List[np.ndarray]) -> None 
    108109        self.dtype = np.dtype('d') 
    … …  
    110111        self.q_input = q_input 
    111112        self.res = np.empty(q_input.nq, q_input.dtype) 
    112         self.kernel = kernel 
    113113        self.dim = '2d' if q_input.is_2d else '1d' 
    114114 
    … …  
    159159        # Generate a closure which calls the form_volume if it exists. 
    160160        form_volume = model_info.form_volume 
    161         self._volume = ((lambda: form_volume(*volume_args)) if form_volume 
    162                         else (lambda: 1.0)) 
     161        self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 
     162                        (lambda: 1.0)) 
    163163 
    164164    def __call__(self, call_details, values, cutoff, magnetic): 
    … …  
    261261    any functions that are not already marked as vectorized. 
    262262    """ 
     263    # Note: must call create_vector_Iq before create_vector_Iqxy 
    263264    _create_vector_Iq(model_info) 
    264     _create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
     265    _create_vector_Iqxy(model_info) 
    265266 
    266267 
    … …  
    280281        model_info.Iq = vector_Iq 
    281282 
     283 
    282284def _create_vector_Iqxy(model_info): 
    283285    """ 
    284286    Define Iqxy as a vector function if it exists, or default it from Iq(). 
    285287    """ 
    286     Iq, Iqxy = model_info.Iq, model_info.Iqxy 
     288    Iqxy = getattr(model_info, 'Iqxy', None) 
    287289    if callable(Iqxy): 
    288290        if not getattr(Iqxy, 'vectorized', False): 
    … …  
    295297            vector_Iqxy.vectorized = True 
    296298            model_info.Iqxy = vector_Iqxy 
    297     elif callable(Iq): 
     299    else: 
    298300        #print("defaulting Iqxy") 
    299301        # Iq is vectorized because create_vector_Iq was already called. 
     302        Iq = model_info.Iq 
    300303        def default_Iqxy(qx, qy, *args): 
    301304            """ 

  • sasmodels/modelinfo.py

    r2d81cfe r108e70e  
    3737 
    3838# assumptions about common parameters exist throughout the code, such as: 
    39 # (1) kernel functions Iq, Iqxy, form_volume, ... don't see them 
     39# (1) kernel functions Iq, Iqxy, Iqac, Iqabc, form_volume, ... don't see them 
    4040# (2) kernel drivers assume scale is par[0] and background is par[1] 
    4141# (3) mixture models drop the background on components and replace the scale 
    … …  
    256256 
    257257    *type* indicates how the parameter will be used.  "volume" parameters 
    258     will be used in all functions.  "orientation" parameters will be used 
    259     in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
    260     *Imagnetic* only.  If *type* is the empty string, the parameter will 
     258    will be used in all functions.  "orientation" parameters are not passed, 
     259    but will be used to convert from *qx*, *qy* to *qa*, *qb*, *qc* in calls to 
     260    *Iqxy* and *Imagnetic*.  If *type* is the empty string, the parameter will 
    261261    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters 
    262262    can automatically be promoted to magnetic parameters, each of which 
    … …  
    386386      with vector parameter p sent as p[]. 
    387387 
    388     * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
    389       function, with vector parameter p sent as p[]. 
    390  
    391388    * *form_volume_parameters* is the list of parameters to the form_volume(...) 
    392389      function, with vector parameter p sent as p[]. 
    … …  
    443440        self.iq_parameters = [p for p in self.kernel_parameters 
    444441                              if p.type not in ('orientation', 'magnetic')] 
    445         # note: orientation no longer sent to Iqxy, so its the same as 
    446         #self.iqxy_parameters = [p for p in self.kernel_parameters 
    447         #                        if p.type != 'magnetic'] 
     442        self.orientation_parameters = [p for p in self.kernel_parameters 
     443                                       if p.type == 'orientation'] 
    448444        self.form_volume_parameters = [p for p in self.kernel_parameters 
    449445                                       if p.type == 'volume'] 
    … …  
    490486                if p.type != 'orientation': 
    491487                    raise TypeError("psi must be an orientation parameter") 
     488            elif p.type == 'orientation': 
     489                raise TypeError("only theta, phi and psi can be orientation parameters") 
    492490        if theta >= 0 and phi >= 0: 
     491            last_par = len(self.kernel_parameters) - 1 
    493492            if phi != theta+1: 
    494493                raise TypeError("phi must follow theta") 
    495494            if psi >= 0 and psi != phi+1: 
    496495                raise TypeError("psi must follow phi") 
     496            #if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): 
     497            #    raise TypeError("orientation parameters must appear at the " 
     498            #                    "end of the parameter table") 
    497499        elif theta >= 0 or phi >= 0 or psi >= 0: 
    498500            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
    … …  
    715717 
    716718 
     719#: Set of variables defined in the model that might contain C code 
     720C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     721 
    717722def _find_source_lines(model_info, kernel_module): 
    718723    # type: (ModelInfo, ModuleType) -> None 
    … …  
    720725    Identify the location of the C source inside the model definition file. 
    721726 
    722     This code runs through the source of the kernel module looking for 
    723     lines that start with 'Iq', 'Iqxy' or 'form_volume'.  Clearly there are 
    724     all sorts of reasons why this might not work (e.g., code commented out 
    725     in a triple-quoted line block, code built using string concatenation, 
    726     or code defined in the branch of an 'if' block), but it should work 
    727     properly in the 95% case, and getting the incorrect line number will 
    728     be harmless. 
    729     """ 
    730     # Check if we need line numbers at all 
    731     if callable(model_info.Iq): 
    732         return None 
    733  
    734     if (model_info.Iq is None 
    735             and model_info.Iqxy is None 
    736             and model_info.Imagnetic is None 
    737             and model_info.form_volume is None): 
     727    This code runs through the source of the kernel module looking for lines 
     728    that contain C code (because it is a c function definition).  Clearly 
     729    there are all sorts of reasons why this might not work (e.g., code 
     730    commented out in a triple-quoted line block, code built using string 
     731    concatenation, code defined in the branch of an 'if' block, code imported 
     732    from another file), but it should work properly in the 95% case, and for 
     733    the remainder, getting the incorrect line number will merely be 
     734    inconvenient. 
     735    """ 
     736    # Only need line numbers if we are creating a C module and the C symbols 
     737    # are defined. 
     738    if (callable(model_info.Iq) 
     739            or not any(hasattr(model_info, s) for s in C_SYMBOLS)): 
    738740        return 
    739741 
    740     # find the defintion lines for the different code blocks 
     742    # load the module source if we can 
    741743    try: 
    742744        source = inspect.getsource(kernel_module) 
    743745    except IOError: 
    744746        return 
    745     for k, v in enumerate(source.split('\n')): 
    746         if v.startswith('Imagnetic'): 
    747             model_info._Imagnetic_line = k+1 
    748         elif v.startswith('Iqxy'): 
    749             model_info._Iqxy_line = k+1 
    750         elif v.startswith('Iq'): 
    751             model_info._Iq_line = k+1 
    752         elif v.startswith('form_volume'): 
    753             model_info._form_volume_line = k+1 
    754  
     747 
     748    # look for symbol at the start of the line 
     749    for lineno, line in enumerate(source.split('\n')): 
     750        for name in C_SYMBOLS: 
     751            if line.startswith(name): 
     752                # Add 1 since some compilers complain about "#line 0" 
     753                model_info.lineno[name] = lineno + 1 
     754                break 
    755755 
    756756def make_model_info(kernel_module): 
    … …  
    761761    Fill in default values for parts of the module that are not provided. 
    762762 
    763     Note: vectorized Iq and Iqxy functions will be created for python 
     763    Note: vectorized Iq and Iqac/Iqabc functions will be created for python 
    764764    models when the model is first called, not when the model is loaded. 
    765765    """ 
    … …  
    790790    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    791791    info.source = getattr(kernel_module, 'source', []) 
     792    info.c_code = getattr(kernel_module, 'c_code', None) 
    792793    # TODO: check the structure of the tests 
    793794    info.tests = getattr(kernel_module, 'tests', []) 
    … …  
    797798    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    798799    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     800    info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore 
     801    info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore 
    799802    info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore 
    800803    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    … …  
    811814    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
    812815 
     816    if callable(info.Iq) and parameters.has_2d: 
     817        raise ValueError("oriented python models not supported") 
     818 
     819    info.lineno = {} 
    813820    _find_source_lines(info, kernel_module) 
    814821 
    … …  
    824831 
    825832    The structure should be mostly static, other than the delayed definition 
    826     of *Iq* and *Iqxy* if they need to be defined. 
     833    of *Iq*, *Iqac* and *Iqabc* if they need to be defined. 
    827834    """ 
    828835    #: Full path to the file defining the kernel, if any. 
    … …  
    906913    structure_factor = None # type: bool 
    907914    #: List of C source files used to define the model.  The source files 
    908     #: should define the *Iq* function, and possibly *Iqxy*, though a default 
    909     #: *Iqxy = Iq(sqrt(qx**2+qy**2)* will be created if no *Iqxy* is provided. 
    910     #: Files containing the most basic functions must appear first in the list, 
    911     #: followed by the files that use those functions.  Form factors are 
    912     #: indicated by providing a :attr:`ER` function. 
     915    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
     916    #: model defines orientation parameters. Files containing the most basic 
     917    #: functions must appear first in the list, followed by the files that 
     918    #: use those functions.  Form factors are indicated by providing 
     919    #: an :attr:`ER` function. 
    913920    source = None           # type: List[str] 
    914921    #: The set of tests that must pass.  The format of the tests is described 
    … …  
    935942    #: See :attr:`ER` for details on the parameters. 
    936943    VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
     944    #: Arbitrary C code containing supporting functions, etc., to be inserted 
     945    #: after everything in source.  This can include Iq and Iqxy functions with 
     946    #: the full function signature, including all parameters. 
     947    c_code = None 
    937948    #: Returns the form volume for python-based models.  Form volume is needed 
    938949    #: for volume normalization in the polydispersity integral.  If no 
    … …  
    955966    #: include the decimal point. See :mod:`generate` for more details. 
    956967    Iq = None               # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    957     #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    958     Iqxy = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     968    #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
     969    Iqac = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     970    #: Returns *I(qa, qb, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
     971    Iqabc = None            # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    959972    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    960973    Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    … …  
    972985    #: Returns a random parameter set for the model 
    973986    random = None           # type: Optional[Callable[[], Dict[str, float]]] 
    974  
    975     # line numbers within the python file for bits of C source, if defined 
    976     # NB: some compilers fail with a "#line 0" directive, so default to 1. 
    977     _Imagnetic_line = 1 
    978     _Iqxy_line = 1 
    979     _Iq_line = 1 
    980     _form_volume_line = 1 
    981  
     987    #: Line numbers for symbols defining C code 
     988    lineno = None           # type: Dict[str, int] 
    982989 
    983990    def __init__(self): 

  • sasmodels/models/_spherepy.py

    ref07e95 r108e70e  
    8888Iq.vectorized = True  # Iq accepts an array of q values 
    8989 
    90 def Iqxy(qx, qy, sld, sld_solvent, radius): 
    91     return Iq(sqrt(qx ** 2 + qy ** 2), sld, sld_solvent, radius) 
    92 Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values 
    93  
    9490def sesans(z, sld, sld_solvent, radius): 
    9591    """ 

  • sasmodels/models/barbell.c

    rbecded3 r108e70e  
    2323    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    2424    double total = 0.0; 
    25     for (int i = 0; i < 76; i++){ 
    26         const double t = Gauss76Z[i]*zm + zb; 
     25    for (int i = 0; i < GAUSS_N; i++){ 
     26        const double t = GAUSS_Z[i]*zm + zb; 
    2727        const double radical = 1.0 - t*t; 
    2828        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    2929        const double Fq = cos(m*t + b) * radical * bj; 
    30         total += Gauss76Wt[i] * Fq; 
     30        total += GAUSS_W[i] * Fq; 
    3131    } 
    3232    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    7373    const double zb = M_PI_4; 
    7474    double total = 0.0; 
    75     for (int i = 0; i < 76; i++){ 
    76         const double alpha= Gauss76Z[i]*zm + zb; 
     75    for (int i = 0; i < GAUSS_N; i++){ 
     76        const double alpha= GAUSS_Z[i]*zm + zb; 
    7777        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    7878        SINCOS(alpha, sin_alpha, cos_alpha); 
    7979        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    80         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     80        total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
    8181    } 
    8282    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9090 
    9191static double 
    92 Iqxy(double qab, double qc, 
     92Iqac(double qab, double qc, 
    9393    double sld, double solvent_sld, 
    9494    double radius_bell, double radius, double length) 

  • sasmodels/models/bcc_paracrystal.c

    rea60e08 r108e70e  
    8181 
    8282    double outer_sum = 0.0; 
    83     for(int i=0; i<150; i++) { 
     83    for(int i=0; i<GAUSS_N; i++) { 
    8484        double inner_sum = 0.0; 
    85         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     85        const double theta = GAUSS_Z[i]*theta_m + theta_b; 
    8686        double sin_theta, cos_theta; 
    8787        SINCOS(theta, sin_theta, cos_theta); 
    8888        const double qc = q*cos_theta; 
    8989        const double qab = q*sin_theta; 
    90         for(int j=0;j<150;j++) { 
    91             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     90        for(int j=0;j<GAUSS_N;j++) { 
     91            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    9292            double sin_phi, cos_phi; 
    9393            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    9595            const double qb = qab*sin_phi; 
    9696            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
    97             inner_sum += Gauss150Wt[j] * form; 
     97            inner_sum += GAUSS_W[j] * form; 
    9898        } 
    9999        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    100         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     100        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    101101    } 
    102102    outer_sum *= theta_m; 
    … …  
    107107 
    108108 
    109 static double Iqxy(double qa, double qb, double qc, 
     109static double Iqabc(double qa, double qb, double qc, 
    110110    double dnn, double d_factor, double radius, 
    111111    double sld, double solvent_sld) 

  • sasmodels/models/capped_cylinder.c

    rbecded3 r108e70e  
    3030    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    32     for (int i=0; i<76 ;i++) { 
    33         const double t = Gauss76Z[i]*zm + zb; 
     32    for (int i=0; i<GAUSS_N; i++) { 
     33        const double t = GAUSS_Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    3535        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    37         total += Gauss76Wt[i] * Fq; 
     37        total += GAUSS_W[i] * Fq; 
    3838    } 
    3939    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9595    const double zb = M_PI_4; 
    9696    double total = 0.0; 
    97     for (int i=0; i<76 ;i++) { 
    98         const double theta = Gauss76Z[i]*zm + zb; 
     97    for (int i=0; i<GAUSS_N ;i++) { 
     98        const double theta = GAUSS_Z[i]*zm + zb; 
    9999        double sin_theta, cos_theta; // slots to hold sincos function output 
    100100        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    103103        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
    104104        // scale by sin_theta for spherical coord integration 
    105         total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
     105        total += GAUSS_W[i] * Aq * Aq * sin_theta; 
    106106    } 
    107107    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    115115 
    116116static double 
    117 Iqxy(double qab, double qc, 
     117Iqac(double qab, double qc, 
    118118    double sld, double solvent_sld, double radius, 
    119119    double radius_cap, double length) 

  • sasmodels/models/core_shell_bicelle.c

    rbecded3 r108e70e  
    5252 
    5353    double total = 0.0; 
    54     for(int i=0;i<N_POINTS_76;i++) { 
    55         double theta = (Gauss76Z[i] + 1.0)*uplim; 
     54    for(int i=0;i<GAUSS_N;i++) { 
     55        double theta = (GAUSS_Z[i] + 1.0)*uplim; 
    5656        double sin_theta, cos_theta; // slots to hold sincos function output 
    5757        SINCOS(theta, sin_theta, cos_theta); 
    5858        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
    5959                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
    60         total += Gauss76Wt[i]*fq*fq*sin_theta; 
     60        total += GAUSS_W[i]*fq*fq*sin_theta; 
    6161    } 
    6262 
    … …  
    6767 
    6868static double 
    69 Iqxy(double qab, double qc, 
     69Iqac(double qab, double qc, 
    7070    double radius, 
    7171    double thick_rim, 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r82592da r108e70e  
    3737    //initialize integral 
    3838    double outer_sum = 0.0; 
    39     for(int i=0;i<76;i++) { 
     39    for(int i=0;i<GAUSS_N;i++) { 
    4040        //setup inner integral over the ellipsoidal cross-section 
    4141        //const double va = 0.0; 
    4242        //const double vb = 1.0; 
    43         //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    44         const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
     43        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
    4545        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    4646        const double qab = q*sin_theta; 
    … …  
    4949        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
    5050        double inner_sum=0.0; 
    51         for(int j=0;j<76;j++) { 
     51        for(int j=0;j<GAUSS_N;j++) { 
    5252            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    5353            // inner integral limits 
    5454            //const double vaj=0.0; 
    5555            //const double vbj=M_PI; 
    56             //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
     56            //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    5858            const double rr = sqrt(r2A - r2B*cos(phi)); 
    5959            const double be1 = sas_2J1x_x(rr*qab); 
    … …  
    6161            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
    6262 
    63             inner_sum += Gauss76Wt[j] * fq * fq; 
     63            inner_sum += GAUSS_W[j] * fq * fq; 
    6464        } 
    6565        //now calculate outer integral 
    66         outer_sum += Gauss76Wt[i] * inner_sum; 
     66        outer_sum += GAUSS_W[i] * inner_sum; 
    6767    } 
    6868 
    … …  
    7171 
    7272static double 
    73 Iqxy(double qa, double qb, double qc, 
     73Iqabc(double qa, double qb, double qc, 
    7474    double r_minor, 
    7575    double x_core, 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r82592da r108e70e  
    77        double length) 
    88{ 
    9     return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length +  
     9    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 
    1010                 square(r_minor)*x_core*2.0*thick_face  ); 
    1111} 
    … …  
    4747    //initialize integral 
    4848    double outer_sum = 0.0; 
    49     for(int i=0;i<76;i++) { 
     49    for(int i=0;i<GAUSS_N;i++) { 
    5050        //setup inner integral over the ellipsoidal cross-section 
    5151        // since we generate these lots of times, why not store them somewhere? 
    52         //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
     52        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     53        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 
    5454        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    5555        double inner_sum=0; 
    … …  
    5858        si1 = sas_sinx_x(sinarg1); 
    5959        si2 = sas_sinx_x(sinarg2); 
    60         for(int j=0;j<76;j++) { 
     60        for(int j=0;j<GAUSS_N;j++) { 
    6161            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    62             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    63             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
     62            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     63            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    6464            const double rr = sqrt(r2A - r2B*cos(beta)); 
    6565            double besarg1 = q*rr*sin_alpha; 
    … …  
    6767            be1 = sas_2J1x_x(besarg1); 
    6868            be2 = sas_2J1x_x(besarg2); 
    69             inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
     69            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
    7070                                              dr2*si1*be2 + 
    7171                                              dr3*si2*be1); 
    7272        } 
    7373        //now calculate outer integral 
    74         outer_sum += Gauss76Wt[i] * inner_sum; 
     74        outer_sum += GAUSS_W[i] * inner_sum; 
    7575    } 
    7676 
    … …  
    7979 
    8080static double 
    81 Iqxy(double qa, double qb, double qc, 
     81Iqabc(double qa, double qb, double qc, 
    8282          double r_minor, 
    8383          double x_core, 
    … …  
    114114    return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 
    115115} 
    116  

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r110f69c r108e70e  
    149149    ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    150150    ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"], 
     151    ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"], 
    151152    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    152153    ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
    153154    ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"], 
    154     ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"] 
    155155    ] 
    156156 

  • sasmodels/models/core_shell_cylinder.c

    rbecded3 r108e70e  
    3030    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3131    double total = 0.0; 
    32     for (int i=0; i<76 ;i++) { 
     32    for (int i=0; i<GAUSS_N ;i++) { 
    3333        // translate a point in [-1,1] to a point in [0, pi/2] 
    34         //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     34        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 
    3535        double sin_theta, cos_theta; 
    36         const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     36        const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4; 
    3737        SINCOS(theta, sin_theta,  cos_theta); 
    3838        const double qab = q*sin_theta; 
    … …  
    4040        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    4141            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    42         total += Gauss76Wt[i] * fq * fq * sin_theta; 
     42        total += GAUSS_W[i] * fq * fq * sin_theta; 
    4343    } 
    4444    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4848 
    4949 
    50 double Iqxy(double qab, double qc, 
     50double Iqac(double qab, double qc, 
    5151    double core_sld, 
    5252    double shell_sld, 

  • sasmodels/models/core_shell_ellipsoid.c

    rbecded3 r108e70e  
    5959    const double b = 0.5; 
    6060    double total = 0.0;     //initialize intergral 
    61     for(int i=0;i<76;i++) { 
    62         const double cos_theta = Gauss76Z[i]*m + b; 
     61    for(int i=0;i<GAUSS_N;i++) { 
     62        const double cos_theta = GAUSS_Z[i]*m + b; 
    6363        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    6464        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 
    … …  
    6666            equat_shell, polar_shell, 
    6767            sld_core_shell, sld_shell_solvent); 
    68         total += Gauss76Wt[i] * fq * fq; 
     68        total += GAUSS_W[i] * fq * fq; 
    6969    } 
    7070    total *= m; 
    … …  
    7575 
    7676static double 
    77 Iqxy(double qab, double qc, 
     77Iqac(double qab, double qc, 
    7878    double radius_equat_core, 
    7979    double x_core, 

  • sasmodels/models/core_shell_parallelepiped.c

    rc69d6d6 re077231  
     1// Set OVERLAPPING to 1 in order to fill in the edges of the box, with 
     2// c endcaps and b overlapping a.  With the proper choice of parameters, 
     3// (setting rim slds to sld, core sld to solvent, rim thickness to thickness 
     4// and subtracting 2*thickness from length, this should match the hollow 
     5// rectangular prism.)  Set it to 0 for the documented behaviour. 
     6#define OVERLAPPING 0 
    17static double 
    28form_volume(double length_a, double length_b, double length_c, 
    39    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    410{ 
    5     //return length_a * length_b * length_c; 
    6     return length_a * length_b * length_c + 
    7            2.0 * thick_rim_a * length_b * length_c + 
    8            2.0 * thick_rim_b * length_a * length_c + 
    9            2.0 * thick_rim_c * length_a * length_b; 
     11    return 
     12#if OVERLAPPING 
     13        // Hollow rectangular prism only includes the volume of the shell 
     14        // so uncomment the next line when comparing.  Solid rectangular 
     15        // prism, or parallelepiped want filled cores, so comment when 
     16        // comparing. 
     17        //-length_a * length_b * length_c + 
     18        (length_a + 2.0*thick_rim_a) * 
     19        (length_b + 2.0*thick_rim_b) * 
     20        (length_c + 2.0*thick_rim_c); 
     21#else 
     22        length_a * length_b * length_c + 
     23        2.0 * thick_rim_a * length_b * length_c + 
     24        2.0 * length_a * thick_rim_b * length_c + 
     25        2.0 * length_a * length_b * thick_rim_c; 
     26#endif 
    1027} 
    1128 
    … …  
    2441    double thick_rim_c) 
    2542{ 
    26     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
     43    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
    2744    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    28     //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     45    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker) 
     46    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK) 
    2947 
    30     const double mu = 0.5 * q * length_b; 
     48    const double half_q = 0.5*q; 
    3149 
    32     //calculate volume before rescaling (in original code, but not used) 
    33     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
    34     //double vol = length_a * length_b * length_c + 
    35     //       2.0 * thick_rim_a * length_b * length_c + 
    36     //       2.0 * thick_rim_b * length_a * length_c + 
    37     //       2.0 * thick_rim_c * length_a * length_b; 
     50    const double tA = length_a + 2.0*thick_rim_a; 
     51    const double tB = length_b + 2.0*thick_rim_b; 
     52    const double tC = length_c + 2.0*thick_rim_c; 
    3853 
    39     // Scale sides by B 
    40     const double a_scaled = length_a / length_b; 
    41     const double c_scaled = length_c / length_b; 
    42  
    43     double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
    44     double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
    45     double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 
    46  
    47     double Vin = length_a * length_b * length_c; 
    48     //double Vot = (length_a * length_b * length_c + 
    49     //            2.0 * thick_rim_a * length_b * length_c + 
    50     //            2.0 * length_a * thick_rim_b * length_c + 
    51     //            2.0 * length_a * length_b * thick_rim_c); 
    52     double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    53     double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    54     double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present 
    55  
    56     // Scale factors (note that drC is not used later) 
    57     const double drho0 = (core_sld-solvent_sld); 
    58     const double drhoA = (arim_sld-solvent_sld); 
    59     const double drhoB = (brim_sld-solvent_sld); 
    60     const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
    61  
    62  
    63     // Precompute scale factors for combining cross terms from the shape 
    64     const double scale23 = drhoA*V1; 
    65     const double scale14 = drhoB*V2; 
    66     const double scale24 = drhoC*V3; 
    67     const double scale11 = drho0*Vin; 
    68     const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 
     54    // Scale factors 
     55    const double dr0 = (core_sld-solvent_sld); 
     56    const double drA = (arim_sld-solvent_sld); 
     57    const double drB = (brim_sld-solvent_sld); 
     58    const double drC = (crim_sld-solvent_sld); 
    6959 
    7060    // outer integral (with gauss points), integration limits = 0, 1 
    71     double outer_total = 0; //initialize integral 
     61    double outer_sum = 0; //initialize integral 
     62    for( int i=0; i<GAUSS_N; i++) { 
     63        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     64        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    7265 
    73     for( int i=0; i<76; i++) { 
    74         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    75         double mu_proj = mu * sqrt(1.0-sigma*sigma); 
     66        // inner integral (with gauss points), integration limits = 0, pi/2 
     67        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
     68        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
     69        double inner_sum = 0.0; 
     70        for(int j=0; j<GAUSS_N; j++) { 
     71            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     72            double sin_beta, cos_beta; 
     73            SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     74            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
     75            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
     76            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta); 
     77            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta); 
    7678 
    77         // inner integral (with gauss points), integration limits = 0, 1 
    78         double inner_total = 0.0; 
    79         double inner_total_crim = 0.0; 
    80         for(int j=0; j<76; j++) { 
    81             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    82             double sin_uu, cos_uu; 
    83             SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    84             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    85             const double si2 = sas_sinx_x(mu_proj * cos_uu ); 
    86             const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
    87             const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
     79#if OVERLAPPING 
     80            const double f = dr0*siA*siB*siC 
     81                + drA*(siAt-siA)*siB*siC 
     82                + drB*siAt*(siBt-siB)*siC 
     83                + drC*siAt*siBt*(siCt-siC); 
     84#else 
     85            const double f = dr0*siA*siB*siC 
     86                + drA*(siAt-siA)*siB*siC 
     87                + drB*siA*(siBt-siB)*siC 
     88                + drC*siA*siB*(siCt-siC); 
     89#endif 
    8890 
    89             // Expression in libCylinder.c (neither drC nor Vot are used) 
    90             const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
    91             const double form_crim = scale11*si1*si2; 
    92  
    93             //  correct FF : sum of square of phase factors 
    94             inner_total += Gauss76Wt[j] * form * form; 
    95             inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 
     91            inner_sum += GAUSS_W[j] * f * f; 
    9692        } 
    97         inner_total *= 0.5; 
    98         inner_total_crim *= 0.5; 
     93        inner_sum *= 0.5; 
    9994        // now sum up the outer integral 
    100         const double si = sas_sinx_x(mu * c_scaled * sigma); 
    101         const double si_crim = sas_sinx_x(mu * tc * sigma); 
    102         outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 
     95        outer_sum += GAUSS_W[i] * inner_sum; 
    10396    } 
    104     outer_total *= 0.5; 
     97    outer_sum *= 0.5; 
    10598 
    10699    //convert from [1e-12 A-1] to [cm-1] 
    107     return 1.0e-4 * outer_total; 
     100    return 1.0e-4 * outer_sum; 
    108101} 
    109102 
    110103static double 
    111 Iqxy(double qa, double qb, double qc, 
     104Iqabc(double qa, double qb, double qc, 
    112105    double core_sld, 
    113106    double arim_sld, 
    … …  
    128121    const double drC = crim_sld-solvent_sld; 
    129122 
    130     double Vin = length_a * length_b * length_c; 
    131     double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc) 
    132     double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    133     double V3 = 2.0 * length_a * length_b * thick_rim_c; 
    134     // As for the 1D case, Vot is not used 
    135     //double Vot = (length_a * length_b * length_c + 
    136     //              2.0 * thick_rim_a * length_b * length_c + 
    137     //              2.0 * length_a * thick_rim_b * length_c + 
    138     //              2.0 * length_a * length_b * thick_rim_c); 
     123    const double tA = length_a + 2.0*thick_rim_a; 
     124    const double tB = length_b + 2.0*thick_rim_b; 
     125    const double tC = length_c + 2.0*thick_rim_c; 
     126    const double siA = length_a*sas_sinx_x(0.5*length_a*qa); 
     127    const double siB = length_b*sas_sinx_x(0.5*length_b*qb); 
     128    const double siC = length_c*sas_sinx_x(0.5*length_c*qc); 
     129    const double siAt = tA*sas_sinx_x(0.5*tA*qa); 
     130    const double siBt = tB*sas_sinx_x(0.5*tB*qb); 
     131    const double siCt = tC*sas_sinx_x(0.5*tC*qc); 
    139132 
    140     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    141     // the scaling by B. 
    142     double ta = length_a + 2.0*thick_rim_a; 
    143     double tb = length_b + 2.0*thick_rim_b; 
    144     double tc = length_c + 2.0*thick_rim_c; 
    145     //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    146     double siA = sas_sinx_x(0.5*length_a*qa); 
    147     double siB = sas_sinx_x(0.5*length_b*qb); 
    148     double siC = sas_sinx_x(0.5*length_c*qc); 
    149     double siAt = sas_sinx_x(0.5*ta*qa); 
    150     double siBt = sas_sinx_x(0.5*tb*qb); 
    151     double siCt = sas_sinx_x(0.5*tc*qc); 
    152  
    153  
    154     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    155     // in the 1D code, but should be checked! 
    156     double f = ( dr0*siA*siB*siC*Vin 
    157                + drA*(siAt-siA)*siB*siC*V1 
    158                + drB*siA*(siBt-siB)*siC*V2 
    159                + drC*siA*siB*(siCt-siC)*V3); 
     133#if OVERLAPPING 
     134    const double f = dr0*siA*siB*siC 
     135        + drA*(siAt-siA)*siB*siC 
     136        + drB*siAt*(siBt-siB)*siC 
     137        + drC*siAt*siBt*(siCt-siC); 
     138#else 
     139    const double f = dr0*siA*siB*siC 
     140        + drA*(siAt-siA)*siB*siC 
     141        + drB*siA*(siBt-siB)*siC 
     142        + drC*siA*siB*(siCt-siC); 
     143#endif 
    160144 
    161145    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r2d81cfe r97be877  
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    66The thickness and the scattering length density of the shell or 
    7 "rim" can be different on each (pair) of faces. However at this time the 1D 
    8 calculation does **NOT** actually calculate a c face rim despite the presence 
    9 of the parameter. Some other aspects of the 1D calculation may be wrong. 
    10  
    11 .. note:: 
    12    This model was originally ported from NIST IGOR macros. However, it is not 
    13    yet fully understood by the SasView developers and is currently under review. 
     7"rim" can be different on each (pair) of faces. 
    148 
    159The form factor is normalized by the particle volume $V$ such that 
    … …  
    2115where $\langle \ldots \rangle$ is an average over all possible orientations 
    2216of the rectangular solid. 
    23  
    2417 
    2518The function calculated is the form factor of the rectangular solid below. 
    … …  
    4134    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    4235 
    43 **meaning that there are "gaps" at the corners of the solid.**  Again note that 
    44 $t_C = 0$ currently. 
     36**meaning that there are "gaps" at the corners of the solid.** 
    4537 
    4638The intensity calculated follows the :ref:`parallelepiped` model, with the 
    4739core-shell intensity being calculated as the square of the sum of the 
    48 amplitudes of the core and shell, in the same manner as a core-shell model. 
    49  
    50 .. math:: 
    51  
    52     F_{a}(Q,\alpha,\beta)= 
    53     \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta) 
    54                }{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta} 
    55     - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta) 
    56            }{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right] 
    57     \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta) 
    58                }{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right] 
    59     \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta) 
    60                }{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right] 
    61  
    62 .. note:: 
    63  
    64     Why does t_B not appear in the above equation? 
    65     For the calculation of the form factor to be valid, the sides of the solid 
    66     MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
    67     If this inequality is not satisfied, the model will not report an error, 
    68     but the calculation will not be correct and thus the result wrong. 
     40amplitudes of the core and the slabs on the edges. 
     41 
     42the scattering amplitude is computed for a particular orientation of the 
     43core-shell parallelepiped with respect to the scattering vector and then 
     44averaged over all possible orientations, where $\alpha$ is the angle between 
     45the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
     46the angle between projection of the particle in the $xy$ detector plane 
     47and the $y$ axis. 
     48 
     49.. math:: 
     50 
     51    F(Q) 
     52    &= (\rho_\text{core}-\rho_\text{solvent}) 
     53       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
     54    &+ (\rho_\text{A}-\rho_\text{solvent}) 
     55        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     56    &+ (\rho_\text{B}-\rho_\text{solvent}) 
     57        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
     58    &+ (\rho_\text{C}-\rho_\text{solvent}) 
     59        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right] 
     60 
     61with 
     62 
     63.. math:: 
     64 
     65    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     66 
     67and 
     68 
     69.. math:: 
     70 
     71    Q_A &= \sin\alpha \sin\beta \\ 
     72    Q_B &= \sin\alpha \cos\beta \\ 
     73    Q_C &= \cos\alpha 
     74 
     75 
     76where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
     77are the scattering length of the parallelepiped core, and the rectangular 
     78slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
     79is the scattering length of the solvent. 
    6980 
    7081FITTING NOTES 
     82~~~~~~~~~~~~~ 
     83 
    7184If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    72 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    73 However, **no interparticle interference effects are included in this 
    74 calculation.** 
     85value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
     86**no interparticle interference effects are included in this calculation.** 
    7587 
    7688There are many parameters in this model. Hold as many fixed as possible with 
    7789known values, or you will certainly end up at a solution that is unphysical. 
    7890 
    79 Constraints must be applied during fitting to ensure that the inequality 
    80 $A < B < C$ is not violated. The calculation will not report an error, 
    81 but the results will not be correct. 
    82  
    8391The returned value is in units of |cm^-1|, on absolute scale. 
    8492 
    8593NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    8694based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    87 and length $(C+2t_C)$ values, after appropriately 
    88 sorting the three dimensions to give an oblate or prolate particle, to give an 
    89 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     95and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
     96to give an oblate or prolate particle, to give an effective radius, 
     97for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    9098 
    9199For 2d data the orientation of the particle is required, described using 
    92 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    93 of the calculation and angular dispersions see :ref:`orientation` . 
     100angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     101details of the calculation and angular dispersions see :ref:`orientation`. 
    94102The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    95103$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
     104 
     105For 2d, constraints must be applied during fitting to ensure that the 
     106inequality $A < B < C$ is not violated, and hence the correct definition 
     107of angles is preserved. The calculation will not report an error, 
     108but the results may be not correct. 
    96109 
    97110.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    114127    Equations (1), (13-14). (in German) 
    115128.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    116    lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
     129   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
    117130   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    118131   =26379>`_ 
    … …  
    123136* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    124137* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
    125 * **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017 
    126 * **Currently Under review by:** Paul Butler 
     138* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 
     139* Cross-checked against hollow rectangular prism and rectangular prism for 
     140  equal thickness overlapping sides, and by Monte Carlo sampling of points 
     141  within the shape for non-uniform, non-overlapping sides. 
    127142""" 
    128143 
    … …  
    175190        Return equivalent radius (ER) 
    176191    """ 
    177  
    178     # surface average radius (rough approximation) 
    179     surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi) 
    180  
    181     height = length_c + 2.0*thick_rim_c 
    182  
    183     ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad)) 
    184     return 0.5 * (ddd) ** (1. / 3.) 
     192    from .parallelepiped import ER as ER_p 
     193 
     194    a = length_a + 2*thick_rim_a 
     195    b = length_b + 2*thick_rim_b 
     196    c = length_c + 2*thick_rim_c 
     197    return ER_p(a, b, c) 
    185198 
    186199# VR defaults to 1.0 
    … …  
    216229            psi_pd=10, psi_pd_n=1) 
    217230 
    218 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
     231# rkh 7/4/17 add random unit test for 2d, note make all params different, 
     232# 2d values not tested against other codes or models 
    219233if 0:  # pak: model rewrite; need to update tests 
    220234    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 

  • sasmodels/models/cylinder.c

    rbecded3 r108e70e  
    2121 
    2222    double total = 0.0; 
    23     for (int i=0; i<76 ;i++) { 
    24         const double theta = Gauss76Z[i]*zm + zb; 
     23    for (int i=0; i<GAUSS_N ;i++) { 
     24        const double theta = GAUSS_Z[i]*zm + zb; 
    2525        double sin_theta, cos_theta; // slots to hold sincos function output 
    2626        // theta (theta,phi) the projection of the cylinder on the detector plane 
    2727        SINCOS(theta , sin_theta, cos_theta); 
    2828        const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
    29         total += Gauss76Wt[i] * form * form * sin_theta; 
     29        total += GAUSS_W[i] * form * form * sin_theta; 
    3030    } 
    3131    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4545 
    4646static double 
    47 Iqxy(double qab, double qc, 
     47Iqac(double qab, double qc, 
    4848    double sld, 
    4949    double solvent_sld, 

  • sasmodels/models/ellipsoid.c

    rbecded3 r108e70e  
    2222 
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    24     // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     24    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    2525    const double zm = 0.5; 
    2626    const double zb = 0.5; 
    2727    double total = 0.0; 
    28     for (int i=0;i<76;i++) { 
    29         const double u = Gauss76Z[i]*zm + zb; 
     28    for (int i=0;i<GAUSS_N;i++) { 
     29        const double u = GAUSS_Z[i]*zm + zb; 
    3030        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
    3131        const double f = sas_3j1x_x(q*r); 
    32         total += Gauss76Wt[i] * f * f; 
     32        total += GAUSS_W[i] * f * f; 
    3333    } 
    3434    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3939 
    4040static double 
    41 Iqxy(double qab, double qc, 
     41Iqac(double qab, double qc, 
    4242    double sld, 
    4343    double sld_solvent, 

  • sasmodels/models/elliptical_cylinder.c

    r82592da r108e70e  
    2222    //initialize integral 
    2323    double outer_sum = 0.0; 
    24     for(int i=0;i<76;i++) { 
     24    for(int i=0;i<GAUSS_N;i++) { 
    2525        //setup inner integral over the ellipsoidal cross-section 
    26         const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     26        const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    2727        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    2828        //const double arg = radius_minor*sin_val; 
    2929        double inner_sum=0; 
    30         for(int j=0;j<76;j++) { 
    31             //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017 
    32             const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     30        for(int j=0;j<GAUSS_N;j++) { 
     31            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3332            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3433            const double be = sas_2J1x_x(q*r); 
    35             inner_sum += Gauss76Wt[j] * be * be; 
     34            inner_sum += GAUSS_W[j] * be * be; 
    3635        } 
    3736        //now calculate the value of the inner integral 
    … …  
    4039        //now calculate outer integral 
    4140        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    42         outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
     41        outer_sum += GAUSS_W[i] * inner_sum * si * si; 
    4342    } 
    4443    outer_sum *= 0.5*(vb-va); 
    … …  
    5554 
    5655static double 
    57 Iqxy(double qa, double qb, double qc, 
     56Iqabc(double qa, double qb, double qc, 
    5857     double radius_minor, double r_ratio, double length, 
    5958     double sld, double solvent_sld) 

  • sasmodels/models/elliptical_cylinder.py

    r2d81cfe r2d81cfe  
    121121# pylint: enable=bad-whitespace, line-too-long 
    122122 
    123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", 
    124           "elliptical_cylinder.c"] 
     123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    125124 
    126125demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 

  • sasmodels/models/fcc_paracrystal.c

    rf728001 r108e70e  
    5353 
    5454    double outer_sum = 0.0; 
    55     for(int i=0; i<150; i++) { 
     55    for(int i=0; i<GAUSS_N; i++) { 
    5656        double inner_sum = 0.0; 
    57         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     57        const double theta = GAUSS_Z[i]*theta_m + theta_b; 
    5858        double sin_theta, cos_theta; 
    5959        SINCOS(theta, sin_theta, cos_theta); 
    6060        const double qc = q*cos_theta; 
    6161        const double qab = q*sin_theta; 
    62         for(int j=0;j<150;j++) { 
    63             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     62        for(int j=0;j<GAUSS_N;j++) { 
     63            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    6464            double sin_phi, cos_phi; 
    6565            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6767            const double qb = qab*sin_phi; 
    6868            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
    69             inner_sum += Gauss150Wt[j] * form; 
     69            inner_sum += GAUSS_W[j] * form; 
    7070        } 
    7171        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    72         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     72        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    7373    } 
    7474    outer_sum *= theta_m; 
    … …  
    8080 
    8181 
    82 static double Iqxy(double qa, double qb, double qc, 
     82static double Iqabc(double qa, double qb, double qc, 
    8383    double dnn, double d_factor, double radius, 
    8484    double sld, double solvent_sld) 

  • sasmodels/models/flexible_cylinder_elliptical.c

    r592343f r74768cb  
    1717    double sum=0.0; 
    1818 
    19     for(int i=0;i<N_POINTS_76;i++) { 
    20         const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4; 
     19    for(int i=0;i<GAUSS_N;i++) { 
     20        const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4; 
    2121        double sn, cn; 
    2222        SINCOS(zi, sn, cn); 
    2323        const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 
    2424        const double yyy = sas_2J1x_x(arg); 
    25         sum += Gauss76Wt[i] * yyy * yyy; 
     25        sum += GAUSS_W[i] * yyy * yyy; 
    2626    } 
    2727    sum *= 0.5; 

  • sasmodels/models/hollow_cylinder.c

    rbecded3 r108e70e  
    3838 
    3939    double summ = 0.0;            //initialize intergral 
    40     for (int i=0;i<76;i++) { 
    41         const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     40    for (int i=0;i<GAUSS_N;i++) { 
     41        const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 
    4242        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    4343        const double form = _fq(q*sin_theta, q*cos_theta, 
    4444                                radius, thickness, length); 
    45         summ += Gauss76Wt[i] * form * form; 
     45        summ += GAUSS_W[i] * form * form; 
    4646    } 
    4747 
    … …  
    5252 
    5353static double 
    54 Iqxy(double qab, double qc, 
     54Iqac(double qab, double qc, 
    5555    double radius, double thickness, double length, 
    5656    double sld, double solvent_sld) 

  • sasmodels/models/hollow_rectangular_prism.c

    r1e7b0db0 r108e70e  
    11double 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,  
     2double Iq(double q, double sld, double solvent_sld, double length_a, 
    33          double b2a_ratio, double c2a_ratio, double thickness); 
    44 
    … …  
    3737    const double v2a = 0.0; 
    3838    const double v2b = M_PI_2;  //phi integration limits 
    39      
     39 
    4040    double outer_sum = 0.0; 
    41     for(int i=0; i<76; i++) { 
     41    for(int i=0; i<GAUSS_N; i++) { 
    4242 
    43         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     43        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    4444        double sin_theta, cos_theta; 
    4545        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    4949 
    5050        double inner_sum = 0.0; 
    51         for(int j=0; j<76; j++) { 
     51        for(int j=0; j<GAUSS_N; j++) { 
    5252 
    53             const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     53            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    5454            double sin_phi, cos_phi; 
    5555            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6666            const double AP2 = vol_core * termA2 * termB2 * termC2; 
    6767 
    68             inner_sum += Gauss76Wt[j] * square(AP1-AP2); 
     68            inner_sum += GAUSS_W[j] * square(AP1-AP2); 
    6969        } 
    7070        inner_sum *= 0.5 * (v2b-v2a); 
    7171 
    72         outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 
     72        outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
    7373    } 
    7474    outer_sum *= 0.5*(v1b-v1a); 
    … …  
    8484    return 1.0e-4 * delrho * delrho * form; 
    8585} 
     86 
     87double Iqabc(double qa, double qb, double qc, 
     88    double sld, 
     89    double solvent_sld, 
     90    double length_a, 
     91    double b2a_ratio, 
     92    double c2a_ratio, 
     93    double thickness) 
     94{ 
     95    const double length_b = length_a * b2a_ratio; 
     96    const double length_c = length_a * c2a_ratio; 
     97    const double a_half = 0.5 * length_a; 
     98    const double b_half = 0.5 * length_b; 
     99    const double c_half = 0.5 * length_c; 
     100    const double vol_total = length_a * length_b * length_c; 
     101    const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness); 
     102 
     103    // Amplitude AP from eqn. (13) 
     104 
     105    const double termA1 = sas_sinx_x(qa * a_half); 
     106    const double termA2 = sas_sinx_x(qa * (a_half-thickness)); 
     107 
     108    const double termB1 = sas_sinx_x(qb * b_half); 
     109    const double termB2 = sas_sinx_x(qb * (b_half-thickness)); 
     110 
     111    const double termC1 = sas_sinx_x(qc * c_half); 
     112    const double termC2 = sas_sinx_x(qc * (c_half-thickness)); 
     113 
     114    const double AP1 = vol_total * termA1 * termB1 * termC1; 
     115    const double AP2 = vol_core * termA2 * termB2 * termC2; 
     116 
     117    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
     118    const double delrho = sld - solvent_sld; 
     119 
     120    // Convert from [1e-12 A-1] to [cm-1] 
     121    return 1.0e-4 * square(delrho * (AP1-AP2)); 
     122} 

  • sasmodels/models/hollow_rectangular_prism.py

    r2d81cfe r2d81cfe  
    55This model provides the form factor, $P(q)$, for a hollow rectangular 
    66parallelepiped with a wall of thickness $\Delta$. 
    7 It computes only the 1D scattering, not the 2D. 
     7 
    88 
    99Definition 
    … …  
    6666(which is unitless). 
    6767 
    68 **The 2D scattering intensity is not computed by this model.** 
     68For 2d data the orientation of the particle is required, described using 
     69angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     70of the calculation and angular dispersions see :ref:`orientation` . 
     71The angle $\Psi$ is the rotational angle around the long *C* axis. For example, 
     72$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector. 
     73 
     74For 2d, constraints must be applied during fitting to ensure that the inequality 
     75$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 
     76but the results may be not correct. 
     77 
     78.. figure:: img/parallelepiped_angle_definition.png 
     79 
     80    Definition of the angles for oriented hollow rectangular prism. 
     81    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 
     82    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the prism. 
     83    The neutron or X-ray beam is along the $z$ axis. 
     84 
     85.. figure:: img/parallelepiped_angle_projection.png 
     86 
     87    Examples of the angles for oriented hollow rectangular prisms against the 
     88    detector plane. 
    6989 
    7090 
    … …  
    113133              ["thickness", "Ang", 1, [0, inf], "volume", 
    114134               "Thickness of parallelepiped"], 
     135              ["theta", "degrees", 0, [-360, 360], "orientation", 
     136               "c axis to beam angle"], 
     137              ["phi", "degrees", 0, [-360, 360], "orientation", 
     138               "rotation about beam"], 
     139              ["psi", "degrees", 0, [-360, 360], "orientation", 
     140               "rotation about c axis"], 
    115141             ] 
    116142 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rab2aea8 r74768cb  
    11double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a,  
     2double Iq(double q, double sld, double solvent_sld, double length_a, 
    33          double b2a_ratio, double c2a_ratio); 
    44 
    … …  
    2929    const double v2a = 0.0; 
    3030    const double v2b = M_PI_2;  //phi integration limits 
    31      
     31 
    3232    double outer_sum = 0.0; 
    33     for(int i=0; i<76; i++) { 
    34         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     33    for(int i=0; i<GAUSS_N; i++) { 
     34        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    3535 
    3636        double sin_theta, cos_theta; 
    … …  
    4444 
    4545        double inner_sum = 0.0; 
    46         for(int j=0; j<76; j++) { 
    47             const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     46        for(int j=0; j<GAUSS_N; j++) { 
     47            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    4848 
    4949            double sin_phi, cos_phi; 
    … …  
    6262                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 
    6363 
    64             inner_sum += Gauss76Wt[j] * square(AL+AT); 
     64            inner_sum += GAUSS_W[j] * square(AL+AT); 
    6565        } 
    6666 
    6767        inner_sum *= 0.5 * (v2b-v2a); 
    68         outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
     68        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    6969    } 
    7070 

  • sasmodels/models/lib/gauss150.c

    r994d77f r99b84ec  
    1 /* 
    2  *  GaussWeights.c 
    3  *  SANSAnalysis 
    4  * 
    5  *  Created by Andrew Jackson on 4/23/07. 
    6  *  Copyright 2007 __MyCompanyName__. All rights reserved. 
    7  * 
    8  */ 
     1// Created by Andrew Jackson on 4/23/07 
    92 
    10 // Gaussians 
    11 constant double Gauss150Z[150]={ 
     3 #ifdef GAUSS_N 
     4 # undef GAUSS_N 
     5 # undef GAUSS_Z 
     6 # undef GAUSS_W 
     7 #endif 
     8 #define GAUSS_N 150 
     9 #define GAUSS_Z Gauss150Z 
     10 #define GAUSS_W Gauss150Wt 
     11 
     12 
     13// Note: using array size 152 rather than 150 so that it is a multiple of 4. 
     14// Some OpenCL devices prefer that vectors start and end on nice boundaries. 
     15constant double Gauss150Z[152]={ 
    1216        -0.9998723404457334, 
    1317        -0.9993274305065947, 
    … …  
    159163        0.9983473449340834, 
    160164        0.9993274305065947, 
    161         0.9998723404457334 
     165        0.9998723404457334, 
     166        0., // zero padding is ignored 
     167        0.  // zero padding is ignored 
    162168}; 
    163169 
    164 constant double Gauss150Wt[150]={ 
     170constant double Gauss150Wt[152]={ 
    165171        0.0003276086705538, 
    166172        0.0007624720924706, 
    … …  
    312318        0.0011976474864367, 
    313319        0.0007624720924706, 
    314         0.0003276086705538 
     320        0.0003276086705538, 
     321        0., // zero padding is ignored 
     322        0.  // zero padding is ignored 
    315323}; 

  • sasmodels/models/lib/gauss20.c

    r994d77f r99b84ec  
    1 /* 
    2  *  GaussWeights.c 
    3  *  SANSAnalysis 
    4  * 
    5  *  Created by Andrew Jackson on 4/23/07. 
    6  *  Copyright 2007 __MyCompanyName__. All rights reserved. 
    7  * 
    8  */ 
     1// Created by Andrew Jackson on 4/23/07 
     2 
     3 #ifdef GAUSS_N 
     4 # undef GAUSS_N 
     5 # undef GAUSS_Z 
     6 # undef GAUSS_W 
     7 #endif 
     8 #define GAUSS_N 20 
     9 #define GAUSS_Z Gauss20Z 
     10 #define GAUSS_W Gauss20Wt 
    911 
    1012// Gaussians 

  • sasmodels/models/lib/gauss76.c

    r66d119f r99b84ec  
    1 /* 
    2  *  GaussWeights.c 
    3  *  SANSAnalysis 
    4  * 
    5  *  Created by Andrew Jackson on 4/23/07. 
    6  *  Copyright 2007 __MyCompanyName__. All rights reserved. 
    7  * 
    8  */ 
    9 #define N_POINTS_76 76 
     1// Created by Andrew Jackson on 4/23/07 
     2 
     3 #ifdef GAUSS_N 
     4 # undef GAUSS_N 
     5 # undef GAUSS_Z 
     6 # undef GAUSS_W 
     7 #endif 
     8 #define GAUSS_N 76 
     9 #define GAUSS_Z Gauss76Z 
     10 #define GAUSS_W Gauss76Wt 
    1011 
    1112// Gaussians 
    12 constant double Gauss76Wt[N_POINTS_76]={ 
     13constant double Gauss76Wt[76]={ 
    1314        .00126779163408536,             //0 
    1415        .00294910295364247, 
    … …  
    8990}; 
    9091 
    91 constant double Gauss76Z[N_POINTS_76]={ 
     92constant double Gauss76Z[76]={ 
    9293        -.999505948362153,              //0 
    9394        -.997397786355355, 

  • sasmodels/models/line.py

    r2d81cfe r108e70e  
    5757Iq.vectorized = True # Iq accepts an array of q values 
    5858 
     59 
    5960def Iqxy(qx, qy, *args): 
    6061    """ 
    … …  
    6970 
    7071Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values 
     72 
     73# uncomment the following to test Iqxy in C models 
     74#del Iq, Iqxy 
     75#c_code = """ 
     76#static double Iq(double q, double b, double m) { return m*q+b; } 
     77#static double Iqxy(double qx, double qy, double b, double m) 
     78#{ return (m*qx+b)*(m*qy+b); } 
     79#""" 
    7180 
    7281def random(): 

  • sasmodels/models/parallelepiped.c

    r9b7b23f r108e70e  
    2323    double outer_total = 0; //initialize integral 
    2424 
    25     for( int i=0; i<76; i++) { 
    26         const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     25    for( int i=0; i<GAUSS_N; i++) { 
     26        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    2727        const double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    2828 
    … …  
    3030        // corresponding to angles from 0 to pi/2. 
    3131        double inner_total = 0.0; 
    32         for(int j=0; j<76; j++) { 
    33             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     32        for(int j=0; j<GAUSS_N; j++) { 
     33            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    3434            double sin_uu, cos_uu; 
    3535            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    3636            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    3737            const double si2 = sas_sinx_x(mu_proj * cos_uu); 
    38             inner_total += Gauss76Wt[j] * square(si1 * si2); 
     38            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
    4040        inner_total *= 0.5; 
    4141 
    4242        const double si = sas_sinx_x(mu * c_scaled * sigma); 
    43         outer_total += Gauss76Wt[i] * inner_total * si * si; 
     43        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4444    } 
    4545    outer_total *= 0.5; 
    … …  
    5353 
    5454static double 
    55 Iqxy(double qa, double qb, double qc, 
     55Iqabc(double qa, double qb, double qc, 
    5656    double sld, 
    5757    double solvent_sld, 

  • sasmodels/models/pringle.c

    r1e7b0db0 r74768cb  
    2929    double sumC = 0.0;          // initialize integral 
    3030    double r; 
    31     for (int i=0; i < 76; i++) { 
    32         r = Gauss76Z[i]*zm + zb; 
     31    for (int i=0; i < GAUSS_N; i++) { 
     32        r = GAUSS_Z[i]*zm + zb; 
    3333 
    3434        const double qrs = r*q_sin_psi; 
    3535        const double qrrc = r*r*q_cos_psi; 
    3636 
    37         double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
     37        double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
    3838        double S, C; 
    3939        SINCOS(alpha*qrrc, S, C); 
    … …  
    8686 
    8787    double sum = 0.0; 
    88     for (int i = 0; i < 76; i++) { 
    89         double psi = Gauss76Z[i]*zm + zb; 
     88    for (int i = 0; i < GAUSS_N; i++) { 
     89        double psi = GAUSS_Z[i]*zm + zb; 
    9090        double sin_psi, cos_psi; 
    9191        SINCOS(psi, sin_psi, cos_psi); 
    … …  
    9393        double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0)); 
    9494        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term; 
    95         sum += Gauss76Wt[i] * pringle_kernel; 
     95        sum += GAUSS_W[i] * pringle_kernel; 
    9696    } 
    9797 

  • sasmodels/models/rectangular_prism.c

    r1e7b0db0 r108e70e  
    11double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a,  
     2double Iq(double q, double sld, double solvent_sld, double length_a, 
    33          double b2a_ratio, double c2a_ratio); 
    44 
    … …  
    2626    const double v2a = 0.0; 
    2727    const double v2b = M_PI_2;  //phi integration limits 
    28      
     28 
    2929    double outer_sum = 0.0; 
    30     for(int i=0; i<76; i++) { 
    31         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     30    for(int i=0; i<GAUSS_N; i++) { 
     31        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    3232        double sin_theta, cos_theta; 
    3333        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    3636 
    3737        double inner_sum = 0.0; 
    38         for(int j=0; j<76; j++) { 
    39             double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     38        for(int j=0; j<GAUSS_N; j++) { 
     39            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    4040            double sin_phi, cos_phi; 
    4141            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    4545            const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 
    4646            const double AP = termA * termB * termC; 
    47             inner_sum += Gauss76Wt[j] * AP * AP; 
     47            inner_sum += GAUSS_W[j] * AP * AP; 
    4848        } 
    4949        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
    50         outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
     50        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    5151    } 
    5252 
    5353    double answer = 0.5*(v1b-v1a)*outer_sum; 
    5454 
    55     // Normalize by Pi (Eqn. 16).  
    56     // The term (ABC)^2 does not appear because it was introduced before on  
     55    // Normalize by Pi (Eqn. 16). 
     56    // The term (ABC)^2 does not appear because it was introduced before on 
    5757    // the definitions of termA, termB, termC. 
    58     // The factor 2 appears because the theta integral has been defined between  
     58    // The factor 2 appears because the theta integral has been defined between 
    5959    // 0 and pi/2, instead of 0 to pi. 
    6060    answer /= M_PI_2; //Form factor P(q) 
    … …  
    6464    answer *= square((sld-solvent_sld)*volume); 
    6565 
    66     // Convert from [1e-12 A-1] to [cm-1]  
     66    // Convert from [1e-12 A-1] to [cm-1] 
    6767    answer *= 1.0e-4; 
    6868 
    6969    return answer; 
    7070} 
     71 
     72 
     73double Iqabc(double qa, double qb, double qc, 
     74    double sld, 
     75    double solvent_sld, 
     76    double length_a, 
     77    double b2a_ratio, 
     78    double c2a_ratio) 
     79{ 
     80    const double length_b = length_a * b2a_ratio; 
     81    const double length_c = length_a * c2a_ratio; 
     82    const double a_half = 0.5 * length_a; 
     83    const double b_half = 0.5 * length_b; 
     84    const double c_half = 0.5 * length_c; 
     85    const double volume = length_a * length_b * length_c; 
     86 
     87    // Amplitude AP from eqn. (13) 
     88 
     89    const double termA = sas_sinx_x(qa * a_half); 
     90    const double termB = sas_sinx_x(qb * b_half); 
     91    const double termC = sas_sinx_x(qc * c_half); 
     92 
     93    const double AP = termA * termB * termC; 
     94 
     95    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
     96    const double delrho = sld - solvent_sld; 
     97 
     98    // Convert from [1e-12 A-1] to [cm-1] 
     99    return 1.0e-4 * square(volume * delrho * AP); 
     100} 

  • sasmodels/models/rectangular_prism.py

    r2d81cfe r2d81cfe  
    1212the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity 
    1313to *a* will generate a distribution of cubes of different sizes). 
    14 Note also that, contrary to :ref:`parallelepiped`, it does not compute 
    15 the 2D scattering. 
    1614 
    1715 
    … …  
    2624that reference), with $\theta$ corresponding to $\alpha$ in that paper, 
    2725and not to the usual convention used for example in the 
    28 :ref:`parallelepiped` model. As the present model does not compute 
    29 the 2D scattering, this has no further consequences. 
     26:ref:`parallelepiped` model. 
    3027 
    3128In this model the scattering from a massive parallelepiped with an 
    … …  
    6562units) *scale* represents the volume fraction (which is unitless). 
    6663 
    67 **The 2D scattering intensity is not computed by this model.** 
     64For 2d data the orientation of the particle is required, described using 
     65angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     66of the calculation and angular dispersions see :ref:`orientation` . 
     67The angle $\Psi$ is the rotational angle around the long *C* axis. For example, 
     68$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector. 
     69 
     70For 2d, constraints must be applied during fitting to ensure that the inequality 
     71$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 
     72but the results may be not correct. 
     73 
     74.. figure:: img/parallelepiped_angle_definition.png 
     75 
     76    Definition of the angles for oriented core-shell parallelepipeds. 
     77    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 
     78    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder. 
     79    The neutron or X-ray beam is along the $z$ axis. 
     80 
     81.. figure:: img/parallelepiped_angle_projection.png 
     82 
     83    Examples of the angles for oriented rectangular prisms against the 
     84    detector plane. 
     85 
    6886 
    6987 
    … …  
    108126              ["c2a_ratio", "", 1, [0, inf], "volume", 
    109127               "Ratio sides c/a"], 
     128              ["theta", "degrees", 0, [-360, 360], "orientation", 
     129               "c axis to beam angle"], 
     130              ["phi", "degrees", 0, [-360, 360], "orientation", 
     131               "rotation about beam"], 
     132              ["psi", "degrees", 0, [-360, 360], "orientation", 
     133               "rotation about c axis"], 
    110134             ] 
    111135 

  • sasmodels/models/sc_paracrystal.c

    rf728001 r108e70e  
    5454 
    5555    double outer_sum = 0.0; 
    56     for(int i=0; i<150; i++) { 
     56    for(int i=0; i<GAUSS_N; i++) { 
    5757        double inner_sum = 0.0; 
    58         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     58        const double theta = GAUSS_Z[i]*theta_m + theta_b; 
    5959        double sin_theta, cos_theta; 
    6060        SINCOS(theta, sin_theta, cos_theta); 
    6161        const double qc = q*cos_theta; 
    6262        const double qab = q*sin_theta; 
    63         for(int j=0;j<150;j++) { 
    64             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     63        for(int j=0;j<GAUSS_N;j++) { 
     64            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    6565            double sin_phi, cos_phi; 
    6666            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6868            const double qb = qab*sin_phi; 
    6969            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
    70             inner_sum += Gauss150Wt[j] * form; 
     70            inner_sum += GAUSS_W[j] * form; 
    7171        } 
    7272        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    73         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     73        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    7474    } 
    7575    outer_sum *= theta_m; 
    … …  
    8282 
    8383static double 
    84 Iqxy(double qa, double qb, double qc, 
     84Iqabc(double qa, double qb, double qc, 
    8585    double dnn, double d_factor, double radius, 
    8686    double sld, double solvent_sld) 

  • sasmodels/models/stacked_disks.c

    rbecded3 r108e70e  
    8181    double halfheight = 0.5*thick_core; 
    8282 
    83     for(int i=0; i<N_POINTS_76; i++) { 
    84         double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
     83    for(int i=0; i<GAUSS_N; i++) { 
     84        double zi = (GAUSS_Z[i] + 1.0)*M_PI_4; 
    8585        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8686        SINCOS(zi, sin_alpha, cos_alpha); 
    … …  
    9595                           solvent_sld, 
    9696                           d); 
    97         summ += Gauss76Wt[i] * yyy * sin_alpha; 
     97        summ += GAUSS_W[i] * yyy * sin_alpha; 
    9898    } 
    9999 
    … …  
    142142 
    143143static double 
    144 Iqxy(double qab, double qc, 
     144Iqac(double qab, double qc, 
    145145    double thick_core, 
    146146    double thick_layer, 

  • sasmodels/models/triaxial_ellipsoid.c

    rbecded3 r108e70e  
    2121    const double zb = M_PI_4; 
    2222    double outer = 0.0; 
    23     for (int i=0;i<76;i++) { 
    24         //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2; 
    25         const double phi = Gauss76Z[i]*zm + zb; 
     23    for (int i=0;i<GAUSS_N;i++) { 
     24        //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2; 
     25        const double phi = GAUSS_Z[i]*zm + zb; 
    2626        const double pa_sinsq_phi = pa*square(sin(phi)); 
    2727 
    … …  
    2929        const double um = 0.5; 
    3030        const double ub = 0.5; 
    31         for (int j=0;j<76;j++) { 
     31        for (int j=0;j<GAUSS_N;j++) { 
    3232            // translate a point in [-1,1] to a point in [0, 1] 
    33             const double usq = square(Gauss76Z[j]*um + ub); 
     33            const double usq = square(GAUSS_Z[j]*um + ub); 
    3434            const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 
    3535            const double fq = sas_3j1x_x(q*r); 
    36             inner += Gauss76Wt[j] * fq * fq; 
     36            inner += GAUSS_W[j] * fq * fq; 
    3737        } 
    38         outer += Gauss76Wt[i] * inner;  // correcting for dx later 
     38        outer += GAUSS_W[i] * inner;  // correcting for dx later 
    3939    } 
    4040    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    … …  
    4646 
    4747static double 
    48 Iqxy(double qa, double qb, double qc, 
     48Iqabc(double qa, double qb, double qc, 
    4949    double sld, 
    5050    double sld_solvent, 

  • sasmodels/special.py

    re65c3ba rdf69efa  
    33................. 
    44 
    5 The C code follows the C99 standard, with the usual math functions, 
    6 as defined in 
    7 `OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_. 
    8 This includes the following: 
     5This following standard C99 math functions are available: 
    96 
    107    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: 
    118        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ 
    129 
    13     exp, log, pow(x,y), expm1, sqrt: 
    14         Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$. 
    15         The function expm1(x) is accurate across all $x$, including $x$ 
    16         very close to zero. 
     10    exp, log, pow(x,y), expm1, log1p, sqrt, cbrt: 
     11        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$, 
     12        $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x) 
     13        are accurate across all $x$, including $x$ very close to zero. 
    1714 
    1815    sin, cos, tan, asin, acos, atan: 
    … …  
    2926        quadrants II and IV when $x$ and $y$ have opposite sign. 
    3027 
    31     fmin(x,y), fmax(x,y), trunc, rint: 
     28    fabs(x), fmin(x,y), fmax(x,y), trunc, rint: 
    3229        Floating point functions.  rint(x) returns the nearest integer. 
    3330 
    … …  
    3532        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that 
    3633        you cannot use :code:`x == NAN` to test for NaN values since that 
    37         will always return false.  NAN does not equal NAN! 
     34        will always return false.  NAN does not equal NAN!  The alternative, 
     35        :code:`x != x` may fail if the compiler optimizes the test away. 
    3836 
    3937    INFINITY: 
    … …  
    8987        for forcing a constant to stay double precision. 
    9088 
    91 The following special functions and scattering calculations are defined in 
    92 `sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_. 
     89The following special functions and scattering calculations are defined. 
    9390These functions have been tuned to be fast and numerically stable down 
    9491to $q=0$ even in single precision.  In some cases they work around bugs 
    … …  
    184181 
    185182 
    186     Gauss76Z[i], Gauss76Wt[i]: 
     183    gauss76.n, gauss76.z[i], gauss76.w[i]: 
    187184        Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively, 
    188185        computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$. 
    189  
    190         Similar arrays are available in :code:`gauss20.c` for 20-point 
    191         quadrature and in :code:`gauss150.c` for 150-point quadrature. 
    192  
     186        When translating the model to C, include 'lib/gauss76.c' in the source 
     187        and use :code:`GAUSS_N`, :code:`GAUSS_Z`, and :code:`GAUSS_W`. 
     188 
     189        Similar arrays are available in :code:`gauss20` for 20-point quadrature 
     190        and :code:`gauss150.c` for 150-point quadrature. By using 
     191        :code:`import gauss76 as gauss` it is easy to change the number of 
     192        points in the integration. 
    193193""" 
    194194# pylint: disable=unused-import 
    … …  
    200200 
    201201# C99 standard math library functions 
    202 from numpy import exp, log, power as pow, expm1, sqrt 
     202from numpy import exp, log, power as pow, expm1, log1p, sqrt, cbrt 
    203203from numpy import sin, cos, tan, arcsin as asin, arccos as acos, arctan as atan 
    204204from numpy import sinh, cosh, tanh, arcsinh as asinh, arccosh as acosh, arctanh as atanh 
    205205from numpy import arctan2 as atan2 
    206 from numpy import fmin, fmax, trunc, rint 
    207 from numpy import NAN, inf as INFINITY 
    208  
     206from numpy import fabs, fmin, fmax, trunc, rint 
     207from numpy import pi, nan, inf 
    209208from scipy.special import gamma as sas_gamma 
    210209from scipy.special import erf as sas_erf 
    … …  
    218217# C99 standard math constants 
    219218M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E = np.pi, np.pi/2, np.pi/4, np.sqrt(0.5), np.e 
     219NAN = nan 
     220INFINITY = inf 
    220221 
    221222# non-standard constants 
    … …  
    226227    """return sin(x), cos(x)""" 
    227228    return sin(x), cos(x) 
     229sincos = SINCOS 
    228230 
    229231def square(x): 
    … …  
    294296 
    295297# Gaussians 
    296  
    297 Gauss20Wt = np.array([ 
    298     .0176140071391521, 
    299     .0406014298003869, 
    300     .0626720483341091, 
    301     .0832767415767047, 
    302     .10193011981724, 
    303     .118194531961518, 
    304     .131688638449177, 
    305     .142096109318382, 
    306     .149172986472604, 
    307     .152753387130726, 
    308     .152753387130726, 
    309     .149172986472604, 
    310     .142096109318382, 
    311     .131688638449177, 
    312     .118194531961518, 
    313     .10193011981724, 
    314     .0832767415767047, 
    315     .0626720483341091, 
    316     .0406014298003869, 
    317     .0176140071391521 
    318 ]) 
    319  
    320 Gauss20Z = np.array([ 
    321     -.993128599185095, 
    322     -.963971927277914, 
    323     -.912234428251326, 
    324     -.839116971822219, 
    325     -.746331906460151, 
    326     -.636053680726515, 
    327     -.510867001950827, 
    328     -.37370608871542, 
    329     -.227785851141645, 
    330     -.076526521133497, 
    331     .0765265211334973, 
    332     .227785851141645, 
    333     .37370608871542, 
    334     .510867001950827, 
    335     .636053680726515, 
    336     .746331906460151, 
    337     .839116971822219, 
    338     .912234428251326, 
    339     .963971927277914, 
    340     .993128599185095 
    341 ]) 
    342  
    343 Gauss76Wt = np.array([ 
    344     .00126779163408536,         #0 
    345     .00294910295364247, 
    346     .00462793522803742, 
    347     .00629918049732845, 
    348     .00795984747723973, 
    349     .00960710541471375, 
    350     .0112381685696677, 
    351     .0128502838475101, 
    352     .0144407317482767, 
    353     .0160068299122486, 
    354     .0175459372914742,          #10 
    355     .0190554584671906, 
    356     .020532847967908, 
    357     .0219756145344162, 
    358     .0233813253070112, 
    359     .0247476099206597, 
    360     .026072164497986, 
    361     .0273527555318275, 
    362     .028587223650054, 
    363     .029773487255905, 
    364     .0309095460374916,          #20 
    365     .0319934843404216, 
    366     .0330234743977917, 
    367     .0339977794120564, 
    368     .0349147564835508, 
    369     .0357728593807139, 
    370     .0365706411473296, 
    371     .0373067565423816, 
    372     .0379799643084053, 
    373     .0385891292645067, 
    374     .0391332242205184,          #30 
    375     .0396113317090621, 
    376     .0400226455325968, 
    377     .040366472122844, 
    378     .0406422317102947, 
    379     .0408494593018285, 
    380     .040987805464794, 
    381     .0410570369162294, 
    382     .0410570369162294, 
    383     .040987805464794, 
    384     .0408494593018285,          #40 
    385     .0406422317102947, 
    386     .040366472122844, 
    387     .0400226455325968, 
    388     .0396113317090621, 
    389     .0391332242205184, 
    390     .0385891292645067, 
    391     .0379799643084053, 
    392     .0373067565423816, 
    393     .0365706411473296, 
    394     .0357728593807139,          #50 
    395     .0349147564835508, 
    396     .0339977794120564, 
    397     .0330234743977917, 
    398     .0319934843404216, 
    399     .0309095460374916, 
    400     .029773487255905, 
    401     .028587223650054, 
    402     .0273527555318275, 
    403     .026072164497986, 
    404     .0247476099206597,          #60 
    405     .0233813253070112, 
    406     .0219756145344162, 
    407     .020532847967908, 
    408     .0190554584671906, 
    409     .0175459372914742, 
    410     .0160068299122486, 
    411     .0144407317482767, 
    412     .0128502838475101, 
    413     .0112381685696677, 
    414     .00960710541471375,         #70 
    415     .00795984747723973, 
    416     .00629918049732845, 
    417     .00462793522803742, 
    418     .00294910295364247, 
    419     .00126779163408536          #75 (indexed from 0) 
    420 ]) 
    421  
    422 Gauss76Z = np.array([ 
    423     -.999505948362153,          #0 
    424     -.997397786355355, 
    425     -.993608772723527, 
    426     -.988144453359837, 
    427     -.981013938975656, 
    428     -.972229228520377, 
    429     -.961805126758768, 
    430     -.949759207710896, 
    431     -.936111781934811, 
    432     -.92088586125215, 
    433     -.904107119545567,          #10 
    434     -.885803849292083, 
    435     -.866006913771982, 
    436     -.844749694983342, 
    437     -.822068037328975, 
    438     -.7980001871612, 
    439     -.77258672828181, 
    440     -.74587051350361, 
    441     -.717896592387704, 
    442     -.688712135277641, 
    443     -.658366353758143,          #20 
    444     -.626910417672267, 
    445     -.594397368836793, 
    446     -.560882031601237, 
    447     -.526420920401243, 
    448     -.491072144462194, 
    449     -.454895309813726, 
    450     -.417951418780327, 
    451     -.380302767117504, 
    452     -.342012838966962, 
    453     -.303146199807908,          #30 
    454     -.263768387584994, 
    455     -.223945802196474, 
    456     -.183745593528914, 
    457     -.143235548227268, 
    458     -.102483975391227, 
    459     -.0615595913906112, 
    460     -.0205314039939986, 
    461     .0205314039939986, 
    462     .0615595913906112, 
    463     .102483975391227,                   #40 
    464     .143235548227268, 
    465     .183745593528914, 
    466     .223945802196474, 
    467     .263768387584994, 
    468     .303146199807908, 
    469     .342012838966962, 
    470     .380302767117504, 
    471     .417951418780327, 
    472     .454895309813726, 
    473     .491072144462194,           #50 
    474     .526420920401243, 
    475     .560882031601237, 
    476     .594397368836793, 
    477     .626910417672267, 
    478     .658366353758143, 
    479     .688712135277641, 
    480     .717896592387704, 
    481     .74587051350361, 
    482     .77258672828181, 
    483     .7980001871612,     #60 
    484     .822068037328975, 
    485     .844749694983342, 
    486     .866006913771982, 
    487     .885803849292083, 
    488     .904107119545567, 
    489     .92088586125215, 
    490     .936111781934811, 
    491     .949759207710896, 
    492     .961805126758768, 
    493     .972229228520377,           #70 
    494     .981013938975656, 
    495     .988144453359837, 
    496     .993608772723527, 
    497     .997397786355355, 
    498     .999505948362153            #75 
    499 ]) 
    500  
    501 Gauss150Z = np.array([ 
    502     -0.9998723404457334, 
    503     -0.9993274305065947, 
    504     -0.9983473449340834, 
    505     -0.9969322929775997, 
    506     -0.9950828645255290, 
    507     -0.9927998590434373, 
    508     -0.9900842691660192, 
    509     -0.9869372772712794, 
    510     -0.9833602541697529, 
    511     -0.9793547582425894, 
    512     -0.9749225346595943, 
    513     -0.9700655145738374, 
    514     -0.9647858142586956, 
    515     -0.9590857341746905, 
    516     -0.9529677579610971, 
    517     -0.9464345513503147, 
    518     -0.9394889610042837, 
    519     -0.9321340132728527, 
    520     -0.9243729128743136, 
    521     -0.9162090414984952, 
    522     -0.9076459563329236, 
    523     -0.8986873885126239, 
    524     -0.8893372414942055, 
    525     -0.8795995893549102, 
    526     -0.8694786750173527, 
    527     -0.8589789084007133, 
    528     -0.8481048644991847, 
    529     -0.8368612813885015, 
    530     -0.8252530581614230, 
    531     -0.8132852527930605, 
    532     -0.8009630799369827, 
    533     -0.7882919086530552, 
    534     -0.7752772600680049, 
    535     -0.7619248049697269, 
    536     -0.7482403613363824, 
    537     -0.7342298918013638, 
    538     -0.7198995010552305, 
    539     -0.7052554331857488, 
    540     -0.6903040689571928, 
    541     -0.6750519230300931, 
    542     -0.6595056411226444, 
    543     -0.6436719971150083, 
    544     -0.6275578900977726, 
    545     -0.6111703413658551, 
    546     -0.5945164913591590, 
    547     -0.5776035965513142, 
    548     -0.5604390262878617, 
    549     -0.5430302595752546, 
    550     -0.5253848818220803, 
    551     -0.5075105815339176, 
    552     -0.4894151469632753, 
    553     -0.4711064627160663, 
    554     -0.4525925063160997, 
    555     -0.4338813447290861, 
    556     -0.4149811308476706, 
    557     -0.3959000999390257, 
    558     -0.3766465660565522, 
    559     -0.3572289184172501, 
    560     -0.3376556177463400, 
    561     -0.3179351925907259, 
    562     -0.2980762356029071, 
    563     -0.2780873997969574, 
    564     -0.2579773947782034, 
    565     -0.2377549829482451, 
    566     -0.2174289756869712, 
    567     -0.1970082295132342, 
    568     -0.1765016422258567, 
    569     -0.1559181490266516, 
    570     -0.1352667186271445, 
    571     -0.1145563493406956, 
    572     -0.0937960651617229, 
    573     -0.0729949118337358, 
    574     -0.0521619529078925, 
    575     -0.0313062657937972, 
    576     -0.0104369378042598, 
    577     0.0104369378042598, 
    578     0.0313062657937972, 
    579     0.0521619529078925, 
    580     0.0729949118337358, 
    581     0.0937960651617229, 
    582     0.1145563493406956, 
    583     0.1352667186271445, 
    584     0.1559181490266516, 
    585     0.1765016422258567, 
    586     0.1970082295132342, 
    587     0.2174289756869712, 
    588     0.2377549829482451, 
    589     0.2579773947782034, 
    590     0.2780873997969574, 
    591     0.2980762356029071, 
    592     0.3179351925907259, 
    593     0.3376556177463400, 
    594     0.3572289184172501, 
    595     0.3766465660565522, 
    596     0.3959000999390257, 
    597     0.4149811308476706, 
    598     0.4338813447290861, 
    599     0.4525925063160997, 
    600     0.4711064627160663, 
    601     0.4894151469632753, 
    602     0.5075105815339176, 
    603     0.5253848818220803, 
    604     0.5430302595752546, 
    605     0.5604390262878617, 
    606     0.5776035965513142, 
    607     0.5945164913591590, 
    608     0.6111703413658551, 
    609     0.6275578900977726, 
    610     0.6436719971150083, 
    611     0.6595056411226444, 
    612     0.6750519230300931, 
    613     0.6903040689571928, 
    614     0.7052554331857488, 
    615     0.7198995010552305, 
    616     0.7342298918013638, 
    617     0.7482403613363824, 
    618     0.7619248049697269, 
    619     0.7752772600680049, 
    620     0.7882919086530552, 
    621     0.8009630799369827, 
    622     0.8132852527930605, 
    623     0.8252530581614230, 
    624     0.8368612813885015, 
    625     0.8481048644991847, 
    626     0.8589789084007133, 
    627     0.8694786750173527, 
    628     0.8795995893549102, 
    629     0.8893372414942055, 
    630     0.8986873885126239, 
    631     0.9076459563329236, 
    632     0.9162090414984952, 
    633     0.9243729128743136, 
    634     0.9321340132728527, 
    635     0.9394889610042837, 
    636     0.9464345513503147, 
    637     0.9529677579610971, 
    638     0.9590857341746905, 
    639     0.9647858142586956, 
    640     0.9700655145738374, 
    641     0.9749225346595943, 
    642     0.9793547582425894, 
    643     0.9833602541697529, 
    644     0.9869372772712794, 
    645     0.9900842691660192, 
    646     0.9927998590434373, 
    647     0.9950828645255290, 
    648     0.9969322929775997, 
    649     0.9983473449340834, 
    650     0.9993274305065947, 
    651     0.9998723404457334 
    652 ]) 
    653  
    654 Gauss150Wt = np.array([ 
    655     0.0003276086705538, 
    656     0.0007624720924706, 
    657     0.0011976474864367, 
    658     0.0016323569986067, 
    659     0.0020663664924131, 
    660     0.0024994789888943, 
    661     0.0029315036836558, 
    662     0.0033622516236779, 
    663     0.0037915348363451, 
    664     0.0042191661429919, 
    665     0.0046449591497966, 
    666     0.0050687282939456, 
    667     0.0054902889094487, 
    668     0.0059094573005900, 
    669     0.0063260508184704, 
    670     0.0067398879387430, 
    671     0.0071507883396855, 
    672     0.0075585729801782, 
    673     0.0079630641773633, 
    674     0.0083640856838475, 
    675     0.0087614627643580, 
    676     0.0091550222717888, 
    677     0.0095445927225849, 
    678     0.0099300043714212, 
    679     0.0103110892851360, 
    680     0.0106876814158841, 
    681     0.0110596166734735, 
    682     0.0114267329968529, 
    683     0.0117888704247183, 
    684     0.0121458711652067, 
    685     0.0124975796646449, 
    686     0.0128438426753249, 
    687     0.0131845093222756, 
    688     0.0135194311690004, 
    689     0.0138484622795371, 
    690     0.0141714592928592, 
    691     0.0144882814685445, 
    692     0.0147987907597169, 
    693     0.0151028518701744, 
    694     0.0154003323133401, 
    695     0.0156911024699895, 
    696     0.0159750356447283, 
    697     0.0162520081211971, 
    698     0.0165218992159766, 
    699     0.0167845913311726, 
    700     0.0170399700056559, 
    701     0.0172879239649355, 
    702     0.0175283451696437, 
    703     0.0177611288626114, 
    704     0.0179861736145128, 
    705     0.0182033813680609, 
    706     0.0184126574807331, 
    707     0.0186139107660094, 
    708     0.0188070535331042, 
    709     0.0189920016251754, 
    710     0.0191686744559934, 
    711     0.0193369950450545, 
    712     0.0194968900511231, 
    713     0.0196482898041878, 
    714     0.0197911283358190, 
    715     0.0199253434079123, 
    716     0.0200508765398072, 
    717     0.0201676730337687, 
    718     0.0202756819988200, 
    719     0.0203748563729175, 
    720     0.0204651529434560, 
    721     0.0205465323660984, 
    722     0.0206189591819181, 
    723     0.0206824018328499, 
    724     0.0207368326754401, 
    725     0.0207822279928917, 
    726     0.0208185680053983, 
    727     0.0208458368787627, 
    728     0.0208640227312962, 
    729     0.0208731176389954, 
    730     0.0208731176389954, 
    731     0.0208640227312962, 
    732     0.0208458368787627, 
    733     0.0208185680053983, 
    734     0.0207822279928917, 
    735     0.0207368326754401, 
    736     0.0206824018328499, 
    737     0.0206189591819181, 
    738     0.0205465323660984, 
    739     0.0204651529434560, 
    740     0.0203748563729175, 
    741     0.0202756819988200, 
    742     0.0201676730337687, 
    743     0.0200508765398072, 
    744     0.0199253434079123, 
    745     0.0197911283358190, 
    746     0.0196482898041878, 
    747     0.0194968900511231, 
    748     0.0193369950450545, 
    749     0.0191686744559934, 
    750     0.0189920016251754, 
    751     0.0188070535331042, 
    752     0.0186139107660094, 
    753     0.0184126574807331, 
    754     0.0182033813680609, 
    755     0.0179861736145128, 
    756     0.0177611288626114, 
    757     0.0175283451696437, 
    758     0.0172879239649355, 
    759     0.0170399700056559, 
    760     0.0167845913311726, 
    761     0.0165218992159766, 
    762     0.0162520081211971, 
    763     0.0159750356447283, 
    764     0.0156911024699895, 
    765     0.0154003323133401, 
    766     0.0151028518701744, 
    767     0.0147987907597169, 
    768     0.0144882814685445, 
    769     0.0141714592928592, 
    770     0.0138484622795371, 
    771     0.0135194311690004, 
    772     0.0131845093222756, 
    773     0.0128438426753249, 
    774     0.0124975796646449, 
    775     0.0121458711652067, 
    776     0.0117888704247183, 
    777     0.0114267329968529, 
    778     0.0110596166734735, 
    779     0.0106876814158841, 
    780     0.0103110892851360, 
    781     0.0099300043714212, 
    782     0.0095445927225849, 
    783     0.0091550222717888, 
    784     0.0087614627643580, 
    785     0.0083640856838475, 
    786     0.0079630641773633, 
    787     0.0075585729801782, 
    788     0.0071507883396855, 
    789     0.0067398879387430, 
    790     0.0063260508184704, 
    791     0.0059094573005900, 
    792     0.0054902889094487, 
    793     0.0050687282939456, 
    794     0.0046449591497966, 
    795     0.0042191661429919, 
    796     0.0037915348363451, 
    797     0.0033622516236779, 
    798     0.0029315036836558, 
    799     0.0024994789888943, 
    800     0.0020663664924131, 
    801     0.0016323569986067, 
    802     0.0011976474864367, 
    803     0.0007624720924706, 
    804     0.0003276086705538 
    805 ]) 
     298class Gauss: 
     299    def __init__(self, w, z): 
     300        self.n = len(w) 
     301        self.w = w 
     302        self.z = z 
     303 
     304gauss20 = Gauss( 
     305    w=np.array([ 
     306        .0176140071391521, 
     307        .0406014298003869, 
     308        .0626720483341091, 
     309        .0832767415767047, 
     310        .10193011981724, 
     311        .118194531961518, 
     312        .131688638449177, 
     313        .142096109318382, 
     314        .149172986472604, 
     315        .152753387130726, 
     316        .152753387130726, 
     317        .149172986472604, 
     318        .142096109318382, 
     319        .131688638449177, 
     320        .118194531961518, 
     321        .10193011981724, 
     322        .0832767415767047, 
     323        .0626720483341091, 
     324        .0406014298003869, 
     325        .0176140071391521 
     326    ]), 
     327    z=np.array([ 
     328        -.993128599185095, 
     329        -.963971927277914, 
     330        -.912234428251326, 
     331        -.839116971822219, 
     332        -.746331906460151, 
     333        -.636053680726515, 
     334        -.510867001950827, 
     335        -.37370608871542, 
     336        -.227785851141645, 
     337        -.076526521133497, 
     338        .0765265211334973, 
     339        .227785851141645, 
     340        .37370608871542, 
     341        .510867001950827, 
     342        .636053680726515, 
     343        .746331906460151, 
     344        .839116971822219, 
     345        .912234428251326, 
     346        .963971927277914, 
     347        .993128599185095 
     348    ]) 
     349) 
     350 
     351gauss76 = Gauss( 
     352    w=np.array([ 
     353        .00126779163408536,             #0 
     354        .00294910295364247, 
     355        .00462793522803742, 
     356        .00629918049732845, 
     357        .00795984747723973, 
     358        .00960710541471375, 
     359        .0112381685696677, 
     360        .0128502838475101, 
     361        .0144407317482767, 
     362        .0160068299122486, 
     363        .0175459372914742,              #10 
     364        .0190554584671906, 
     365        .020532847967908, 
     366        .0219756145344162, 
     367        .0233813253070112, 
     368        .0247476099206597, 
     369        .026072164497986, 
     370        .0273527555318275, 
     371        .028587223650054, 
     372        .029773487255905, 
     373        .0309095460374916,              #20 
     374        .0319934843404216, 
     375        .0330234743977917, 
     376        .0339977794120564, 
     377        .0349147564835508, 
     378        .0357728593807139, 
     379        .0365706411473296, 
     380        .0373067565423816, 
     381        .0379799643084053, 
     382        .0385891292645067, 
     383        .0391332242205184,              #30 
     384        .0396113317090621, 
     385        .0400226455325968, 
     386        .040366472122844, 
     387        .0406422317102947, 
     388        .0408494593018285, 
     389        .040987805464794, 
     390        .0410570369162294, 
     391        .0410570369162294, 
     392        .040987805464794, 
     393        .0408494593018285,              #40 
     394        .0406422317102947, 
     395        .040366472122844, 
     396        .0400226455325968, 
     397        .0396113317090621, 
     398        .0391332242205184, 
     399        .0385891292645067, 
     400        .0379799643084053, 
     401        .0373067565423816, 
     402        .0365706411473296, 
     403        .0357728593807139,              #50 
     404        .0349147564835508, 
     405        .0339977794120564, 
     406        .0330234743977917, 
     407        .0319934843404216, 
     408        .0309095460374916, 
     409        .029773487255905, 
     410        .028587223650054, 
     411        .0273527555318275, 
     412        .026072164497986, 
     413        .0247476099206597,              #60 
     414        .0233813253070112, 
     415        .0219756145344162, 
     416        .020532847967908, 
     417        .0190554584671906, 
     418        .0175459372914742, 
     419        .0160068299122486, 
     420        .0144407317482767, 
     421        .0128502838475101, 
     422        .0112381685696677, 
     423        .00960710541471375,             #70 
     424        .00795984747723973, 
     425        .00629918049732845, 
     426        .00462793522803742, 
     427        .00294910295364247, 
     428        .00126779163408536              #75 (indexed from 0) 
     429    ]), 
     430    z=np.array([ 
     431        -.999505948362153,              #0 
     432        -.997397786355355, 
     433        -.993608772723527, 
     434        -.988144453359837, 
     435        -.981013938975656, 
     436        -.972229228520377, 
     437        -.961805126758768, 
     438        -.949759207710896, 
     439        -.936111781934811, 
     440        -.92088586125215, 
     441        -.904107119545567,              #10 
     442        -.885803849292083, 
     443        -.866006913771982, 
     444        -.844749694983342, 
     445        -.822068037328975, 
     446        -.7980001871612, 
     447        -.77258672828181, 
     448        -.74587051350361, 
     449        -.717896592387704, 
     450        -.688712135277641, 
     451        -.658366353758143,              #20 
     452        -.626910417672267, 
     453        -.594397368836793, 
     454        -.560882031601237, 
     455        -.526420920401243, 
     456        -.491072144462194, 
     457        -.454895309813726, 
     458        -.417951418780327, 
     459        -.380302767117504, 
     460        -.342012838966962, 
     461        -.303146199807908,              #30 
     462        -.263768387584994, 
     463        -.223945802196474, 
     464        -.183745593528914, 
     465        -.143235548227268, 
     466        -.102483975391227, 
     467        -.0615595913906112, 
     468        -.0205314039939986, 
     469        .0205314039939986, 
     470        .0615595913906112, 
     471        .102483975391227,                       #40 
     472        .143235548227268, 
     473        .183745593528914, 
     474        .223945802196474, 
     475        .263768387584994, 
     476        .303146199807908, 
     477        .342012838966962, 
     478        .380302767117504, 
     479        .417951418780327, 
     480        .454895309813726, 
     481        .491072144462194,               #50 
     482        .526420920401243, 
     483        .560882031601237, 
     484        .594397368836793, 
     485        .626910417672267, 
     486        .658366353758143, 
     487        .688712135277641, 
     488        .717896592387704, 
     489        .74587051350361, 
     490        .77258672828181, 
     491        .7980001871612, #60 
     492        .822068037328975, 
     493        .844749694983342, 
     494        .866006913771982, 
     495        .885803849292083, 
     496        .904107119545567, 
     497        .92088586125215, 
     498        .936111781934811, 
     499        .949759207710896, 
     500        .961805126758768, 
     501        .972229228520377,               #70 
     502        .981013938975656, 
     503        .988144453359837, 
     504        .993608772723527, 
     505        .997397786355355, 
     506        .999505948362153                #75 
     507    ]) 
     508) 
     509 
     510gauss150 = Gauss( 
     511    z=np.array([ 
     512        -0.9998723404457334, 
     513        -0.9993274305065947, 
     514        -0.9983473449340834, 
     515        -0.9969322929775997, 
     516        -0.9950828645255290, 
     517        -0.9927998590434373, 
     518        -0.9900842691660192, 
     519        -0.9869372772712794, 
     520        -0.9833602541697529, 
     521        -0.9793547582425894, 
     522        -0.9749225346595943, 
     523        -0.9700655145738374, 
     524        -0.9647858142586956, 
     525        -0.9590857341746905, 
     526        -0.9529677579610971, 
     527        -0.9464345513503147, 
     528        -0.9394889610042837, 
     529        -0.9321340132728527, 
     530        -0.9243729128743136, 
     531        -0.9162090414984952, 
     532        -0.9076459563329236, 
     533        -0.8986873885126239, 
     534        -0.8893372414942055, 
     535        -0.8795995893549102, 
     536        -0.8694786750173527, 
     537        -0.8589789084007133, 
     538        -0.8481048644991847, 
     539        -0.8368612813885015, 
     540        -0.8252530581614230, 
     541        -0.8132852527930605, 
     542        -0.8009630799369827, 
     543        -0.7882919086530552, 
     544        -0.7752772600680049, 
     545        -0.7619248049697269, 
     546        -0.7482403613363824, 
     547        -0.7342298918013638, 
     548        -0.7198995010552305, 
     549        -0.7052554331857488, 
     550        -0.6903040689571928, 
     551        -0.6750519230300931, 
     552        -0.6595056411226444, 
     553        -0.6436719971150083, 
     554        -0.6275578900977726, 
     555        -0.6111703413658551, 
     556        -0.5945164913591590, 
     557        -0.5776035965513142, 
     558        -0.5604390262878617, 
     559        -0.5430302595752546, 
     560        -0.5253848818220803, 
     561        -0.5075105815339176, 
     562        -0.4894151469632753, 
     563        -0.4711064627160663, 
     564        -0.4525925063160997, 
     565        -0.4338813447290861, 
     566        -0.4149811308476706, 
     567        -0.3959000999390257, 
     568        -0.3766465660565522, 
     569        -0.3572289184172501, 
     570        -0.3376556177463400, 
     571        -0.3179351925907259, 
     572        -0.2980762356029071, 
     573        -0.2780873997969574, 
     574        -0.2579773947782034, 
     575        -0.2377549829482451, 
     576        -0.2174289756869712, 
     577        -0.1970082295132342, 
     578        -0.1765016422258567, 
     579        -0.1559181490266516, 
     580        -0.1352667186271445, 
     581        -0.1145563493406956, 
     582        -0.0937960651617229, 
     583        -0.0729949118337358, 
     584        -0.0521619529078925, 
     585        -0.0313062657937972, 
     586        -0.0104369378042598, 
     587        0.0104369378042598, 
     588        0.0313062657937972, 
     589        0.0521619529078925, 
     590        0.0729949118337358, 
     591        0.0937960651617229, 
     592        0.1145563493406956, 
     593        0.1352667186271445, 
     594        0.1559181490266516, 
     595        0.1765016422258567, 
     596        0.1970082295132342, 
     597        0.2174289756869712, 
     598        0.2377549829482451, 
     599        0.2579773947782034, 
     600        0.2780873997969574, 
     601        0.2980762356029071, 
     602        0.3179351925907259, 
     603        0.3376556177463400, 
     604        0.3572289184172501, 
     605        0.3766465660565522, 
     606        0.3959000999390257, 
     607        0.4149811308476706, 
     608        0.4338813447290861, 
     609        0.4525925063160997, 
     610        0.4711064627160663, 
     611        0.4894151469632753, 
     612        0.5075105815339176, 
     613        0.5253848818220803, 
     614        0.5430302595752546, 
     615        0.5604390262878617, 
     616        0.5776035965513142, 
     617        0.5945164913591590, 
     618        0.6111703413658551, 
     619        0.6275578900977726, 
     620        0.6436719971150083, 
     621        0.6595056411226444, 
     622        0.6750519230300931, 
     623        0.6903040689571928, 
     624        0.7052554331857488, 
     625        0.7198995010552305, 
     626        0.7342298918013638, 
     627        0.7482403613363824, 
     628        0.7619248049697269, 
     629        0.7752772600680049, 
     630        0.7882919086530552, 
     631        0.8009630799369827, 
     632        0.8132852527930605, 
     633        0.8252530581614230, 
     634        0.8368612813885015, 
     635        0.8481048644991847, 
     636        0.8589789084007133, 
     637        0.8694786750173527, 
     638        0.8795995893549102, 
     639        0.8893372414942055, 
     640        0.8986873885126239, 
     641        0.9076459563329236, 
     642        0.9162090414984952, 
     643        0.9243729128743136, 
     644        0.9321340132728527, 
     645        0.9394889610042837, 
     646        0.9464345513503147, 
     647        0.9529677579610971, 
     648        0.9590857341746905, 
     649        0.9647858142586956, 
     650        0.9700655145738374, 
     651        0.9749225346595943, 
     652        0.9793547582425894, 
     653        0.9833602541697529, 
     654        0.9869372772712794, 
     655        0.9900842691660192, 
     656        0.9927998590434373, 
     657        0.9950828645255290, 
     658        0.9969322929775997, 
     659        0.9983473449340834, 
     660        0.9993274305065947, 
     661        0.9998723404457334 
     662    ]), 
     663    w=np.array([ 
     664        0.0003276086705538, 
     665        0.0007624720924706, 
     666        0.0011976474864367, 
     667        0.0016323569986067, 
     668        0.0020663664924131, 
     669        0.0024994789888943, 
     670        0.0029315036836558, 
     671        0.0033622516236779, 
     672        0.0037915348363451, 
     673        0.0042191661429919, 
     674        0.0046449591497966, 
     675        0.0050687282939456, 
     676        0.0054902889094487, 
     677        0.0059094573005900, 
     678        0.0063260508184704, 
     679        0.0067398879387430, 
     680        0.0071507883396855, 
     681        0.0075585729801782, 
     682        0.0079630641773633, 
     683        0.0083640856838475, 
     684        0.0087614627643580, 
     685        0.0091550222717888, 
     686        0.0095445927225849, 
     687        0.0099300043714212, 
     688        0.0103110892851360, 
     689        0.0106876814158841, 
     690        0.0110596166734735, 
     691        0.0114267329968529, 
     692        0.0117888704247183, 
     693        0.0121458711652067, 
     694        0.0124975796646449, 
     695        0.0128438426753249, 
     696        0.0131845093222756, 
     697        0.0135194311690004, 
     698        0.0138484622795371, 
     699        0.0141714592928592, 
     700        0.0144882814685445, 
     701        0.0147987907597169, 
     702        0.0151028518701744, 
     703        0.0154003323133401, 
     704        0.0156911024699895, 
     705        0.0159750356447283, 
     706        0.0162520081211971, 
     707        0.0165218992159766, 
     708        0.0167845913311726, 
     709        0.0170399700056559, 
     710        0.0172879239649355, 
     711        0.0175283451696437, 
     712        0.0177611288626114, 
     713        0.0179861736145128, 
     714        0.0182033813680609, 
     715        0.0184126574807331, 
     716        0.0186139107660094, 
     717        0.0188070535331042, 
     718        0.0189920016251754, 
     719        0.0191686744559934, 
     720        0.0193369950450545, 
     721        0.0194968900511231, 
     722        0.0196482898041878, 
     723        0.0197911283358190, 
     724        0.0199253434079123, 
     725        0.0200508765398072, 
     726        0.0201676730337687, 
     727        0.0202756819988200, 
     728        0.0203748563729175, 
     729        0.0204651529434560, 
     730        0.0205465323660984, 
     731        0.0206189591819181, 
     732        0.0206824018328499, 
     733        0.0207368326754401, 
     734        0.0207822279928917, 
     735        0.0208185680053983, 
     736        0.0208458368787627, 
     737        0.0208640227312962, 
     738        0.0208731176389954, 
     739        0.0208731176389954, 
     740        0.0208640227312962, 
     741        0.0208458368787627, 
     742        0.0208185680053983, 
     743        0.0207822279928917, 
     744        0.0207368326754401, 
     745        0.0206824018328499, 
     746        0.0206189591819181, 
     747        0.0205465323660984, 
     748        0.0204651529434560, 
     749        0.0203748563729175, 
     750        0.0202756819988200, 
     751        0.0201676730337687, 
     752        0.0200508765398072, 
     753        0.0199253434079123, 
     754        0.0197911283358190, 
     755        0.0196482898041878, 
     756        0.0194968900511231, 
     757        0.0193369950450545, 
     758        0.0191686744559934, 
     759        0.0189920016251754, 
     760        0.0188070535331042, 
     761        0.0186139107660094, 
     762        0.0184126574807331, 
     763        0.0182033813680609, 
     764        0.0179861736145128, 
     765        0.0177611288626114, 
     766        0.0175283451696437, 
     767        0.0172879239649355, 
     768        0.0170399700056559, 
     769        0.0167845913311726, 
     770        0.0165218992159766, 
     771        0.0162520081211971, 
     772        0.0159750356447283, 
     773        0.0156911024699895, 
     774        0.0154003323133401, 
     775        0.0151028518701744, 
     776        0.0147987907597169, 
     777        0.0144882814685445, 
     778        0.0141714592928592, 
     779        0.0138484622795371, 
     780        0.0135194311690004, 
     781        0.0131845093222756, 
     782        0.0128438426753249, 
     783        0.0124975796646449, 
     784        0.0121458711652067, 
     785        0.0117888704247183, 
     786        0.0114267329968529, 
     787        0.0110596166734735, 
     788        0.0106876814158841, 
     789        0.0103110892851360, 
     790        0.0099300043714212, 
     791        0.0095445927225849, 
     792        0.0091550222717888, 
     793        0.0087614627643580, 
     794        0.0083640856838475, 
     795        0.0079630641773633, 
     796        0.0075585729801782, 
     797        0.0071507883396855, 
     798        0.0067398879387430, 
     799        0.0063260508184704, 
     800        0.0059094573005900, 
     801        0.0054902889094487, 
     802        0.0050687282939456, 
     803        0.0046449591497966, 
     804        0.0042191661429919, 
     805        0.0037915348363451, 
     806        0.0033622516236779, 
     807        0.0029315036836558, 
     808        0.0024994789888943, 
     809        0.0020663664924131, 
     810        0.0016323569986067, 
     811        0.0011976474864367, 
     812        0.0007624720924706, 
     813        0.0003276086705538 
     814    ]) 
     815) 
Note: See TracChangeset for help on using the changeset viewer.