Changes in / [924a119:c1bccff] in sasmodels
- Files:
-
- 37 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/developer/overview.rst
r2a7e20e r3d40839 185 185 jitter applied before particle orientation. 186 186 187 When computing the orientation dispersity integral, the weights for188 the individual points depends on the map projection used to translate jitter189 angles into latitude/longitude. The choice of projection is set by190 *sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for191 equirectangular and *PROJECTION=2* for sinusoidal. The more complicated192 Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh193 for details.194 195 187 For numerical integration within form factors etc. sasmodels is mostly using 196 188 Gaussian quadrature with 20, 76 or 150 points depending on the model. It also … … 207 199 Useful testing routines include: 208 200 209 The *sascomp* utility is used to view and compare models with different 210 parameters and calculation engines. The usual case is to simply plot a 211 model that you are developing:: 212 213 python sascomp path/to/model.py 214 215 Once the obvious problems are addressed, check the numerical precision 216 across a variety of randomly generated inputs:: 217 218 python sascomp -engine=single,double path/to/model.py -sets=10 219 220 You can compare different parameter values for the same or different models. 221 For example when looking along the long axis of a cylinder ($\theta=0$), 222 dispersity in $\theta$ should be equivalent to dispersity in $\phi$ 223 when $\phi=90$:: 224 225 python sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle \\ 226 phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10 227 228 It turns out that they are not because the equirectangular map projection 229 weights the points by $\cos(\theta)$ so $\Delta\theta$ is not identical 230 to $\Delta\phi$. Setting *PROJECTION=2* in :mod:`sasmodels.generate` helps 231 somewhat. Postel would help even more in this case, though leading 232 to distortions elsewhere. See :mod:`sasmodels.compare` for many more details. 233 234 *sascomp -ngauss=n* allows you to set the number of quadrature points used 235 for the 1D integration for any model. For example, a carbon nanotube with 236 length 10 $\mu$\ m and radius 1 nm is not computed correctly at high $q$:: 237 238 python sascomp cylinder length=100000 radius=10 -ngauss=76,500 -double -highq 239 240 Note: ticket 702 gives some forms for long rods and thin disks that may 241 be used in place of the integration, depending on $q$, radius and length; if 242 the cylinder model is updated with these corrections then above call will show 243 no difference. 244 245 *explore/check1d.py* uses sasmodels 1D integration and compares that with a 201 :mod:`asymint` a direct implementation of the surface integral for certain 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 246 209 rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to 247 210 $\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle … … 251 214 similar reasoning.] This should rotate the sample through the entire 252 215 $\theta$-$\phi$ surface according to the pattern that you see in jitter.py when 253 you use 'rectangle' rather than 'gaussian' for its distribution without254 changing the viewing angle. In order to match the 1-D pattern for an arbitrary 255 viewing angle on triaxial shapes, we need to integrate216 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 256 219 over $\theta$, $\phi$ and $\Psi$. 257 220 258 *sascomp -sphere=n* uses the same rectangular distribution as check1d to 259 compute the pattern of the $q_x$-$q_y$ grid. This ought to produce a 260 spherically symmetric pattern. 261 262 *explore/precision.py* investigates the accuracy of individual functions 263 on the different execution platforms. This includes the basic special 264 functions as well as more complex expressions used in scattering. In many 265 cases the OpenCL function in sasmodels will use a polynomial approximation 266 over part of the range to improve accuracy, as shown in:: 267 268 python explore/precision.py 3j1/x 269 270 *explore/realspace.py* allows you to set up a Monte Carlo simulation of your 271 model by sampling random points within and computing the 1D and 2D scattering 272 patterns. This was used to check the core shell parallelepiped example. This 273 is similar to the general sas calculator in sasview, though it uses different 274 code. 275 276 *explore/jitter.py* is for exploring different options for handling 277 orientation and orientation dispersity. It uses *explore/guyou.py* to 278 generate the Guyou projection. 279 280 *explore/asymint.py* is a direct implementation of the 1D integration for 281 a number of different models. It calculates the integral for a particular 282 $q$ using several different integration schemes, including mpmath with 100 283 bits of precision (33 digits), so you can use it to check the target values 284 for the 1D model tests. This is not a general purpose tool; you will need to 285 edit the file to change the model parameters. It does not currently 286 apply the $u=cos(\theta)$ substitution that many models are using 287 internally so the 76-point gaussian quadrature may not match the value 288 produced by the eqivalent model from sasmodels. 289 290 *explore/symint.py* is for rotationally symmetric models (just cylinder for 291 now), so it can compute an entire curve rather than a single $q$ point. It 292 includes code to compare the long cylinder approximation to cylinder. 293 294 *explore/rpa.py* is for checking different implementations of the RPA model 295 against calculations for specific blends. This is a work in (suspended) 296 progress. 297 298 There are a few modules left over in *explore* that can probably be removed. 299 *angular_pd.py* was an early version of *jitter.py*. *sc.py* and *sc.c* 300 was used to explore different calculations for paracrystal models; it has 301 been absorbed into *asymint.py*. *transform_angles.py* translates orientation 302 parameters from the SasView 3.x definition to sasmodels. 303 304 *explore/angles.py* generates code for the view and jitter transformations. 305 Keep this around since it may be needed if we add new projections. 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 306 245 307 246 Testing -
doc/guide/orientation/orientation.rst
r5fb0634 r82592da 4 4 ================== 5 5 6 With two dimensional small angle diffraction data sasmodelswill calculate6 With two dimensional small angle diffraction data SasView will calculate 7 7 scattering from oriented particles, applicable for example to shear flow 8 8 or orientation in a magnetic field. 9 9 10 10 In general we first need to define the reference orientation 11 of the particle's $a$-$b$-$c$ axes with respect to the incoming 12 neutron or X-ray beam. This is done using three angles: $\theta$ and $\phi$ 13 define the orientation of the $c$-axis of the particle, and angle $\Psi$ is 14 defined as the orientation of the major axis of the particle cross section 15 with respect to its starting position along the beam direction (or 16 equivalently, as rotation about the $c$ axis). There is an unavoidable 17 ambiguity when $c$ is aligned with $z$ in that $\phi$ and $\Psi$ both 18 serve to rotate the particle about $c$, but this symmetry is destroyed 19 when $\theta$ is not a multiple of 180. 20 21 The figures below are for an elliptical cross section cylinder, but may 22 be applied analogously to other shapes of particle. 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. 23 18 24 19 .. note:: … … 34 29 35 30 Definition of angles for oriented elliptical cylinder, where axis_ratio 36 b/a is shown >1 .Note that rotation $\theta$, initially in the $x$-$z$31 b/a is shown >1, Note that rotation $\theta$, initially in the $x$-$z$ 37 32 plane, is carried out first, then rotation $\phi$ about the $z$-axis, 38 33 finally rotation $\Psi$ is around the axis of the cylinder. The neutron 39 or X-ray beam is along the $ -z$ axis.34 or X-ray beam is along the $z$ axis. 40 35 41 36 .. figure:: … … 45 40 with $\Psi$ = 0. 46 41 47 Having established the mean direction of the particle (the view) we can then 48 apply angular orientation distributions (jitter). This is done by a numerical 49 integration over a range of angles in a similar way to particle size 50 dispersity. The orientation dispersity is defined with respect to the 51 $a$-$b$-$c$ axes of the particle, with roll angle $\Psi$ about the $c$-axis, 52 yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the 53 $a$-axis. 54 55 More formally, starting with axes $a$-$b$-$c$ of the particle aligned 56 with axes $x$-$y$-$z$ of the laboratory frame, the orientation dispersity 57 is applied first, using the 58 `Tait-Bryan <https://en.wikipedia.org/wiki/Euler_angles#Conventions_2>`_ 59 $x$-$y'$-$z''$ convention with angles $\Delta\phi$-$\Delta\theta$-$\Delta\Psi$. 60 The reference orientation then follows, using the 61 `Euler angles <https://en.wikipedia.org/wiki/Euler_angles#Conventions>`_ 62 $z$-$y'$-$z''$ with angles $\phi$-$\theta$-$\Psi$. This is implemented 63 using rotation matrices as 64 65 .. math:: 66 67 R = R_z(\phi)\, R_y(\theta)\, R_z(\Psi)\, 68 R_x(\Delta\phi)\, R_y(\Delta\theta)\, R_z(\Delta\Psi) 69 70 To transform detector $(q_x, q_y)$ values into $(q_a, q_b, q_c)$ for the 71 shape in its canonical orientation, use 72 73 .. math:: 74 75 [q_a, q_b, q_c]^T = R^{-1} \, [q_x, q_y, 0]^T 76 77 78 The inverse rotation is easily calculated by rotating the opposite directions 79 in the reverse order, so 80 81 .. math:: 82 83 R^{-1} = R_z(-\Delta\Psi)\, R_y(-\Delta\theta)\, R_x(-\Delta\phi)\, 84 R_z(-\Psi)\, R_y(-\theta)\, R_z(-\phi) 85 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. 86 47 87 48 The $\theta$ and $\phi$ orientation parameters for the cylinder only appear 88 when fitting 2d data. On introducing "Orientational Distribution" in the89 angles, "distribution of theta" and "distribution of phi" parameters will49 when fitting 2d data. On introducing "Orientational Distribution" in 50 the angles, "distribution of theta" and "distribution of phi" parameters will 90 51 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ 91 of the cylinder, which correspond to the $b$ and $a$ axes of the cylinder 92 cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and 93 $X$ axes of the instrument.) The third orientation distribution, in $\Psi$, 94 is about the $c$ axis of the particle. Some experimentation may be required 95 to understand the 2d patterns fully. A number of different shapes of 96 distribution are available, as described for size dispersity, see 97 :ref:`polydispersityhelp`. 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` . 98 58 99 Given that the angular dispersion distribution is defined in cartesian space, 100 over a cube defined by 101 102 .. math:: 103 104 [-\Delta \theta, \Delta \theta] \times 105 [-\Delta \phi, \Delta \phi] \times 106 [-\Delta \Psi, \Delta \Psi] 107 108 but the orientation is defined over a sphere, we are left with a 109 `map projection <https://en.wikipedia.org/wiki/List_of_map_projections>`_ 110 problem, with different tradeoffs depending on how values in $\Delta\theta$ 111 and $\Delta\phi$ are translated into latitude/longitude on the sphere. 112 113 Sasmodels is using the 114 `equirectangular projection <https://en.wikipedia.org/wiki/Equirectangular_projection>`_. 115 In this projection, square patches in angular dispersity become wedge-shaped 116 patches on the sphere. To correct for the changing point density, there is a 117 scale factor of $\sin(\Delta\theta)$ that applies to each point in the 118 integral. This is not enough, though. Consider a shape which is tumbling 119 freely around the $b$ axis, with $\Delta\theta$ uniform in $[-180, 180]$. At 120 $\pm 90$, all points in $\Delta\phi$ map to the pole, so the jitter will have 121 a distinct angular preference. If the spin axis is along the beam (which 122 will be the case for $\theta=90$ and $\Psi=90$) the scattering pattern 123 should be circularly symmetric, but it will go to zero at $q_x = 0$ due to the 124 $\sin(\Delta\theta)$ correction. This problem does not appear for a shape 125 that is tumbling freely around the $a$ axis, with $\Delta\phi$ uniform in 126 $[-180, 180]$, so swap the $a$ and $b$ axes so $\Delta\theta < \Delta\phi$ 127 and adjust $\Psi$ by 90. This works with the current sasmodels shapes due to 128 symmetry. 129 130 Alternative projections were considered. 131 The `sinusoidal projection <https://en.wikipedia.org/wiki/Sinusoidal_projection>`_ 132 works by scaling $\Delta\phi$ as $\Delta\theta$ increases, and dropping those 133 points outside $[-180, 180]$. The distortions are a little less for middle 134 ranges of $\Delta\theta$, but they are still severe for large $\Delta\theta$ 135 and the model is much harder to explain. 136 The `azimuthal equidistance projection <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>`_ 137 also improves on the equirectangular projection by extending the range of 138 reasonable values for the $\Delta\theta$ range, with $\Delta\phi$ forming a 139 wedge that cuts to the opposite side of the sphere rather than cutting to the 140 pole. This projection has the nice property that distance from the center are 141 preserved, and that $\Delta\theta$ and $\Delta\phi$ act the same. 142 The `azimuthal equal area projection <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_ 143 is like the azimuthal equidistance projection, but it preserves area instead 144 of distance. It also has the same behaviour for $\Delta\theta$ and $\Delta\phi$. 145 The `Guyou projection <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>`_ 146 has an excellent balance with reasonable distortion in both $\Delta\theta$ 147 and $\Delta\phi$, as well as preserving small patches. However, it requires 148 considerably more computational overhead, and we have not yet derived the 149 formula for the distortion correction, measuring the degree of stretch at 150 the point $(\Delta\theta, \Delta\phi)$ on the map. 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. 151 63 152 64 .. note:: 153 Note that the form factors for oriented particles are performing 154 numerical integrations over one or more variables, so care should be 155 taken, especially with very large particles or more extreme aspect 156 ratios. In such cases results may not be accurate, particularly at very 157 high Q, unless the model has been specifically coded to use limiting 158 forms of the scattering equations. 65 Note that the form factors for oriented particles are also performing 66 numerical integrations over one or more variables, so care should be taken, 67 especially with very large particles or more extreme aspect ratios. In such 68 cases results may not be accurate, particularly at very high Q, unless the model 69 has been specifically coded to use limiting forms of the scattering equations. 70 71 For best numerical results keep the $\theta$ distribution narrower than the $\phi$ 72 distribution. Thus for asymmetric particles, such as elliptical_cylinder, you may 73 need to reorder the sizes of the three axes to acheive the desired result. 74 This is due to the issues of mapping a rectangular distribution onto the 75 surface of a sphere. 159 76 160 For best numerical results keep the $\theta$ distribution narrower than 161 the $\phi$ distribution. Thus for asymmetric particles, such as 162 elliptical_cylinder, you may need to reorder the sizes of the three axes 163 to acheive the desired result. This is due to the issues of mapping a 164 rectanglar distribution onto the surface of a sphere. 165 166 Users can experiment with the values of *Npts* and *Nsigs*, the number of steps 167 used in the integration and the range spanned in number of standard deviations. 168 The standard deviation is entered in units of degrees. For a "rectangular" 169 distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations. 170 The new "uniform" distribution avoids this by letting you directly specify the 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 171 82 half width. 172 83 173 The angular distributions may be truncated outside of the range -180 to +180174 degrees, so beware of using saying a broad Gaussian distribution with large 175 value of *Nsigs*, as the array of *Npts* may be truncated to many fewer 176 points than would give a good integration,as well as becoming rather 177 meaningless. (At some point in the future the actual dispersion arrays may be 178 made available to the userfor inspection.)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.) 179 90 180 91 Some more detailed technical notes are provided in the developer section of 181 92 this manual :ref:`orientation_developer` . 182 93 183 This definition of orientation is new to SasView 4.2. In earlier versions,184 the orientation distribution appeared as a distribution of view angles.185 This led to strange effects when $c$ was aligned with $z$, where changes186 to the $\phi$ angle served only to rotate the shape about $c$, rather than187 having a consistent interpretation as the pitch of the shape relative to188 the flow field defining the reference orientation. Prior to SasView 4.1,189 the reference orientation was defined using a Tait-Bryan convention, making190 it difficult to control. Now, rotation in $\theta$ modifies the spacings191 in the refraction pattern, and rotation in $\phi$ rotates it in the detector192 plane.193 194 195 94 *Document History* 196 95 197 96 | 2017-11-06 Richard Heenan 198 | 2017-12-20 Paul Kienzle -
doc/guide/plugin.rst
r7e6bc45e rd0dc9a3 292 292 **Note: The order of the parameters in the definition will be the order of the 293 293 parameters in the user interface and the order of the parameters in Iq(), 294 Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background* 295 **parameters are implicit to all models, so they do not need to be included 296 in the parameter table.** 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.** 297 296 298 297 - **"name"** is the name of the parameter shown on the FitPage. … … 363 362 scattered intensity. 364 363 365 - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(), 366 and have polydispersity loops generated automatically. 367 368 - "orientation" parameters are not passed, but instead are combined with 369 orientation dispersity to translate *qx* and *qy* to *qa*, *qb* and *qc*. 370 These parameters should appear at the end of the table with the specific 371 names *theta*, *phi* and for asymmetric shapes *psi*, in that order. 364 - "volume" parameters are passed to Iq(), Iqxy(), and form_volume(), and 365 have polydispersity loops generated automatically. 366 367 - "orientation" parameters are only passed to Iqxy(), and have angular 368 dispersion. 372 369 373 370 Some models will have integer parameters, such as number of pearls in the … … 422 419 That is, the individual models do not need to include polydispersity 423 420 calculations, but instead rely on numerical integration to compute the 424 appropriately smeared pattern. 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. 425 424 426 425 Python Models … … 469 468 barbell). If I(q; pars) is NaN for any $q$, then those parameters will be 470 469 ignored, 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, ...)*. 471 474 472 475 Models should define *form_volume(par1, par2, ...)* where the parameter … … 494 497 } 495 498 499 *Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*, 500 and it includes orientation parameters. 501 496 502 *form_volume* defines the volume of the shape. As in python models, it 497 503 includes only the volume parameters. 498 504 505 *Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and 506 *form_volume* will default to 1.0. 507 499 508 **source=['fn.c', ...]** includes the listed C source files in the 500 program before *Iq* and *form_volume* are defined. This allows you to 501 extend the library of C functions available to your model. 502 503 *c_code* includes arbitrary C code into your kernel, which can be 504 handy for defining helper functions for *Iq* and *form_volume*. Note that 505 you can put the full function definition for *Iq* and *form_volume* 506 (include function declaration) into *c_code* as well, or put them into an 507 external C file and add that file to the list of sources. 509 program before *Iq* and *Iqxy* are defined. This allows you to extend the 510 library of C functions available to your model. 508 511 509 512 Models are defined using double precision declarations for the … … 529 532 530 533 #define INVALID(v) (v.bell_radius < v.radius) 531 532 The INVALID define can go into *Iq*, or *c_code*, or an external C file533 listed in *source*.534 535 Oriented Shapes536 ...............537 538 If the scattering is dependent on the orientation of the shape, then you539 will need to include *orientation* parameters *theta*, *phi* and *psi*540 at the end of the parameter table. As described in the section541 :ref:`orientation`, the individual $(q_x, q_y)$ points on the detector will542 be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its543 canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the544 laboratory frame and beam travelling along $-z$.545 546 The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where547 *par1*, etc. are the parameters to the model. If the shape is rotationally548 symmetric about *c* then *psi* is not needed, and the model is called549 as *Iqac(qab, qc, par1, par2, ...)*. In either case, the orientation550 parameters are not included in the function call.551 552 For 1D oriented shapes, an integral over all angles is usually needed for553 the *Iq* function. Given symmetry and the substitution $u = \cos(\alpha)$,554 $du = -\sin(\alpha)\,d\alpha$ this becomes555 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\,du564 565 for566 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 u572 573 Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the574 numerical integration is then::575 576 double outer_sum = 0.0;577 for (int i = 0; i < GAUSS_N; i++) {578 const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5;579 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);580 const double qc = cos_alpha * q;581 double inner_sum = 0.0;582 for (int j = 0; j < GAUSS_N; j++) {583 const double beta = M_PI_4 * GAUSS_Z[j] + M_PI_4;584 double sin_beta, cos_beta;585 SINCOS(beta, sin_beta, cos_beta);586 const double qa = sin_alpha * sin_beta * q;587 const double qb = sin_alpha * cos_beta * q;588 const double form = Fq(qa, qb, qc, ...);589 inner_sum += GAUSS_W[j] * form * form;590 }591 outer_sum += GAUSS_W[i] * inner_sum;592 }593 outer_sum *= 0.25; // = 8/(4 pi) * outer_sum * (pi/2) / 4594 595 The *z* values for the Gauss-Legendre integration extends from -1 to 1, so596 the double sum of *w[i]w[j]* explains the factor of 4. Correcting for the597 average *dz[i]dz[j]* gives $(1-0) \cdot (\pi/2-0) = \pi/2$. The $8/(4 \pi)$598 factor comes from the integral over the quadrant. With less symmetry (eg.,599 in the bcc and fcc paracrystal models), then an integral over the entire600 sphere may be necessary.601 602 For simpler models which are rotationally symmetric a single integral603 suffices:604 605 .. math::606 607 I(q) &= \frac{1}{\pi}\int_{-\pi/2}^{\pi/2}608 F(q_{ab}, q_c)^2 \sin(\alpha)\,d\alpha/\pi \\609 &= \frac{2}{\pi} \int_0^1 F^2\,du610 611 for612 613 .. math::614 615 q_{ab} &= q \sin(\alpha) = q \sqrt{1 - u^2} \\616 q_c &= q \cos(\alpha) = q u617 618 619 with integration loop::620 621 double sum = 0.0;622 for (int i = 0; i < GAUSS_N; i++) {623 const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5;624 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);625 const double qab = sin_alpha * q;626 const double qc = cos_alpha * q;627 const double form = Fq(qab, qc, ...);628 sum += GAUSS_W[j] * form * form;629 }630 sum *= 0.5; // = 2/pi * sum * (pi/2) / 2631 632 Magnetism633 .........634 635 Magnetism is supported automatically for all shapes by modifying the636 effective SLD of particle according to the Halpern-Johnson vector637 describing the interaction between neutron spin and magnetic field. All638 parameters marked as type *sld* in the parameter table are treated as639 possibly magnetic particles with magnitude *M0* and direction640 *mtheta* and *mphi*. Polarization parameters are also provided641 automatically for magnetic models to set the spin state of the measurement.642 643 For more complicated systems where magnetism is not uniform throughout644 the individual particles, you will need to write your own models.645 You should not mark the nuclear sld as type *sld*, but instead leave646 them unmarked and provide your own magnetism and polarization parameters.647 For 2D measurements you will need $(q_x, q_y)$ values for the measurement648 to compute the proper magnetism and orientation, which you can implement649 using *Iqxy(qx, qy, par1, par2, ...)*.650 534 651 535 Special Functions … … 912 796 show a 50x improvement or more over the equivalent pure python model. 913 797 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 914 807 915 808 .. _Form_Factors: … … 1113 1006 variable name *Rg* for example because $R_g$ is the right name for the model 1114 1007 parameter then ignore the lint errors. Also, ignore *missing-docstring* 1115 for standard model functions *Iq*, *Iq ac*, etc.1008 for standard model functions *Iq*, *Iqxy*, etc. 1116 1009 1117 1010 We will have delinting sessions at the SasView Code Camps, where we can -
explore/jitter.py
r8cfb486 rff10479 165 165 # constants in kernel_iq.c 166 166 'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance', 167 'azimuthal_equal_area',168 167 ] 169 168 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', -
explore/precision.py
r2a7e20e r2a602c7 345 345 ) 346 346 add_function( 347 name="(1/2 -sin(x)/x+(1-cos(x))/x^2)/x",347 name="(1/2+(1-cos(x))/x^2-sin(x)/x)/x", 348 348 mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x, 349 349 np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x, … … 609 609 names = ", ".join(sorted(ALL_FUNCTIONS)) 610 610 print("""\ 611 usage: precision.py [-f/a/r] [-x<range>] "name"...611 usage: precision.py [-f/a/r] [-x<range>] name... 612 612 where 613 613 -f indicates that the function value should be plotted, … … 620 620 zoom indicates linear stepping in [1000, 1010] 621 621 neg indicates linear stepping in [-100.1, 100.1] 622 and name is "all " or one of:622 and name is "all [first]" or one of: 623 623 """+names) 624 624 sys.exit(1) -
explore/realspace.py
r8fb2a94 ra1c32c2 44 44 """ 45 45 return Rx(phi)*Ry(theta)*Rz(psi) 46 47 def apply_view(points, view):48 """49 Apply the view transform (theta, phi, psi) to a set of points.50 51 Points are stored in a 3 x n numpy array.52 53 View angles are in degrees.54 """55 theta, phi, psi = view56 return np.asarray((Rz(phi)*Ry(theta)*Rz(psi))*np.matrix(points.T)).T57 58 59 def invert_view(qx, qy, view):60 """61 Return (qa, qb, qc) for the (theta, phi, psi) view angle at detector62 pixel (qx, qy).63 64 View angles are in degrees.65 """66 theta, phi, psi = view67 q = np.vstack((qx.flatten(), qy.flatten(), 0*qx.flatten()))68 return np.asarray((Rz(-psi)*Ry(-theta)*Rz(-phi))*np.matrix(q))69 70 46 71 47 class Shape: … … 243 219 values = self.value.repeat(points.shape[0]) 244 220 return values, self._adjust(points) 245 246 NUMBA = False247 if NUMBA:248 from numba import njit249 @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. + 0j254 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)**2257 return Iq258 259 def calc_Iqxy(qx, qy, rho, points, volume=1.0, view=(0, 0, 0)):260 qx, qy = np.broadcast_arrays(qx, qy)261 qa, qb, qc = invert_view(qx, qy, view)262 rho, volume = np.broadcast_arrays(rho, volume)263 values = rho*volume264 x, y, z = points.T265 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))))**2271 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)273 221 274 222 def _calc_Pr_nonuniform(r, rho, points): … … 291 239 print("processing %d of %d"%(k, len(rho)-1)) 292 240 Pr = extended_Pr[1:-1] 293 return Pr 294 295 def _calc_Pr_uniform(r, rho, points): 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 296 249 # Make Pr a little be bigger than necessary so that only distances 297 250 # min < d < max end up in Pr 298 dr, n_max = r[0],len(r)251 n_max = len(r) 299 252 extended_Pr = np.zeros(n_max+1, 'd') 300 253 t0 = time.clock() 301 254 t_next = t0 + 3 255 row_zero = np.zeros(len(rho), 'i') 302 256 for k, rho_k in enumerate(rho[:-1]): 303 257 distances = sqrt(np.sum((points[k] - points[k+1:])**2, axis=1)) 304 258 weights = rho_k * rho[k+1:] 305 index = np.minimum(np.asarray(distances/ dr, 'i'), n_max)259 index = np.minimum(np.asarray(distances/r[0], 'i'), n_max) 306 260 # Note: indices may be duplicated, so "Pr[index] += w" will not work!! 307 261 extended_Pr += np.bincount(index, weights, n_max+1) … … 315 269 # we are only accumulating the upper triangular distances. 316 270 #Pr = Pr * 2 / len(rho)**2 317 return Pr 271 return Pr / Pr.max() 318 272 319 273 # Can get an additional 2x by going to C. Cuda/OpenCL will allow even … … 337 291 } 338 292 """ 339 340 if NUMBA:341 @njit("f8[:](f8[:], f8[:], f8[:,:])")342 def _calc_Pr_uniform_jit(r, rho, points):343 dr = r[0]344 n_max = len(r)345 Pr = np.zeros_like(r)346 for j in range(len(rho) - 1):347 x, y, z = points[j, 0], points[j, 1], points[j, 2]348 for k in range(j+1, len(rho)):349 distance = sqrt((x - points[k, 0])**2350 + (y - points[k, 1])**2351 + (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 because356 # we are only accumulating the upper triangular distances.357 #Pr = Pr * 2 / len(rho)**2358 return Pr359 360 361 def calc_Pr(r, rho, points):362 # P(r) with uniform steps in r is 3x faster; check if we are uniform363 # before continuing364 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 373 293 374 294 def j0(x): … … 413 333 plt.legend() 414 334 415 def plot_calc_2d(qx, qy, Iqxy, theory=None):416 import matplotlib.pyplot as plt417 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 429 335 def plot_points(rho, points): 430 336 import mpl_toolkits.mplot3d … … 437 343 pass 438 344 n = len(points) 439 #print("len points", n) 440 index = np.random.choice(n, size=500) if n > 500 else slice(None, None) 345 index = np.random.choice(n, size=1000) if n > 1000 else slice(None, None) 441 346 ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index]) 442 347 #low, high = points.min(axis=0), points.max(axis=0) 443 348 #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]]) 444 349 ax.autoscale(True) 445 446 def check_shape(shape, fn=None):447 rho_solvent = 0448 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 None457 458 import pylab459 #plot_points(rho, points); pylab.figure()460 plot_calc(r, Pr, q, Iq, theory=theory)461 pylab.show()462 463 def check_shape_2d(shape, fn=None, view=(0, 0, 0)):464 rho_solvent = 0465 nq, qmax = 100, 1.0466 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 None477 Iqxy += 0.001 * Iqxy.max()478 if theory is not None:479 theory += 0.001 * theory.max()480 481 import pylab482 #plot_points(rho, points); pylab.figure()483 plot_calc_2d(qx, qy, Iqxy, theory=theory)484 pylab.show()485 486 def sas_sinx_x(x):487 with np.errstate(all='ignore'):488 retvalue = sin(x)/x489 retvalue[x == 0.] = 1.490 return retvalue491 350 492 351 def sas_2J1x_x(x): … … 514 373 Iq = Iq/Iq[0] 515 374 return Iq 516 517 def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)):518 qa, qb, qc = invert_view(qx, qy, view)519 qab = np.sqrt(qa**2 + qb**2)520 Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2)521 Iq = Fq**2522 return Iq.reshape(qx.shape)523 375 524 376 def sphere_Iq(q, radius): … … 563 415 return Iq/Iq[0] 564 416 565 def csbox_Iqxy(qx, qy, a, b, c, da, db, dc, slda, sldb, sldc, sld_core, view=(0,0,0)): 566 qa, qb, qc = invert_view(qx, qy, view) 567 568 sld_solvent = 0 569 overlapping = False 570 dr0 = sld_core - sld_solvent 571 drA, drB, drC = slda-sld_solvent, sldb-sld_solvent, sldc-sld_solvent 572 tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc 573 siA = a*sas_sinx_x(a*qa/2) 574 siB = b*sas_sinx_x(b*qb/2) 575 siC = c*sas_sinx_x(c*qc/2) 576 siAt = tA*sas_sinx_x(tA*qa/2) 577 siBt = tB*sas_sinx_x(tB*qb/2) 578 siCt = tC*sas_sinx_x(tC*qc/2) 579 Fq = (dr0*siA*siB*siC 580 + drA*(siAt-siA)*siB*siC 581 + drB*siA*(siBt-siB)*siC 582 + drC*siA*siB*(siCt-siC)) 583 Iq = Fq**2 584 return Iq.reshape(qx.shape) 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() 585 431 586 432 def check_cylinder(radius=25, length=125, rho=2.): … … 588 434 fn = lambda q: cylinder_Iq(q, radius, length) 589 435 check_shape(shape, fn) 590 591 def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)):592 shape = EllipticalCylinder(radius, radius, length, rho)593 fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view)594 check_shape_2d(shape, fn, view=view)595 596 def check_cylinder_2d_lattice(radius=25, length=125, rho=2.,597 view=(0, 0, 0)):598 nx, dx = 1, 2*radius599 ny, dy = 30, 2*radius600 nz, dz = 30, length601 dx, dy, dz = 2*dx, 2*dy, 2*dz602 def center(*args):603 sigma = 0.333604 space = 2605 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)617 436 618 437 def check_sphere(radius=125, rho=2): … … 630 449 side_c2 = copy(side_c).shift(0, 0, -c-dc) 631 450 shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 632 def fn(q): 633 return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 634 #check_shape(shape, fn) 635 636 view = (20, 30, 40) 637 def fn_xy(qx, qy): 638 return csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 639 slda, sldb, sldc, sld_core, view=view) 640 check_shape_2d(shape, fn_xy, view=view) 451 fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 452 check_shape(shape, fn) 641 453 642 454 if __name__ == "__main__": 643 455 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))646 456 #check_sphere() 647 457 #check_csbox() -
sasmodels/compare.py
r2a7e20e r2d81cfe 92 92 -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh 93 93 -neval=1 sets the number of evals for more accurate timing 94 -ngauss=0 overrides the number of points in the 1-D gaussian quadrature95 94 96 95 === precision options === -
sasmodels/details.py
r108e70e r2d81cfe 258 258 # type: (...) -> Sequence[np.ndarray] 259 259 """ 260 **Deprecated** Theta weights will be computed in the kernel wrapper if261 they are needed.262 263 260 If there is a theta parameter, update the weights of that parameter so that 264 261 the cosine weighting required for polar integration is preserved. … … 275 272 Returns updated weights vectors 276 273 """ 274 # TODO: explain in a comment why scale and background are missing 277 275 # Apparently the parameters.theta_offset similarly skips scale and 278 276 # and background, so the indexing works out, but they are still shipped … … 281 279 index = parameters.theta_offset 282 280 theta = dispersity[index] 281 # TODO: modify the dispersity vector to avoid the theta=-90,90,270,... 283 282 theta_weight = abs(cos(radians(theta))) 284 283 weights = tuple(theta_weight*w if k == index else w -
sasmodels/generate.py
r108e70e r2d81cfe 7 7 particular dimensions averaged over all orientations. 8 8 9 *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc 10 for a rotationally symmetric form with particular dimensions. 11 qab, qc are determined from shape orientation and scattering angles. 12 This call is used if the shape has orientation parameters theta and phi. 13 14 *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc 15 for a form with particular dimensions. qa, qb, qc are determined from 16 shape orientation and scattering angles. This call is used if the shape 17 has orientation parameters theta, phi and psi. 18 19 *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy. Use this 20 to create an arbitrary 2D theory function, needed for q-dependent 21 background functions and for models with non-uniform magnetism. 9 *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form 10 with particular dimensions for a single orientation. 11 12 *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the 13 polarized neutron spin states (up-up, up-down, down-up, down-down) for 14 a form with particular dimensions for a single orientation. 22 15 23 16 *form_volume(p1, p2, ...)* returns the volume of the form with particular … … 38 31 scale and background parameters for each model. 39 32 40 C code should be stylized C-99 functions written for OpenCL. All functions 41 need prototype declarations even if the are defined before they are used. 42 Although OpenCL supports *#include* preprocessor directives, the list of 43 includes should be given as part of the metadata in the kernel module 44 definition. The included files should be listed using a path relative to the 45 kernel module, or if using "lib/file.c" if it is one of the standard includes 46 provided with the sasmodels source. The includes need to be listed in order 47 so that functions are defined before they are used. 33 *Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 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. 48 42 49 43 Floating point values should be declared as *double*. For single precision … … 113 107 present, the volume ratio is 1. 114 108 115 *form_volume*, *Iq*, *Iq ac*, *Iqabc* are strings containing116 the C source code for the body of the volume, Iq, and Iqacfunctions109 *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the 110 C source code for the body of the volume, Iq, and Iqxy functions 117 111 respectively. These can also be defined in the last source file. 118 112 119 *Iq* , *Iqac*, *Iqabc* also be instead be python functions defining the113 *Iq* and *Iqxy* also be instead be python functions defining the 120 114 kernel. If they are marked as *Iq.vectorized = True* then the 121 115 kernel is passed the entire *q* vector at once, otherwise it is … … 174 168 from zlib import crc32 175 169 from inspect import currentframe, getframeinfo 176 import logging177 170 178 171 import numpy as np # type: ignore … … 188 181 pass 189 182 # pylint: enable=unused-import 190 191 logger = logging.getLogger(__name__)192 183 193 184 # jitter projection to use in the kernel code. See explore/jitter.py … … 635 626 636 627 """ 637 def _gen_fn( model_info, name, pars):638 # type: ( ModelInfo, str, List[Parameter]) -> str628 def _gen_fn(name, pars, body, filename, line): 629 # type: (str, List[Parameter], str, str, int) -> str 639 630 """ 640 631 Generate a function given pars and body. … … 648 639 """ 649 640 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.filename652 # Note: if symbol is defined strangely in the module then default it to 1653 lineno = model_info.lineno.get(name, 1)654 641 return _FN_TEMPLATE % { 655 642 'name': name, 'pars': par_decl, 'body': body, 656 'filename': filename.replace('\\', '\\\\'), 'line': line no,643 'filename': filename.replace('\\', '\\\\'), 'line': line, 657 644 } 658 645 … … 669 656 670 657 # type in IQXY pattern could be single, float, double, long double, ... 671 _IQXY_PATTERN = re.compile( r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]",658 _IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)", 672 659 flags=re.MULTILINE) 673 def find_xy_mode(source):660 def _have_Iqxy(sources): 674 661 # type: (List[str]) -> bool 675 662 """ 676 Return t he xy mode as qa, qac, qabc orqxy.663 Return true if any file defines Iqxy. 677 664 678 665 Note this is not a C parser, and so can be easily confused by 679 666 non-standard syntax. Also, it will incorrectly identify the following 680 as having 2D models::667 as having Iqxy:: 681 668 682 669 /* 683 double Iq ac(qab, qc, ...) { ... fill this in later ... }670 double Iqxy(qx, qy, ...) { ... fill this in later ... } 684 671 */ 685 672 686 If you want to comment out the function, use // on the front of the 687 line:: 688 689 /* 690 // double Iqac(qab, qc, ...) { ... fill this in later ... } 691 */ 692 693 """ 694 for code in source: 695 m = _IQXY_PATTERN.search(code) 696 if m is not None: 697 return m.group('mode') 698 return 'qa' 699 700 701 def _add_source(source, code, path, lineno=1): 673 If you want to comment out an Iqxy function, use // on the front of the 674 line instead. 675 """ 676 for _path, code in sources: 677 if _IQXY_PATTERN.search(code): 678 return True 679 return False 680 681 682 def _add_source(source, code, path): 702 683 """ 703 684 Add a file to the list of source code chunks, tagged with path and line. 704 685 """ 705 686 path = path.replace('\\', '\\\\') 706 source.append('#line %d "%s"' % (lineno, path))687 source.append('#line 1 "%s"' % path) 707 688 source.append(code) 708 689 … … 735 716 user_code = [(f, open(f).read()) for f in model_sources(model_info)] 736 717 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 737 723 # Build initial sources 738 724 source = [] … … 741 727 _add_source(source, code, path) 742 728 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 747 729 # Make parameters for q, qx, qy so that we can use them in declarations 748 q, qx, qy, qab, qa, qb, qc \ 749 = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 730 q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 750 731 # Generate form_volume function, etc. from body only 751 732 if isinstance(model_info.form_volume, str): 752 733 pars = partable.form_volume_parameters 753 source.append(_gen_fn(model_info, 'form_volume', pars)) 734 source.append(_gen_fn('form_volume', pars, model_info.form_volume, 735 model_info.filename, model_info._form_volume_line)) 754 736 if isinstance(model_info.Iq, str): 755 737 pars = [q] + partable.iq_parameters 756 source.append(_gen_fn(model_info, 'Iq', pars)) 738 source.append(_gen_fn('Iq', pars, model_info.Iq, 739 model_info.filename, model_info._Iq_line)) 757 740 if isinstance(model_info.Iqxy, str): 758 pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters 759 source.append(_gen_fn(model_info, 'Iqxy', pars)) 760 if isinstance(model_info.Iqac, str): 761 pars = [qab, qc] + partable.iq_parameters 762 source.append(_gen_fn(model_info, 'Iqac', pars)) 763 if isinstance(model_info.Iqabc, str): 764 pars = [qa, qb, qc] + partable.iq_parameters 765 source.append(_gen_fn(model_info, 'Iqabc', pars)) 766 767 # What kind of 2D model do we need? Is it consistent with the parameters? 768 xy_mode = find_xy_mode(source) 769 if xy_mode == 'qabc' and not partable.is_asymmetric: 770 raise ValueError("asymmetric oriented models need to define Iqabc") 771 elif xy_mode == 'qac' and partable.is_asymmetric: 772 raise ValueError("symmetric oriented models need to define Iqac") 773 elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'): 774 raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode) 775 elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'): 776 if xy_mode == 'qxy': 777 logger.warn("oriented shapes should define Iqac or Iqabc") 778 else: 779 raise ValueError("Expected function Iqac or Iqabc for oriented shape") 741 pars = [qx, qy] + partable.iqxy_parameters 742 source.append(_gen_fn('Iqxy', pars, model_info.Iqxy, 743 model_info.filename, model_info._Iqxy_line)) 780 744 781 745 # Define the parameter table … … 803 767 if xy_mode == 'qabc': 804 768 pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 805 call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iq abc(%s)" % pars769 call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 806 770 clear_iqxy = "#undef CALL_IQ_ABC" 807 771 elif xy_mode == 'qac': 808 772 pars = ",".join(["_qa", "_qc"] + model_refs) 809 call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iq ac(%s)" % pars773 call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 810 774 clear_iqxy = "#undef CALL_IQ_AC" 811 el if xy_mode == 'qa':775 else: # xy_mode == 'qa' 812 776 pars = ",".join(["_qa"] + model_refs) 813 777 call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 814 778 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)" % pars819 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 827 779 828 780 magpars = [k-2 for k, p in enumerate(partable.call_parameters) -
sasmodels/kernel_header.c
r108e70e r8698a0d 150 150 inline double cube(double x) { return x*x*x; } 151 151 inline 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 angles154 155 // To rotate from the canonical position to theta, phi, psi, first rotate by156 // psi about the major axis, oriented along z, which is a rotation in the157 // detector plane xy. Next rotate by theta about the y axis, aligning the major158 // 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 -psi161 // The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit162 // vector in the q direction.163 // To change between counterclockwise and clockwise rotation, change the164 // sign of phi and psi.165 166 #if 1167 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017168 #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 #else175 // SasView 3.x definition of orientation176 #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 #endif183 184 #if 1185 #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 #else203 // SasView 3.x definition of orientation204 #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 31 31 // CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 32 32 // 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 models34 33 // INVALID(table) : test if the current point is feesible to calculate. This 35 34 // will be defined in the kernel definition file. … … 470 469 #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 471 470 #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 472 #elif defined(CALL_IQ_XY) 473 // direct call to qx,qy calculator 474 double qx, qy; 475 #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 476 #define BUILD_ROTATION() do {} while(0) 477 #define APPLY_ROTATION() do {} while(0) 478 #define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table) 479 #endif 480 481 // Define APPLY_PROJECTION depending on model symmetries. We do this outside 482 // the previous if block so that we don't need to repeat the identical 483 // logic in the IQ_AC and IQ_ABC branches. This will become more important 484 // if we implement more projections, or more complicated projections. 485 #if defined(CALL_IQ) || defined(CALL_IQ_A) // no orientation 471 #endif 472 473 // Doing jitter projection code outside the previous if block so that we don't 474 // need to repeat the identical logic in the IQ_AC and IQ_ABC branches. This 475 // will become more important if we implement more projections, or more 476 // complicated projections. 477 #if defined(CALL_IQ) || defined(CALL_IQ_A) 486 478 #define APPLY_PROJECTION() const double weight=weight0 487 #elif defined(CALL_IQ_XY) // pass orientation to the model 488 // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc 489 // Need to plug the values for the orientation angles back into parameter 490 // table in case they were overridden by the orientation offset. This 491 // means that orientation dispersity will not work for these models, but 492 // it was broken anyway, so no matter. Still want to provide Iqxy in case 493 // the user model wants full control of orientation/magnetism. 494 #if defined(HAVE_PSI) 495 const double theta = values[details->theta_par+2]; 496 const double phi = values[details->theta_par+3]; 497 const double psi = values[details->theta_par+4]; 498 double weight; 499 #define APPLY_PROJECTION() do { \ 500 local_values.table.theta = theta; \ 501 local_values.table.phi = phi; \ 502 local_values.table.psi = psi; \ 503 weight=weight0; \ 504 } while (0) 505 #elif defined(HAVE_THETA) 506 const double theta = values[details->theta_par+2]; 507 const double phi = values[details->theta_par+3]; 508 double weight; 509 #define APPLY_PROJECTION() do { \ 510 local_values.table.theta = theta; \ 511 local_values.table.phi = phi; \ 512 weight=weight0; \ 513 } while (0) 514 #else 515 #define APPLY_PROJECTION() const double weight=weight0 516 #endif 517 #else // apply jitter and view before calling the model 479 #else // !spherosymmetric projection 518 480 // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 519 481 const double theta = values[details->theta_par+2]; … … 526 488 // we go through the mesh. 527 489 double dtheta, dphi, weight; 528 #if PROJECTION == 1 // equirectangular490 #if PROJECTION == 1 529 491 #define APPLY_PROJECTION() do { \ 530 492 dtheta = local_values.table.theta; \ … … 532 494 weight = fabs(cos(dtheta*M_PI_180)) * weight0; \ 533 495 } while (0) 534 #elif PROJECTION == 2 // sinusoidal496 #elif PROJECTION == 2 535 497 #define APPLY_PROJECTION() do { \ 536 498 dtheta = local_values.table.theta; \ … … 542 504 } while (0) 543 505 #endif 544 #endif // done defining APPLY_PROJECTION506 #endif // !spherosymmetric projection 545 507 546 508 // ** define looping macros ** -
sasmodels/kernelpy.py
r108e70e r2d81cfe 26 26 # pylint: enable=unused-import 27 27 28 logger = logging.getLogger(__name__)29 30 28 class PyModel(KernelModel): 31 29 """ … … 33 31 """ 34 32 def __init__(self, model_info): 35 # Make sure Iq isavailable and vectorized33 # Make sure Iq and Iqxy are available and vectorized 36 34 _create_default_functions(model_info) 37 35 self.info = model_info 38 36 self.dtype = np.dtype('d') 39 logg er.info("load python model " + self.info.name)37 logging.info("load python model " + self.info.name) 40 38 41 39 def make_kernel(self, q_vectors): 42 40 q_input = PyInput(q_vectors, dtype=F64) 43 return PyKernel(self.info, q_input) 41 kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 42 return PyKernel(kernel, self.info, q_input) 44 43 45 44 def release(self): … … 90 89 Callable SAS kernel. 91 90 92 *kernel* is the kernel object to call.91 *kernel* is the DllKernel object to call. 93 92 94 93 *model_info* is the module information … … 105 104 Call :meth:`release` when done with the kernel instance. 106 105 """ 107 def __init__(self, model_info, q_input):106 def __init__(self, kernel, model_info, q_input): 108 107 # type: (callable, ModelInfo, List[np.ndarray]) -> None 109 108 self.dtype = np.dtype('d') … … 111 110 self.q_input = q_input 112 111 self.res = np.empty(q_input.nq, q_input.dtype) 112 self.kernel = kernel 113 113 self.dim = '2d' if q_input.is_2d else '1d' 114 114 … … 159 159 # Generate a closure which calls the form_volume if it exists. 160 160 form_volume = model_info.form_volume 161 self._volume = ((lambda: form_volume(*volume_args)) if form_volume else162 (lambda: 1.0))161 self._volume = ((lambda: form_volume(*volume_args)) if form_volume 162 else (lambda: 1.0)) 163 163 164 164 def __call__(self, call_details, values, cutoff, magnetic): … … 261 261 any functions that are not already marked as vectorized. 262 262 """ 263 # Note: must call create_vector_Iq before create_vector_Iqxy264 263 _create_vector_Iq(model_info) 265 _create_vector_Iqxy(model_info) 264 _create_vector_Iqxy(model_info) # call create_vector_Iq() first 266 265 267 266 … … 281 280 model_info.Iq = vector_Iq 282 281 283 284 282 def _create_vector_Iqxy(model_info): 285 283 """ 286 284 Define Iqxy as a vector function if it exists, or default it from Iq(). 287 285 """ 288 Iq xy = getattr(model_info, 'Iqxy', None)286 Iq, Iqxy = model_info.Iq, model_info.Iqxy 289 287 if callable(Iqxy): 290 288 if not getattr(Iqxy, 'vectorized', False): … … 297 295 vector_Iqxy.vectorized = True 298 296 model_info.Iqxy = vector_Iqxy 299 el se:297 elif callable(Iq): 300 298 #print("defaulting Iqxy") 301 299 # Iq is vectorized because create_vector_Iq was already called. 302 Iq = model_info.Iq303 300 def default_Iqxy(qx, qy, *args): 304 301 """ -
sasmodels/modelinfo.py
r108e70e r2d81cfe 37 37 38 38 # assumptions about common parameters exist throughout the code, such as: 39 # (1) kernel functions Iq, Iqxy, Iqac, Iqabc,form_volume, ... don't see them39 # (1) kernel functions Iq, Iqxy, form_volume, ... don't see them 40 40 # (2) kernel drivers assume scale is par[0] and background is par[1] 41 41 # (3) mixture models drop the background on components and replace the scale … … 256 256 257 257 *type* indicates how the parameter will be used. "volume" parameters 258 will be used in all functions. "orientation" parameters are not passed,259 but will be used to convert from *qx*, *qy* to *qa*, *qb*, *qc* in calls to260 *I qxy* and *Imagnetic*. If *type* is the empty string, the parameter will258 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 261 261 be used in all of *Iq*, *Iqxy* and *Imagnetic*. "sld" parameters 262 262 can automatically be promoted to magnetic parameters, each of which … … 386 386 with vector parameter p sent as p[]. 387 387 388 * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 389 function, with vector parameter p sent as p[]. 390 388 391 * *form_volume_parameters* is the list of parameters to the form_volume(...) 389 392 function, with vector parameter p sent as p[]. … … 440 443 self.iq_parameters = [p for p in self.kernel_parameters 441 444 if p.type not in ('orientation', 'magnetic')] 442 self.orientation_parameters = [p for p in self.kernel_parameters 443 if p.type == 'orientation'] 445 # note: orientation no longer sent to Iqxy, so its the same as 446 #self.iqxy_parameters = [p for p in self.kernel_parameters 447 # if p.type != 'magnetic'] 444 448 self.form_volume_parameters = [p for p in self.kernel_parameters 445 449 if p.type == 'volume'] … … 486 490 if p.type != 'orientation': 487 491 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")490 492 if theta >= 0 and phi >= 0: 491 last_par = len(self.kernel_parameters) - 1492 493 if phi != theta+1: 493 494 raise TypeError("phi must follow theta") 494 495 if psi >= 0 and psi != phi+1: 495 496 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")499 497 elif theta >= 0 or phi >= 0 or psi >= 0: 500 498 raise TypeError("oriented shapes must have both theta and phi and maybe psi") … … 717 715 718 716 719 #: Set of variables defined in the model that might contain C code720 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code']721 722 717 def _find_source_lines(model_info, kernel_module): 723 718 # type: (ModelInfo, ModuleType) -> None … … 725 720 Identify the location of the C source inside the model definition file. 726 721 727 This code runs through the source of the kernel module looking for lines 728 that contain C code (because it is a c function definition). Clearly 729 there are all sorts of reasons why this might not work (e.g., code 730 commented out in a triple-quoted line block, code built using string 731 concatenation, code defined in the branch of an 'if' block, code imported 732 from another file), but it should work properly in the 95% case, and for 733 the remainder, getting the incorrect line number will merely be 734 inconvenient. 735 """ 736 # Only need line numbers if we are creating a C module and the C symbols 737 # are defined. 738 if (callable(model_info.Iq) 739 or not any(hasattr(model_info, s) for s in C_SYMBOLS)): 722 This code runs through the source of the kernel module looking for 723 lines that start with 'Iq', 'Iqxy' or 'form_volume'. Clearly there are 724 all sorts of reasons why this might not work (e.g., code commented out 725 in a triple-quoted line block, code built using string concatenation, 726 or code defined in the branch of an 'if' block), but it should work 727 properly in the 95% case, and getting the incorrect line number will 728 be harmless. 729 """ 730 # Check if we need line numbers at all 731 if callable(model_info.Iq): 732 return None 733 734 if (model_info.Iq is None 735 and model_info.Iqxy is None 736 and model_info.Imagnetic is None 737 and model_info.form_volume is None): 740 738 return 741 739 742 # load the module source if we can740 # find the defintion lines for the different code blocks 743 741 try: 744 742 source = inspect.getsource(kernel_module) 745 743 except IOError: 746 744 return 747 748 # look for symbol at the start of the line 749 for lineno, line in enumerate(source.split('\n')): 750 for name in C_SYMBOLS: 751 if line.startswith(name): 752 # Add 1 since some compilers complain about "#line 0" 753 model_info.lineno[name] = lineno + 1 754 break 745 for k, v in enumerate(source.split('\n')): 746 if v.startswith('Imagnetic'): 747 model_info._Imagnetic_line = k+1 748 elif v.startswith('Iqxy'): 749 model_info._Iqxy_line = k+1 750 elif v.startswith('Iq'): 751 model_info._Iq_line = k+1 752 elif v.startswith('form_volume'): 753 model_info._form_volume_line = k+1 754 755 755 756 756 def make_model_info(kernel_module): … … 761 761 Fill in default values for parts of the module that are not provided. 762 762 763 Note: vectorized Iq and Iq ac/Iqabcfunctions will be created for python763 Note: vectorized Iq and Iqxy functions will be created for python 764 764 models when the model is first called, not when the model is loaded. 765 765 """ … … 790 790 info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 791 791 info.source = getattr(kernel_module, 'source', []) 792 info.c_code = getattr(kernel_module, 'c_code', None)793 792 # TODO: check the structure of the tests 794 793 info.tests = getattr(kernel_module, 'tests', []) … … 798 797 info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 799 798 info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 800 info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore801 info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore802 799 info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore 803 800 info.profile = getattr(kernel_module, 'profile', None) # type: ignore … … 814 811 info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 815 812 816 if callable(info.Iq) and parameters.has_2d:817 raise ValueError("oriented python models not supported")818 819 info.lineno = {}820 813 _find_source_lines(info, kernel_module) 821 814 … … 831 824 832 825 The structure should be mostly static, other than the delayed definition 833 of *Iq* , *Iqac* and *Iqabc* if they need to be defined.826 of *Iq* and *Iqxy* if they need to be defined. 834 827 """ 835 828 #: Full path to the file defining the kernel, if any. … … 913 906 structure_factor = None # type: bool 914 907 #: List of C source files used to define the model. The source files 915 #: should define the *Iq* function, and possibly *Iq ac* or *Iqabc* if the916 #: model defines orientation parameters. Files containing the most basic917 #: functions must appear first in the list, followed by the files that918 #: use those functions. Form factors are indicated by providing919 #: an:attr:`ER` function.908 #: should define the *Iq* function, and possibly *Iqxy*, though a default 909 #: *Iqxy = Iq(sqrt(qx**2+qy**2)* will be created if no *Iqxy* is provided. 910 #: Files containing the most basic functions must appear first in the list, 911 #: followed by the files that use those functions. Form factors are 912 #: indicated by providing a :attr:`ER` function. 920 913 source = None # type: List[str] 921 914 #: The set of tests that must pass. The format of the tests is described … … 942 935 #: See :attr:`ER` for details on the parameters. 943 936 VR = None # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 944 #: Arbitrary C code containing supporting functions, etc., to be inserted945 #: after everything in source. This can include Iq and Iqxy functions with946 #: the full function signature, including all parameters.947 c_code = None948 937 #: Returns the form volume for python-based models. Form volume is needed 949 938 #: for volume normalization in the polydispersity integral. If no … … 966 955 #: include the decimal point. See :mod:`generate` for more details. 967 956 Iq = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 968 #: Returns *I(qab, qc, a, b, ...)*. The interface follows :attr:`Iq`. 969 Iqac = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 970 #: Returns *I(qa, qb, qc, a, b, ...)*. The interface follows :attr:`Iq`. 971 Iqabc = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 957 #: Returns *I(qx, qy, a, b, ...)*. The interface follows :attr:`Iq`. 958 Iqxy = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 972 959 #: Returns *I(qx, qy, a, b, ...)*. The interface follows :attr:`Iq`. 973 960 Imagnetic = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] … … 985 972 #: Returns a random parameter set for the model 986 973 random = None # type: Optional[Callable[[], Dict[str, float]]] 987 #: Line numbers for symbols defining C code 988 lineno = None # type: Dict[str, int] 974 975 # line numbers within the python file for bits of C source, if defined 976 # NB: some compilers fail with a "#line 0" directive, so default to 1. 977 _Imagnetic_line = 1 978 _Iqxy_line = 1 979 _Iq_line = 1 980 _form_volume_line = 1 981 989 982 990 983 def __init__(self): -
sasmodels/models/_spherepy.py
r108e70e ref07e95 88 88 Iq.vectorized = True # Iq accepts an array of q values 89 89 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 90 94 def sesans(z, sld, sld_solvent, radius): 91 95 """ -
sasmodels/models/barbell.c
r108e70e r74768cb 90 90 91 91 static double 92 Iq ac(double qab, double qc,92 Iqxy(double qab, double qc, 93 93 double sld, double solvent_sld, 94 94 double radius_bell, double radius, double length) -
sasmodels/models/bcc_paracrystal.c
r108e70e r74768cb 107 107 108 108 109 static double Iq abc(double qa, double qb, double qc,109 static double Iqxy(double qa, double qb, double qc, 110 110 double dnn, double d_factor, double radius, 111 111 double sld, double solvent_sld) -
sasmodels/models/capped_cylinder.c
r108e70e r74768cb 115 115 116 116 static double 117 Iq ac(double qab, double qc,117 Iqxy(double qab, double qc, 118 118 double sld, double solvent_sld, double radius, 119 119 double radius_cap, double length) -
sasmodels/models/core_shell_bicelle.c
r108e70e r74768cb 67 67 68 68 static double 69 Iq ac(double qab, double qc,69 Iqxy(double qab, double qc, 70 70 double radius, 71 71 double thick_rim, -
sasmodels/models/core_shell_bicelle_elliptical.c
r108e70e r74768cb 71 71 72 72 static double 73 Iq abc(double qa, double qb, double qc,73 Iqxy(double qa, double qb, double qc, 74 74 double r_minor, 75 75 double x_core, -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
r108e70e r74768cb 79 79 80 80 static double 81 Iq abc(double qa, double qb, double qc,81 Iqxy(double qa, double qb, double qc, 82 82 double r_minor, 83 83 double x_core, … … 114 114 return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 115 115 } 116 -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py
r108e70e r110f69c 149 149 ["sld_rim", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Cylinder rim scattering length density"], 150 150 ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld", "Solvent scattering length density"], 151 ["sigma", "Ang", 0, [0, inf], "", "interfacial roughness"],152 151 ["theta", "degrees", 90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 153 152 ["phi", "degrees", 0, [-360, 360], "orientation", "rotation about beam"], 154 153 ["psi", "degrees", 0, [-360, 360], "orientation", "rotation about cylinder axis"], 154 ["sigma", "Ang", 0, [0, inf], "", "interfacial roughness"] 155 155 ] 156 156 -
sasmodels/models/core_shell_cylinder.c
r108e70e r74768cb 48 48 49 49 50 double Iq ac(double qab, double qc,50 double Iqxy(double qab, double qc, 51 51 double core_sld, 52 52 double shell_sld, -
sasmodels/models/core_shell_ellipsoid.c
r108e70e r74768cb 75 75 76 76 static double 77 Iq ac(double qab, double qc,77 Iqxy(double qab, double qc, 78 78 double radius_equat_core, 79 79 double x_core, -
sasmodels/models/core_shell_parallelepiped.c
re077231 r4493288 43 43 // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 44 44 // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 45 // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)46 // Code rewritten ; cross checked against hollow rectangular prism and realspace(PAK)45 // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 46 // Code rewritten (PAK) 47 47 48 48 const double half_q = 0.5*q; … … 102 102 103 103 static double 104 Iq abc(double qa, double qb, double qc,104 Iqxy(double qa, double qb, double qc, 105 105 double core_sld, 106 106 double arim_sld, … … 121 121 const double drC = crim_sld-solvent_sld; 122 122 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. 123 125 const double tA = length_a + 2.0*thick_rim_a; 124 126 const double tB = length_b + 2.0*thick_rim_b; -
sasmodels/models/core_shell_parallelepiped.py
r97be877 r10ee838 136 136 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 137 137 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 138 * **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 139 * Cross-checked against hollow rectangular prism and rectangular prism for 140 equal thickness overlapping sides, and by Monte Carlo sampling of points 141 within the shape for non-uniform, non-overlapping sides. 138 * **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017 139 * **Currently Under review by:** Paul Butler 142 140 """ 143 141 -
sasmodels/models/cylinder.c
r108e70e r74768cb 45 45 46 46 static double 47 Iq ac(double qab, double qc,47 Iqxy(double qab, double qc, 48 48 double sld, 49 49 double solvent_sld, -
sasmodels/models/ellipsoid.c
r108e70e r74768cb 39 39 40 40 static double 41 Iq ac(double qab, double qc,41 Iqxy(double qab, double qc, 42 42 double sld, 43 43 double sld_solvent, -
sasmodels/models/elliptical_cylinder.c
r108e70e r74768cb 54 54 55 55 static double 56 Iq abc(double qa, double qb, double qc,56 Iqxy(double qa, double qb, double qc, 57 57 double radius_minor, double r_ratio, double length, 58 58 double sld, double solvent_sld) -
sasmodels/models/fcc_paracrystal.c
r108e70e r74768cb 80 80 81 81 82 static double Iq abc(double qa, double qb, double qc,82 static double Iqxy(double qa, double qb, double qc, 83 83 double dnn, double d_factor, double radius, 84 84 double sld, double solvent_sld) -
sasmodels/models/hollow_cylinder.c
r108e70e r74768cb 52 52 53 53 static double 54 Iq ac(double qab, double qc,54 Iqxy(double qab, double qc, 55 55 double radius, double thickness, double length, 56 56 double sld, double solvent_sld) -
sasmodels/models/hollow_rectangular_prism.c
r108e70e r74768cb 85 85 } 86 86 87 double Iq abc(double qa, double qb, double qc,87 double Iqxy(double qa, double qb, double qc, 88 88 double sld, 89 89 double solvent_sld, -
sasmodels/models/line.py
r108e70e r2d81cfe 57 57 Iq.vectorized = True # Iq accepts an array of q values 58 58 59 60 59 def Iqxy(qx, qy, *args): 61 60 """ … … 70 69 71 70 Iqxy.vectorized = True # Iqxy accepts an array of qx qy values 72 73 # uncomment the following to test Iqxy in C models74 #del Iq, Iqxy75 #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 #"""80 71 81 72 def random(): -
sasmodels/models/parallelepiped.c
r108e70e r74768cb 53 53 54 54 static double 55 Iq abc(double qa, double qb, double qc,55 Iqxy(double qa, double qb, double qc, 56 56 double sld, 57 57 double solvent_sld, -
sasmodels/models/rectangular_prism.c
r108e70e r74768cb 71 71 72 72 73 double Iq abc(double qa, double qb, double qc,73 double Iqxy(double qa, double qb, double qc, 74 74 double sld, 75 75 double solvent_sld, -
sasmodels/models/sc_paracrystal.c
r108e70e r74768cb 82 82 83 83 static double 84 Iq abc(double qa, double qb, double qc,84 Iqxy(double qa, double qb, double qc, 85 85 double dnn, double d_factor, double radius, 86 86 double sld, double solvent_sld) -
sasmodels/models/stacked_disks.c
r108e70e r74768cb 142 142 143 143 static double 144 Iq ac(double qab, double qc,144 Iqxy(double qab, double qc, 145 145 double thick_core, 146 146 double thick_layer, -
sasmodels/models/triaxial_ellipsoid.c
r108e70e r74768cb 46 46 47 47 static double 48 Iq abc(double qa, double qb, double qc,48 Iqxy(double qa, double qb, double qc, 49 49 double sld, 50 50 double sld_solvent,
Note: See TracChangeset
for help on using the changeset viewer.