Changes in / [c1bccff:924a119] in sasmodels


Ignore:
Files:
37 edited

Legend:

Unmodified
Added
Removed

  • doc/developer/overview.rst

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

  • doc/guide/orientation/orientation.rst

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

  • doc/guide/plugin.rst

    rd0dc9a3 r7e6bc45e  
    292292**Note: The order of the parameters in the definition will be the order of the 
    293293parameters in the user interface and the order of the parameters in Iq(), 
    294 Iqxy() and form_volume(). And** *scale* **and** *background* **parameters are 
    295 implicit to all models, so they do not need to be included in the parameter table.** 
     294Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background* 
     295**parameters are implicit to all models, so they do not need to be included 
     296in the parameter table.** 
    296297 
    297298- **"name"** is the name of the parameter shown on the FitPage. 
    … …  
    362363    scattered intensity. 
    363364 
    364   - "volume" parameters are passed to Iq(), Iqxy(), and form_volume(), and 
    365     have polydispersity loops generated automatically. 
    366  
    367   - "orientation" parameters are only passed to Iqxy(), and have angular 
    368     dispersion. 
     365  - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(), 
     366    and have polydispersity loops generated automatically. 
     367 
     368  - "orientation" parameters are not passed, but instead are combined with 
     369    orientation dispersity to translate *qx* and *qy* to *qa*, *qb* and *qc*. 
     370    These parameters should appear at the end of the table with the specific 
     371    names *theta*, *phi* and for asymmetric shapes *psi*, in that order. 
    369372 
    370373Some models will have integer parameters, such as number of pearls in the 
    … …  
    419422That is, the individual models do not need to include polydispersity 
    420423calculations, but instead rely on numerical integration to compute the 
    421 appropriately smeared pattern.   Angular dispersion values over polar angle 
    422 $\theta$ requires an additional $\cos \theta$ weighting due to decreased 
    423 arc length for the equatorial angle $\phi$ with increasing latitude. 
     424appropriately smeared pattern. 
    424425 
    425426Python Models 
    … …  
    468469barbell).  If I(q; pars) is NaN for any $q$, then those parameters will be 
    469470ignored, and not included in the calculation of the weighted polydispersity. 
    470  
    471 Similar to *Iq*, you can define *Iqxy(qx, qy, par1, par2, ...)* where the 
    472 parameter list includes any orientation parameters.  If *Iqxy* is not defined, 
    473 then it will default to *Iqxy = Iq(sqrt(qx**2+qy**2), par1, par2, ...)*. 
    474471 
    475472Models should define *form_volume(par1, par2, ...)* where the parameter 
    … …  
    497494    } 
    498495 
    499 *Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*, 
    500 and it includes orientation parameters. 
    501  
    502496*form_volume* defines the volume of the shape. As in python models, it 
    503497includes only the volume parameters. 
    504498 
    505 *Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and 
    506 *form_volume* will default to 1.0. 
    507  
    508499**source=['fn.c', ...]** includes the listed C source files in the 
    509 program before *Iq* and *Iqxy* are defined. This allows you to extend the 
    510 library of C functions available to your model. 
     500program before *Iq* and *form_volume* are defined. This allows you to 
     501extend the library of C functions available to your model. 
     502 
     503*c_code* includes arbitrary C code into your kernel, which can be 
     504handy for defining helper functions for *Iq* and *form_volume*. Note that 
     505you can put the full function definition for *Iq* and *form_volume* 
     506(include function declaration) into *c_code* as well, or put them into an 
     507external C file and add that file to the list of sources. 
    511508 
    512509Models are defined using double precision declarations for the 
    … …  
    532529 
    533530    #define INVALID(v) (v.bell_radius < v.radius) 
     531 
     532The INVALID define can go into *Iq*, or *c_code*, or an external C file 
     533listed in *source*. 
     534 
     535Oriented Shapes 
     536............... 
     537 
     538If the scattering is dependent on the orientation of the shape, then you 
     539will need to include *orientation* parameters *theta*, *phi* and *psi* 
     540at the end of the parameter table.  As described in the section 
     541:ref:`orientation`, the individual $(q_x, q_y)$ points on the detector will 
     542be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its 
     543canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the 
     544laboratory frame and beam travelling along $-z$. 
     545 
     546The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
     547*par1*, etc. are the parameters to the model.  If the shape is rotationally 
     548symmetric about *c* then *psi* is not needed, and the model is called 
     549as *Iqac(qab, qc, par1, par2, ...)*.  In either case, the orientation 
     550parameters are not included in the function call. 
     551 
     552For 1D oriented shapes, an integral over all angles is usually needed for 
     553the *Iq* function. Given symmetry and the substitution $u = \cos(\alpha)$, 
     554$du = -\sin(\alpha)\,d\alpha$ this becomes 
     555 
     556.. math:: 
     557 
     558    I(q) &= \frac{1}{4\pi} \int_{-\pi/2}^{pi/2} \int_{-pi}^{pi} 
     559            F(q_a, q_b, q_c)^2 \sin(\alpha)\,d\beta\,d\alpha \\ 
     560        &= \frac{8}{4\pi} \int_{0}^{pi/2} \int_{0}^{\pi/2} 
     561            F^2 \sin(\alpha)\,d\beta\,d\alpha \\ 
     562        &= \frac{8}{4\pi} \int_1^0 \int_{0}^{\pi/2} - F^2 \,d\beta\,du \\ 
     563        &= \frac{8}{4\pi} \int_0^1 \int_{0}^{\pi/2} F^2 \,d\beta\,du 
     564 
     565for 
     566 
     567.. math:: 
     568 
     569    q_a &= q \sin(\alpha)\sin(\beta) = q \sqrt{1-u^2} \sin(\beta) \\ 
     570    q_b &= q \sin(\alpha)\cos(\beta) = q \sqrt{1-u^2} \cos(\beta) \\ 
     571    q_c &= q \cos(\alpha) = q u 
     572 
     573Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the 
     574numerical integration is then:: 
     575 
     576    double outer_sum = 0.0; 
     577    for (int i = 0; i < GAUSS_N; i++) { 
     578        const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5; 
     579        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
     580        const double qc = cos_alpha * q; 
     581        double inner_sum = 0.0; 
     582        for (int j = 0; j < GAUSS_N; j++) { 
     583            const double beta = M_PI_4 * GAUSS_Z[j] + M_PI_4; 
     584            double sin_beta, cos_beta; 
     585            SINCOS(beta, sin_beta, cos_beta); 
     586            const double qa = sin_alpha * sin_beta * q; 
     587            const double qb = sin_alpha * cos_beta * q; 
     588            const double form = Fq(qa, qb, qc, ...); 
     589            inner_sum += GAUSS_W[j] * form * form; 
     590        } 
     591        outer_sum += GAUSS_W[i] * inner_sum; 
     592    } 
     593    outer_sum *= 0.25; // = 8/(4 pi) * outer_sum * (pi/2) / 4 
     594 
     595The *z* values for the Gauss-Legendre integration extends from -1 to 1, so 
     596the double sum of *w[i]w[j]* explains the factor of 4.  Correcting for the 
     597average *dz[i]dz[j]* gives $(1-0) \cdot (\pi/2-0) = \pi/2$.  The $8/(4 \pi)$ 
     598factor comes from the integral over the quadrant.  With less symmetry (eg., 
     599in the bcc and fcc paracrystal models), then an integral over the entire 
     600sphere may be necessary. 
     601 
     602For simpler models which are rotationally symmetric a single integral 
     603suffices: 
     604 
     605.. math:: 
     606 
     607    I(q) &= \frac{1}{\pi}\int_{-\pi/2}^{\pi/2} 
     608            F(q_{ab}, q_c)^2 \sin(\alpha)\,d\alpha/\pi \\ 
     609        &= \frac{2}{\pi} \int_0^1 F^2\,du 
     610 
     611for 
     612 
     613.. math:: 
     614 
     615    q_{ab} &= q \sin(\alpha) = q \sqrt{1 - u^2} \\ 
     616    q_c &= q \cos(\alpha) = q u 
     617 
     618 
     619with integration loop:: 
     620 
     621    double sum = 0.0; 
     622    for (int i = 0; i < GAUSS_N; i++) { 
     623        const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5; 
     624        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
     625        const double qab = sin_alpha * q; 
     626        const double qc = cos_alpha * q; 
     627        const double form = Fq(qab, qc, ...); 
     628        sum += GAUSS_W[j] * form * form; 
     629    } 
     630    sum *= 0.5; // = 2/pi * sum * (pi/2) / 2 
     631 
     632Magnetism 
     633......... 
     634 
     635Magnetism is supported automatically for all shapes by modifying the 
     636effective SLD of particle according to the Halpern-Johnson vector 
     637describing the interaction between neutron spin and magnetic field.  All 
     638parameters marked as type *sld* in the parameter table are treated as 
     639possibly magnetic particles with magnitude *M0* and direction 
     640*mtheta* and *mphi*.  Polarization parameters are also provided 
     641automatically for magnetic models to set the spin state of the measurement. 
     642 
     643For more complicated systems where magnetism is not uniform throughout 
     644the individual particles, you will need to write your own models. 
     645You should not mark the nuclear sld as type *sld*, but instead leave 
     646them unmarked and provide your own magnetism and polarization parameters. 
     647For 2D measurements you will need $(q_x, q_y)$ values for the measurement 
     648to compute the proper magnetism and orientation, which you can implement 
     649using *Iqxy(qx, qy, par1, par2, ...)*. 
    534650 
    535651Special Functions 
    … …  
    796912show a 50x improvement or more over the equivalent pure python model. 
    797913 
    798 External C Models 
    799 ................. 
    800  
    801 External C models are very much like embedded C models, except that 
    802 *Iq*, *Iqxy* and *form_volume* are defined in an external source file 
    803 loaded using the *source=[...]* statement. You need to supply the function 
    804 declarations for each of these that you need instead of building them 
    805 automatically from the parameter table. 
    806  
    807914 
    808915.. _Form_Factors: 
    … …  
    10061113variable name *Rg* for example because $R_g$ is the right name for the model 
    10071114parameter then ignore the lint errors.  Also, ignore *missing-docstring* 
    1008 for standard model functions *Iq*, *Iqxy*, etc. 
     1115for standard model functions *Iq*, *Iqac*, etc. 
    10091116 
    10101117We will have delinting sessions at the SasView Code Camps, where we can 

  • explore/jitter.py

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

  • explore/precision.py

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

  • explore/realspace.py

    ra1c32c2 r8fb2a94  
    4444    """ 
    4545    return Rx(phi)*Ry(theta)*Rz(psi) 
     46 
     47def 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 
     59def 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 
    4670 
    4771class Shape: 
    … …  
    219243        values = self.value.repeat(points.shape[0]) 
    220244        return values, self._adjust(points) 
     245 
     246NUMBA = False 
     247if 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 
     259def 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) 
    221273 
    222274def _calc_Pr_nonuniform(r, rho, points): 
    … …  
    239291            print("processing %d of %d"%(k, len(rho)-1)) 
    240292    Pr = extended_Pr[1:-1] 
    241     return Pr / Pr.max() 
    242  
    243 def 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  
     293    return Pr 
     294 
     295def _calc_Pr_uniform(r, rho, points): 
    249296    # Make Pr a little be bigger than necessary so that only distances 
    250297    # min < d < max end up in Pr 
    251     n_max = len(r) 
     298    dr, n_max = r[0], len(r) 
    252299    extended_Pr = np.zeros(n_max+1, 'd') 
    253300    t0 = time.clock() 
    254301    t_next = t0 + 3 
    255     row_zero = np.zeros(len(rho), 'i') 
    256302    for k, rho_k in enumerate(rho[:-1]): 
    257303        distances = sqrt(np.sum((points[k] - points[k+1:])**2, axis=1)) 
    258304        weights = rho_k * rho[k+1:] 
    259         index = np.minimum(np.asarray(distances/r[0], 'i'), n_max) 
     305        index = np.minimum(np.asarray(distances/dr, 'i'), n_max) 
    260306        # Note: indices may be duplicated, so "Pr[index] += w" will not work!! 
    261307        extended_Pr += np.bincount(index, weights, n_max+1) 
    … …  
    269315    # we are only accumulating the upper triangular distances. 
    270316    #Pr = Pr * 2 / len(rho)**2 
    271     return Pr / Pr.max() 
     317    return Pr 
    272318 
    273319    # Can get an additional 2x by going to C.  Cuda/OpenCL will allow even 
    … …  
    291337} 
    292338""" 
     339 
     340if 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 
     361def 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 
    293373 
    294374def j0(x): 
    … …  
    333413        plt.legend() 
    334414 
     415def 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 
    335429def plot_points(rho, points): 
    336430    import mpl_toolkits.mplot3d 
    … …  
    343437        pass 
    344438    n = len(points) 
    345     index = np.random.choice(n, size=1000) if n > 1000 else slice(None, None) 
     439    #print("len points", n) 
     440    index = np.random.choice(n, size=500) if n > 500 else slice(None, None) 
    346441    ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index]) 
    347442    #low, high = points.min(axis=0), points.max(axis=0) 
    348443    #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]]) 
    349444    ax.autoscale(True) 
     445 
     446def 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 
     463def 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 
     486def sas_sinx_x(x): 
     487    with np.errstate(all='ignore'): 
     488        retvalue = sin(x)/x 
     489    retvalue[x == 0.] = 1. 
     490    return retvalue 
    350491 
    351492def sas_2J1x_x(x): 
    … …  
    373514    Iq = Iq/Iq[0] 
    374515    return Iq 
     516 
     517def 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) 
    375523 
    376524def sphere_Iq(q, radius): 
    … …  
    415563    return Iq/Iq[0] 
    416564 
    417 def 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() 
     565def 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) 
    431585 
    432586def check_cylinder(radius=25, length=125, rho=2.): 
    … …  
    434588    fn = lambda q: cylinder_Iq(q, radius, length) 
    435589    check_shape(shape, fn) 
     590 
     591def 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 
     596def 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) 
    436617 
    437618def check_sphere(radius=125, rho=2): 
    … …  
    449630    side_c2 = copy(side_c).shift(0, 0, -c-dc) 
    450631    shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
    451     fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
    452     check_shape(shape, fn) 
     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) 
    453641 
    454642if __name__ == "__main__": 
    455643    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)) 
    456646    #check_sphere() 
    457647    #check_csbox() 

  • sasmodels/compare.py

    r2d81cfe r2a7e20e  
    9292    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh 
    9393    -neval=1 sets the number of evals for more accurate timing 
     94    -ngauss=0 overrides the number of points in the 1-D gaussian quadrature 
    9495 
    9596    === precision options === 

  • sasmodels/details.py

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

  • sasmodels/generate.py

    r2d81cfe r108e70e  
    77    particular dimensions averaged over all orientations. 
    88 
    9     *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form 
    10     with particular dimensions for a single orientation. 
    11  
    12     *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the 
    13     polarized neutron spin states (up-up, up-down, down-up, down-down) for 
    14     a form with particular dimensions for a single orientation. 
     9    *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc 
     10    for a rotationally symmetric form with particular dimensions. 
     11    qab, qc are determined from shape orientation and scattering angles. 
     12    This call is used if the shape has orientation parameters theta and phi. 
     13 
     14    *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc 
     15    for a form with particular dimensions.  qa, qb, qc are determined from 
     16    shape orientation and scattering angles. This call is used if the shape 
     17    has orientation parameters theta, phi and psi. 
     18 
     19    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy.  Use this 
     20    to create an arbitrary 2D theory function, needed for q-dependent 
     21    background functions and for models with non-uniform magnetism. 
    1522 
    1623    *form_volume(p1, p2, ...)* returns the volume of the form with particular 
    … …  
    3138scale and background parameters for each model. 
    3239 
    33 *Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 
    34 functions written for OpenCL.  All functions need prototype declarations 
    35 even if the are defined before they are used.  OpenCL does not support 
    36 *#include* preprocessor directives, so instead the list of includes needs 
    37 to be given as part of the metadata in the kernel module definition. 
    38 The included files should be listed using a path relative to the kernel 
    39 module, or if using "lib/file.c" if it is one of the standard includes 
    40 provided with the sasmodels source.  The includes need to be listed in 
    41 order so that functions are defined before they are used. 
     40C code should be stylized C-99 functions written for OpenCL. All functions 
     41need prototype declarations even if the are defined before they are used. 
     42Although OpenCL supports *#include* preprocessor directives, the list of 
     43includes should be given as part of the metadata in the kernel module 
     44definition. The included files should be listed using a path relative to the 
     45kernel module, or if using "lib/file.c" if it is one of the standard includes 
     46provided with the sasmodels source. The includes need to be listed in order 
     47so that functions are defined before they are used. 
    4248 
    4349Floating point values should be declared as *double*.  For single precision 
    … …  
    107113    present, the volume ratio is 1. 
    108114 
    109     *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the 
    110     C source code for the body of the volume, Iq, and Iqxy functions 
     115    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
     116    the C source code for the body of the volume, Iq, and Iqac functions 
    111117    respectively.  These can also be defined in the last source file. 
    112118 
    113     *Iq* and *Iqxy* also be instead be python functions defining the 
     119    *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the 
    114120    kernel.  If they are marked as *Iq.vectorized = True* then the 
    115121    kernel is passed the entire *q* vector at once, otherwise it is 
    … …  
    168174from zlib import crc32 
    169175from inspect import currentframe, getframeinfo 
     176import logging 
    170177 
    171178import numpy as np  # type: ignore 
    … …  
    181188    pass 
    182189# pylint: enable=unused-import 
     190 
     191logger = logging.getLogger(__name__) 
    183192 
    184193# jitter projection to use in the kernel code.  See explore/jitter.py 
    … …  
    626635 
    627636""" 
    628 def _gen_fn(name, pars, body, filename, line): 
    629     # type: (str, List[Parameter], str, str, int) -> str 
     637def _gen_fn(model_info, name, pars): 
     638    # type: (ModelInfo, str, List[Parameter]) -> str 
    630639    """ 
    631640    Generate a function given pars and body. 
    … …  
    639648    """ 
    640649    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) 
    641654    return _FN_TEMPLATE % { 
    642655        'name': name, 'pars': par_decl, 'body': body, 
    643         'filename': filename.replace('\\', '\\\\'), 'line': line, 
     656        'filename': filename.replace('\\', '\\\\'), 'line': lineno, 
    644657    } 
    645658 
    … …  
    656669 
    657670# type in IQXY pattern could be single, float, double, long double, ... 
    658 _IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)", 
     671_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
    659672                           flags=re.MULTILINE) 
    660 def _have_Iqxy(sources): 
     673def find_xy_mode(source): 
    661674    # type: (List[str]) -> bool 
    662675    """ 
    663     Return true if any file defines Iqxy. 
     676    Return the xy mode as qa, qac, qabc or qxy. 
    664677 
    665678    Note this is not a C parser, and so can be easily confused by 
    666679    non-standard syntax.  Also, it will incorrectly identify the following 
    667     as having Iqxy:: 
     680    as having 2D models:: 
    668681 
    669682        /* 
    670         double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     683        double Iqac(qab, qc, ...) { ... fill this in later ... } 
    671684        */ 
    672685 
    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  
    682 def _add_source(source, code, path): 
     686    If you want to comment out the function, use // on the front of the 
     687    line:: 
     688 
     689        /* 
     690        // double Iqac(qab, qc, ...) { ... fill this in later ... } 
     691        */ 
     692 
     693    """ 
     694    for code in source: 
     695        m = _IQXY_PATTERN.search(code) 
     696        if m is not None: 
     697            return m.group('mode') 
     698    return 'qa' 
     699 
     700 
     701def _add_source(source, code, path, lineno=1): 
    683702    """ 
    684703    Add a file to the list of source code chunks, tagged with path and line. 
    685704    """ 
    686705    path = path.replace('\\', '\\\\') 
    687     source.append('#line 1 "%s"' % path) 
     706    source.append('#line %d "%s"' % (lineno, path)) 
    688707    source.append(code) 
    689708 
    … …  
    716735    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    717736 
    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  
    723737    # Build initial sources 
    724738    source = [] 
    … …  
    727741        _add_source(source, code, path) 
    728742 
     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 
    729747    # Make parameters for q, qx, qy so that we can use them in declarations 
    730     q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     748    q, qx, qy, qab, qa, qb, qc \ 
     749        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
    731750    # Generate form_volume function, etc. from body only 
    732751    if isinstance(model_info.form_volume, str): 
    733752        pars = partable.form_volume_parameters 
    734         source.append(_gen_fn('form_volume', pars, model_info.form_volume, 
    735                               model_info.filename, model_info._form_volume_line)) 
     753        source.append(_gen_fn(model_info, 'form_volume', pars)) 
    736754    if isinstance(model_info.Iq, str): 
    737755        pars = [q] + partable.iq_parameters 
    738         source.append(_gen_fn('Iq', pars, model_info.Iq, 
    739                               model_info.filename, model_info._Iq_line)) 
     756        source.append(_gen_fn(model_info, 'Iq', pars)) 
    740757    if isinstance(model_info.Iqxy, str): 
    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)) 
     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") 
    744780 
    745781    # Define the parameter table 
    … …  
    767803    if xy_mode == 'qabc': 
    768804        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    769         call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
     805        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars 
    770806        clear_iqxy = "#undef CALL_IQ_ABC" 
    771807    elif xy_mode == 'qac': 
    772808        pars = ",".join(["_qa", "_qc"] + model_refs) 
    773         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
     809        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    774810        clear_iqxy = "#undef CALL_IQ_AC" 
    775     else:  # xy_mode == 'qa' 
     811    elif xy_mode == 'qa': 
    776812        pars = ",".join(["_qa"] + model_refs) 
    777813        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    778814        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 
    779827 
    780828    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 

  • sasmodels/kernel_header.c

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

  • sasmodels/kernel_iq.c

    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 
    3334//  INVALID(table) : test if the current point is feesible to calculate.  This 
    3435//      will be defined in the kernel definition file. 
    … …  
    469470  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    470471  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
    471 #endif 
    472  
    473 // Doing jitter projection code outside the previous if block so that we don't 
    474 // need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This 
    475 // will become more important if we implement more projections, or more 
    476 // complicated projections. 
    477 #if defined(CALL_IQ) || defined(CALL_IQ_A) 
     472#elif defined(CALL_IQ_XY) 
     473  // direct call to qx,qy calculator 
     474  double qx, qy; 
     475  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     476  #define BUILD_ROTATION() do {} while(0) 
     477  #define APPLY_ROTATION() do {} while(0) 
     478  #define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table) 
     479#endif 
     480 
     481// Define APPLY_PROJECTION depending on model symmetries. We do this outside 
     482// the previous if block so that we don't need to repeat the identical 
     483// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
     484// if we implement more projections, or more complicated projections. 
     485#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation 
    478486  #define APPLY_PROJECTION() const double weight=weight0 
    479 #else // !spherosymmetric projection 
     487#elif defined(CALL_IQ_XY) // pass orientation to the model 
     488  // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc 
     489  // Need to plug the values for the orientation angles back into parameter 
     490  // table in case they were overridden by the orientation offset.  This 
     491  // means that orientation dispersity will not work for these models, but 
     492  // it was broken anyway, so no matter.  Still want to provide Iqxy in case 
     493  // the user model wants full control of orientation/magnetism. 
     494  #if defined(HAVE_PSI) 
     495    const double theta = values[details->theta_par+2]; 
     496    const double phi = values[details->theta_par+3]; 
     497    const double psi = values[details->theta_par+4]; 
     498    double weight; 
     499    #define APPLY_PROJECTION() do { \ 
     500      local_values.table.theta = theta; \ 
     501      local_values.table.phi = phi; \ 
     502      local_values.table.psi = psi; \ 
     503      weight=weight0; \ 
     504    } while (0) 
     505  #elif defined(HAVE_THETA) 
     506    const double theta = values[details->theta_par+2]; 
     507    const double phi = values[details->theta_par+3]; 
     508    double weight; 
     509    #define APPLY_PROJECTION() do { \ 
     510      local_values.table.theta = theta; \ 
     511      local_values.table.phi = phi; \ 
     512      weight=weight0; \ 
     513    } while (0) 
     514  #else 
     515    #define APPLY_PROJECTION() const double weight=weight0 
     516  #endif 
     517#else // apply jitter and view before calling the model 
    480518  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 
    481519  const double theta = values[details->theta_par+2]; 
    … …  
    488526  // we go through the mesh. 
    489527  double dtheta, dphi, weight; 
    490   #if PROJECTION == 1 
     528  #if PROJECTION == 1 // equirectangular 
    491529    #define APPLY_PROJECTION() do { \ 
    492530      dtheta = local_values.table.theta; \ 
    … …  
    494532      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \ 
    495533    } while (0) 
    496   #elif PROJECTION == 2 
     534  #elif PROJECTION == 2 // sinusoidal 
    497535    #define APPLY_PROJECTION() do { \ 
    498536      dtheta = local_values.table.theta; \ 
    … …  
    504542    } while (0) 
    505543  #endif 
    506 #endif // !spherosymmetric projection 
     544#endif // done defining APPLY_PROJECTION 
    507545 
    508546// ** define looping macros ** 

  • sasmodels/kernelpy.py

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

  • sasmodels/modelinfo.py

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

  • sasmodels/models/_spherepy.py

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

  • sasmodels/models/barbell.c

    r74768cb r108e70e  
    9090 
    9191static double 
    92 Iqxy(double qab, double qc, 
     92Iqac(double qab, double qc, 
    9393    double sld, double solvent_sld, 
    9494    double radius_bell, double radius, double length) 

  • sasmodels/models/bcc_paracrystal.c

    r74768cb r108e70e  
    107107 
    108108 
    109 static double Iqxy(double qa, double qb, double qc, 
     109static double Iqabc(double qa, double qb, double qc, 
    110110    double dnn, double d_factor, double radius, 
    111111    double sld, double solvent_sld) 

  • sasmodels/models/capped_cylinder.c

    r74768cb r108e70e  
    115115 
    116116static double 
    117 Iqxy(double qab, double qc, 
     117Iqac(double qab, double qc, 
    118118    double sld, double solvent_sld, double radius, 
    119119    double radius_cap, double length) 

  • sasmodels/models/core_shell_bicelle.c

    r74768cb r108e70e  
    6767 
    6868static double 
    69 Iqxy(double qab, double qc, 
     69Iqac(double qab, double qc, 
    7070    double radius, 
    7171    double thick_rim, 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r74768cb r108e70e  
    7171 
    7272static double 
    73 Iqxy(double qa, double qb, double qc, 
     73Iqabc(double qa, double qb, double qc, 
    7474    double r_minor, 
    7575    double x_core, 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r74768cb r108e70e  
    7979 
    8080static double 
    81 Iqxy(double qa, double qb, double qc, 
     81Iqabc(double qa, double qb, double qc, 
    8282          double r_minor, 
    8383          double x_core, 
    … …  
    114114    return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 
    115115} 
    116  

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

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

  • sasmodels/models/core_shell_cylinder.c

    r74768cb r108e70e  
    4848 
    4949 
    50 double Iqxy(double qab, double qc, 
     50double Iqac(double qab, double qc, 
    5151    double core_sld, 
    5252    double shell_sld, 

  • sasmodels/models/core_shell_ellipsoid.c

    r74768cb r108e70e  
    7575 
    7676static double 
    77 Iqxy(double qab, double qc, 
     77Iqac(double qab, double qc, 
    7878    double radius_equat_core, 
    7979    double x_core, 

  • sasmodels/models/core_shell_parallelepiped.c

    r4493288 re077231  
    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 Singhs thesis now (Dirk Honecker) 
    46     // Code rewritten (PAK) 
     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) 
    4747 
    4848    const double half_q = 0.5*q; 
    … …  
    102102 
    103103static double 
    104 Iqxy(double qa, double qb, double qc, 
     104Iqabc(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. 
    125123    const double tA = length_a + 2.0*thick_rim_a; 
    126124    const double tB = length_b + 2.0*thick_rim_b; 

  • sasmodels/models/core_shell_parallelepiped.py

    r10ee838 r97be877  
    136136* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    137137* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
    138 * **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017 
    139 * **Currently Under review by:** Paul Butler 
     138* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 
     139* Cross-checked against hollow rectangular prism and rectangular prism for 
     140  equal thickness overlapping sides, and by Monte Carlo sampling of points 
     141  within the shape for non-uniform, non-overlapping sides. 
    140142""" 
    141143 

  • sasmodels/models/cylinder.c

    r74768cb r108e70e  
    4545 
    4646static double 
    47 Iqxy(double qab, double qc, 
     47Iqac(double qab, double qc, 
    4848    double sld, 
    4949    double solvent_sld, 

  • sasmodels/models/ellipsoid.c

    r74768cb r108e70e  
    3939 
    4040static double 
    41 Iqxy(double qab, double qc, 
     41Iqac(double qab, double qc, 
    4242    double sld, 
    4343    double sld_solvent, 

  • sasmodels/models/elliptical_cylinder.c

    r74768cb r108e70e  
    5454 
    5555static double 
    56 Iqxy(double qa, double qb, double qc, 
     56Iqabc(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

    r74768cb r108e70e  
    8080 
    8181 
    82 static double Iqxy(double qa, double qb, double qc, 
     82static double Iqabc(double qa, double qb, double qc, 
    8383    double dnn, double d_factor, double radius, 
    8484    double sld, double solvent_sld) 

  • sasmodels/models/hollow_cylinder.c

    r74768cb r108e70e  
    5252 
    5353static double 
    54 Iqxy(double qab, double qc, 
     54Iqac(double qab, double qc, 
    5555    double radius, double thickness, double length, 
    5656    double sld, double solvent_sld) 

  • sasmodels/models/hollow_rectangular_prism.c

    r74768cb r108e70e  
    8585} 
    8686 
    87 double Iqxy(double qa, double qb, double qc, 
     87double Iqabc(double qa, double qb, double qc, 
    8888    double sld, 
    8989    double solvent_sld, 

  • sasmodels/models/line.py

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

  • sasmodels/models/parallelepiped.c

    r74768cb r108e70e  
    5353 
    5454static double 
    55 Iqxy(double qa, double qb, double qc, 
     55Iqabc(double qa, double qb, double qc, 
    5656    double sld, 
    5757    double solvent_sld, 

  • sasmodels/models/rectangular_prism.c

    r74768cb r108e70e  
    7171 
    7272 
    73 double Iqxy(double qa, double qb, double qc, 
     73double Iqabc(double qa, double qb, double qc, 
    7474    double sld, 
    7575    double solvent_sld, 

  • sasmodels/models/sc_paracrystal.c

    r74768cb r108e70e  
    8282 
    8383static double 
    84 Iqxy(double qa, double qb, double qc, 
     84Iqabc(double qa, double qb, double qc, 
    8585    double dnn, double d_factor, double radius, 
    8686    double sld, double solvent_sld) 

  • sasmodels/models/stacked_disks.c

    r74768cb r108e70e  
    142142 
    143143static double 
    144 Iqxy(double qab, double qc, 
     144Iqac(double qab, double qc, 
    145145    double thick_core, 
    146146    double thick_layer, 

  • sasmodels/models/triaxial_ellipsoid.c

    r74768cb r108e70e  
    4646 
    4747static double 
    48 Iqxy(double qa, double qb, double qc, 
     48Iqabc(double qa, double qb, double qc, 
    4949    double sld, 
    5050    double sld_solvent, 
Note: See TracChangeset for help on using the changeset viewer.