Changes in / [72e41a0:92d330fd] in sasmodels


Ignore:
Files:
2 deleted
50 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    re9ed2de r8a5f021  
    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

    r2d09df1 rce8c388  
    66    env: 
    77    - PY=2.7 
    8     - DEPLOY=True 
    98  - os: linux 
    109    env: 
    … …  
    5251  on: 
    5352    branch: master 
    54     condition: $DEPLOY = True 
    5553notifications: 
    5654  slack: 

  • appveyor.yml

    r1258e32 rd810d96  
    6767    - cmd: conda install --yes --quiet obvious-ci 
    6868    - cmd: conda install --yes --quiet numpy toolchain scipy cython cffi 
    69     #- cmd: conda install --yes --channel conda-forge pyopencl 
     69    - cmd: conda install --yes --channel conda-forge pyopencl 
    7070    - cmd: pip install bumps unittest-xml-reporting tinycc 
    7171 

  • doc/developer/overview.rst

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

  • doc/guide/orientation/orientation.rst

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

  • doc/guide/plugin.rst

    r7e6bc45e r3048ec6  
    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 Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background* 
    295 **parameters are implicit to all models, so they do not need to be included 
    296 in the parameter table.** 
     294Iqxy() and form_volume(). And** *scale* **and** *background* **parameters are 
     295implicit to all models, so they do not need to be included in the parameter table.** 
    297296 
    298297- **"name"** is the name of the parameter shown on the FitPage. 
    … …  
    363362    scattered intensity. 
    364363 
    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. 
     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. 
    372369 
    373370Some models will have integer parameters, such as number of pearls in the 
    … …  
    422419That is, the individual models do not need to include polydispersity 
    423420calculations, but instead rely on numerical integration to compute the 
    424 appropriately smeared pattern. 
     421appropriately smeared pattern.   Angular dispersion values over polar angle 
     422$\theta$ requires an additional $\cos \theta$ weighting due to decreased 
     423arc length for the equatorial angle $\phi$ with increasing latitude. 
    425424 
    426425Python Models 
    … …  
    469468barbell).  If I(q; pars) is NaN for any $q$, then those parameters will be 
    470469ignored, and not included in the calculation of the weighted polydispersity. 
     470 
     471Similar to *Iq*, you can define *Iqxy(qx, qy, par1, par2, ...)* where the 
     472parameter list includes any orientation parameters.  If *Iqxy* is not defined, 
     473then it will default to *Iqxy = Iq(sqrt(qx**2+qy**2), par1, par2, ...)*. 
    471474 
    472475Models should define *form_volume(par1, par2, ...)* where the parameter 
    … …  
    494497    } 
    495498 
     499*Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*, 
     500and it includes orientation parameters. 
     501 
    496502*form_volume* defines the volume of the shape. As in python models, it 
    497503includes only the volume parameters. 
    498504 
     505*Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and 
     506*form_volume* will default to 1.0. 
     507 
    499508**source=['fn.c', ...]** includes the listed C source files in the 
    500 program before *Iq* and *form_volume* are defined. This allows you to 
    501 extend the library of C functions available to your model. 
    502  
    503 *c_code* includes arbitrary C code into your kernel, which can be 
    504 handy for defining helper functions for *Iq* and *form_volume*. Note that 
    505 you 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 
    507 external C file and add that file to the list of sources. 
     509program before *Iq* and *Iqxy* are defined. This allows you to extend the 
     510library of C functions available to your model. 
    508511 
    509512Models are defined using double precision declarations for the 
    … …  
    530533    #define INVALID(v) (v.bell_radius < v.radius) 
    531534 
    532 The INVALID define can go into *Iq*, or *c_code*, or an external C file 
    533 listed in *source*. 
    534  
    535 Oriented Shapes 
    536 ............... 
    537  
    538 If the scattering is dependent on the orientation of the shape, then you 
    539 will need to include *orientation* parameters *theta*, *phi* and *psi* 
    540 at 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 
    542 be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its 
    543 canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the 
    544 laboratory frame and beam travelling along $-z$. 
    545  
    546 The 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 
    548 symmetric about *c* then *psi* is not needed, and the model is called 
    549 as *Iqac(qab, qc, par1, par2, ...)*.  In either case, the orientation 
    550 parameters are not included in the function call. 
    551  
    552 For 1D oriented shapes, an integral over all angles is usually needed for 
    553 the *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  
    565 for 
    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  
    573 Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the 
    574 numerical 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  
    595 The *z* values for the Gauss-Legendre integration extends from -1 to 1, so 
    596 the double sum of *w[i]w[j]* explains the factor of 4.  Correcting for the 
    597 average *dz[i]dz[j]* gives $(1-0) \cdot (\pi/2-0) = \pi/2$.  The $8/(4 \pi)$ 
    598 factor comes from the integral over the quadrant.  With less symmetry (eg., 
    599 in the bcc and fcc paracrystal models), then an integral over the entire 
    600 sphere may be necessary. 
    601  
    602 For simpler models which are rotationally symmetric a single integral 
    603 suffices: 
    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  
    611 for 
    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  
    619 with 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  
    632 Magnetism 
    633 ......... 
    634  
    635 Magnetism is supported automatically for all shapes by modifying the 
    636 effective SLD of particle according to the Halpern-Johnson vector 
    637 describing the interaction between neutron spin and magnetic field.  All 
    638 parameters marked as type *sld* in the parameter table are treated as 
    639 possibly magnetic particles with magnitude *M0* and direction 
    640 *mtheta* and *mphi*.  Polarization parameters are also provided 
    641 automatically for magnetic models to set the spin state of the measurement. 
    642  
    643 For more complicated systems where magnetism is not uniform throughout 
    644 the individual particles, you will need to write your own models. 
    645 You should not mark the nuclear sld as type *sld*, but instead leave 
    646 them unmarked and provide your own magnetism and polarization parameters. 
    647 For 2D measurements you will need $(q_x, q_y)$ values for the measurement 
    648 to compute the proper magnetism and orientation, which you can implement 
    649 using *Iqxy(qx, qy, par1, par2, ...)*. 
    650  
    651535Special Functions 
    652536................. 
    … …  
    659543    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: 
    660544        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ 
    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. 
     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. 
    665549    sin, cos, tan, asin, acos, atan: 
    666550        Trigonometry functions and inverses, operating on radians. 
    … …  
    673557        atan(y/x) would return a value in quadrant I. Similarly for 
    674558        quadrants II and IV when $x$ and $y$ have opposite sign. 
    675     fabs(x), fmin(x,y), fmax(x,y), trunc, rint: 
     559    fmin(x,y), fmax(x,y), trunc, rint: 
    676560        Floating point functions.  rint(x) returns the nearest integer. 
    677561    NAN: 
    678562        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that 
    679563        you cannot use :code:`x == NAN` to test for NaN values since that 
    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. 
     564        will always return false.  NAN does not equal NAN! 
    682565    INFINITY: 
    683566        $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x) 
    … …  
    851734        Similar arrays are available in :code:`gauss20.c` for 20-point 
    852735        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. 
    856736 
    857737        :code:`source = ["lib/gauss76.c", ...]` 
    … …  
    911791can be a model that runs 1000x faster on a good card.  Even your laptop may 
    912792show a 50x improvement or more over the equivalent pure python model. 
     793 
     794External C Models 
     795................. 
     796 
     797External C models are very much like embedded C models, except that 
     798*Iq*, *Iqxy* and *form_volume* are defined in an external source file 
     799loaded using the *source=[...]* statement. You need to supply the function 
     800declarations for each of these that you need instead of building them 
     801automatically from the parameter table. 
    913802 
    914803 
    … …  
    11131002variable name *Rg* for example because $R_g$ is the right name for the model 
    11141003parameter then ignore the lint errors.  Also, ignore *missing-docstring* 
    1115 for standard model functions *Iq*, *Iqac*, etc. 
     1004for standard model functions *Iq*, *Iqxy*, etc. 
    11161005 
    11171006We will have delinting sessions at the SasView Code Camps, where we can 

  • explore/asymint.py

    ra1c32c2 r1820208  
    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(a*qa/2) 
    89         siB = env.sas_sinx_x(b*qb/2) 
    90         siC = env.sas_sinx_x(c*qc/2) 
     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) 
    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  
    97 def 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) 
    12895    return norm, Fq 
    12996 
    … …  
    217184    NORM, KERNEL = make_parallelepiped(A, B, C) 
    218185    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
    219 elif 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) 
    237186elif shape == 'paracrystal': 
    238187    LATTICE = 'bcc' 
    … …  
    393342    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    394343    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)) 
    397344    #gridded_2d(Q, n=2**8+1) 
    398345    gridded_2d(Q, n=2**10+1) 
    399     #gridded_2d(Q, n=2**12+1) 
     346    #gridded_2d(Q, n=2**13+1) 
    400347    #gridded_2d(Q, n=2**15+1) 
    401     if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
    402         # adaptive forms on models for which the calculations are fast enough 
     348    if shape != 'paracrystal':  # adaptive forms are too slow! 
    403349        print("dblquad", *scipy_dblquad_2d(Q)) 
    404350        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 

  • explore/jitter.py

    r8cfb486 rff10479  
    165165    # constants in kernel_iq.c 
    166166    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance', 
    167     'azimuthal_equal_area', 
    168167] 
    169168def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 

  • explore/precision.py

    r2a7e20e r2a602c7  
    345345) 
    346346add_function( 
    347     name="(1/2-sin(x)/x+(1-cos(x))/x^2)/x", 
     347    name="(1/2+(1-cos(x))/x^2-sin(x)/x)/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" or one of: 
     622and name is "all [first]" or one of: 
    623623    """+names) 
    624624    sys.exit(1) 

  • sasmodels/compare.py

    r2a7e20e r2d81cfe  
    4242from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4343from .direct_model import DirectModel, get_mesh 
    44 from .generate import FLOAT_RE, set_integration_size 
     44from .generate import FLOAT_RE 
    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 
    9594 
    9695    === precision options === 
    … …  
    694693        data = empty_data2D(q, resolution=res) 
    695694        data.accuracy = opts['accuracy'] 
    696         set_beam_stop(data, qmin) 
     695        set_beam_stop(data, 0.0004) 
    697696        index = ~data.mask 
    698697    else: 
    … …  
    707706    return data, index 
    708707 
    709 def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
     708def make_engine(model_info, data, dtype, cutoff): 
    710709    # type: (ModelInfo, Data, str, float) -> Calculator 
    711710    """ 
    … …  
    715714    than OpenCL. 
    716715    """ 
    717     if ngauss: 
    718         set_integration_size(model_info, ngauss) 
    719  
    720716    if dtype is None or not dtype.endswith('!'): 
    721717        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    … …  
    958954    'poly', 'mono', 'cutoff=', 
    959955    'magnetic', 'nonmagnetic', 
    960     'accuracy=', 'ngauss=', 
     956    'accuracy=', 
    961957    'neval=',  # for timing... 
    962958 
    … …  
    10931089        'show_weights' : False, 
    10941090        'sphere'    : 0, 
    1095         'ngauss'    : '0', 
    10961091    } 
    10971092    for arg in flags: 
    … …  
    11201115        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
    11211116        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
    1122         elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:] 
    11231117        elif arg.startswith('-random='): 
    11241118            opts['seed'] = int(arg[8:]) 
    … …  
    11751169 
    11761170    comparison = any(PAR_SPLIT in v for v in values) 
    1177  
    11781171    if PAR_SPLIT in name: 
    11791172        names = name.split(PAR_SPLIT, 2) 
    … …  
    11881181        return None 
    11891182 
    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  
    11961183    if PAR_SPLIT in opts['engine']: 
    11971184        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    … …  
    12121199        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12131200 
    1214     base = make_engine(model_info[0], data, opts['engine'][0], 
    1215                        opts['cutoff'][0], opts['ngauss'][0]) 
     1201    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
    12161202    if comparison: 
    1217         comp = make_engine(model_info[1], data, opts['engine'][1], 
    1218                            opts['cutoff'][1], opts['ngauss'][1]) 
     1203        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
    12191204    else: 
    12201205        comp = None 
    … …  
    12891274        if model_info != model_info2: 
    12901275            pars2 = randomize_pars(model_info2, pars2) 
    1291             limit_dimensions(model_info2, pars2, maxdim) 
     1276            limit_dimensions(model_info, pars2, maxdim) 
    12921277            # Share values for parameters with the same name 
    12931278            for k, v in pars.items(): 

  • sasmodels/details.py

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

  • sasmodels/generate.py

    r108e70e r2d81cfe  
    77    particular dimensions averaged over all orientations. 
    88 
    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. 
     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. 
    2215 
    2316    *form_volume(p1, p2, ...)* returns the volume of the form with particular 
    … …  
    3831scale and background parameters for each model. 
    3932 
    40 C code should be stylized C-99 functions written for OpenCL. All functions 
    41 need prototype declarations even if the are defined before they are used. 
    42 Although OpenCL supports *#include* preprocessor directives, the list of 
    43 includes should be given as part of the metadata in the kernel module 
    44 definition. The included files should be listed using a path relative to the 
    45 kernel module, or if using "lib/file.c" if it is one of the standard includes 
    46 provided with the sasmodels source. The includes need to be listed in order 
    47 so that functions are defined before they are used. 
     33*Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 
     34functions written for OpenCL.  All functions need prototype declarations 
     35even if the are defined before they are used.  OpenCL does not support 
     36*#include* preprocessor directives, so instead the list of includes needs 
     37to be given as part of the metadata in the kernel module definition. 
     38The included files should be listed using a path relative to the kernel 
     39module, or if using "lib/file.c" if it is one of the standard includes 
     40provided with the sasmodels source.  The includes need to be listed in 
     41order so that functions are defined before they are used. 
    4842 
    4943Floating point values should be declared as *double*.  For single precision 
    … …  
    113107    present, the volume ratio is 1. 
    114108 
    115     *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
    116     the C source code for the body of the volume, Iq, and Iqac functions 
     109    *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the 
     110    C source code for the body of the volume, Iq, and Iqxy functions 
    117111    respectively.  These can also be defined in the last source file. 
    118112 
    119     *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the 
     113    *Iq* and *Iqxy* also be instead be python functions defining the 
    120114    kernel.  If they are marked as *Iq.vectorized = True* then the 
    121115    kernel is passed the entire *q* vector at once, otherwise it is 
    … …  
    174168from zlib import crc32 
    175169from inspect import currentframe, getframeinfo 
    176 import logging 
    177170 
    178171import numpy as np  # type: ignore 
    … …  
    188181    pass 
    189182# pylint: enable=unused-import 
    190  
    191 logger = logging.getLogger(__name__) 
    192183 
    193184# jitter projection to use in the kernel code.  See explore/jitter.py 
    … …  
    279270""" 
    280271 
    281  
    282 def 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] 
    299272 
    300273def format_units(units): 
    … …  
    635608 
    636609""" 
    637 def _gen_fn(model_info, name, pars): 
    638     # type: (ModelInfo, str, List[Parameter]) -> str 
     610def _gen_fn(name, pars, body, filename, line): 
     611    # type: (str, List[Parameter], str, str, int) -> str 
    639612    """ 
    640613    Generate a function given pars and body. 
    … …  
    648621    """ 
    649622    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) 
    654623    return _FN_TEMPLATE % { 
    655624        'name': name, 'pars': par_decl, 'body': body, 
    656         'filename': filename.replace('\\', '\\\\'), 'line': lineno, 
     625        'filename': filename.replace('\\', '\\\\'), 'line': line, 
    657626    } 
    658627 
    … …  
    669638 
    670639# type in IQXY pattern could be single, float, double, long double, ... 
    671 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     640_IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)", 
    672641                           flags=re.MULTILINE) 
    673 def find_xy_mode(source): 
     642def _have_Iqxy(sources): 
    674643    # type: (List[str]) -> bool 
    675644    """ 
    676     Return the xy mode as qa, qac, qabc or qxy. 
     645    Return true if any file defines Iqxy. 
    677646 
    678647    Note this is not a C parser, and so can be easily confused by 
    679648    non-standard syntax.  Also, it will incorrectly identify the following 
    680     as having 2D models:: 
     649    as having Iqxy:: 
    681650 
    682651        /* 
    683         double Iqac(qab, qc, ...) { ... fill this in later ... } 
     652        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
    684653        */ 
    685654 
    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  
    701 def _add_source(source, code, path, lineno=1): 
     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 
     664def _add_source(source, code, path): 
    702665    """ 
    703666    Add a file to the list of source code chunks, tagged with path and line. 
    704667    """ 
    705668    path = path.replace('\\', '\\\\') 
    706     source.append('#line %d "%s"' % (lineno, path)) 
     669    source.append('#line 1 "%s"' % path) 
    707670    source.append(code) 
    708671 
    … …  
    735698    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    736699 
     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 
    737705    # Build initial sources 
    738706    source = [] 
    … …  
    741709        _add_source(source, code, path) 
    742710 
    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  
    747711    # Make parameters for q, qx, qy so that we can use them in declarations 
    748     q, qx, qy, qab, qa, qb, qc \ 
    749         = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     712    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
    750713    # Generate form_volume function, etc. from body only 
    751714    if isinstance(model_info.form_volume, str): 
    752715        pars = partable.form_volume_parameters 
    753         source.append(_gen_fn(model_info, 'form_volume', pars)) 
     716        source.append(_gen_fn('form_volume', pars, model_info.form_volume, 
     717                              model_info.filename, model_info._form_volume_line)) 
    754718    if isinstance(model_info.Iq, str): 
    755719        pars = [q] + partable.iq_parameters 
    756         source.append(_gen_fn(model_info, 'Iq', pars)) 
     720        source.append(_gen_fn('Iq', pars, model_info.Iq, 
     721                              model_info.filename, model_info._Iq_line)) 
    757722    if isinstance(model_info.Iqxy, str): 
    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") 
     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)) 
    780726 
    781727    # Define the parameter table 
    … …  
    803749    if xy_mode == 'qabc': 
    804750        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    805         call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars 
     751        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
    806752        clear_iqxy = "#undef CALL_IQ_ABC" 
    807753    elif xy_mode == 'qac': 
    808754        pars = ",".join(["_qa", "_qc"] + model_refs) 
    809         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
     755        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
    810756        clear_iqxy = "#undef CALL_IQ_AC" 
    811     elif xy_mode == 'qa': 
     757    else:  # xy_mode == 'qa' 
    812758        pars = ",".join(["_qa"] + model_refs) 
    813759        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    814760        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  
    827761 
    828762    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 

  • sasmodels/kernel_header.c

    r108e70e r8698a0d  
    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

    rec8d4ac r6aee3ab  
    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 
    3433//  INVALID(table) : test if the current point is feesible to calculate.  This 
    3534//      will be defined in the kernel definition file. 
    … …  
    174173static double 
    175174qac_apply( 
    176     QACRotation *rotation, 
     175    QACRotation rotation, 
    177176    double qx, double qy, 
    178177    double *qa_out, double *qc_out) 
    179178{ 
    180     const double dqc = rotation->R31*qx + rotation->R32*qy; 
     179    const double dqc = rotation.R31*qx + rotation.R32*qy; 
    181180    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    182181    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
    … …  
    247246static double 
    248247qabc_apply( 
    249     QABCRotation *rotation, 
     248    QABCRotation rotation, 
    250249    double qx, double qy, 
    251250    double *qa_out, double *qb_out, double *qc_out) 
    252251{ 
    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; 
     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; 
    256255} 
    257256 
    … …  
    454453  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code. 
    455454  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi); 
    456   #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc) 
     455  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 
    457456  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) 
    458457 
    … …  
    468467  local_values.table.psi = 0.; 
    469468  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi) 
    470   #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
     469  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
    471470  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
    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 
     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) 
    486478  #define APPLY_PROJECTION() const double weight=weight0 
    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 
     479#else // !spherosymmetric projection 
    518480  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 
    519481  const double theta = values[details->theta_par+2]; 
    … …  
    526488  // we go through the mesh. 
    527489  double dtheta, dphi, weight; 
    528   #if PROJECTION == 1 // equirectangular 
     490  #if PROJECTION == 1 
    529491    #define APPLY_PROJECTION() do { \ 
    530492      dtheta = local_values.table.theta; \ 
    … …  
    532494      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \ 
    533495    } while (0) 
    534   #elif PROJECTION == 2 // sinusoidal 
     496  #elif PROJECTION == 2 
    535497    #define APPLY_PROJECTION() do { \ 
    536498      dtheta = local_values.table.theta; \ 
    … …  
    542504    } while (0) 
    543505  #endif 
    544 #endif // done defining APPLY_PROJECTION 
     506#endif // !spherosymmetric projection 
    545507 
    546508// ** define looping macros ** 

  • sasmodels/kernelpy.py

    r108e70e r2d81cfe  
    2626# pylint: enable=unused-import 
    2727 
    28 logger = logging.getLogger(__name__) 
    29  
    3028class PyModel(KernelModel): 
    3129    """ 
    … …  
    3331    """ 
    3432    def __init__(self, model_info): 
    35         # Make sure Iq is available and vectorized 
     33        # Make sure Iq and Iqxy are available and vectorized 
    3634        _create_default_functions(model_info) 
    3735        self.info = model_info 
    3836        self.dtype = np.dtype('d') 
    39         logger.info("load python model " + self.info.name) 
     37        logging.info("load python model " + self.info.name) 
    4038 
    4139    def make_kernel(self, q_vectors): 
    4240        q_input = PyInput(q_vectors, dtype=F64) 
    43         return PyKernel(self.info, q_input) 
     41        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
     42        return PyKernel(kernel, self.info, q_input) 
    4443 
    4544    def release(self): 
    … …  
    9089    Callable SAS kernel. 
    9190 
    92     *kernel* is the kernel object to call. 
     91    *kernel* is the DllKernel object to call. 
    9392 
    9493    *model_info* is the module information 
    … …  
    105104    Call :meth:`release` when done with the kernel instance. 
    106105    """ 
    107     def __init__(self, model_info, q_input): 
     106    def __init__(self, kernel, model_info, q_input): 
    108107        # type: (callable, ModelInfo, List[np.ndarray]) -> None 
    109108        self.dtype = np.dtype('d') 
    … …  
    111110        self.q_input = q_input 
    112111        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 else 
    162                         (lambda: 1.0)) 
     161        self._volume = ((lambda: form_volume(*volume_args)) if form_volume 
     162                        else (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 
    264263    _create_vector_Iq(model_info) 
    265     _create_vector_Iqxy(model_info) 
     264    _create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
    266265 
    267266 
    … …  
    281280        model_info.Iq = vector_Iq 
    282281 
    283  
    284282def _create_vector_Iqxy(model_info): 
    285283    """ 
    286284    Define Iqxy as a vector function if it exists, or default it from Iq(). 
    287285    """ 
    288     Iqxy = getattr(model_info, 'Iqxy', None) 
     286    Iq, Iqxy = model_info.Iq, model_info.Iqxy 
    289287    if callable(Iqxy): 
    290288        if not getattr(Iqxy, 'vectorized', False): 
    … …  
    297295            vector_Iqxy.vectorized = True 
    298296            model_info.Iqxy = vector_Iqxy 
    299     else: 
     297    elif callable(Iq): 
    300298        #print("defaulting Iqxy") 
    301299        # Iq is vectorized because create_vector_Iq was already called. 
    302         Iq = model_info.Iq 
    303300        def default_Iqxy(qx, qy, *args): 
    304301            """ 

  • sasmodels/modelinfo.py

    r108e70e r2d81cfe  
    3737 
    3838# assumptions about common parameters exist throughout the code, such as: 
    39 # (1) kernel functions Iq, Iqxy, Iqac, Iqabc, form_volume, ... don't see them 
     39# (1) kernel functions Iq, Iqxy, 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 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 
     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 
    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 
    388391    * *form_volume_parameters* is the list of parameters to the form_volume(...) 
    389392      function, with vector parameter p sent as p[]. 
    … …  
    440443        self.iq_parameters = [p for p in self.kernel_parameters 
    441444                              if p.type not in ('orientation', 'magnetic')] 
    442         self.orientation_parameters = [p for p in self.kernel_parameters 
    443                                        if p.type == 'orientation'] 
     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'] 
    444448        self.form_volume_parameters = [p for p in self.kernel_parameters 
    445449                                       if p.type == 'volume'] 
    … …  
    486490                if p.type != 'orientation': 
    487491                    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") 
    490492        if theta >= 0 and phi >= 0: 
    491             last_par = len(self.kernel_parameters) - 1 
    492493            if phi != theta+1: 
    493494                raise TypeError("phi must follow theta") 
    494495            if psi >= 0 and psi != phi+1: 
    495496                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") 
    499497        elif theta >= 0 or phi >= 0 or psi >= 0: 
    500498            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
    … …  
    717715 
    718716 
    719 #: Set of variables defined in the model that might contain C code 
    720 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
    721  
    722717def _find_source_lines(model_info, kernel_module): 
    723718    # type: (ModelInfo, ModuleType) -> None 
    … …  
    725720    Identify the location of the C source inside the model definition file. 
    726721 
    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)): 
     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): 
    740738        return 
    741739 
    742     # load the module source if we can 
     740    # find the defintion lines for the different code blocks 
    743741    try: 
    744742        source = inspect.getsource(kernel_module) 
    745743    except IOError: 
    746744        return 
    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 
     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 
    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 Iqac/Iqabc functions will be created for python 
     763    Note: vectorized Iq and Iqxy 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) 
    793792    # TODO: check the structure of the tests 
    794793    info.tests = getattr(kernel_module, 'tests', []) 
    … …  
    798797    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    799798    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 
    802799    info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore 
    803800    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    … …  
    814811    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
    815812 
    816     if callable(info.Iq) and parameters.has_2d: 
    817         raise ValueError("oriented python models not supported") 
    818  
    819     info.lineno = {} 
    820813    _find_source_lines(info, kernel_module) 
    821814 
    … …  
    831824 
    832825    The structure should be mostly static, other than the delayed definition 
    833     of *Iq*, *Iqac* and *Iqabc* if they need to be defined. 
     826    of *Iq* and *Iqxy* if they need to be defined. 
    834827    """ 
    835828    #: Full path to the file defining the kernel, if any. 
    … …  
    913906    structure_factor = None # type: bool 
    914907    #: List of C source files used to define the model.  The source files 
    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. 
     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. 
    920913    source = None           # type: List[str] 
    921914    #: The set of tests that must pass.  The format of the tests is described 
    … …  
    942935    #: See :attr:`ER` for details on the parameters. 
    943936    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 
    948937    #: Returns the form volume for python-based models.  Form volume is needed 
    949938    #: for volume normalization in the polydispersity integral.  If no 
    … …  
    966955    #: include the decimal point. See :mod:`generate` for more details. 
    967956    Iq = 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]] 
     957    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
     958    Iqxy = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    972959    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    973960    Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    … …  
    985972    #: Returns a random parameter set for the model 
    986973    random = None           # type: Optional[Callable[[], Dict[str, float]]] 
    987     #: Line numbers for symbols defining C code 
    988     lineno = None           # type: Dict[str, int] 
     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 
    989982 
    990983    def __init__(self): 

  • sasmodels/models/_spherepy.py

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

  • sasmodels/models/barbell.c

    r108e70e rbecded3  
    2323    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    2424    double total = 0.0; 
    25     for (int i = 0; i < GAUSS_N; i++){ 
    26         const double t = GAUSS_Z[i]*zm + zb; 
     25    for (int i = 0; i < 76; i++){ 
     26        const double t = Gauss76Z[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 += GAUSS_W[i] * Fq; 
     30        total += Gauss76Wt[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 < GAUSS_N; i++){ 
    76         const double alpha= GAUSS_Z[i]*zm + zb; 
     75    for (int i = 0; i < 76; i++){ 
     76        const double alpha= Gauss76Z[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 += GAUSS_W[i] * Aq * Aq * sin_alpha; 
     80        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8181    } 
    8282    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9090 
    9191static double 
    92 Iqac(double qab, double qc, 
     92Iqxy(double qab, double qc, 
    9393    double sld, double solvent_sld, 
    9494    double radius_bell, double radius, double length) 

  • sasmodels/models/bcc_paracrystal.c

    r108e70e rea60e08  
    8181 
    8282    double outer_sum = 0.0; 
    83     for(int i=0; i<GAUSS_N; i++) { 
     83    for(int i=0; i<150; i++) { 
    8484        double inner_sum = 0.0; 
    85         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     85        const double theta = Gauss150Z[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<GAUSS_N;j++) { 
    91             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     90        for(int j=0;j<150;j++) { 
     91            const double phi = Gauss150Z[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 += GAUSS_W[j] * form; 
     97            inner_sum += Gauss150Wt[j] * form; 
    9898        } 
    9999        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    100         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     100        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    101101    } 
    102102    outer_sum *= theta_m; 
    … …  
    107107 
    108108 
    109 static double Iqabc(double qa, double qb, double qc, 
     109static double Iqxy(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

    r108e70e rbecded3  
    3030    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    32     for (int i=0; i<GAUSS_N; i++) { 
    33         const double t = GAUSS_Z[i]*zm + zb; 
     32    for (int i=0; i<76 ;i++) { 
     33        const double t = Gauss76Z[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 += GAUSS_W[i] * Fq; 
     37        total += Gauss76Wt[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<GAUSS_N ;i++) { 
    98         const double theta = GAUSS_Z[i]*zm + zb; 
     97    for (int i=0; i<76 ;i++) { 
     98        const double theta = Gauss76Z[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 += GAUSS_W[i] * Aq * Aq * sin_theta; 
     105        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    106106    } 
    107107    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    115115 
    116116static double 
    117 Iqac(double qab, double qc, 
     117Iqxy(double qab, double qc, 
    118118    double sld, double solvent_sld, double radius, 
    119119    double radius_cap, double length) 

  • sasmodels/models/core_shell_bicelle.c

    r108e70e rbecded3  
    5252 
    5353    double total = 0.0; 
    54     for(int i=0;i<GAUSS_N;i++) { 
    55         double theta = (GAUSS_Z[i] + 1.0)*uplim; 
     54    for(int i=0;i<N_POINTS_76;i++) { 
     55        double theta = (Gauss76Z[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 += GAUSS_W[i]*fq*fq*sin_theta; 
     60        total += Gauss76Wt[i]*fq*fq*sin_theta; 
    6161    } 
    6262 
    … …  
    6767 
    6868static double 
    69 Iqac(double qab, double qc, 
     69Iqxy(double qab, double qc, 
    7070    double radius, 
    7171    double thick_rim, 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r108e70e r82592da  
    3737    //initialize integral 
    3838    double outer_sum = 0.0; 
    39     for(int i=0;i<GAUSS_N;i++) { 
     39    for(int i=0;i<76;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 = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    44         const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.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; 
    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<GAUSS_N;j++) { 
     51        for(int j=0;j<76;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 = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 
     56            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( Gauss76Z[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 += GAUSS_W[j] * fq * fq; 
     63            inner_sum += Gauss76Wt[j] * fq * fq; 
    6464        } 
    6565        //now calculate outer integral 
    66         outer_sum += GAUSS_W[i] * inner_sum; 
     66        outer_sum += Gauss76Wt[i] * inner_sum; 
    6767    } 
    6868 
    … …  
    7171 
    7272static double 
    73 Iqabc(double qa, double qb, double qc, 
     73Iqxy(double qa, double qb, double qc, 
    7474    double r_minor, 
    7575    double x_core, 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r108e70e r82592da  
    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<GAUSS_N;i++) { 
     49    for(int i=0;i<76;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 = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 
     52        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     53        const double cos_alpha = ( Gauss76Z[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<GAUSS_N;j++) { 
     60        for(int j=0;j<76;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 = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    63             const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
     62            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     63            const double beta = ( Gauss76Z[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 += GAUSS_W[j] *square(dr1*si1*be1 + 
     69            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
    7070                                              dr2*si1*be2 + 
    7171                                              dr3*si2*be1); 
    7272        } 
    7373        //now calculate outer integral 
    74         outer_sum += GAUSS_W[i] * inner_sum; 
     74        outer_sum += Gauss76Wt[i] * inner_sum; 
    7575    } 
    7676 
    … …  
    7979 
    8080static double 
    81 Iqabc(double qa, double qb, double qc, 
     81Iqxy(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

    r108e70e r110f69c  
    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"], 
    152151    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    153152    ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
    154153    ["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

    r108e70e rbecded3  
    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<GAUSS_N ;i++) { 
     32    for (int i=0; i<76 ;i++) { 
    3333        // translate a point in [-1,1] to a point in [0, pi/2] 
    34         //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 
     34        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    3535        double sin_theta, cos_theta; 
    36         const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4; 
     36        const double theta = Gauss76Z[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 += GAUSS_W[i] * fq * fq * sin_theta; 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4343    } 
    4444    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4848 
    4949 
    50 double Iqac(double qab, double qc, 
     50double Iqxy(double qab, double qc, 
    5151    double core_sld, 
    5252    double shell_sld, 

  • sasmodels/models/core_shell_ellipsoid.c

    r108e70e rbecded3  
    5959    const double b = 0.5; 
    6060    double total = 0.0;     //initialize intergral 
    61     for(int i=0;i<GAUSS_N;i++) { 
    62         const double cos_theta = GAUSS_Z[i]*m + b; 
     61    for(int i=0;i<76;i++) { 
     62        const double cos_theta = Gauss76Z[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 += GAUSS_W[i] * fq * fq; 
     68        total += Gauss76Wt[i] * fq * fq; 
    6969    } 
    7070    total *= m; 
    … …  
    7575 
    7676static double 
    77 Iqac(double qab, double qc, 
     77Iqxy(double qab, double qc, 
    7878    double radius_equat_core, 
    7979    double x_core, 

  • sasmodels/models/core_shell_parallelepiped.c

    re077231 rc69d6d6  
    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 
    71static double 
    82form_volume(double length_a, double length_b, double length_c, 
    93    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    104{ 
    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 
     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; 
    2710} 
    2811 
    … …  
    4124    double thick_rim_c) 
    4225{ 
    43     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
     26    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    4427    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    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) 
     28    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    4729 
    48     const double half_q = 0.5*q; 
     30    const double mu = 0.5 * q * length_b; 
    4931 
    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; 
     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; 
    5338 
    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); 
     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; 
    5969 
    6070    // outer integral (with gauss points), integration limits = 0, 1 
    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); 
     71    double outer_total = 0; //initialize integral 
    6572 
    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); 
     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); 
    7876 
    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 
     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); 
    9088 
    91             inner_sum += GAUSS_W[j] * f * f; 
     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; 
    9296        } 
    93         inner_sum *= 0.5; 
     97        inner_total *= 0.5; 
     98        inner_total_crim *= 0.5; 
    9499        // now sum up the outer integral 
    95         outer_sum += GAUSS_W[i] * inner_sum; 
     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); 
    96103    } 
    97     outer_sum *= 0.5; 
     104    outer_total *= 0.5; 
    98105 
    99106    //convert from [1e-12 A-1] to [cm-1] 
    100     return 1.0e-4 * outer_sum; 
     107    return 1.0e-4 * outer_total; 
    101108} 
    102109 
    103110static double 
    104 Iqabc(double qa, double qb, double qc, 
     111Iqxy(double qa, double qb, double qc, 
    105112    double core_sld, 
    106113    double arim_sld, 
    … …  
    121128    const double drC = crim_sld-solvent_sld; 
    122129 
    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); 
     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); 
    132139 
    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 
     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); 
    144160 
    145161    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r97be877 r2d81cfe  
    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. 
     7"rim" can be different on each (pair) of faces. However at this time the 1D 
     8calculation does **NOT** actually calculate a c face rim despite the presence 
     9of 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. 
    814 
    915The form factor is normalized by the particle volume $V$ such that 
    … …  
    1521where $\langle \ldots \rangle$ is an average over all possible orientations 
    1622of the rectangular solid. 
     23 
    1724 
    1825The function calculated is the form factor of the rectangular solid below. 
    … …  
    3441    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    3542 
    36 **meaning that there are "gaps" at the corners of the solid.** 
     43**meaning that there are "gaps" at the corners of the solid.**  Again note that 
     44$t_C = 0$ currently. 
    3745 
    3846The intensity calculated follows the :ref:`parallelepiped` model, with the 
    3947core-shell intensity being calculated as the square of the sum of the 
    40 amplitudes of the core and the slabs on the edges. 
    41  
    42 the scattering amplitude is computed for a particular orientation of the 
    43 core-shell parallelepiped with respect to the scattering vector and then 
    44 averaged over all possible orientations, where $\alpha$ is the angle between 
    45 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
    46 the angle between projection of the particle in the $xy$ detector plane 
    47 and the $y$ axis. 
     48amplitudes of the core and shell, in the same manner as a core-shell model. 
    4849 
    4950.. math:: 
    5051 
    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  
    61 with 
    62  
    63 .. math:: 
    64  
    65     S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
    66  
    67 and 
    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  
    76 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
    77 are the scattering length of the parallelepiped core, and the rectangular 
    78 slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    79 is the scattering length of the solvent. 
     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. 
    8069 
    8170FITTING NOTES 
    82 ~~~~~~~~~~~~~ 
    83  
    8471If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    85 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
    86 **no interparticle interference effects are included in this calculation.** 
     72value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
     73However, **no interparticle interference effects are included in this 
     74calculation.** 
    8775 
    8876There are many parameters in this model. Hold as many fixed as possible with 
    8977known values, or you will certainly end up at a solution that is unphysical. 
    9078 
     79Constraints must be applied during fitting to ensure that the inequality 
     80$A < B < C$ is not violated. The calculation will not report an error, 
     81but the results will not be correct. 
     82 
    9183The returned value is in units of |cm^-1|, on absolute scale. 
    9284 
    9385NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    9486based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    95 and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
    96 to give an oblate or prolate particle, to give an effective radius, 
    97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     87and length $(C+2t_C)$ values, after appropriately 
     88sorting the three dimensions to give an oblate or prolate particle, to give an 
     89effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    9890 
    9991For 2d data the orientation of the particle is required, described using 
    100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
    101 details of the calculation and angular dispersions see :ref:`orientation`. 
     92angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     93of the calculation and angular dispersions see :ref:`orientation` . 
    10294The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    10395$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    104  
    105 For 2d, constraints must be applied during fitting to ensure that the 
    106 inequality $A < B < C$ is not violated, and hence the correct definition 
    107 of angles is preserved. The calculation will not report an error, 
    108 but the results may be not correct. 
    10996 
    11097.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    127114    Equations (1), (13-14). (in German) 
    128115.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    129    lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
     116   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
    130117   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    131118   =26379>`_ 
    … …  
    136123* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    137124* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
    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. 
     125* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017 
     126* **Currently Under review by:** Paul Butler 
    142127""" 
    143128 
    … …  
    190175        Return equivalent radius (ER) 
    191176    """ 
    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) 
     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.) 
    198185 
    199186# VR defaults to 1.0 
    … …  
    229216            psi_pd=10, psi_pd_n=1) 
    230217 
    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 
     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 
    233219if 0:  # pak: model rewrite; need to update tests 
    234220    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 

  • sasmodels/models/cylinder.c

    r108e70e rbecded3  
    2121 
    2222    double total = 0.0; 
    23     for (int i=0; i<GAUSS_N ;i++) { 
    24         const double theta = GAUSS_Z[i]*zm + zb; 
     23    for (int i=0; i<76 ;i++) { 
     24        const double theta = Gauss76Z[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 += GAUSS_W[i] * form * form * sin_theta; 
     29        total += Gauss76Wt[i] * form * form * sin_theta; 
    3030    } 
    3131    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4545 
    4646static double 
    47 Iqac(double qab, double qc, 
     47Iqxy(double qab, double qc, 
    4848    double sld, 
    4949    double solvent_sld, 

  • sasmodels/models/ellipsoid.c

    r108e70e rbecded3  
    2222 
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    24     // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     24    // const double u = Gauss76Z[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<GAUSS_N;i++) { 
    29         const double u = GAUSS_Z[i]*zm + zb; 
     28    for (int i=0;i<76;i++) { 
     29        const double u = Gauss76Z[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 += GAUSS_W[i] * f * f; 
     32        total += Gauss76Wt[i] * f * f; 
    3333    } 
    3434    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3939 
    4040static double 
    41 Iqac(double qab, double qc, 
     41Iqxy(double qab, double qc, 
    4242    double sld, 
    4343    double sld_solvent, 

  • sasmodels/models/elliptical_cylinder.c

    r108e70e r82592da  
    2222    //initialize integral 
    2323    double outer_sum = 0.0; 
    24     for(int i=0;i<GAUSS_N;i++) { 
     24    for(int i=0;i<76;i++) { 
    2525        //setup inner integral over the ellipsoidal cross-section 
    26         const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     26        const double cos_val = ( Gauss76Z[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<GAUSS_N;j++) { 
    31             const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.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; 
    3233            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3334            const double be = sas_2J1x_x(q*r); 
    34             inner_sum += GAUSS_W[j] * be * be; 
     35            inner_sum += Gauss76Wt[j] * be * be; 
    3536        } 
    3637        //now calculate the value of the inner integral 
    … …  
    3940        //now calculate outer integral 
    4041        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    41         outer_sum += GAUSS_W[i] * inner_sum * si * si; 
     42        outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
    4243    } 
    4344    outer_sum *= 0.5*(vb-va); 
    … …  
    5455 
    5556static double 
    56 Iqabc(double qa, double qb, double qc, 
     57Iqxy(double qa, double qb, double qc, 
    5758     double radius_minor, double r_ratio, double length, 
    5859     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", "elliptical_cylinder.c"] 
     123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", 
     124          "elliptical_cylinder.c"] 
    124125 
    125126demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 

  • sasmodels/models/fcc_paracrystal.c

    r108e70e rf728001  
    5353 
    5454    double outer_sum = 0.0; 
    55     for(int i=0; i<GAUSS_N; i++) { 
     55    for(int i=0; i<150; i++) { 
    5656        double inner_sum = 0.0; 
    57         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     57        const double theta = Gauss150Z[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<GAUSS_N;j++) { 
    63             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     62        for(int j=0;j<150;j++) { 
     63            const double phi = Gauss150Z[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 += GAUSS_W[j] * form; 
     69            inner_sum += Gauss150Wt[j] * form; 
    7070        } 
    7171        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    72         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     72        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    7373    } 
    7474    outer_sum *= theta_m; 
    … …  
    8080 
    8181 
    82 static double Iqabc(double qa, double qb, double qc, 
     82static double Iqxy(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

    r74768cb r592343f  
    1717    double sum=0.0; 
    1818 
    19     for(int i=0;i<GAUSS_N;i++) { 
    20         const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4; 
     19    for(int i=0;i<N_POINTS_76;i++) { 
     20        const double zi = ( Gauss76Z[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 += GAUSS_W[i] * yyy * yyy; 
     25        sum += Gauss76Wt[i] * yyy * yyy; 
    2626    } 
    2727    sum *= 0.5; 

  • sasmodels/models/hollow_cylinder.c

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

  • sasmodels/models/hollow_rectangular_prism.c

    r108e70e r1e7b0db0  
    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     
     40    double outer_sum = 0.0; 
     41    for(int i=0; i<76; i++) { 
    3942 
    40     double outer_sum = 0.0; 
    41     for(int i=0; i<GAUSS_N; i++) { 
    42  
    43         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     43        const double theta = 0.5 * ( Gauss76Z[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<GAUSS_N; j++) { 
     51        for(int j=0; j<76; j++) { 
    5252 
    53             const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     53            const double phi = 0.5 * ( Gauss76Z[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 += GAUSS_W[j] * square(AP1-AP2); 
     68            inner_sum += Gauss76Wt[j] * square(AP1-AP2); 
    6969        } 
    7070        inner_sum *= 0.5 * (v2b-v2a); 
    7171 
    72         outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
     72        outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 
    7373    } 
    7474    outer_sum *= 0.5*(v1b-v1a); 
    … …  
    8484    return 1.0e-4 * delrho * delrho * form; 
    8585} 
    86  
    87 double 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  
     7It computes only the 1D scattering, not the 2D. 
    88 
    99Definition 
    … …  
    6666(which is unitless). 
    6767 
    68 For 2d data the orientation of the particle is required, described using 
    69 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    70 of the calculation and angular dispersions see :ref:`orientation` . 
    71 The 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  
    74 For 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, 
    76 but 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. 
     68**The 2D scattering intensity is not computed by this model.** 
    8969 
    9070 
    … …  
    133113              ["thickness", "Ang", 1, [0, inf], "volume", 
    134114               "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"], 
    141115             ] 
    142116 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r74768cb rab2aea8  
    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<GAUSS_N; i++) { 
    34         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     33    for(int i=0; i<76; i++) { 
     34        const double theta = 0.5 * ( Gauss76Z[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<GAUSS_N; j++) { 
    47             const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     46        for(int j=0; j<76; j++) { 
     47            const double phi = 0.5 * ( Gauss76Z[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 += GAUSS_W[j] * square(AL+AT); 
     64            inner_sum += Gauss76Wt[j] * square(AL+AT); 
    6565        } 
    6666 
    6767        inner_sum *= 0.5 * (v2b-v2a); 
    68         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     68        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
    6969    } 
    7070 

  • sasmodels/models/lib/gauss150.c

    r99b84ec r994d77f  
    1 // Created by Andrew Jackson on 4/23/07 
     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 */ 
    29 
    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. 
    15 constant double Gauss150Z[152]={ 
     10// Gaussians 
     11constant double Gauss150Z[150]={ 
    1612        -0.9998723404457334, 
    1713        -0.9993274305065947, 
    … …  
    163159        0.9983473449340834, 
    164160        0.9993274305065947, 
    165         0.9998723404457334, 
    166         0., // zero padding is ignored 
    167         0.  // zero padding is ignored 
     161        0.9998723404457334 
    168162}; 
    169163 
    170 constant double Gauss150Wt[152]={ 
     164constant double Gauss150Wt[150]={ 
    171165        0.0003276086705538, 
    172166        0.0007624720924706, 
    … …  
    318312        0.0011976474864367, 
    319313        0.0007624720924706, 
    320         0.0003276086705538, 
    321         0., // zero padding is ignored 
    322         0.  // zero padding is ignored 
     314        0.0003276086705538 
    323315}; 

  • sasmodels/models/lib/gauss20.c

    r99b84ec r994d77f  
    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 
     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 */ 
    119 
    1210// Gaussians 

  • sasmodels/models/lib/gauss76.c

    r99b84ec r66d119f  
    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 
     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 
    1110 
    1211// Gaussians 
    13 constant double Gauss76Wt[76]={ 
     12constant double Gauss76Wt[N_POINTS_76]={ 
    1413        .00126779163408536,             //0 
    1514        .00294910295364247, 
    … …  
    9089}; 
    9190 
    92 constant double Gauss76Z[76]={ 
     91constant double Gauss76Z[N_POINTS_76]={ 
    9392        -.999505948362153,              //0 
    9493        -.997397786355355, 

  • sasmodels/models/line.py

    r108e70e r2d81cfe  
    5757Iq.vectorized = True # Iq accepts an array of q values 
    5858 
    59  
    6059def Iqxy(qx, qy, *args): 
    6160    """ 
    … …  
    7069 
    7170Iqxy.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 #""" 
    8071 
    8172def random(): 

  • sasmodels/models/parallelepiped.c

    r108e70e r9b7b23f  
    2323    double outer_total = 0; //initialize integral 
    2424 
    25     for( int i=0; i<GAUSS_N; i++) { 
    26         const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     25    for( int i=0; i<76; i++) { 
     26        const double sigma = 0.5 * ( Gauss76Z[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<GAUSS_N; j++) { 
    33             const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     32        for(int j=0; j<76; j++) { 
     33            const double uu = 0.5 * ( Gauss76Z[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 += GAUSS_W[j] * square(si1 * si2); 
     38            inner_total += Gauss76Wt[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 += GAUSS_W[i] * inner_total * si * si; 
     43        outer_total += Gauss76Wt[i] * inner_total * si * si; 
    4444    } 
    4545    outer_total *= 0.5; 
    … …  
    5353 
    5454static double 
    55 Iqabc(double qa, double qb, double qc, 
     55Iqxy(double qa, double qb, double qc, 
    5656    double sld, 
    5757    double solvent_sld, 

  • sasmodels/models/pringle.c

    r74768cb r1e7b0db0  
    2929    double sumC = 0.0;          // initialize integral 
    3030    double r; 
    31     for (int i=0; i < GAUSS_N; i++) { 
    32         r = GAUSS_Z[i]*zm + zb; 
     31    for (int i=0; i < 76; i++) { 
     32        r = Gauss76Z[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 = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
     37        double y = Gauss76Wt[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 < GAUSS_N; i++) { 
    89         double psi = GAUSS_Z[i]*zm + zb; 
     88    for (int i = 0; i < 76; i++) { 
     89        double psi = Gauss76Z[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 += GAUSS_W[i] * pringle_kernel; 
     95        sum += Gauss76Wt[i] * pringle_kernel; 
    9696    } 
    9797 

  • sasmodels/models/rectangular_prism.c

    r108e70e r1e7b0db0  
    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<GAUSS_N; i++) { 
    31         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     30    for(int i=0; i<76; i++) { 
     31        const double theta = 0.5 * ( Gauss76Z[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<GAUSS_N; j++) { 
    39             double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     38        for(int j=0; j<76; j++) { 
     39            double phi = 0.5 * ( Gauss76Z[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 += GAUSS_W[j] * AP * AP; 
     47            inner_sum += Gauss76Wt[j] * AP * AP; 
    4848        } 
    4949        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
    50         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     50        outer_sum += Gauss76Wt[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  
    73 double 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). 
     14Note also that, contrary to :ref:`parallelepiped`, it does not compute 
     15the 2D scattering. 
    1416 
    1517 
    … …  
    2426that reference), with $\theta$ corresponding to $\alpha$ in that paper, 
    2527and not to the usual convention used for example in the 
    26 :ref:`parallelepiped` model. 
     28:ref:`parallelepiped` model. As the present model does not compute 
     29the 2D scattering, this has no further consequences. 
    2730 
    2831In this model the scattering from a massive parallelepiped with an 
    … …  
    6265units) *scale* represents the volume fraction (which is unitless). 
    6366 
    64 For 2d data the orientation of the particle is required, described using 
    65 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    66 of the calculation and angular dispersions see :ref:`orientation` . 
    67 The 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  
    70 For 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, 
    72 but 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  
     67**The 2D scattering intensity is not computed by this model.** 
    8668 
    8769 
    … …  
    126108              ["c2a_ratio", "", 1, [0, inf], "volume", 
    127109               "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"], 
    134110             ] 
    135111 

  • sasmodels/models/sc_paracrystal.c

    r108e70e rf728001  
    5454 
    5555    double outer_sum = 0.0; 
    56     for(int i=0; i<GAUSS_N; i++) { 
     56    for(int i=0; i<150; i++) { 
    5757        double inner_sum = 0.0; 
    58         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     58        const double theta = Gauss150Z[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<GAUSS_N;j++) { 
    64             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     63        for(int j=0;j<150;j++) { 
     64            const double phi = Gauss150Z[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 += GAUSS_W[j] * form; 
     70            inner_sum += Gauss150Wt[j] * form; 
    7171        } 
    7272        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    73         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     73        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    7474    } 
    7575    outer_sum *= theta_m; 
    … …  
    8282 
    8383static double 
    84 Iqabc(double qa, double qb, double qc, 
     84Iqxy(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

    r108e70e rbecded3  
    8181    double halfheight = 0.5*thick_core; 
    8282 
    83     for(int i=0; i<GAUSS_N; i++) { 
    84         double zi = (GAUSS_Z[i] + 1.0)*M_PI_4; 
     83    for(int i=0; i<N_POINTS_76; i++) { 
     84        double zi = (Gauss76Z[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 += GAUSS_W[i] * yyy * sin_alpha; 
     97        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    9898    } 
    9999 
    … …  
    142142 
    143143static double 
    144 Iqac(double qab, double qc, 
     144Iqxy(double qab, double qc, 
    145145    double thick_core, 
    146146    double thick_layer, 

  • sasmodels/models/triaxial_ellipsoid.c

    r108e70e rbecded3  
    2121    const double zb = M_PI_4; 
    2222    double outer = 0.0; 
    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; 
     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; 
    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<GAUSS_N;j++) { 
     31        for (int j=0;j<76;j++) { 
    3232            // translate a point in [-1,1] to a point in [0, 1] 
    33             const double usq = square(GAUSS_Z[j]*um + ub); 
     33            const double usq = square(Gauss76Z[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 += GAUSS_W[j] * fq * fq; 
     36            inner += Gauss76Wt[j] * fq * fq; 
    3737        } 
    38         outer += GAUSS_W[i] * inner;  // correcting for dx later 
     38        outer += Gauss76Wt[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 Iqabc(double qa, double qb, double qc, 
     48Iqxy(double qa, double qb, double qc, 
    4949    double sld, 
    5050    double sld_solvent, 

  • sasmodels/special.py

    rdf69efa re65c3ba  
    33................. 
    44 
    5 This following standard C99 math functions are available: 
     5The C code follows the C99 standard, with the usual math functions, 
     6as defined in 
     7`OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_. 
     8This includes the following: 
    69 
    710    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: 
    811        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ 
    912 
    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. 
     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. 
    1417 
    1518    sin, cos, tan, asin, acos, atan: 
    … …  
    2629        quadrants II and IV when $x$ and $y$ have opposite sign. 
    2730 
    28     fabs(x), fmin(x,y), fmax(x,y), trunc, rint: 
     31    fmin(x,y), fmax(x,y), trunc, rint: 
    2932        Floating point functions.  rint(x) returns the nearest integer. 
    3033 
    … …  
    3235        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that 
    3336        you cannot use :code:`x == NAN` to test for NaN values since that 
    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. 
     37        will always return false.  NAN does not equal NAN! 
    3638 
    3739    INFINITY: 
    … …  
    8789        for forcing a constant to stay double precision. 
    8890 
    89 The following special functions and scattering calculations are defined. 
     91The following special functions and scattering calculations are defined in 
     92`sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_. 
    9093These functions have been tuned to be fast and numerically stable down 
    9194to $q=0$ even in single precision.  In some cases they work around bugs 
    … …  
    181184 
    182185 
    183     gauss76.n, gauss76.z[i], gauss76.w[i]: 
     186    Gauss76Z[i], Gauss76Wt[i]: 
    184187        Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively, 
    185188        computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$. 
    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. 
     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 
    193193""" 
    194194# pylint: disable=unused-import 
    … …  
    200200 
    201201# C99 standard math library functions 
    202 from numpy import exp, log, power as pow, expm1, log1p, sqrt, cbrt 
     202from numpy import exp, log, power as pow, expm1, sqrt 
    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 fabs, fmin, fmax, trunc, rint 
    207 from numpy import pi, nan, inf 
     206from numpy import fmin, fmax, trunc, rint 
     207from numpy import NAN, inf as INFINITY 
     208 
    208209from scipy.special import gamma as sas_gamma 
    209210from scipy.special import erf as sas_erf 
    … …  
    217218# C99 standard math constants 
    218219M_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 
    219 NAN = nan 
    220 INFINITY = inf 
    221220 
    222221# non-standard constants 
    … …  
    227226    """return sin(x), cos(x)""" 
    228227    return sin(x), cos(x) 
    229 sincos = SINCOS 
    230228 
    231229def square(x): 
    … …  
    296294 
    297295# Gaussians 
    298 class Gauss: 
    299     def __init__(self, w, z): 
    300         self.n = len(w) 
    301         self.w = w 
    302         self.z = z 
    303  
    304 gauss20 = 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  
    351 gauss76 = 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  
    510 gauss150 = 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 ) 
     296 
     297Gauss20Wt = 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 
     320Gauss20Z = 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 
     343Gauss76Wt = 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 
     422Gauss76Z = 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 
     501Gauss150Z = 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 
     654Gauss150Wt = 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]) 
Note: See TracChangeset for help on using the changeset viewer.