source: sasmodels/sasmodels/models/barbell.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, 9 months ago

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

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