Changes in / [924a119:c1bccff] in sasmodels


Ignore:
Files:
37 edited

Legend:

Unmodified
Added
Removed
  • 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 rd0dc9a3  
    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 
     
    529532 
    530533    #define INVALID(v) (v.bell_radius < v.radius) 
    531  
    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, ...)*. 
    650534 
    651535Special Functions 
     
    912796show a 50x improvement or more over the equivalent pure python model. 
    913797 
     798External C Models 
     799................. 
     800 
     801External C models are very much like embedded C models, except that 
     802*Iq*, *Iqxy* and *form_volume* are defined in an external source file 
     803loaded using the *source=[...]* statement. You need to supply the function 
     804declarations for each of these that you need instead of building them 
     805automatically from the parameter table. 
     806 
    914807 
    915808.. _Form_Factors: 
     
    11131006variable name *Rg* for example because $R_g$ is the right name for the model 
    11141007parameter then ignore the lint errors.  Also, ignore *missing-docstring* 
    1115 for standard model functions *Iq*, *Iqac*, etc. 
     1008for standard model functions *Iq*, *Iqxy*, etc. 
    11161009 
    11171010We will have delinting sessions at the SasView Code Camps, where we can 
  • 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) 
  • explore/realspace.py

    r8fb2a94 ra1c32c2  
    4444    """ 
    4545    return Rx(phi)*Ry(theta)*Rz(psi) 
    46  
    47 def apply_view(points, view): 
    48     """ 
    49     Apply the view transform (theta, phi, psi) to a set of points. 
    50  
    51     Points are stored in a 3 x n numpy array. 
    52  
    53     View angles are in degrees. 
    54     """ 
    55     theta, phi, psi = view 
    56     return np.asarray((Rz(phi)*Ry(theta)*Rz(psi))*np.matrix(points.T)).T 
    57  
    58  
    59 def invert_view(qx, qy, view): 
    60     """ 
    61     Return (qa, qb, qc) for the (theta, phi, psi) view angle at detector 
    62     pixel (qx, qy). 
    63  
    64     View angles are in degrees. 
    65     """ 
    66     theta, phi, psi = view 
    67     q = np.vstack((qx.flatten(), qy.flatten(), 0*qx.flatten())) 
    68     return np.asarray((Rz(-psi)*Ry(-theta)*Rz(-phi))*np.matrix(q)) 
    69  
    7046 
    7147class Shape: 
     
    243219        values = self.value.repeat(points.shape[0]) 
    244220        return values, self._adjust(points) 
    245  
    246 NUMBA = False 
    247 if NUMBA: 
    248     from numba import njit 
    249     @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])") 
    250     def _Iqxy(values, x, y, z, qa, qb, qc): 
    251         Iq = np.zeros_like(qa) 
    252         for j in range(len(Iq)): 
    253             total = 0. + 0j 
    254             for k in range(len(Iq)): 
    255                 total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k])) 
    256             Iq[j] = abs(total)**2 
    257         return Iq 
    258  
    259 def calc_Iqxy(qx, qy, rho, points, volume=1.0, view=(0, 0, 0)): 
    260     qx, qy = np.broadcast_arrays(qx, qy) 
    261     qa, qb, qc = invert_view(qx, qy, view) 
    262     rho, volume = np.broadcast_arrays(rho, volume) 
    263     values = rho*volume 
    264     x, y, z = points.T 
    265  
    266     # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r) 
    267     if NUMBA: 
    268         Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 
    269     else: 
    270         Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 
    271               for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 
    272     return np.asarray(Iq).reshape(qx.shape) / np.sum(volume) 
    273221 
    274222def _calc_Pr_nonuniform(r, rho, points): 
     
    291239            print("processing %d of %d"%(k, len(rho)-1)) 
    292240    Pr = extended_Pr[1:-1] 
    293     return Pr 
    294  
    295 def _calc_Pr_uniform(r, rho, points): 
     241    return Pr / Pr.max() 
     242 
     243def calc_Pr(r, rho, points): 
     244    # P(r) with uniform steps in r is 3x faster; check if we are uniform 
     245    # before continuing 
     246    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 
     247        return _calc_Pr_nonuniform(r, rho, points) 
     248 
    296249    # Make Pr a little be bigger than necessary so that only distances 
    297250    # min < d < max end up in Pr 
    298     dr, n_max = r[0], len(r) 
     251    n_max = len(r) 
    299252    extended_Pr = np.zeros(n_max+1, 'd') 
    300253    t0 = time.clock() 
    301254    t_next = t0 + 3 
     255    row_zero = np.zeros(len(rho), 'i') 
    302256    for k, rho_k in enumerate(rho[:-1]): 
    303257        distances = sqrt(np.sum((points[k] - points[k+1:])**2, axis=1)) 
    304258        weights = rho_k * rho[k+1:] 
    305         index = np.minimum(np.asarray(distances/dr, 'i'), n_max) 
     259        index = np.minimum(np.asarray(distances/r[0], 'i'), n_max) 
    306260        # Note: indices may be duplicated, so "Pr[index] += w" will not work!! 
    307261        extended_Pr += np.bincount(index, weights, n_max+1) 
     
    315269    # we are only accumulating the upper triangular distances. 
    316270    #Pr = Pr * 2 / len(rho)**2 
    317     return Pr 
     271    return Pr / Pr.max() 
    318272 
    319273    # Can get an additional 2x by going to C.  Cuda/OpenCL will allow even 
     
    337291} 
    338292""" 
    339  
    340 if NUMBA: 
    341     @njit("f8[:](f8[:], f8[:], f8[:,:])") 
    342     def _calc_Pr_uniform_jit(r, rho, points): 
    343         dr = r[0] 
    344         n_max = len(r) 
    345         Pr = np.zeros_like(r) 
    346         for j in range(len(rho) - 1): 
    347             x, y, z = points[j, 0], points[j, 1], points[j, 2] 
    348             for k in range(j+1, len(rho)): 
    349                 distance = sqrt((x - points[k, 0])**2 
    350                                 + (y - points[k, 1])**2 
    351                                 + (z - points[k, 2])**2) 
    352                 index = int(distance/dr) 
    353                 if index < n_max: 
    354                     Pr[index] += rho[j] * rho[k] 
    355         # Make Pr independent of sampling density.  The factor of 2 comes because 
    356         # we are only accumulating the upper triangular distances. 
    357         #Pr = Pr * 2 / len(rho)**2 
    358         return Pr 
    359  
    360  
    361 def calc_Pr(r, rho, points): 
    362     # P(r) with uniform steps in r is 3x faster; check if we are uniform 
    363     # before continuing 
    364     if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 
    365         Pr = _calc_Pr_nonuniform(r, rho, points) 
    366     else: 
    367         if NUMBA: 
    368             Pr = _calc_Pr_uniform_jit(r, rho, points) 
    369         else: 
    370             Pr = _calc_Pr_uniform(r, rho, points) 
    371     return Pr / Pr.max() 
    372  
    373293 
    374294def j0(x): 
     
    413333        plt.legend() 
    414334 
    415 def plot_calc_2d(qx, qy, Iqxy, theory=None): 
    416     import matplotlib.pyplot as plt 
    417     qx, qy = bin_edges(qx), bin_edges(qy) 
    418     #qx, qy = np.meshgrid(qx, qy) 
    419     if theory is not None: 
    420         plt.subplot(121) 
    421     plt.pcolormesh(qx, qy, np.log10(Iqxy)) 
    422     plt.xlabel('qx (1/A)') 
    423     plt.ylabel('qy (1/A)') 
    424     if theory is not None: 
    425         plt.subplot(122) 
    426         plt.pcolormesh(qx, qy, np.log10(theory)) 
    427         plt.xlabel('qx (1/A)') 
    428  
    429335def plot_points(rho, points): 
    430336    import mpl_toolkits.mplot3d 
     
    437343        pass 
    438344    n = len(points) 
    439     #print("len points", n) 
    440     index = np.random.choice(n, size=500) if n > 500 else slice(None, None) 
     345    index = np.random.choice(n, size=1000) if n > 1000 else slice(None, None) 
    441346    ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index]) 
    442347    #low, high = points.min(axis=0), points.max(axis=0) 
    443348    #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]]) 
    444349    ax.autoscale(True) 
    445  
    446 def check_shape(shape, fn=None): 
    447     rho_solvent = 0 
    448     q = np.logspace(-3, 0, 200) 
    449     r = shape.r_bins(q, r_step=0.01) 
    450     sampling_density = 6*5000 / shape.volume() 
    451     rho, points = shape.sample(sampling_density) 
    452     t0 = time.time() 
    453     Pr = calc_Pr(r, rho-rho_solvent, points) 
    454     print("calc Pr time", time.time() - t0) 
    455     Iq = calc_Iq(q, r, Pr) 
    456     theory = (q, fn(q)) if fn is not None else None 
    457  
    458     import pylab 
    459     #plot_points(rho, points); pylab.figure() 
    460     plot_calc(r, Pr, q, Iq, theory=theory) 
    461     pylab.show() 
    462  
    463 def check_shape_2d(shape, fn=None, view=(0, 0, 0)): 
    464     rho_solvent = 0 
    465     nq, qmax = 100, 1.0 
    466     qx = np.linspace(0.0, qmax, nq) 
    467     qy = np.linspace(0.0, qmax, nq) 
    468     Qx, Qy = np.meshgrid(qx, qy) 
    469     sampling_density = 50000 / shape.volume() 
    470     #t0 = time.time() 
    471     rho, points = shape.sample(sampling_density) 
    472     #print("sample time", time.time() - t0) 
    473     t0 = time.time() 
    474     Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
    475     print("calc time", time.time() - t0) 
    476     theory = fn(Qx, Qy) if fn is not None else None 
    477     Iqxy += 0.001 * Iqxy.max() 
    478     if theory is not None: 
    479         theory += 0.001 * theory.max() 
    480  
    481     import pylab 
    482     #plot_points(rho, points); pylab.figure() 
    483     plot_calc_2d(qx, qy, Iqxy, theory=theory) 
    484     pylab.show() 
    485  
    486 def sas_sinx_x(x): 
    487     with np.errstate(all='ignore'): 
    488         retvalue = sin(x)/x 
    489     retvalue[x == 0.] = 1. 
    490     return retvalue 
    491350 
    492351def sas_2J1x_x(x): 
     
    514373    Iq = Iq/Iq[0] 
    515374    return Iq 
    516  
    517 def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)): 
    518     qa, qb, qc = invert_view(qx, qy, view) 
    519     qab = np.sqrt(qa**2 + qb**2) 
    520     Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2) 
    521     Iq = Fq**2 
    522     return Iq.reshape(qx.shape) 
    523375 
    524376def sphere_Iq(q, radius): 
     
    563415    return Iq/Iq[0] 
    564416 
    565 def csbox_Iqxy(qx, qy, a, b, c, da, db, dc, slda, sldb, sldc, sld_core, view=(0,0,0)): 
    566     qa, qb, qc = invert_view(qx, qy, view) 
    567  
    568     sld_solvent = 0 
    569     overlapping = False 
    570     dr0 = sld_core - sld_solvent 
    571     drA, drB, drC = slda-sld_solvent, sldb-sld_solvent, sldc-sld_solvent 
    572     tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc 
    573     siA = a*sas_sinx_x(a*qa/2) 
    574     siB = b*sas_sinx_x(b*qb/2) 
    575     siC = c*sas_sinx_x(c*qc/2) 
    576     siAt = tA*sas_sinx_x(tA*qa/2) 
    577     siBt = tB*sas_sinx_x(tB*qb/2) 
    578     siCt = tC*sas_sinx_x(tC*qc/2) 
    579     Fq = (dr0*siA*siB*siC 
    580           + drA*(siAt-siA)*siB*siC 
    581           + drB*siA*(siBt-siB)*siC 
    582           + drC*siA*siB*(siCt-siC)) 
    583     Iq = Fq**2 
    584     return Iq.reshape(qx.shape) 
     417def check_shape(shape, fn=None): 
     418    rho_solvent = 0 
     419    q = np.logspace(-3, 0, 200) 
     420    r = shape.r_bins(q, r_step=0.01) 
     421    sampling_density = 15000 / shape.volume() 
     422    rho, points = shape.sample(sampling_density) 
     423    Pr = calc_Pr(r, rho-rho_solvent, points) 
     424    Iq = calc_Iq(q, r, Pr) 
     425    theory = (q, fn(q)) if fn is not None else None 
     426 
     427    import pylab 
     428    #plot_points(rho, points); pylab.figure() 
     429    plot_calc(r, Pr, q, Iq, theory=theory) 
     430    pylab.show() 
    585431 
    586432def check_cylinder(radius=25, length=125, rho=2.): 
     
    588434    fn = lambda q: cylinder_Iq(q, radius, length) 
    589435    check_shape(shape, fn) 
    590  
    591 def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)): 
    592     shape = EllipticalCylinder(radius, radius, length, rho) 
    593     fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    594     check_shape_2d(shape, fn, view=view) 
    595  
    596 def check_cylinder_2d_lattice(radius=25, length=125, rho=2., 
    597                               view=(0, 0, 0)): 
    598     nx, dx = 1, 2*radius 
    599     ny, dy = 30, 2*radius 
    600     nz, dz = 30, length 
    601     dx, dy, dz = 2*dx, 2*dy, 2*dz 
    602     def center(*args): 
    603         sigma = 0.333 
    604         space = 2 
    605         return [(space*n+np.random.randn()*sigma)*x for n, x in args] 
    606     shapes = [EllipticalCylinder(radius, radius, length, rho, 
    607                                  #center=(ix*dx, iy*dy, iz*dz) 
    608                                  orientation=np.random.randn(3)*0, 
    609                                  center=center((ix, dx), (iy, dy), (iz, dz)) 
    610                                 ) 
    611               for ix in range(nx) 
    612               for iy in range(ny) 
    613               for iz in range(nz)] 
    614     shape = Composite(shapes) 
    615     fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    616     check_shape_2d(shape, fn, view=view) 
    617436 
    618437def check_sphere(radius=125, rho=2): 
     
    630449    side_c2 = copy(side_c).shift(0, 0, -c-dc) 
    631450    shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
    632     def fn(q): 
    633         return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
    634     #check_shape(shape, fn) 
    635  
    636     view = (20, 30, 40) 
    637     def fn_xy(qx, qy): 
    638         return csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 
    639                           slda, sldb, sldc, sld_core, view=view) 
    640     check_shape_2d(shape, fn_xy, view=view) 
     451    fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
     452    check_shape(shape, fn) 
    641453 
    642454if __name__ == "__main__": 
    643455    check_cylinder(radius=10, length=40) 
    644     #check_cylinder_2d(radius=10, length=40, view=(90,30,0)) 
    645     #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0)) 
    646456    #check_sphere() 
    647457    #check_csbox() 
  • sasmodels/compare.py

    r2a7e20e r2d81cfe  
    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 === 
  • 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 
     
    635626 
    636627""" 
    637 def _gen_fn(model_info, name, pars): 
    638     # type: (ModelInfo, str, List[Parameter]) -> str 
     628def _gen_fn(name, pars, body, filename, line): 
     629    # type: (str, List[Parameter], str, str, int) -> str 
    639630    """ 
    640631    Generate a function given pars and body. 
     
    648639    """ 
    649640    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) 
    654641    return _FN_TEMPLATE % { 
    655642        'name': name, 'pars': par_decl, 'body': body, 
    656         'filename': filename.replace('\\', '\\\\'), 'line': lineno, 
     643        'filename': filename.replace('\\', '\\\\'), 'line': line, 
    657644    } 
    658645 
     
    669656 
    670657# 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*[(]", 
     658_IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)", 
    672659                           flags=re.MULTILINE) 
    673 def find_xy_mode(source): 
     660def _have_Iqxy(sources): 
    674661    # type: (List[str]) -> bool 
    675662    """ 
    676     Return the xy mode as qa, qac, qabc or qxy. 
     663    Return true if any file defines Iqxy. 
    677664 
    678665    Note this is not a C parser, and so can be easily confused by 
    679666    non-standard syntax.  Also, it will incorrectly identify the following 
    680     as having 2D models:: 
     667    as having Iqxy:: 
    681668 
    682669        /* 
    683         double Iqac(qab, qc, ...) { ... fill this in later ... } 
     670        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
    684671        */ 
    685672 
    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): 
     673    If you want to comment out an Iqxy function, use // on the front of the 
     674    line instead. 
     675    """ 
     676    for _path, code in sources: 
     677        if _IQXY_PATTERN.search(code): 
     678            return True 
     679    return False 
     680 
     681 
     682def _add_source(source, code, path): 
    702683    """ 
    703684    Add a file to the list of source code chunks, tagged with path and line. 
    704685    """ 
    705686    path = path.replace('\\', '\\\\') 
    706     source.append('#line %d "%s"' % (lineno, path)) 
     687    source.append('#line 1 "%s"' % path) 
    707688    source.append(code) 
    708689 
     
    735716    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    736717 
     718    # What kind of 2D model do we need? 
     719    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str) 
     720               else 'qac' if not partable.is_asymmetric 
     721               else 'qabc') 
     722 
    737723    # Build initial sources 
    738724    source = [] 
     
    741727        _add_source(source, code, path) 
    742728 
    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  
    747729    # 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()] 
     730    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
    750731    # Generate form_volume function, etc. from body only 
    751732    if isinstance(model_info.form_volume, str): 
    752733        pars = partable.form_volume_parameters 
    753         source.append(_gen_fn(model_info, 'form_volume', pars)) 
     734        source.append(_gen_fn('form_volume', pars, model_info.form_volume, 
     735                              model_info.filename, model_info._form_volume_line)) 
    754736    if isinstance(model_info.Iq, str): 
    755737        pars = [q] + partable.iq_parameters 
    756         source.append(_gen_fn(model_info, 'Iq', pars)) 
     738        source.append(_gen_fn('Iq', pars, model_info.Iq, 
     739                              model_info.filename, model_info._Iq_line)) 
    757740    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") 
     741        pars = [qx, qy] + partable.iqxy_parameters 
     742        source.append(_gen_fn('Iqxy', pars, model_info.Iqxy, 
     743                              model_info.filename, model_info._Iqxy_line)) 
    780744 
    781745    # Define the parameter table 
     
    803767    if xy_mode == 'qabc': 
    804768        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    805         call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars 
     769        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
    806770        clear_iqxy = "#undef CALL_IQ_ABC" 
    807771    elif xy_mode == 'qac': 
    808772        pars = ",".join(["_qa", "_qc"] + model_refs) 
    809         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
     773        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
    810774        clear_iqxy = "#undef CALL_IQ_AC" 
    811     elif xy_mode == 'qa': 
     775    else:  # xy_mode == 'qa' 
    812776        pars = ",".join(["_qa"] + model_refs) 
    813777        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    814778        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  
    827779 
    828780    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 rec8d4ac  
    3131//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
    3232//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
    33 //  CALL_IQ_XY(qx, qy, table) : call the Iqxy function for arbitrary models 
    3433//  INVALID(table) : test if the current point is feesible to calculate.  This 
    3534//      will be defined in the kernel definition file. 
     
    470469  #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 r74768cb  
    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 r74768cb  
    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 r74768cb  
    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 r74768cb  
    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 r74768cb  
    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 r74768cb  
    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 r74768cb  
    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 r74768cb  
    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 r4493288  
    4343    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
    4444    // 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) 
     45    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     46    // Code rewritten (PAK) 
    4747 
    4848    const double half_q = 0.5*q; 
     
    102102 
    103103static double 
    104 Iqabc(double qa, double qb, double qc, 
     104Iqxy(double qa, double qb, double qc, 
    105105    double core_sld, 
    106106    double arim_sld, 
     
    121121    const double drC = crim_sld-solvent_sld; 
    122122 
     123    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
     124    // the scaling by B. 
    123125    const double tA = length_a + 2.0*thick_rim_a; 
    124126    const double tB = length_b + 2.0*thick_rim_b; 
  • sasmodels/models/core_shell_parallelepiped.py

    r97be877 r10ee838  
    136136* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    137137* **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. 
     138* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017 
     139* **Currently Under review by:** Paul Butler 
    142140""" 
    143141 
  • sasmodels/models/cylinder.c

    r108e70e r74768cb  
    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 r74768cb  
    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 r74768cb  
    5454 
    5555static double 
    56 Iqabc(double qa, double qb, double qc, 
     56Iqxy(double qa, double qb, double qc, 
    5757     double radius_minor, double r_ratio, double length, 
    5858     double sld, double solvent_sld) 
  • sasmodels/models/fcc_paracrystal.c

    r108e70e r74768cb  
    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/hollow_cylinder.c

    r108e70e r74768cb  
    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 r74768cb  
    8585} 
    8686 
    87 double Iqabc(double qa, double qb, double qc, 
     87double Iqxy(double qa, double qb, double qc, 
    8888    double sld, 
    8989    double solvent_sld, 
  • 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 r74768cb  
    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/rectangular_prism.c

    r108e70e r74768cb  
    7171 
    7272 
    73 double Iqabc(double qa, double qb, double qc, 
     73double Iqxy(double qa, double qb, double qc, 
    7474    double sld, 
    7575    double solvent_sld, 
  • sasmodels/models/sc_paracrystal.c

    r108e70e r74768cb  
    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 r74768cb  
    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 r74768cb  
    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, 
Note: See TracChangeset for help on using the changeset viewer.