source: sasmodels/sasmodels/models/parallelepiped.py @ 3e428ec

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 3e428ec was 3e428ec, checked in by Doucet, Mathieu <doucetm@…>, 8 years ago

Cleaning up

  • Property mode set to 100644
File size: 5.6 KB
Line 
1# parallelepiped model
2# Note: model title and parameter table are inserted automatically
3r"""
4The form factor is normalized by the particle volume.
5
6For information about polarised and magnetic scattering, click here_.
7
8Definition
9----------
10
11This model provides the form factor, *P(q)*, for a rectangular parallelepiped
12(below) where the form factor is normalized by the volume of the
13parallelepiped. If you need to apply polydispersity, see also the
14RectangularPrismModel_.
15
16The calculated form factor is:
17
18.. math::
19
20    P(Q) = {\text{scale} \over V} F^2(Q) + \text{background}
21
22where the volume *V* = *A B C* and the averaging < > is applied over all
23orientations for 1D.
24
25.. image:: img/parallelepiped.jpg
26
27*Figure. Parallelepiped with the corresponding Definition of sides.
28
29The edge of the solid must satisfy the condition that** *A* < *B* < *C*.
30Then, assuming *a* = *A* / *B* < 1, *b* = *B* / *B* = 1, and
31*c* = *C* / *B* > 1, the form factor is
32
33.. math::
34
35    P(q) = \frac{\textstyle{scale}}{V}\int_0^1 \phi(\mu \sqrt{1-\sigma^2},a)
36    [S(\mu c \sigma/2)]^2 d\sigma
37
38with
39
40.. math::
41
42    \phi(\mu,a) = \int_0^1 \{S[\frac{\mu}{2}\cos(\frac{\pi}{2}u)]
43    S[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)]\}^2 du
44
45    S(x) = \frac{\sin x}{x}
46
47    \mu = qB
48
49and the contrast is defined as
50
51.. math::
52
53    \Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent}
54
55The scattering intensity per unit volume is returned in units of |cm^-1|;
56ie, *I(q)* = |phi| *P(q)*\ .
57
58NB: The 2nd virial coefficient of the parallelpiped is calculated based on
59the averaged effective radius (= sqrt(*short_a* \* *short_b* / |pi|)) and
60length(= *long_c*) values, and used as the effective radius for
61*S(Q)* when *P(Q)* \* *S(Q)* is applied.
62
63To provide easy access to the orientation of the parallelepiped, we define
64three angles |theta|, |phi| and |bigpsi|. The definition of |theta| and |phi|
65is the same as for the cylinder model (see also figures below).
66The angle |bigpsi| is the rotational angle around the *long_c* axis against
67the *q* plane. For example, |bigpsi| = 0 when the *short_b* axis is parallel
68to the *x*-axis of the detector.
69
70
71.. _parallelepiped-orientation:
72
73.. figure:: img/orientation.jpg
74
75    Definition of the angles for oriented parallelepipeds.
76
77.. figure:: img/orientation2.jpg
78
79    Examples of the angles for oriented parallelepipeds against the detector plane.
80
81
82Validation
83----------
84
85Validation of the code was done by comparing the output of the 1D calculation
86to the angular average of the output of a 2D calculation over all possible
87angles. The Figure below shows the comparison where the solid dot refers to
88averaged 2D while the line represents the result of the 1D calculation (for
89the averaging, 76, 180, 76 points are taken for the angles of |theta|, |phi|,
90and |psi| respectively).
91
92.. _parallelepiped-compare:
93
94.. figure:: img/parallelepiped_compare.jpg
95
96*Figure. Comparison between 1D and averaged 2D.*
97
98This model reimplements the form factor calculations implemented in a c-library
99provided by the NIST Center for Neutron Research (Kline, 2006).
100
101"""
102
103from numpy import pi, inf, sqrt
104
105name = "parallelepiped"
106title = "Rectangular parallelepiped with uniform scattering length density."
107description = """
108     P(q)= scale/V*integral from 0 to 1 of ...
109           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
110
111            phi(mu,a) = integral from 0 to 1 of ..
112        (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
113            S(x) = sin(x)/x
114        mu = q*B
115"""
116category = "shape:parallelpiped"
117
118#             ["name", "units", default, [lower, upper], "type","description"],
119parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "",
120               "Parallelepiped scattering length density"],
121              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
122               "Solvent scattering length density"],
123              ["a_side", "Ang", 35, [0, inf], "volume",
124               "Shorter side of the parallelepiped"],
125              ["b_side", "Ang", 75, [0, inf], "volume",
126               "Second side of the parallelepiped"],
127              ["c_side", "Ang", 400, [0, inf], "volume",
128               "Larger side of the parallelepiped"],
129              ["theta", "degrees", 60, [-inf, inf], "orientation",
130               "In plane angle"],
131              ["phi", "degrees", 60, [-inf, inf], "orientation",
132               "Out of plane angle"],
133              ["psi", "degrees", 60, [-inf, inf], "orientation",
134               "Rotation angle around its own c axis against q plane"],
135             ]
136
137source = ["lib/J1.c", "lib/gauss76.c", "parallelepiped.c"]
138
139def ER(a_side, b_side, c_side):
140
141    # surface average radius (rough approximation)
142    surf_rad = sqrt(a_side * b_side / pi)
143
144    # DiamCyl recoded here (to check and possibly put in a library?)
145    a = surf_rad
146    b = 0.5 * c_side
147    t1 = a * a * b
148    t2 = 1.0 + (b / a) * (1.0 + a / b / 2.0) * (1.0 + pi * a / b / 2.0)
149    ddd = 3.0 * t1 * t2
150
151    return 0.5 * (ddd) ** (1. / 3.)
152
153# parameters for demo
154demo = dict(scale=1, background=0,
155            sld=6.3e-6, solvent_sld=1.0e-6,
156            a_side=35, b_side=75, c_side=400,
157            theta=45, phi=30, psi=15,
158            a_side_pd=0.1, a_side_pd_n=10,
159            b_side_pd=0.1, b_side_pd_n=1,
160            c_side_pd=0.1, c_side_pd_n=10,
161            theta_pd=10, theta_pd_n=1,
162            phi_pd=10, phi_pd_n=1,
163            psi_pd=10, psi_pd_n=10)
164
165# For testing against the old sasview models, include the converted parameter
166# names and the target sasview model name.
167oldname = 'ParallelepipedModel'
168oldpars = dict(theta='parallel_theta', phi='parallel_phi', psi='parallel_psi',
169               a_side='short_a', b_side='short_b', c_side='long_c',
170               sld='sldPipe', solvent_sld='sldSolv')
171
Note: See TracBrowser for help on using the repository browser.