Changes in / [1fbadb2:5400ad4] in sasmodels
- Location:
- sasmodels/models
- Files:
-
- 4 edited
-
core_shell_parallelepiped.c (modified) (2 diffs)
-
core_shell_parallelepiped.py (modified) (7 diffs)
-
parallelepiped.c (modified) (2 diffs)
-
parallelepiped.py (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_shell_parallelepiped.c
re077231 rdbf1a60 59 59 60 60 // outer integral (with gauss points), integration limits = 0, 1 61 // substitute d_cos_alpha for sin_alpha d_alpha 61 62 double outer_sum = 0; //initialize integral 62 63 for( int i=0; i<GAUSS_N; i++) { 63 64 const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 64 65 const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 65 66 // inner integral (with gauss points), integration limits = 0, pi/267 66 const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 68 67 const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 68 69 // inner integral (with gauss points), integration limits = 0, 1 70 // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta) 69 71 double inner_sum = 0.0; 70 72 for(int j=0; j<GAUSS_N; j++) { 71 const double beta= 0.5 * ( GAUSS_Z[j] + 1.0 );73 const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 72 74 double sin_beta, cos_beta; 73 SINCOS(M_PI_2* beta, sin_beta, cos_beta);75 SINCOS(M_PI_2*u, sin_beta, cos_beta); 74 76 const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 75 77 const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); … … 91 93 inner_sum += GAUSS_W[j] * f * f; 92 94 } 95 // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 93 96 inner_sum *= 0.5; 94 97 // now sum up the outer integral 95 98 outer_sum += GAUSS_W[i] * inner_sum; 96 99 } 100 // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 97 101 outer_sum *= 0.5; 98 102 -
sasmodels/models/core_shell_parallelepiped.py
r97be877 r5bc6d21 11 11 .. math:: 12 12 13 I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 13 I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle 14 + \text{background} 14 15 15 16 where $\langle \ldots \rangle$ is an average over all possible orientations 16 of the rectangular solid .17 18 The function calculated is the form factor of the rectangular solid below. 17 of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be 18 pulled out of the form factor term due to the multiple slds in the model. 19 19 20 The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 20 21 $A < B < C$. 21 22 22 .. image:: img/core_shell_parallelepiped_geometry.jpg 23 .. figure:: img/parallelepiped_geometry.jpg 24 25 Core of the core shell parallelepiped with the corresponding definition 26 of sides. 27 23 28 24 29 There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 25 30 (on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 26 $(=t_C)$ faces. The projection in the $AB$ plane is then 27 28 .. image:: img/core_shell_parallelepiped_projection.jpg 29 30 The volume of the solid is 31 $(=t_C)$ faces. The projection in the $AB$ plane is 32 33 .. figure:: img/core_shell_parallelepiped_projection.jpg 34 35 AB cut through the core-shell parallelipiped showing the cross secion of 36 four of the six shell slabs. As can be seen, this model leaves **"gaps"** 37 at the corners of the solid. 38 39 40 The total volume of the solid is thus given as 31 41 32 42 .. math:: 33 43 34 44 V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 35 36 **meaning that there are "gaps" at the corners of the solid.**37 45 38 46 The intensity calculated follows the :ref:`parallelepiped` model, with the 39 47 core-shell intensity being calculated as the square of the sum of the 40 amplitudes of the core and the slabs on the edges. 41 42 the scattering amplitude is computed for a particular orientation of the 43 core-shell parallelepiped with respect to the scattering vector and then 44 averaged over all possible orientations, where $\alpha$ is the angle between 45 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 46 the angle between projection of the particle in the $xy$ detector plane 47 and the $y$ axis. 48 49 .. math:: 50 51 F(Q) 48 amplitudes of the core and the slabs on the edges. The scattering amplitude is 49 computed for a particular orientation of the core-shell parallelepiped with 50 respect to the scattering vector and then averaged over all possible 51 orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis 52 of the parallelepiped, and $\beta$ is the angle between the projection of the 53 particle in the $xy$ detector plane and the $y$ axis. 54 55 .. math:: 56 57 P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha 58 \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta} 59 60 and 61 62 .. math:: 63 64 F(q,\alpha,\beta) 52 65 &= (\rho_\text{core}-\rho_\text{solvent}) 53 66 S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 54 67 &+ (\rho_\text{A}-\rho_\text{solvent}) 55 \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\68 \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ 56 69 &+ (\rho_\text{B}-\rho_\text{solvent}) 57 70 S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ … … 63 76 .. math:: 64 77 65 S(Q , L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} QL}78 S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L} 66 79 67 80 and … … 69 82 .. math:: 70 83 71 Q_A &= \sin\alpha \sin\beta \\72 Q_B &= \sin\alpha \cos\beta \\73 Q_C &= \cos\alpha84 Q_A &= q \sin\alpha \sin\beta \\ 85 Q_B &= q \sin\alpha \cos\beta \\ 86 Q_C &= q \cos\alpha 74 87 75 88 76 89 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 77 are the scattering length of the parallelepiped core, and the rectangular90 are the scattering lengths of the parallelepiped core, and the rectangular 78 91 slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 79 92 is the scattering length of the solvent. 93 94 .. note:: 95 96 the code actually implements two substitutions: $d(cos\alpha)$ is 97 substituted for -$sin\alpha \ d\alpha$ (note that in the 98 :ref:`parallelepiped` code this is explicitly implemented with 99 $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that 100 $du = \pi/2 \ d\beta$. Thus both integrals go from 0 to 1 rather than 0 101 to $\pi/2$. 80 102 81 103 FITTING NOTES … … 94 116 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 95 117 and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 96 to give an oblate or prolate particle, to give an effective radius ,97 for $S( Q)$ when $P(Q) * S(Q)$ is applied.118 to give an oblate or prolate particle, to give an effective radius 119 for $S(q)$ when $P(q) * S(q)$ is applied. 98 120 99 121 For 2d data the orientation of the particle is required, described using 100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below , for further122 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below. For further 101 123 details of the calculation and angular dispersions see :ref:`orientation`. 102 124 The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 103 125 $\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 104 126 105 For 2d, constraints must be applied during fitting to ensure that the106 inequality $A < B < C$ is not violated, and hence the correct definition107 of angles is preserved. The calculation will not report an error,108 but the results may be not correct.127 .. note:: For 2d, constraints must be applied during fitting to ensure that the 128 inequality $A < B < C$ is not violated, and hence the correct definition 129 of angles is preserved. The calculation will not report an error, 130 but the results may be not correct. 109 131 110 132 .. figure:: img/parallelepiped_angle_definition.png … … 113 135 Note that rotation $\theta$, initially in the $xz$ plane, is carried 114 136 out first, then rotation $\phi$ about the $z$ axis, finally rotation 115 $\Psi$ is now around the axis of the cylinder. The neutron or X-ray137 $\Psi$ is now around the axis of the particle. The neutron or X-ray 116 138 beam is along the $z$ axis. 117 139 … … 120 142 Examples of the angles for oriented core-shell parallelepipeds against the 121 143 detector plane. 144 145 146 Validation 147 ---------- 148 149 Cross-checked against hollow rectangular prism and rectangular prism for equal 150 thickness overlapping sides, and by Monte Carlo sampling of points within the 151 shape for non-uniform, non-overlapping sides. 152 122 153 123 154 References … … 135 166 136 167 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 137 * **Converted to sasmodels by:** Miguel Gonzale s**Date:** February 26, 2016168 * **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016 138 169 * **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 139 * Cross-checked against hollow rectangular prism and rectangular prism for140 equal thickness overlapping sides, and by Monte Carlo sampling of points141 within the shape for non-uniform, non-overlapping sides.142 170 """ 143 171 -
sasmodels/models/parallelepiped.c
r108e70e rdbf1a60 38 38 inner_total += GAUSS_W[j] * square(si1 * si2); 39 39 } 40 // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 40 41 inner_total *= 0.5; 41 42 … … 43 44 outer_total += GAUSS_W[i] * inner_total * si * si; 44 45 } 46 // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 45 47 outer_total *= 0.5; 46 48 -
sasmodels/models/parallelepiped.py
ref07e95 r5bc6d21 2 2 # Note: model title and parameter table are inserted automatically 3 3 r""" 4 The form factor is normalized by the particle volume.5 For information about polarised and magnetic scattering, see6 the :ref:`magnetism` documentation.7 8 4 Definition 9 5 ---------- 10 6 11 7 This model calculates the scattering from a rectangular parallelepiped 12 (\:numref:`parallelepiped-image`\). 13 If you need to apply polydispersity, see also :ref:`rectangular-prism`. 8 (:numref:`parallelepiped-image`). 9 If you need to apply polydispersity, see also :ref:`rectangular-prism`. For 10 information about polarised and magnetic scattering, see 11 the :ref:`magnetism` documentation. 14 12 15 13 .. _parallelepiped-image: … … 26 24 error, or fixing of some dimensions at expected values, may help. 27 25 28 The 1D scattering intensity $I(q)$ is calculated as: 26 The form factor is normalized by the particle volume and the 1D scattering 27 intensity $I(q)$ is then calculated as: 29 28 30 29 .. Comment by Miguel Gonzalez: … … 39 38 40 39 I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 41 \left< P(q, \alpha ) \right> + \text{background}40 \left< P(q, \alpha, \beta) \right> + \text{background} 42 41 43 42 where the volume $V = A B C$, the contrast is defined as 44 $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, 45 $P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 46 at an angle $\alpha$ (angle between the long axis C and $\vec q$), 47 and the averaging $\left<\ldots\right>$ is applied over all orientations. 43 $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ 44 is the form factor corresponding to a parallelepiped oriented 45 at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$ 46 ( the angle between the projection of the particle in the $xy$ detector plane 47 and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 48 orientations. 48 49 49 50 Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 50 form factor is given by (Mittelbach and Porod, 1961 )51 form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 51 52 52 53 .. math:: … … 66 67 \mu &= qB 67 68 68 The scattering intensity per unit volume is returned in units of |cm^-1|. 69 where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 70 applied. 69 71 70 72 NB: The 2nd virial coefficient of the parallelepiped is calculated based on … … 120 122 .. math:: 121 123 122 P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2 123 \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2 124 \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2 124 P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1} 125 {2}qA\cos\alpha)}\right]^2 126 \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1} 127 {2}qB\cos\beta)}\right]^2 128 \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1} 129 {2}qC\cos\gamma)}\right]^2 125 130 126 131 with … … 160 165 ---------- 161 166 162 P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 163 164 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854167 .. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*, 168 14 (1961) 185-211 169 .. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 165 170 166 171 Authorship and Verification
Note: See TracChangeset
for help on using the changeset viewer.
