Changeset 4493288 in sasmodels


Ignore:
Timestamp:
Nov 29, 2017 11:25:39 PM (7 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:
10ee838
Parents:
0e55afe
Message:

core_shell_parallelepiped: simplify the calculation and update the docs

Location:
sasmodels/models
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/core_shell_parallelepiped.c

    rfc0b7aa r4493288  
    4343    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
    4444    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    45     //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     45    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     46    // Code rewritten (PAK) 
    4647 
    47     const double mu = 0.5 * q * length_b; 
     48    const double half_q = 0.5*q; 
    4849 
    49     // Scale sides by B 
    50     const double a_over_b = length_a / length_b; 
    51     const double c_over_b = length_c / length_b; 
     50    const double tA = length_a + 2.0*thick_rim_a; 
     51    const double tB = length_b + 2.0*thick_rim_b; 
     52    const double tC = length_c + 2.0*thick_rim_c; 
    5253 
    53     double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b; 
    54     double tB_over_b = 1+ 2.0*thick_rim_b/length_b; 
    55     double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b; 
    56  
    57     double Vin = length_a * length_b * length_c; 
    58 #if OVERLAPPING 
    59     const double capA_area = length_b*length_c; 
    60     const double capB_area = (length_a+2.*thick_rim_a)*length_c; 
    61     const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b); 
    62 #else 
    63     const double capA_area = length_b*length_c; 
    64     const double capB_area = length_a*length_c; 
    65     const double capC_area = length_a*length_b; 
    66 #endif 
    67     const double Va = length_a * capA_area; 
    68     const double Vb = length_b * capB_area; 
    69     const double Vc = length_c * capC_area; 
    70     const double Vat = Va + 2.0 * thick_rim_a * capA_area; 
    71     const double Vbt = Vb + 2.0 * thick_rim_b * capB_area; 
    72     const double Vct = Vc + 2.0 * thick_rim_c * capC_area; 
    73  
    74     // Scale factors (note that drC is not used later) 
     54    // Scale factors 
    7555    const double dr0 = (core_sld-solvent_sld); 
    7656    const double drA = (arim_sld-solvent_sld); 
    … …  
    8161    double outer_sum = 0; //initialize integral 
    8262    for( int i=0; i<76; i++) { 
    83         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    84         double mu_proj = mu * sqrt(1.0-sigma*sigma); 
     63        const double cos_alpha = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     64        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    8565 
    8666        // inner integral (with gauss points), integration limits = 0, pi/2 
    87         const double siC = sas_sinx_x(mu * sigma * c_over_b); 
    88         const double siCt = sas_sinx_x(mu * sigma * tC_over_b); 
     67        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
     68        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
    8969        double inner_sum = 0.0; 
    9070        for(int j=0; j<76; j++) { 
    91             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    92             double sin_uu, cos_uu; 
    93             SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    94             const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b); 
    95             const double siB = sas_sinx_x(mu_proj * cos_uu ); 
    96             const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b); 
    97             const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b); 
     71            const double beta = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     72            double sin_beta, cos_beta; 
     73            SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     74            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
     75            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
     76            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta); 
     77            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta); 
    9878 
    9979#if OVERLAPPING 
    100             const double f = dr0*Vin*siA*siB*siC 
    101                 + drA*(Vat*siAt-Va*siA)*siB*siC 
    102                 + drB*siAt*(Vbt*siBt-Vb*siB)*siC 
    103                 + drC*siAt*siBt*(Vct*siCt-Vc*siC); 
     80            const double f = dr0*siA*siB*siC 
     81                + drA*(siAt-siA)*siB*siC 
     82                + drB*siAt*(siBt-siB)*siC 
     83                + drC*siAt*siBt*(siCt-siC); 
    10484#else 
    105             const double f = dr0*Vin*siA*siB*siC 
    106                 + drA*(Vat*siAt-Va*siA)*siB*siC 
    107                 + drB*siA*(Vbt*siBt-Vb*siB)*siC 
    108                 + drC*siA*siB*(Vct*siCt-Vc*siC); 
     85            const double f = dr0*siA*siB*siC 
     86                + drA*(siAt-siA)*siB*siC 
     87                + drB*siA*(siBt-siB)*siC 
     88                + drC*siA*siB*(siCt-siC); 
    10989#endif 
    11090 
    … …  
    141121    const double drC = crim_sld-solvent_sld; 
    142122 
    143     double Vin = length_a * length_b * length_c; 
    144 #if OVERLAPPING 
    145     const double capA_area = length_b*length_c; 
    146     const double capB_area = (length_a+2.*thick_rim_a)*length_c; 
    147     const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b); 
    148 #else 
    149     const double capA_area = length_b*length_c; 
    150     const double capB_area = length_a*length_c; 
    151     const double capC_area = length_a*length_b; 
    152 #endif 
    153     const double Va = length_a * capA_area; 
    154     const double Vb = length_b * capB_area; 
    155     const double Vc = length_c * capC_area; 
    156     const double Vat = Va + 2.0 * thick_rim_a * capA_area; 
    157     const double Vbt = Vb + 2.0 * thick_rim_b * capB_area; 
    158     const double Vct = Vc + 2.0 * thick_rim_c * capC_area; 
    159  
    160123    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    161124    // the scaling by B. 
    … …  
    163126    const double tB = length_b + 2.0*thick_rim_b; 
    164127    const double tC = length_c + 2.0*thick_rim_c; 
    165     const double siA = sas_sinx_x(0.5*length_a*qa); 
    166     const double siB = sas_sinx_x(0.5*length_b*qb); 
    167     const double siC = sas_sinx_x(0.5*length_c*qc); 
    168     const double siAt = sas_sinx_x(0.5*tA*qa); 
    169     const double siBt = sas_sinx_x(0.5*tB*qb); 
    170     const double siCt = sas_sinx_x(0.5*tC*qc); 
     128    const double siA = length_a*sas_sinx_x(0.5*length_a*qa); 
     129    const double siB = length_b*sas_sinx_x(0.5*length_b*qb); 
     130    const double siC = length_c*sas_sinx_x(0.5*length_c*qc); 
     131    const double siAt = tA*sas_sinx_x(0.5*tA*qa); 
     132    const double siBt = tB*sas_sinx_x(0.5*tB*qb); 
     133    const double siCt = tC*sas_sinx_x(0.5*tC*qc); 
    171134 
    172135#if OVERLAPPING 
    173     const double f = dr0*Vin*siA*siB*siC 
    174         + drA*(Vat*siAt-Va*siA)*siB*siC 
    175         + drB*siAt*(Vbt*siBt-Vb*siB)*siC 
    176         + drC*siAt*siBt*(Vct*siCt-Vc*siC); 
     136    const double f = dr0*siA*siB*siC 
     137        + drA*(siAt-siA)*siB*siC 
     138        + drB*siAt*(siBt-siB)*siC 
     139        + drC*siAt*siBt*(siCt-siC); 
    177140#else 
    178     const double f = dr0*Vin*siA*siB*siC 
    179         + drA*(Vat*siAt-Va*siA)*siB*siC 
    180         + drB*siA*(Vbt*siBt-Vb*siB)*siC 
    181         + drC*siA*siB*(Vct*siCt-Vc*siC); 
     141    const double f = dr0*siA*siB*siC 
     142        + drA*(siAt-siA)*siB*siC 
     143        + drB*siA*(siBt-siB)*siC 
     144        + drC*siA*siB*(siCt-siC); 
    182145#endif 
    183146 

  • sasmodels/models/core_shell_parallelepiped.py

    r0e55afe r4493288  
    4040amplitudes of the core and the slabs on the edges. 
    4141 
    42 the scattering amplitude is computed for a particular orientation of the core-shell 
    43 parallelepiped with respect to the scattering vector and then averaged over all 
    44 possible orientations, where $\alpha$ is the angle between the $z$ axis and the longest axis $C$ 
    45 of the parallelepiped, $\beta$ is the angle between projection of the particle in the $xy$ detector plane and the $y$ axis. 
    46  
    47 .. math:: 
    48     \begin{align*} 
    49     F(Q)&=A B C (\rho_\text{core}-\rho_\text{solvent})  S(A \sin\alpha \sin\beta)S(B \sin\alpha \cos\beta)S(C \cos\alpha) \\ 
    50     &+ 2t_A B C (\rho_\text{A}-\rho_\text{solvent})  \left[S((A+t_A) \sin\alpha \sin\beta)-S(A \sin\alpha \sin\beta)\right] S(B \sin\alpha \cos\beta) S(C \cos\alpha)\\ 
    51     &+ 2 A t_B C (\rho_\text{B}-\rho_\text{solvent})  S(A \sin\alpha \sin\beta) \left[S((B+t_B) \sin\alpha \cos\beta)-S(B \sin\alpha \cos\beta)\right] S(C \cos\alpha)\\ 
    52     &+ 2 A B t_C (\rho_\text{C}-\rho_\text{solvent}) S(A \sin\alpha \sin\beta) S(B \sin\alpha \cos\beta) \left[S((C+t_C) \cos\alpha)-S(C \cos\alpha)\right] 
    53     \end{align*} 
     42the scattering amplitude is computed for a particular orientation of the 
     43core-shell parallelepiped with respect to the scattering vector and then 
     44averaged over all possible orientations, where $\alpha$ is the angle between 
     45the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
     46the angle between projection of the particle in the $xy$ detector plane 
     47and the $y$ axis. 
     48 
     49.. math:: 
     50 
     51    F(Q) 
     52    &= (\rho_\text{core}-\rho_\text{solvent}) 
     53       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
     54    &+ (\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) \\ 
     56    &+ (\rho_\text{B}-\rho_\text{solvent}) 
     57        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
     58    &+ (\rho_\text{C}-\rho_\text{solvent}) 
     59        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right] 
    5460 
    5561with 
    … …  
    5763.. math:: 
    5864 
    59     S(x) = \frac{\sin \tfrac{1}{2}Q x}{\tfrac{1}{2}Q x} 
    60  
    61 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ are 
    62 the scattering length of the parallelepiped core, and the rectangular slabs of 
    63 thickness $t_A$, $t_B$ and $t_C$, respectively. 
    64 $\rho_\text{solvent}$ is the scattering length of the solvent. 
     65    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     66 
     67and 
     68 
     69.. math:: 
     70 
     71    Q_A &= \sin\alpha \sin\beta \\ 
     72    Q_B &= \sin\alpha \cos\beta \\ 
     73    Q_C &= \cos\alpha 
     74 
     75 
     76where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
     77are the scattering length of the parallelepiped core, and the rectangular 
     78slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
     79is the scattering length of the solvent. 
    6580 
    6681FITTING NOTES 
     82~~~~~~~~~~~~~ 
     83 
    6784If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    68 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    69 However, **no interparticle interference effects are included in this 
    70 calculation.** 
     85value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
     86**no interparticle interference effects are included in this calculation.** 
    7187 
    7288There are many parameters in this model. Hold as many fixed as possible with 
    … …  
    7793NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    7894based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    79 and length $(C+2t_C)$ values, after appropriately 
    80 sorting the three dimensions to give an oblate or prolate particle, to give an 
    81 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     95and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
     96to give an oblate or prolate particle, to give an effective radius, 
     97for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    8298 
    8399For 2d data the orientation of the particle is required, described using 
    84 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    85 of the calculation and angular dispersions see :ref:`orientation` . 
     100angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     101details of the calculation and angular dispersions see :ref:`orientation`. 
    86102The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    87103$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    88104 
    89 For 2d, constraints must be applied during fitting to ensure that the inequality 
    90 $A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 
     105For 2d, constraints must be applied during fitting to ensure that the 
     106inequality $A < B < C$ is not violated, and hence the correct definition 
     107of angles is preserved. The calculation will not report an error, 
    91108but the results may be not correct. 
    92109 
    … …  
    173190 
    174191    # surface average radius (rough approximation) 
    175     surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi) 
     192    surf_rad = sqrt((length_a + 2.0*thick_rim_a) 
     193                    * (length_b + 2.0*thick_rim_b) / pi) 
    176194 
    177195    height = length_c + 2.0*thick_rim_c 
    178196 
    179     ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad)) 
     197    ddd = (0.75 * surf_rad * (2 * surf_rad * height 
     198           + (height + surf_rad) * (height + pi * surf_rad))) 
    180199    return 0.5 * (ddd) ** (1. / 3.) 
    181200 
    … …  
    212231            psi_pd=10, psi_pd_n=1) 
    213232 
    214 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
     233# rkh 7/4/17 add random unit test for 2d, note make all params different, 
     234# 2d values not tested against other codes or models 
    215235if 0:  # pak: model rewrite; need to update tests 
    216236    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
Note: See TracChangeset for help on using the changeset viewer.