Changeset dbf1a60 in sasmodels


Ignore:
Timestamp:
Mar 11, 2018 12:29:41 PM (6 years ago)
Author:
butler
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
9616dfe
Parents:
367886f
Message:

Add comments to c code and clean documentatin

Added comments to c code in both parallelepiped and core shell
parallelepiped noting the change of integration varialbes in the
computation. Cleaned up and final corrections to the core shell
documentation and did some cleaning of parallelipiped. In particular
tried to bring a bit more consistency between the docs.

addresses #896

Location:
sasmodels/models
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • 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

    r367886f rdbf1a60  
    1111.. math:: 
    1212 
    13     I(q) = \text{scale}\frac{\langle P(q,\alpha,\beta) \rangle}{V}  
     13    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle  
    1414    + \text{background} 
    1515 
    1616where $\langle \ldots \rangle$ is an average over all possible orientations 
    17 of the rectangular solid. 
    18  
    19 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 
    2020The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 
    2121$A < B < C$. 
    … …  
    2929There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
    3030(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 then 
     31$(=t_C)$ faces. The projection in the $AB$ plane is 
    3232 
    3333.. figure:: img/core_shell_parallelepiped_projection.jpg 
    3434 
    3535   AB cut through the core-shell parllelipiped showing the cross secion of 
    36    four of the six shell slabs 
    37  
    38 The volume of the solid is 
     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 
    3941 
    4042.. math:: 
    4143 
    4244    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    43  
    44 **meaning that there are "gaps" at the corners of the solid.** 
    4545 
    4646The intensity calculated follows the :ref:`parallelepiped` model, with the 
    4747core-shell intensity being calculated as the square of the sum of the 
    48 amplitudes of the core and the slabs on the edges. 
    49  
    50 the scattering amplitude is computed for a particular orientation of the 
    51 core-shell parallelepiped with respect to the scattering vector and then 
    52 averaged over all possible orientations, where $\alpha$ is the angle between 
    53 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
    54 the angle between projection of the particle in the $xy$ detector plane 
    55 and the $y$ axis. 
    56  
    57 .. math:: 
    58  
    59     P(q)=\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ cos\alpha 
    60     \ d\alpha \ d\beta 
     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} 
    6159 
    6260and 
    … …  
    6462.. math:: 
    6563 
    66     F(q) 
     64    F(q,\alpha,\beta) 
    6765    &= (\rho_\text{core}-\rho_\text{solvent}) 
    6866       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    … …  
    7876.. math:: 
    7977 
    80     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} 
    8179 
    8280and 
    … …  
    9391slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    9492is 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$. 
    95102 
    96103FITTING NOTES 
    … …  
    118125$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    119126 
    120 For 2d, constraints must be applied during fitting to ensure that the 
    121 inequality $A < B < C$ is not violated, and hence the correct definition 
    122 of angles is preserved. The calculation will not report an error, 
    123 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. 
    124131 
    125132.. figure:: img/parallelepiped_angle_definition.png 

  • 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

    r5bc373b rdbf1a60  
    3939 
    4040    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 
    41            \left< P(q, \alpha) \right> + \text{background} 
     41           \left< P(q, \alpha, \beta) \right> + \text{background} 
    4242 
    4343where 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. 
     44$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ 
     45is the form factor corresponding to a parallelepiped oriented 
     46at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$ 
     47( the angle between the projection of the particle in the $xy$ detector plane 
     48and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 
     49orientations. 
    4850 
    4951Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    50 form factor is given by (Mittelbach and Porod, 1961) 
     52form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 
    5153 
    5254.. math:: 
    … …  
    6668    \mu &= qB 
    6769 
    68 The scattering intensity per unit volume is returned in units of |cm^-1|. 
     70where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 
     71applied. 
    6972 
    7073NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    … …  
    120123.. math:: 
    121124 
    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 
     125    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1} 
     126                   {2}qA\cos\alpha)}\right]^2 
     127                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1} 
     128                   {2}qB\cos\beta)}\right]^2 
     129                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1} 
     130                   {2}qC\cos\gamma)}\right]^2 
    125131 
    126132with 
    … …  
    160166---------- 
    161167 
    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 
     168.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*, 
     169   14 (1961) 185-211 
     170.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    165171 
    166172Authorship and Verification 
Note: See TracChangeset for help on using the changeset viewer.