Changes in / [c64a68e:33969b6] 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
rdbf1a60 re077231 59 59 60 60 // outer integral (with gauss points), integration limits = 0, 1 61 // substitute d_cos_alpha for sin_alpha d_alpha62 61 double outer_sum = 0; //initialize integral 63 62 for( int i=0; i<GAUSS_N; i++) { 64 63 const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 65 64 const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 65 66 // inner integral (with gauss points), integration limits = 0, pi/2 66 67 const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 67 68 const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 68 69 // inner integral (with gauss points), integration limits = 0, 170 // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta)71 69 double inner_sum = 0.0; 72 70 for(int j=0; j<GAUSS_N; j++) { 73 const double u= 0.5 * ( GAUSS_Z[j] + 1.0 );71 const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 74 72 double sin_beta, cos_beta; 75 SINCOS(M_PI_2* u, sin_beta, cos_beta);73 SINCOS(M_PI_2*beta, sin_beta, cos_beta); 76 74 const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 77 75 const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); … … 93 91 inner_sum += GAUSS_W[j] * f * f; 94 92 } 95 // now complete change of inner integration variable (1-0)/(1-(-1))= 0.596 93 inner_sum *= 0.5; 97 94 // now sum up the outer integral 98 95 outer_sum += GAUSS_W[i] * inner_sum; 99 96 } 100 // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5101 97 outer_sum *= 0.5; 102 98 -
sasmodels/models/core_shell_parallelepiped.py
r5bc6d21 r97be877 11 11 .. math:: 12 12 13 I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle 14 + \text{background} 13 I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 15 14 16 15 where $\langle \ldots \rangle$ is an average over all possible orientations 17 of the rectangular solid , and the usual $\Delta \rho^2 \ V^2$ term cannot be18 pulled out of the form factor term due to the multiple slds in the model. 19 16 of the rectangular solid. 17 18 The function calculated is the form factor of the rectangular solid below. 20 19 The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 21 20 $A < B < C$. 22 21 23 .. figure:: img/parallelepiped_geometry.jpg 24 25 Core of the core shell parallelepiped with the corresponding definition 26 of sides. 27 22 .. image:: img/core_shell_parallelepiped_geometry.jpg 28 23 29 24 There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 30 25 (on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 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 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 41 31 42 32 .. math:: 43 33 44 34 V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 35 36 **meaning that there are "gaps" at the corners of the solid.** 45 37 46 38 The intensity calculated follows the :ref:`parallelepiped` model, with the 47 39 core-shell intensity being calculated as the square of the sum of the 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) 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) 65 52 &= (\rho_\text{core}-\rho_\text{solvent}) 66 53 S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 67 54 &+ (\rho_\text{A}-\rho_\text{solvent}) 68 \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\55 \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 69 56 &+ (\rho_\text{B}-\rho_\text{solvent}) 70 57 S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ … … 76 63 .. math:: 77 64 78 S(Q _X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_XL}65 S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 79 66 80 67 and … … 82 69 .. math:: 83 70 84 Q_A &= q\sin\alpha \sin\beta \\85 Q_B &= q\sin\alpha \cos\beta \\86 Q_C &= q\cos\alpha71 Q_A &= \sin\alpha \sin\beta \\ 72 Q_B &= \sin\alpha \cos\beta \\ 73 Q_C &= \cos\alpha 87 74 88 75 89 76 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 90 are the scattering length sof the parallelepiped core, and the rectangular77 are the scattering length of the parallelepiped core, and the rectangular 91 78 slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 92 79 is the scattering length of the solvent. 93 94 .. note::95 96 the code actually implements two substitutions: $d(cos\alpha)$ is97 substituted for -$sin\alpha \ d\alpha$ (note that in the98 :ref:`parallelepiped` code this is explicitly implemented with99 $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that100 $du = \pi/2 \ d\beta$. Thus both integrals go from 0 to 1 rather than 0101 to $\pi/2$.102 80 103 81 FITTING NOTES … … 116 94 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 117 95 and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 118 to give an oblate or prolate particle, to give an effective radius 119 for $S( q)$ when $P(q) * S(q)$ is applied.96 to give an oblate or prolate particle, to give an effective radius, 97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 120 98 121 99 For 2d data the orientation of the particle is required, described using 122 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below . For further100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 123 101 details of the calculation and angular dispersions see :ref:`orientation`. 124 102 The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 125 103 $\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 126 104 127 .. note::For 2d, constraints must be applied during fitting to ensure that the128 129 130 105 For 2d, constraints must be applied during fitting to ensure that the 106 inequality $A < B < C$ is not violated, and hence the correct definition 107 of angles is preserved. The calculation will not report an error, 108 but the results may be not correct. 131 109 132 110 .. figure:: img/parallelepiped_angle_definition.png … … 135 113 Note that rotation $\theta$, initially in the $xz$ plane, is carried 136 114 out first, then rotation $\phi$ about the $z$ axis, finally rotation 137 $\Psi$ is now around the axis of the particle. The neutron or X-ray115 $\Psi$ is now around the axis of the cylinder. The neutron or X-ray 138 116 beam is along the $z$ axis. 139 117 … … 142 120 Examples of the angles for oriented core-shell parallelepipeds against the 143 121 detector plane. 144 145 146 Validation147 ----------148 149 Cross-checked against hollow rectangular prism and rectangular prism for equal150 thickness overlapping sides, and by Monte Carlo sampling of points within the151 shape for non-uniform, non-overlapping sides.152 153 122 154 123 References … … 166 135 167 136 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 168 * **Converted to sasmodels by:** Miguel Gonzale z**Date:** February 26, 2016137 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 169 138 * **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 139 * Cross-checked against hollow rectangular prism and rectangular prism for 140 equal thickness overlapping sides, and by Monte Carlo sampling of points 141 within the shape for non-uniform, non-overlapping sides. 170 142 """ 171 143 -
sasmodels/models/parallelepiped.c
rdbf1a60 r108e70e 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.541 40 inner_total *= 0.5; 42 41 … … 44 43 outer_total += GAUSS_W[i] * inner_total * si * si; 45 44 } 46 // now complete change of outer integration variable (1-0)/(1-(-1))= 0.547 45 outer_total *= 0.5; 48 46 -
sasmodels/models/parallelepiped.py
rb343226 ref07e95 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, see 6 the :ref:`magnetism` documentation. 7 4 8 Definition 5 9 ---------- 6 10 7 This model calculates the scattering from a rectangular parallelepiped 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. 11 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`. 12 14 13 15 .. _parallelepiped-image: … … 24 26 error, or fixing of some dimensions at expected values, may help. 25 27 26 The form factor is normalized by the particle volume and the 1D scattering 27 intensity $I(q)$ is then calculated as: 28 The 1D scattering intensity $I(q)$ is calculated as: 28 29 29 30 .. Comment by Miguel Gonzalez: … … 38 39 39 40 I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 40 \left< P(q, \alpha , \beta) \right> + \text{background}41 \left< P(q, \alpha) \right> + \text{background} 41 42 42 43 where the volume $V = A B C$, the contrast is defined as 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. 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. 49 48 50 49 Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 51 form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_)50 form factor is given by (Mittelbach and Porod, 1961) 52 51 53 52 .. math:: … … 67 66 \mu &= qB 68 67 69 where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 70 applied. 68 The scattering intensity per unit volume is returned in units of |cm^-1|. 71 69 72 70 NB: The 2nd virial coefficient of the parallelepiped is calculated based on … … 122 120 .. math:: 123 121 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 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 130 125 131 126 with … … 165 160 ---------- 166 161 167 .. [#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-854162 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-854 170 165 171 166 Authorship and Verification
Note: See TracChangeset
for help on using the changeset viewer.
