Changeset 72e41a0 in sasmodels for doc


Ignore:
Timestamp:
Jan 17, 2018 4:14:03 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
c11d09f
Parents:
92d330fd (diff), 1258e32 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into boltzmann

Location:
doc
Files:
2 added
4 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

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

  • doc/guide/pd/polydispersity.rst

    reda8b30 r92d330fd  
    4242calculations are generally more robust with more data points or more angles. 
    4343 
    44 The following five distribution functions are provided: 
     44The following distribution functions are provided: 
    4545 
    4646*  *Rectangular Distribution* 
     47*  *Uniform Distribution* 
    4748*  *Gaussian Distribution* 
    4849*  *Lognormal Distribution* 
    4950*  *Schulz Distribution* 
    5051*  *Array Distribution* 
     52*  *Boltzmann Distribution* 
    5153 
    5254These are all implemented as *number-average* distributions. 
    … …  
    8587    Rectangular distribution. 
    8688 
     89 
     90 
     91Uniform Distribution 
     92^^^^^^^^^^^^^^^^^^^^^^^^ 
     93 
     94The Uniform Distribution is defined as 
     95 
     96    .. math:: 
     97 
     98        f(x) = \frac{1}{\text{Norm}} 
     99        \begin{cases} 
     100          1 & \text{for } |x - \bar x| \leq \sigma \\ 
     101          0 & \text{for } |x - \bar x| > \sigma 
     102        \end{cases} 
     103 
     104    where $\bar x$ is the mean of the distribution, $\sigma$ is the half-width, and 
     105    *Norm* is a normalization factor which is determined during the numerical 
     106    calculation. 
     107 
     108    Note that the polydispersity is given by 
     109 
     110    .. math:: \text{PD} = \sigma / \bar x 
     111 
     112    .. figure:: pd_uniform.jpg 
     113 
     114        Uniform distribution. 
     115 
     116The value $N_\sigma$ is ignored for this distribution. 
     117 
    87118.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    88119 
    … …  
    183214^^^^^^^^^^^^^^^^^^ 
    184215 
    185 This user-definable distribution should be given as as a simple ASCII text 
     216This user-definable distribution should be given as a simple ASCII text 
    186217file where the array is defined by two columns of numbers: $x$ and $f(x)$. 
    187218The $f(x)$ will be normalized to 1 during the computation. 
    … …  
    202233given for the model will have no affect, and will be ignored when computing 
    203234the average.  This means that any parameter with an array distribution will 
    204 not be fittable. 
     235not be fitable. 
     236 
     237.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     238 
     239Boltzmann Distribution 
     240^^^^^^^^^^^^^^^^^^^^^^ 
     241 
     242The Boltzmann Distribution is defined as 
     243 
     244.. math:: 
     245 
     246    f(x) = \frac{1}{\text{Norm}} 
     247           \exp\left(-\frac{ | x - \bar x | }{\sigma}\right) 
     248 
     249where $\bar x$ is the mean of the distribution and *Norm* is a normalization 
     250factor which is determined during the numerical calculation. 
     251The width is defined as 
     252 
     253.. math:: \sigma=\frac{k T}{E} 
     254 
     255which is the inverse Boltzmann factor, 
     256where $k$ is the Boltzmann constant, $T$ the temperature in Kelvin and $E$ a 
     257characteristic energy per particle. 
     258 
     259.. figure:: pd_boltzmann.jpg 
     260 
     261    Boltzmann distribution. 
    205262 
    206263.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
Note: See TracChangeset for help on using the changeset viewer.