Changes in / [139c528:a8bd500] in sasmodels


Ignore:
Location:
sasmodels
Files:
12 added
2 deleted
3 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/convert.py

    r34d6cab r44bd2be  
    132132        _remove_pd(oldpars, 'num_pearls', name) 
    133133        _remove_pd(oldpars, 'thick_string', name) 
     134    elif name == 'core_shell_parallelepiped': 
     135        _remove_pd(oldpars, 'rimA', name) 
     136        _remove_pd(oldpars, 'rimB', name) 
     137        _remove_pd(oldpars, 'rimC', name) 
    134138    elif name == 'rpa': 
    135139        # convert scattering lengths from femtometers to centimeters 

  • sasmodels/models/parallelepiped.c

    rc5b7d07 rdeb7ee0  
    99    double argA,argB,argC,tmp1,tmp2,tmp3; 
    1010    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    11     argA = a*ala/2.0; 
    12     argB = b*alb/2.0; 
    13     argC = c*alc/2.0; 
     11    argA = 0.5*a*ala; 
     12    argB = 0.5*b*alb; 
     13    argC = 0.5*c*alc; 
    1414    if(argA==0.0) { 
    1515        tmp1 = 1.0; 
    … …  
    3131} 
    3232 
     33 
    3334double form_volume(double a_side, double b_side, double c_side) 
    3435{ 
    3536    return a_side * b_side * c_side; 
    3637} 
     38 
    3739 
    3840double Iq(double q, 
    … …  
    4345    double c_side) 
    4446{ 
    45     double tmp1, tmp2, yyy; 
     47    double tmp1, tmp2; 
    4648     
    4749    double mu = q * b_side; 
    … …  
    5052    double a_scaled = a_side / b_side; 
    5153    double c_scaled = c_side / b_side; 
    52      
    53     // outer integral (with 76 gauss points), integration limits = 0, 1 
     54         
     55    //Order of integration 
     56    int nordi=76;                                
     57    int nordj=76; 
     58 
     59    // outer integral (with nordi gauss points), integration limits = 0, 1 
    5460    double summ = 0; //initialize integral 
    5561 
    56     for( int i=0; i<76; i++) { 
     62    for( int i=0; i<nordi; i++) { 
    5763                 
    58         // inner integral (with 76 gauss points), integration limits = 0, 1 
     64        // inner integral (with nordj gauss points), integration limits = 0, 1 
    5965         
    6066        double summj = 0.0; 
    61         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );              
     67            double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );          
    6268                 
    63         for(int j=0; j<76; j++) { 
     69            for(int j=0; j<nordj; j++) { 
    6470 
    6571            double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    6672            double mudum = mu * sqrt(1.0-sigma*sigma); 
    67             double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu); 
    68             double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu); 
     73                double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu); 
     74                double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu); 
    6975            if(arg1==0.0) { 
    7076                tmp1 = 1.0; 
    … …  
    7783                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 
    7884            } 
    79             yyy = Gauss76Wt[j] * tmp1 * tmp2; 
    8085 
    81             summj += yyy; 
     86            summj += Gauss76Wt[j] * tmp1 * tmp2; 
    8287        } 
    8388                 
    … …  
    9297        } 
    9398                 
    94         // sum of outer integral 
    95         yyy = Gauss76Wt[i] * answer; 
    96         summ += yyy; 
     99            // sum of outer integral 
     100        summ += Gauss76Wt[i] * answer; 
    97101         
    98102    }    
    99103    
    100104    const double vd = (sld-solvent_sld) * form_volume(a_side, b_side, c_side); 
     105     
     106    // convert from [1e-12 A-1] to [cm-1] and 0.5 factor for outer integral 
    101107    return 1.0e-4 * 0.5 * vd * vd * summ; 
    102108     

  • sasmodels/models/parallelepiped.py

    reb69cce rdeb7ee0  
    77---------- 
    88 
    9 This model provides the form factor, $P(q)$, for a rectangular parallelepiped 
    10 (below) where the form factor is normalized by the volume of the 
    11 parallelepiped. If you need to apply polydispersity, see also 
    12 rectangular_prism_. 
    13  
    14 The calculated form factor is: 
    15  
    16 .. math:: 
    17  
    18     P(q) = \frac{\text{scale}}{V} F^2(q) + \text{background} 
    19  
    20 where the volume $V = A B C$ and the averaging $\left<\ldots\right>$ is 
    21 applied over all orientations for 1D. 
     9| This model calculates the scattering from a rectangular parallelepiped (:num:`Figure #parallelepiped-image`). 
     10| If you need to apply polydispersity, see also :ref:`rectangular-prism`. 
     11 
     12.. _parallelepiped-image: 
    2213 
    2314.. figure:: img/parallelepiped.jpg 
    … …  
    2516   Parallelepiped with the corresponding definition of sides. 
    2617 
    27 The edge of the solid must satisfy the condition that $A < B < C$. 
    28 Then, assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    29 form factor is 
    30  
    31 .. math:: 
    32  
    33     P(q) = \frac{\text{scale}}{V}\int_0^1 
    34         \phi\left(\mu \sqrt{1-\sigma^2},a\right) 
     18.. note:: 
     19 
     20   The edge of the solid must satisfy the condition that $A < B < C$. 
     21   This requirement is not enforced in the model, so it is up to the 
     22   user to check this during the analysis. 
     23 
     24The 1D scattering intensity $I(q)$ is calculated as: 
     25 
     26.. Comment by Miguel Gonzalez: 
     27   I am modifying the original text because I find the notation a little bit 
     28   confusing. I think that in most textbooks/papers, the notation P(Q) is 
     29   used for the form factor (adim, P(Q=0)=1), although F(q) seems also to 
     30   be used. But here (as for many other models), P(q) is used to represent 
     31   the scattering intensity (in cm-1 normally). It would be good to agree on 
     32   a common notation. 
     33 
     34.. math:: 
     35 
     36    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 \left< P(q, \alpha) \right> 
     37    + \text{background} 
     38 
     39where the volume $V = A B C$, the contrast is defined as 
     40:math:`\Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent}`, 
     41$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 
     42at an angle $\alpha$ (angle between the long axis C and :math:`\vec q`), 
     43and the averaging $\left<\ldots\right>$ is applied over all orientations. 
     44 
     45Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
     46form factor is given by (Mittelbach and Porod, 1961) 
     47 
     48.. math:: 
     49 
     50    P(q, \alpha) = \int_0^1 \phi_Q\left(\mu \sqrt{1-\sigma^2},a\right) 
    3551        \left[S(\mu c \sigma/2)\right]^2 d\sigma 
    3652 
    … …  
    3955.. math:: 
    4056 
    41     \phi(\mu,a) = \int_0^1 
    42         \left\{S\left[\frac{\mu}{2}\cos(\frac{\pi}{2}u)\right] 
    43                S\left[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)\right] 
     57    \phi_Q(\mu,a) = \int_0^1 
     58        \left\{S\left[\frac{\mu}{2}\cos\left(\frac{\pi}{2}u\right)\right] 
     59               S\left[\frac{\mu a}{2}\sin\left(\frac{\pi}{2}u\right)\right] 
    4460               \right\}^2 du 
    4561 
    … …  
    4864    \mu = qB 
    4965 
    50 and the contrast is defined as 
    51  
    52 .. math:: 
    53  
    54     \Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent} 
    55  
    56 The scattering intensity per unit volume is returned in units of |cm^-1|; 
    57 i.e., $I(q) = \phi P(q)$. 
    58  
    59 NB: The 2nd virial coefficient of the parallelpiped is calculated based on 
     66 
     67The scattering intensity per unit volume is returned in units of |cm^-1|. 
     68 
     69NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    6070the averaged effective radius $(=\sqrt{A B / \pi})$ and 
    6171length $(= C)$ values, and used as the effective radius for 
    … …  
    6575three angles $\theta$, $\phi$ and $\Psi$. The definition of $\theta$ and 
    6676$\phi$ is the same as for the cylinder model (see also figures below). 
    67 The angle $\Psi$ is the rotational angle around the $C$ axis against 
    68 the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel 
    69 to the $x$-axis of the detector. 
    70  
     77 
     78.. Comment by Miguel Gonzalez: 
     79   The following text has been commented because I think there are two 
     80   mistakes. Psi is the rotational angle around C (but I cannot understand 
     81   what it means against the q plane) and psi=0 corresponds to a||x and b||y. 
     82 
     83   The angle $\Psi$ is the rotational angle around the $C$ axis against 
     84   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel 
     85   to the $x$-axis of the detector. 
     86 
     87The angle $\Psi$ is the rotational angle around the $C$ axis. 
     88For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis 
     89oriented parallel to the y-axis of the detector with $A$ along the z-axis. 
     90For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated 
     91$\theta$ degrees around $z$ and $\phi$ degrees around $y$, 
     92before doing a final rotation of $\Psi$ degrees around the resulting $C$ to 
     93obtain the final orientation of the parallelepiped. 
     94For example, for $\theta = 0$ and $\phi = 90$, we have that $\Psi = 0$ 
     95corresponds to $A$ along $x$ and $B$ along $y$, 
     96while for $\theta = 90$ and $\phi = 0$, $\Psi = 0$ corresponds to 
     97$A$ along $z$ and $B$ along $x$. 
    7198 
    7299.. _parallelepiped-orientation: 
    73100 
    74 .. figure:: img/orientation.jpg 
     101.. figure:: img/parallelepiped_angles_definition.jpg 
    75102 
    76103    Definition of the angles for oriented parallelepipeds. 
    77104 
    78 .. figure:: img/orientation2.jpg 
     105.. figure:: img/parallelepiped_angles_examples.jpg 
    79106 
    80107    Examples of the angles for oriented parallelepipeds against the detector plane. 
     108 
     109For a given orientation of the parallelepiped, the 2D form factor is calculated as 
     110 
     111.. math:: 
     112 
     113    P(q_x, q_y) = \left[\frac{sin(qA\cos\alpha/2)}{(qA\cos\alpha/2)}\right]^2 
     114                  \left[\frac{sin(qB\cos\beta/2)}{(qB\cos\beta/2)}\right]^2 
     115                  \left[\frac{sin(qC\cos\gamma/2)}{(qC\cos\gamma/2)}\right]^2 
     116 
     117with 
     118 
     119.. math:: 
     120 
     121    \cos\alpha = \hat A \cdot \hat q, \; 
     122    \cos\beta  = \hat B \cdot \hat q, \; 
     123    \cos\gamma = \hat C \cdot \hat q 
     124 
     125and the scattering intensity as: 
     126 
     127.. math:: 
     128 
     129    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y) + \text{background} 
     130 
     131.. Comment by Miguel Gonzalez: 
     132   This reflects the logic of the code, as in parallelepiped.c the call 
     133   to _pkernel returns P(q_x, q_y) and then this is multiplied by V^2 * (delta rho)^2. 
     134   And finally outside parallelepiped.c it will be multiplied by scale, normalized by 
     135   V and the background added. But mathematically it makes more sense to write 
     136   I(q_x, q_y) = \text{scale} V \Delta\rho^2 P(q_x, q_y) + \text{background}, 
     137   with scale being the volume fraction. 
    81138 
    82139 
    … …  
    97154   Comparison between 1D and averaged 2D. 
    98155 
    99 This model reimplements the form factor calculations implemented in a c-library 
     156This model is based on form factor calculations implemented in a c-library 
    100157provided by the NIST Center for Neutron Research (Kline, 2006). 
    101158 
    102159""" 
    103160 
     161import numpy as np 
    104162from numpy import pi, inf, sqrt 
    105163 
    … …  
    107165title = "Rectangular parallelepiped with uniform scattering length density." 
    108166description = """ 
    109      P(q)= scale/V*integral from 0 to 1 of ... 
     167    I(q)= scale*V*(sld - solvent_sld)^2*P(q,alpha)+background 
     168        P(q,alpha) = integral from 0 to 1 of ... 
    110169           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma 
    111  
     170        with 
    112171            phi(mu,a) = integral from 0 to 1 of .. 
    113         (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du 
     172            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du 
    114173            S(x) = sin(x)/x 
    115         mu = q*B 
     174            mu = q*B 
     175        V: Volume of the rectangular parallelepiped 
     176        alpha: angle between the long axis of the  
     177            parallelepiped and the q-vector for 1D 
    116178""" 
    117 category = "shape:parallelpiped" 
     179category = "shape:parallelepiped" 
    118180 
    119181#             ["name", "units", default, [lower, upper], "type","description"], 
    … …  
    139201 
    140202def ER(a_side, b_side, c_side): 
     203    """ 
     204        Return effective radius (ER) for P(q)*S(q) 
     205    """ 
    141206 
    142207    # surface average radius (rough approximation) 
    143208    surf_rad = sqrt(a_side * b_side / pi) 
    144209 
    145     # DiamCyl recoded here (to check and possibly put in a library?) 
    146     a = surf_rad 
    147     b = 0.5 * c_side 
    148     t1 = a * a * b 
    149     t2 = 1.0 + (b / a) * (1.0 + a / b / 2.0) * (1.0 + pi * a / b / 2.0) 
    150     ddd = 3.0 * t1 * t2 
    151  
     210    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    152211    return 0.5 * (ddd) ** (1. / 3.) 
     212 
     213# VR defaults to 1.0 
    153214 
    154215# parameters for demo 
    … …  
    171232               sld='sldPipe', solvent_sld='sldSolv') 
    172233 
     234 
     235qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
     236tests = [[{}, 0.2, 0.17658004974], 
     237         [{}, [0.2], [0.17658004974]], 
     238         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00460296014], 
     239         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00460296014]], 
     240        ] 
     241del qx, qy  # not necessary to delete, but cleaner 
Note: See TracChangeset for help on using the changeset viewer.