[3330bb4] | 1 | # cylinder model |
---|
| 2 | # Note: model title and parameter table are inserted automatically |
---|
| 3 | r""" |
---|
| 4 | |
---|
| 5 | For information about polarised and magnetic scattering, see |
---|
| 6 | the :ref:`magnetism` documentation. |
---|
| 7 | |
---|
| 8 | Definition |
---|
| 9 | ---------- |
---|
| 10 | |
---|
| 11 | The output of the 2D scattering intensity function for oriented cylinders is |
---|
| 12 | given by (Guinier, 1955) |
---|
| 13 | |
---|
| 14 | .. math:: |
---|
| 15 | |
---|
| 16 | P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha).sin(\alpha) + \text{background} |
---|
| 17 | |
---|
| 18 | where |
---|
| 19 | |
---|
| 20 | .. math:: |
---|
| 21 | |
---|
| 22 | F(q,\alpha) = 2 (\Delta \rho) V |
---|
| 23 | \frac{\sin \left(\tfrac12 qL\cos\alpha \right)} |
---|
| 24 | {\tfrac12 qL \cos \alpha} |
---|
| 25 | \frac{J_1 \left(q R \sin \alpha\right)}{q R \sin \alpha} |
---|
| 26 | |
---|
| 27 | and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V =\pi R^2L$ |
---|
| 28 | is the volume of the cylinder, $L$ is the length of the cylinder, $R$ is the |
---|
| 29 | radius of the cylinder, and $\Delta\rho$ (contrast) is the scattering length |
---|
| 30 | density difference between the scatterer and the solvent. $J_1$ is the |
---|
| 31 | first order Bessel function. |
---|
| 32 | |
---|
| 33 | For randomly oriented particles: |
---|
| 34 | |
---|
| 35 | .. math:: |
---|
| 36 | |
---|
| 37 | F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}=\int_{0}^{1}{F^2(q,u)du} |
---|
| 38 | |
---|
| 39 | |
---|
[3fd0499] | 40 | Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with |
---|
| 41 | $sin(\alpha)=\sqrt{1-u^2}$. |
---|
[3330bb4] | 42 | |
---|
| 43 | The output of the 1D scattering intensity function for randomly oriented |
---|
| 44 | cylinders is thus given by |
---|
| 45 | |
---|
| 46 | .. math:: |
---|
| 47 | |
---|
| 48 | P(q) = \frac{\text{scale}}{V} |
---|
| 49 | \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | NB: The 2nd virial coefficient of the cylinder is calculated based on the |
---|
| 53 | radius and length values, and used as the effective radius for $S(q)$ |
---|
| 54 | when $P(q) \cdot S(q)$ is applied. |
---|
| 55 | |
---|
[eda8b30] | 56 | For 2d scattering from oriented cylinders, we define the direction of the |
---|
[3330bb4] | 57 | axis of the cylinder using two angles $\theta$ (note this is not the |
---|
| 58 | same as the scattering angle used in q) and $\phi$. Those angles |
---|
[eda8b30] | 59 | are defined in :numref:`cylinder-angle-definition` , for further details see :ref:`orientation` . |
---|
[3330bb4] | 60 | |
---|
| 61 | .. _cylinder-angle-definition: |
---|
| 62 | |
---|
[9802ab3] | 63 | .. figure:: img/cylinder_angle_definition.png |
---|
[3330bb4] | 64 | |
---|
[eda8b30] | 65 | Angles $\theta$ and $\phi$ orient the cylinder relative |
---|
[2d81cfe] | 66 | to the beam line coordinates, where the beam is along the $z$ axis. Rotation $\theta$, initially |
---|
[eda8b30] | 67 | in the $xz$ plane, is carried out first, then rotation $\phi$ about the $z$ axis. Orientation distributions |
---|
| 68 | are described as rotations about two perpendicular axes $\delta_1$ and $\delta_2$ |
---|
[9802ab3] | 69 | in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes. |
---|
[3330bb4] | 70 | |
---|
[3fd0499] | 71 | .. figure:: img/cylinder_angle_projection.png |
---|
| 72 | |
---|
| 73 | Examples for oriented cylinders. |
---|
| 74 | |
---|
[31df0c9] | 75 | The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. |
---|
[3330bb4] | 76 | |
---|
| 77 | Validation |
---|
| 78 | ---------- |
---|
| 79 | |
---|
| 80 | Validation of the code was done by comparing the output of the 1D model |
---|
| 81 | to the output of the software provided by the NIST (Kline, 2006). |
---|
| 82 | The implementation of the intensity for fully oriented cylinders was done |
---|
| 83 | by averaging over a uniform distribution of orientations using |
---|
| 84 | |
---|
| 85 | .. math:: |
---|
| 86 | |
---|
| 87 | P(q) = \int_0^{\pi/2} d\phi |
---|
| 88 | \int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | where $p(\theta,\phi) = 1$ is the probability distribution for the orientation |
---|
| 92 | and $P_0(q,\theta)$ is the scattering intensity for the fully oriented |
---|
| 93 | system, and then comparing to the 1D result. |
---|
| 94 | |
---|
| 95 | References |
---|
| 96 | ---------- |
---|
| 97 | |
---|
[0507e09] | 98 | .. [#] J. Pedersen, *Adv. Colloid Interface Sci.*, 70 (1997) 171-210 |
---|
| 99 | .. [#] G. Fournet, *Bull. Soc. Fr. Mineral. Cristallogr.*, 74 (1951) 39-113 |
---|
| 100 | .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 |
---|
| 101 | |
---|
| 102 | Source |
---|
| 103 | ------ |
---|
| 104 | |
---|
| 105 | `cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/cylinder.py>`_ |
---|
| 106 | |
---|
| 107 | `cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/cylinder.c>`_ |
---|
| 108 | |
---|
| 109 | Authorship and Verification |
---|
| 110 | ---------------------------- |
---|
| 111 | |
---|
| 112 | * **Author:** |
---|
| 113 | * **Last Modified by:** |
---|
| 114 | * **Last Reviewed by:** |
---|
| 115 | * **Source added by :** Steve King **Date:** March 25, 2019 |
---|
[3330bb4] | 116 | """ |
---|
| 117 | |
---|
| 118 | import numpy as np # type: ignore |
---|
| 119 | from numpy import pi, inf # type: ignore |
---|
| 120 | |
---|
| 121 | name = "cylinder" |
---|
| 122 | title = "Right circular cylinder with uniform scattering length density." |
---|
| 123 | description = """ |
---|
| 124 | f(q,alpha) = 2*(sld - sld_solvent)*V*sin(qLcos(alpha)/2)) |
---|
| 125 | /[qLcos(alpha)/2]*J1(qRsin(alpha))/[qRsin(alpha)] |
---|
| 126 | |
---|
| 127 | P(q,alpha)= scale/V*f(q,alpha)^(2)+background |
---|
| 128 | V: Volume of the cylinder |
---|
| 129 | R: Radius of the cylinder |
---|
| 130 | L: Length of the cylinder |
---|
| 131 | J1: The bessel function |
---|
| 132 | alpha: angle between the axis of the |
---|
| 133 | cylinder and the q-vector for 1D |
---|
| 134 | :the ouput is P(q)=scale/V*integral |
---|
| 135 | from pi/2 to zero of... |
---|
| 136 | f(q,alpha)^(2)*sin(alpha)*dalpha + background |
---|
| 137 | """ |
---|
| 138 | category = "shape:cylinder" |
---|
| 139 | |
---|
| 140 | # [ "name", "units", default, [lower, upper], "type", "description"], |
---|
| 141 | parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", |
---|
| 142 | "Cylinder scattering length density"], |
---|
| 143 | ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", |
---|
| 144 | "Solvent scattering length density"], |
---|
| 145 | ["radius", "Ang", 20, [0, inf], "volume", |
---|
| 146 | "Cylinder radius"], |
---|
| 147 | ["length", "Ang", 400, [0, inf], "volume", |
---|
| 148 | "Cylinder length"], |
---|
[9b79f29] | 149 | ["theta", "degrees", 60, [-360, 360], "orientation", |
---|
| 150 | "cylinder axis to beam angle"], |
---|
[2d81cfe] | 151 | ["phi", "degrees", 60, [-360, 360], "orientation", |
---|
[9b79f29] | 152 | "rotation about beam"], |
---|
[3330bb4] | 153 | ] |
---|
| 154 | |
---|
[2d81cfe] | 155 | source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] |
---|
[71b751d] | 156 | have_Fq = True |
---|
[ee60aa7] | 157 | effective_radius_type = [ |
---|
[99658f6] | 158 | "excluded volume", "equivalent volume sphere", "radius", |
---|
[ee60aa7] | 159 | "half length", "half min dimension", "half max dimension", "half diagonal", |
---|
| 160 | ] |
---|
[3330bb4] | 161 | |
---|
[31df0c9] | 162 | def random(): |
---|
[b297ba9] | 163 | """Return a random parameter set for the model.""" |
---|
[2d81cfe] | 164 | volume = 10**np.random.uniform(5, 12) |
---|
| 165 | length = 10**np.random.uniform(-2, 2)*volume**0.333 |
---|
| 166 | radius = np.sqrt(volume/length/np.pi) |
---|
[31df0c9] | 167 | pars = dict( |
---|
| 168 | #scale=1, |
---|
| 169 | #background=0, |
---|
| 170 | length=length, |
---|
| 171 | radius=radius, |
---|
| 172 | ) |
---|
| 173 | return pars |
---|
| 174 | |
---|
| 175 | |
---|
[3330bb4] | 176 | # parameters for demo |
---|
| 177 | demo = dict(scale=1, background=0, |
---|
| 178 | sld=6, sld_solvent=1, |
---|
| 179 | radius=20, length=300, |
---|
| 180 | theta=60, phi=60, |
---|
| 181 | radius_pd=.2, radius_pd_n=9, |
---|
| 182 | length_pd=.2, length_pd_n=10, |
---|
| 183 | theta_pd=10, theta_pd_n=5, |
---|
| 184 | phi_pd=10, phi_pd_n=5) |
---|
| 185 | |
---|
[304c775] | 186 | # pylint: disable=bad-whitespace, line-too-long |
---|
[3330bb4] | 187 | qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) |
---|
| 188 | # After redefinition of angles, find new tests values. Was 10 10 in old coords |
---|
[2d81cfe] | 189 | tests = [ |
---|
| 190 | [{}, 0.2, 0.042761386790780453], |
---|
| 191 | [{}, [0.2], [0.042761386790780453]], |
---|
| 192 | # new coords |
---|
| 193 | [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], |
---|
| 194 | [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], |
---|
| 195 | # old coords |
---|
| 196 | #[{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], |
---|
| 197 | #[{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], |
---|
| 198 | ] |
---|
[3330bb4] | 199 | del qx, qy # not necessary to delete, but cleaner |
---|
[304c775] | 200 | |
---|
| 201 | # Default radius and length |
---|
[b297ba9] | 202 | def calc_volume(radius, length): |
---|
| 203 | """Return form volume for cylinder.""" |
---|
| 204 | return pi*radius**2*length |
---|
| 205 | def calc_r_effs(radius, length): |
---|
| 206 | """Return effective radii for modes 0-7 of cylinder.""" |
---|
| 207 | return [ |
---|
| 208 | 0., |
---|
| 209 | 0.5*(0.75*radius*(2.0*radius*length |
---|
| 210 | + (radius + length)*(pi*radius + length)))**(1./3.), |
---|
| 211 | (0.75*radius**2*length)**(1./3.), |
---|
| 212 | radius, |
---|
| 213 | length/2., |
---|
| 214 | min(radius, length/2.), |
---|
| 215 | max(radius, length/2.), |
---|
| 216 | np.sqrt(4*radius**2 + length**2)/2., |
---|
| 217 | ] |
---|
| 218 | r_effs = calc_r_effs(parameters[2][2], parameters[3][2]) |
---|
| 219 | cyl_vol = calc_volume(parameters[2][2], parameters[3][2]) |
---|
[304c775] | 220 | tests.extend([ |
---|
[b297ba9] | 221 | ({'radius_effective_mode': 0}, 0.1, None, None, r_effs[0], cyl_vol, 1.0), |
---|
| 222 | ({'radius_effective_mode': 1}, 0.1, None, None, r_effs[1], None, None), |
---|
| 223 | ({'radius_effective_mode': 2}, 0.1, None, None, r_effs[2], None, None), |
---|
| 224 | ({'radius_effective_mode': 3}, 0.1, None, None, r_effs[3], None, None), |
---|
| 225 | ({'radius_effective_mode': 4}, 0.1, None, None, r_effs[4], None, None), |
---|
| 226 | ({'radius_effective_mode': 5}, 0.1, None, None, r_effs[5], None, None), |
---|
| 227 | ({'radius_effective_mode': 6}, 0.1, None, None, r_effs[6], None, None), |
---|
| 228 | ({'radius_effective_mode': 7}, 0.1, None, None, r_effs[7], None, None), |
---|
[304c775] | 229 | ]) |
---|
[b297ba9] | 230 | del r_effs, cyl_vol |
---|
[304c775] | 231 | # pylint: enable=bad-whitespace, line-too-long |
---|
| 232 | |
---|
[2d81cfe] | 233 | # ADDED by: RKH ON: 18Mar2016 renamed sld's etc |
---|