Changeset c64a68e in sasmodels


Ignore:
Timestamp:
May 2, 2018 4:58:15 PM (3 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
96153e4
Parents:
b343226 (diff), 33969b6 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-896

Location:
sasmodels
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/kerneldll.py

    rbf94e6e r33969b6  
    123123    # add openmp support if not running on a mac 
    124124    if sys.platform != "darwin": 
    125         CC.append("-fopenmp") 
     125        # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9) 
     126        # Shut it off for all unix until we can investigate. 
     127        #CC.append("-fopenmp") 
     128        pass 
    126129    def compile_command(source, output): 
    127130        """unix compiler command""" 
  • sasmodels/sasview_model.py

    r7c3fb15 rd533590  
    686686        self._intermediate_results = getattr(calculator, 'results', None) 
    687687        calculator.release() 
    688         self._model.release() 
     688        #self._model.release() 
    689689        return result 
    690690 
  • sasmodels/models/core_shell_parallelepiped.c

    re077231 rdbf1a60  
    5959 
    6060    // outer integral (with gauss points), integration limits = 0, 1 
     61    // substitute d_cos_alpha for sin_alpha d_alpha 
    6162    double outer_sum = 0; //initialize integral 
    6263    for( int i=0; i<GAUSS_N; i++) { 
    6364        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    6465        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    65  
    66         // inner integral (with gauss points), integration limits = 0, pi/2 
    6766        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
    6867        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) 
    6971        double inner_sum = 0.0; 
    7072        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 ); 
    7274            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); 
    7476            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
    7577            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
     
    9193            inner_sum += GAUSS_W[j] * f * f; 
    9294        } 
     95        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    9396        inner_sum *= 0.5; 
    9497        // now sum up the outer integral 
    9598        outer_sum += GAUSS_W[i] * inner_sum; 
    9699    } 
     100    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    97101    outer_sum *= 0.5; 
    98102 
  • sasmodels/models/core_shell_parallelepiped.py

    r97be877 r5bc6d21  
    1111.. math:: 
    1212 
    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} 
    1415 
    1516where $\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. 
     17of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be 
     18pulled out of the form factor term due to the multiple slds in the model. 
     19 
    1920The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 
    2021$A < B < C$. 
    2122 
    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 
    2328 
    2429There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
    2530(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 
     40The total volume of the solid is thus given as 
    3141 
    3242.. math:: 
    3343 
    3444    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    35  
    36 **meaning that there are "gaps" at the corners of the solid.** 
    3745 
    3846The intensity calculated follows the :ref:`parallelepiped` model, with the 
    3947core-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) 
     48amplitudes of the core and the slabs on the edges. The scattering amplitude is 
     49computed for a particular orientation of the core-shell parallelepiped with 
     50respect to the scattering vector and then averaged over all possible 
     51orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis 
     52of the parallelepiped, and $\beta$ is the angle between the projection of the 
     53particle 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 
     60and 
     61 
     62.. math:: 
     63 
     64    F(q,\alpha,\beta) 
    5265    &= (\rho_\text{core}-\rho_\text{solvent}) 
    5366       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    5467    &+ (\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) \\ 
    5669    &+ (\rho_\text{B}-\rho_\text{solvent}) 
    5770        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
     
    6376.. math:: 
    6477 
    65     S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     78    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L} 
    6679 
    6780and 
     
    6982.. math:: 
    7083 
    71     Q_A &= \sin\alpha \sin\beta \\ 
    72     Q_B &= \sin\alpha \cos\beta \\ 
    73     Q_C &= \cos\alpha 
     84    Q_A &= q \sin\alpha \sin\beta \\ 
     85    Q_B &= q \sin\alpha \cos\beta \\ 
     86    Q_C &= q \cos\alpha 
    7487 
    7588 
    7689where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
    77 are the scattering length of the parallelepiped core, and the rectangular 
     90are the scattering lengths of the parallelepiped core, and the rectangular 
    7891slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    7992is 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$. 
    80102 
    81103FITTING NOTES 
     
    94116based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    95117and 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. 
     118to give an oblate or prolate particle, to give an effective radius 
     119for $S(q)$ when $P(q) * S(q)$ is applied. 
    98120 
    99121For 2d data the orientation of the particle is required, described using 
    100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     122angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below. For further 
    101123details of the calculation and angular dispersions see :ref:`orientation`. 
    102124The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    103125$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    104126 
    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. 
     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. 
    109131 
    110132.. figure:: img/parallelepiped_angle_definition.png 
     
    113135    Note that rotation $\theta$, initially in the $xz$ plane, is carried 
    114136    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-ray 
     137    $\Psi$ is now around the axis of the particle. The neutron or X-ray 
    116138    beam is along the $z$ axis. 
    117139 
     
    120142    Examples of the angles for oriented core-shell parallelepipeds against the 
    121143    detector plane. 
     144 
     145 
     146Validation 
     147---------- 
     148 
     149Cross-checked against hollow rectangular prism and rectangular prism for equal 
     150thickness overlapping sides, and by Monte Carlo sampling of points within the 
     151shape for non-uniform, non-overlapping sides. 
     152 
    122153 
    123154References 
     
    135166 
    136167* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    137 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
     168* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016 
    138169* **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. 
    142170""" 
    143171 
  • sasmodels/models/parallelepiped.c

    r108e70e rdbf1a60  
    3838            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
     40        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    4041        inner_total *= 0.5; 
    4142 
     
    4344        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4445    } 
     46    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    4547    outer_total *= 0.5; 
    4648 
  • sasmodels/models/parallelepiped.py

    ref07e95 rb343226  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    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  
    84Definition 
    95---------- 
    106 
    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`. 
     7This model calculates the scattering from a rectangular parallelepiped 
     8(:numref:`parallelepiped-image`). 
     9If you need to apply polydispersity, see also :ref:`rectangular-prism`. For 
     10information about polarised and magnetic scattering, see 
     11the :ref:`magnetism` documentation. 
    1412 
    1513.. _parallelepiped-image: 
     
    2624error, or fixing of some dimensions at expected values, may help. 
    2725 
    28 The 1D scattering intensity $I(q)$ is calculated as: 
     26The form factor is normalized by the particle volume and the 1D scattering 
     27intensity $I(q)$ is then calculated as: 
    2928 
    3029.. Comment by Miguel Gonzalez: 
     
    3938 
    4039    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} 
    4241 
    4342where 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)$ 
     44is the form factor corresponding to a parallelepiped oriented 
     45at 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 
     47and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 
     48orientations. 
    4849 
    4950Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    50 form factor is given by (Mittelbach and Porod, 1961) 
     51form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 
    5152 
    5253.. math:: 
     
    6667    \mu &= qB 
    6768 
    68 The scattering intensity per unit volume is returned in units of |cm^-1|. 
     69where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 
     70applied. 
    6971 
    7072NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
     
    120122.. math:: 
    121123 
    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 
    125130 
    126131with 
     
    160165---------- 
    161166 
    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-854 
     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-854 
    165170 
    166171Authorship and Verification 
Note: See TracChangeset for help on using the changeset viewer.