source: sasmodels/sasmodels/models/capped_cylinder.py @ 0507e09

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0507e09 was 0507e09, checked in by smk78, 5 months ago

Added link to source code to each model. Closes #883

  • Property mode set to 100644
File size: 7.3 KB
Line 
1r"""
2Definitions
3-----------
4
5Calculates the scattering from a cylinder with spherical section end-caps.
6Like :ref:`barbell`, this is a sphereocylinder with end caps that have a
7radius larger than that of the cylinder, but with the center of the end cap
8radius lying within the cylinder. This model simply becomes a convex
9lens when the length of the cylinder $L=0$. See the diagram for the details
10of the geometry and restrictions on parameter values.
11
12.. figure:: img/capped_cylinder_geometry.jpg
13
14    Capped cylinder geometry, where $r$ is *radius*, $R$ is *bell_radius* and
15    $L$ is *length*. Since the end cap radius $R \geq r$ and by definition
16    for this geometry $h < 0$, $h$ is then defined by $r$ and $R$ as
17    $h = - \sqrt{R^2 - r^2}$
18
19The scattered intensity $I(q)$ is calculated as
20
21.. math::
22
23    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right>
24
25where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as
26
27.. math::
28
29    A(q) =&\ \pi r^2L
30        \frac{\sin\left(\tfrac12 qL\cos\alpha\right)}
31            {\tfrac12 qL\cos\alpha}
32        \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\
33        &\ + 4 \pi R^3 \int_{-h/R}^1 dt
34        \cos\left[ q\cos\alpha
35            \left(Rt + h + {\tfrac12} L\right)\right]
36        \times (1-t^2)
37        \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]}
38             {qR\sin\alpha \left(1-t^2\right)^{1/2}}
39
40The $\left<\ldots\right>$ brackets denote an average of the structure over
41all orientations. $\left< A^2(q)\right>$ is then the form factor, $P(q)$.
42The scale factor is equivalent to the volume fraction of cylinders, each of
43volume, $V$. Contrast $\Delta\rho$ is the difference of scattering length
44densities of the cylinder and the surrounding solvent.
45
46The volume of the capped cylinder is (with $h$ as a positive value here)
47
48.. math::
49
50    V = \pi r_c^2 L + \tfrac{2\pi}{3}(R-h)^2(2R + h)
51
52
53and its radius of gyration is
54
55.. math::
56
57    R_g^2 =&\ \left[ \tfrac{12}{5}R^5
58        + R^4\left(6h+\tfrac32 L\right)
59        + R^2\left(4h^2 + L^2 + 4Lh\right)
60        + R^2\left(3Lh^2 + \tfrac32 L^2h\right) \right. \\
61        &\ \left. + \tfrac25 h^5 - \tfrac12 Lh^4 - \tfrac12 L^2h^3
62        + \tfrac14 L^3r^2 + \tfrac32 Lr^4 \right]
63        \left( 4R^3 6R^2h - 2h^3 + 3r^2L \right)^{-1}
64
65
66.. note::
67
68    The requirement that $R \geq r$ is not enforced in the model!
69    It is up to you to restrict this during analysis.
70
71The 2D scattering intensity is calculated similar to the 2D cylinder model.
72
73.. figure:: img/cylinder_angle_definition.png
74
75    Definition of the angles for oriented 2D cylinders.
76
77
78References
79----------
80
81.. [#] H Kaya, *J. Appl. Cryst.*, 37 (2004) 223-230
82.. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda
83   and errata)
84.. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659
85
86Source
87------
88
89`capped_cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/capped_cylinder.py>`_
90
91`capped_cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/capped_cylinder.c>`_
92
93Authorship and Verification
94----------------------------
95
96* **Author:** NIST IGOR/DANSE **Date:** pre 2010
97* **Last Modified by:** Paul Butler **Date:** September 30, 2016
98* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017
99* **Source added by :** Steve King **Date:** March 25, 2019
100"""
101
102import numpy as np
103from numpy import inf, sin, cos, pi
104
105name = "capped_cylinder"
106title = "Right circular cylinder with spherical end caps and uniform SLD"
107description = """That is, a sphereocylinder
108    with end caps that have a radius larger than
109    that of the cylinder and the center of the
110    end cap radius lies within the cylinder.
111    Note: As the length of cylinder -->0,
112    it becomes a Convex Lens.
113    It must be that radius <(=) radius_cap.
114    [Parameters];
115    scale: volume fraction of spheres,
116    background:incoherent background,
117    radius: radius of the cylinder,
118    length: length of the cylinder,
119    radius_cap: radius of the semi-spherical cap,
120    sld: SLD of the capped cylinder,
121    sld_solvent: SLD of the solvent.
122"""
123category = "shape:cylinder"
124# pylint: disable=bad-whitespace, line-too-long
125#             ["name", "units", default, [lower, upper], "type", "description"],
126parameters = [["sld",         "1e-6/Ang^2", 4, [-inf, inf], "sld",    "Cylinder scattering length density"],
127              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",    "Solvent scattering length density"],
128              ["radius",      "Ang",       20, [0, inf],    "volume", "Cylinder radius"],
129
130              # TODO: use an expression for cap radius with fixed bounds.
131              # The current form requires cap radius R bigger than cylinder radius r.
132              # Could instead use R/r in [1,inf], r/R in [0,1], or the angle between
133              # cylinder and cap in [0,90].  The problem is similar for the barbell
134              # model.  Propose r/R in [0,1] in both cases, with the model specifying
135              # cylinder radius in the capped cylinder model and sphere radius in the
136              # barbell model.  This leads to the natural value of zero for no cap
137              # in the capped cylinder, and zero for no bar in the barbell model.  In
138              # both models, one would be a pill.
139              ["radius_cap", "Ang",     20, [0, inf],    "volume", "Cap radius"],
140              ["length",     "Ang",    400, [0, inf],    "volume", "Cylinder length"],
141              ["theta",      "degrees", 60, [-360, 360], "orientation", "cylinder axis to beam angle"],
142              ["phi",        "degrees", 60, [-360, 360], "orientation", "rotation about beam"],
143             ]
144# pylint: enable=bad-whitespace, line-too-long
145
146source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"]
147have_Fq = True
148effective_radius_type = [
149    "equivalent cylinder excluded volume", "equivalent volume sphere",
150    "radius", "half length", "half total length",
151    ]
152
153def random():
154    """Return a random parameter set for the model."""
155    # TODO: increase volume range once problem with bell radius is fixed
156    # The issue is that bell radii of more than about 200 fail at high q
157    volume = 10**np.random.uniform(7, 9)
158    bar_volume = 10**np.random.uniform(-4, -1)*volume
159    bell_volume = volume - bar_volume
160    bell_radius = (bell_volume/6)**0.3333  # approximate
161    min_bar = bar_volume/np.pi/bell_radius**2
162    bar_length = 10**np.random.uniform(0, 3)*min_bar
163    bar_radius = np.sqrt(bar_volume/bar_length/np.pi)
164    if bar_radius > bell_radius:
165        bell_radius, bar_radius = bar_radius, bell_radius
166    pars = dict(
167        #background=0,
168        radius_cap=bell_radius,
169        radius=bar_radius,
170        length=bar_length,
171    )
172    return pars
173
174
175demo = dict(scale=1, background=0,
176            sld=6, sld_solvent=1,
177            radius=260, radius_cap=290, length=290,
178            theta=30, phi=15,
179            radius_pd=.2, radius_pd_n=1,
180            radius_cap_pd=.2, radius_cap_pd_n=1,
181            length_pd=.2, length_pd_n=10,
182            theta_pd=15, theta_pd_n=45,
183            phi_pd=15, phi_pd_n=1)
184q = 0.1
185# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
186qx = q*cos(pi/6.0)
187qy = q*sin(pi/6.0)
188tests = [
189    [{}, 0.075, 26.0698570695],
190    [{'theta':80., 'phi':10.}, (qx, qy), 0.561811990502],
191]
192del qx, qy  # not necessary to delete, but cleaner
Note: See TracBrowser for help on using the repository browser.