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 years ago

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

  • Property mode set to 100644
File size: 7.3 KB
RevLine 
[5d4777d]1r"""
[b0c4271]2Definitions
3-----------
4
[5d4777d]5Calculates the scattering from a cylinder with spherical section end-caps.
[eb69cce]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
[e65a3e7]8radius lying within the cylinder. This model simply becomes a convex
[eb69cce]9lens when the length of the cylinder $L=0$. See the diagram for the details
10of the geometry and restrictions on parameter values.
[5d4777d]11
[eb69cce]12.. figure:: img/capped_cylinder_geometry.jpg
[5d4777d]13
[eb69cce]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}$
[5d4777d]18
[eb69cce]19The scattered intensity $I(q)$ is calculated as
[5d4777d]20
21.. math::
22
[fcb33e4]23    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right>
[5d4777d]24
[fcb33e4]25where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as
[5d4777d]26
27.. math::
28
[eb69cce]29    A(q) =&\ \pi r^2L
[fcb33e4]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} \\
[19dcb933]33        &\ + 4 \pi R^3 \int_{-h/R}^1 dt
[fcb33e4]34        \cos\left[ q\cos\alpha
[19dcb933]35            \left(Rt + h + {\tfrac12} L\right)\right]
36        \times (1-t^2)
[fcb33e4]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}}
[5d4777d]39
[eb69cce]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)$.
[5d4777d]42The scale factor is equivalent to the volume fraction of cylinders, each of
[eb69cce]43volume, $V$. Contrast $\Delta\rho$ is the difference of scattering length
44densities of the cylinder and the surrounding solvent.
[5d4777d]45
46The volume of the capped cylinder is (with $h$ as a positive value here)
47
48.. math::
49
[19dcb933]50    V = \pi r_c^2 L + \tfrac{2\pi}{3}(R-h)^2(2R + h)
[5d4777d]51
52
[eb69cce]53and its radius of gyration is
[5d4777d]54
[19dcb933]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
[5d4777d]65
[19dcb933]66.. note::
[5d4777d]67
[eb69cce]68    The requirement that $R \geq r$ is not enforced in the model!
[19dcb933]69    It is up to you to restrict this during analysis.
[5d4777d]70
71The 2D scattering intensity is calculated similar to the 2D cylinder model.
72
[9802ab3]73.. figure:: img/cylinder_angle_definition.png
[5d4777d]74
75    Definition of the angles for oriented 2D cylinders.
76
77
[eb69cce]78References
79----------
[5d4777d]80
[b0c4271]81.. [#] H Kaya, *J. Appl. Cryst.*, 37 (2004) 223-230
[31df0c9]82.. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda
[b0c4271]83   and errata)
[0507e09]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>`_
[b0c4271]92
93Authorship and Verification
94----------------------------
[5d4777d]95
[b0c4271]96* **Author:** NIST IGOR/DANSE **Date:** pre 2010
97* **Last Modified by:** Paul Butler **Date:** September 30, 2016
[fcb33e4]98* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017
[0507e09]99* **Source added by :** Steve King **Date:** March 25, 2019
[5d4777d]100"""
[2d81cfe]101
102import numpy as np
[0b56f38]103from numpy import inf, sin, cos, pi
[5d4777d]104
105name = "capped_cylinder"
106title = "Right circular cylinder with spherical end caps and uniform SLD"
107description = """That is, a sphereocylinder
[485aee2]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,
[e65a3e7]112    it becomes a Convex Lens.
[2222134]113    It must be that radius <(=) radius_cap.
[485aee2]114    [Parameters];
115    scale: volume fraction of spheres,
116    background:incoherent background,
117    radius: radius of the cylinder,
118    length: length of the cylinder,
[2222134]119    radius_cap: radius of the semi-spherical cap,
[485aee2]120    sld: SLD of the capped cylinder,
[e65a3e7]121    sld_solvent: SLD of the solvent.
[5d4777d]122"""
[a5d0d00]123category = "shape:cylinder"
[dcdf29d]124# pylint: disable=bad-whitespace, line-too-long
[485aee2]125#             ["name", "units", default, [lower, upper], "type", "description"],
[42356c8]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"],
[dcdf29d]128              ["radius",      "Ang",       20, [0, inf],    "volume", "Cylinder radius"],
129
[485aee2]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.
[2222134]139              ["radius_cap", "Ang",     20, [0, inf],    "volume", "Cap radius"],
[dcdf29d]140              ["length",     "Ang",    400, [0, inf],    "volume", "Cylinder length"],
[9b79f29]141              ["theta",      "degrees", 60, [-360, 360], "orientation", "cylinder axis to beam angle"],
142              ["phi",        "degrees", 60, [-360, 360], "orientation", "rotation about beam"],
[485aee2]143             ]
[dcdf29d]144# pylint: enable=bad-whitespace, line-too-long
[485aee2]145
[26141cb]146source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"]
[71b751d]147have_Fq = True
[ee60aa7]148effective_radius_type = [
[b297ba9]149    "equivalent cylinder excluded volume", "equivalent volume sphere",
150    "radius", "half length", "half total length",
[ee60aa7]151    ]
[485aee2]152
[31df0c9]153def random():
[b297ba9]154    """Return a random parameter set for the model."""
[31df0c9]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
[2d81cfe]157    volume = 10**np.random.uniform(7, 9)
158    bar_volume = 10**np.random.uniform(-4, -1)*volume
159    bell_volume = volume - bar_volume
[31df0c9]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
[485aee2]175demo = dict(scale=1, background=0,
[e65a3e7]176            sld=6, sld_solvent=1,
[2222134]177            radius=260, radius_cap=290, length=290,
[485aee2]178            theta=30, phi=15,
179            radius_pd=.2, radius_pd_n=1,
[2222134]180            radius_cap_pd=.2, radius_cap_pd_n=1,
[485aee2]181            length_pd=.2, length_pd_n=10,
182            theta_pd=15, theta_pd_n=45,
183            phi_pd=15, phi_pd_n=1)
[0b56f38]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)
[2d81cfe]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.