source: sasmodels/sasmodels/models/elliptical_cylinder.py @ db1d9d5

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since db1d9d5 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

  • Property mode set to 100644
File size: 6.4 KB
Line 
1# pylint: disable=line-too-long
2r"""
3
4.. figure:: img/elliptical_cylinder_geometry.png
5
6   Elliptical cylinder geometry $a = r_\text{minor}$
7   and $\nu = r_\text{major} / r_\text{minor}$ is the *axis_ratio*.
8
9The function calculated is
10
11.. math::
12
13    I(\vec q)=\frac{1}{V_\text{cyl}}\int{d\psi}\int{d\phi}\int{
14        p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\alpha)d\alpha}
15
16with the functions
17
18.. math::
19
20    F(q,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab}
21
22where
23
24.. math::
25
26    a = qr'\sin(\alpha)
27
28    b = q\frac{L}{2}\cos(\alpha)
29
30    r'=\frac{r_{minor}}{\sqrt{2}}\sqrt{(1+\nu^{2}) + (1-\nu^{2})cos(\psi)}
31
32
33and the angle $\psi$ is defined as the orientation of the major axis of the
34ellipse with respect to the vector $\vec q$. The angle $\alpha$ is the angle
35between the axis of the cylinder and $\vec q$.
36
37
38For 1D scattering, with no preferred orientation, the form factor is averaged over all possible orientations and normalized
39by the particle volume
40
41.. math::
42
43    P(q) = \text{scale}  <F^2> / V
44
45For 2d data the orientation of the particle is required, described using a different set
46of angles as in the diagrams below, for further details of the calculation and angular
47dispersions  see :ref:`orientation` .
48
49
50.. figure:: img/elliptical_cylinder_angle_definition.png
51
52    Note that the angles here are not the same as in the equations for the scattering function.
53    Rotation $\theta$, initially in the $xz$ plane, is carried out first, then
54    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder.
55    The neutron or X-ray beam is along the $z$ axis.
56
57.. figure:: img/elliptical_cylinder_angle_projection.png
58
59    Examples of the angles for oriented elliptical cylinders against the
60    detector plane, with $\Psi$ = 0.
61
62The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
63
64
65NB: The 2nd virial coefficient of the cylinder is calculated based on the
66averaged radius $(=\sqrt{r_\text{minor}^2 * \text{axis ratio}})$ and length
67values, and used as the effective radius for $S(Q)$ when $P(Q)*S(Q)$ is applied.
68
69
70Validation
71----------
72
73Validation of our code was done by comparing the output of the 1D calculation
74to the angular average of the output of the 2D calculation over all possible
75angles.
76
77In the 2D average, more binning in the angle $\phi$ is necessary to get the
78proper result. The following figure shows the results of the averaging by
79varying the number of angular bins.
80
81.. figure:: img/elliptical_cylinder_averaging.png
82
83    The intensities averaged from 2D over different numbers of bins and angles.
84
85References
86----------
87
88.. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, New York, (1987) [see table 3.4]
89.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659
90
91Source
92------
93
94`elliptical_cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/elliptical_cylinder.py>`_
95
96`elliptical_cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/elliptical_cylinder.c>`_
97
98Authorship and Verification
99----------------------------
100
101* **Author:**
102* **Last Modified by:**
103* **Last Reviewed by:**  Richard Heenan - corrected equation in docs **Date:** December 21, 2016
104* **Source added by :** Steve King **Date:** March 25, 2019
105"""
106
107import numpy as np
108from numpy import pi, inf, sqrt, sin, cos
109
110name = "elliptical_cylinder"
111title = "Form factor for an elliptical cylinder."
112description = """
113    Form factor for an elliptical cylinder.
114    See L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray and Neutron Scattering, Plenum, New York, (1987).
115"""
116category = "shape:cylinder"
117
118# pylint: disable=bad-whitespace, line-too-long
119#             ["name", "units", default, [lower, upper], "type","description"],
120parameters = [["radius_minor",     "Ang",        20.0,  [0, inf],    "volume",      "Ellipse minor radius"],
121              ["axis_ratio",   "",          1.5,   [1, inf],    "volume",      "Ratio of major radius over minor radius"],
122              ["length",      "Ang",        400.0, [1, inf],    "volume",      "Length of the cylinder"],
123              ["sld",         "1e-6/Ang^2", 4.0,   [-inf, inf], "sld",         "Cylinder scattering length density"],
124              ["sld_solvent", "1e-6/Ang^2", 1.0,   [-inf, inf], "sld",         "Solvent scattering length density"],
125              ["theta",       "degrees",    90.0,  [-360, 360], "orientation", "cylinder axis to beam angle"],
126              ["phi",         "degrees",    0,     [-360, 360], "orientation", "rotation about beam"],
127              ["psi",         "degrees",    0,     [-360, 360], "orientation", "rotation about cylinder axis"]]
128
129# pylint: enable=bad-whitespace, line-too-long
130
131source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"]
132have_Fq = True
133radius_effective_modes = [
134    "equivalent cylinder excluded volume",
135    "equivalent volume sphere", "average radius", "min radius", "max radius",
136    "equivalent circular cross-section",
137    "half length", "half min dimension", "half max dimension", "half diagonal",
138    ]
139
140demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0,
141            sld=4.0, sld_solvent=1.0, theta=10.0, phi=20, psi=30,
142            theta_pd=10, phi_pd=2, psi_pd=3)
143
144def random():
145    """Return a random parameter set for the model."""
146    # V = pi * radius_major * radius_minor * length;
147    volume = 10**np.random.uniform(3, 9)
148    length = 10**np.random.uniform(1, 3)
149    axis_ratio = 10**np.random.uniform(0, 2)
150    radius_minor = sqrt(volume/length/axis_ratio)
151    volfrac = 10**np.random.uniform(-4, -2)
152    pars = dict(
153        #background=0, sld=0, sld_solvent=1,
154        scale=1e9*volfrac/volume,
155        length=length,
156        radius_minor=radius_minor,
157        axis_ratio=axis_ratio,
158    )
159    return pars
160
161q = 0.1
162# april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct!
163qx = q*cos(pi/6.0)
164qy = q*sin(pi/6.0)
165
166tests = [
167    #[{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024],
168    #[{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1],
169
170    # The SasView test result was 0.00169, with a background of 0.001
171    [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'sld': 4.0, 'length':400.0,
172      'sld_solvent':1.0, 'background':0.0},
173     0.001, 675.504402],
174    #[{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ],
175]
Note: See TracBrowser for help on using the repository browser.