Changeset 1fdb555 in sasmodels


Ignore:
Timestamp:
Oct 11, 2016 1:42:26 PM (8 years ago)
Author:
dirk
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
ef5a314
Parents:
0d6e865 (diff), 483a022 (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' of ​https://github.com/SasView/sasmodels

Location:
sasmodels
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/core_shell_bicelle.py

    r0d6e865 r1fdb555  
    4040.. math:: 
    4141 
    42     I(Q,\alpha) = \frac{\text{scale}}{V} \cdot 
     42    I(Q,\alpha) = \frac{\text{scale}}{V_t} \cdot 
    4343        F(Q,\alpha)^2 + \text{background} 
    4444 

  • sasmodels/models/multilayer_vesicle.py

    r19e7ca7 r041bc75  
    1515 
    1616See the :ref:`core-shell-sphere` model for more documentation. 
     17 
     18The 1D scattering intensity is calculated in the following way (Guinier, 1955) 
     19 
     20.. math:: 
     21 
     22    P(q) = \frac{\text{scale.volfraction}}{V_t} F^2(q) + \text{background} 
     23 
     24where 
     25 
     26.. math:: 
     27 
     28     F(q) = (\rho_{shell}-\rho_{solv}) \sum_{i=1}^{n\_pairs} \left[ 
     29     3V(R_i)\frac{\sin(qR_i)-qR_i\cos(qR_i)}{(qR_i)^3} \\ 
     30      - 3V(R_i+t_s)\frac{\sin(q(R_i+t_s))-q(R_i+t_s)\cos(q(R_i+t_s))}{(q(R_i+t_s))^3} 
     31     \right] 
     32 
     33 
     34where $R_i = r_c + (i-1)(t_s + t_w)$ 
     35    
     36where $V_t$ is the volume of the whole particle, $V(R)$ is the volume of a sphere 
     37of radius $R$, $r_c$ is the radius of the core, $\rho_{shell}$ is the scattering length  
     38density of a shell, $\rho_{solv}$ is the scattering length density of the solvent. 
     39 
    1740 
    1841The 2D scattering intensity is the same as 1D, regardless of the orientation 

  • sasmodels/sasview_model.py

    re4bf271 r64f0a1c  
    125125        if not hasattr(model, 'filename'): 
    126126            model.filename = kernel_module.__file__ 
     127            # For old models, treat .pyc and .py files interchangeably. 
     128            # This is needed because of the Sum|Multi(p1,p2) types of models 
     129            # and the convoluted way in which they are created. 
     130            if model.filename.endswith(".py"): 
     131                logging.info("Loading %s as .pyc", model.filename) 
     132                model.filename = model.filename+'c' 
    127133        if not hasattr(model, 'id'): 
    128134            model.id = splitext(basename(model.filename))[0] 

  • sasmodels/kernel_template.c

    r6a0d6aa r0d6e865  
    264264      #if defined(IQXY_HAS_THETA) 
    265265        // Force a nominal value for the spherical correction even when 
    266         // theta is +90/-90 so that there are no divide by zero problems. 
    267         // For cos(theta) fixed at 90, we effectively multiply top and bottom 
     266        // theta is +0/180 so that there are no divide by zero problems. 
     267        // For sin(theta) fixed at 0 and 180, we effectively multiply top and bottom 
    268268        // by 1e-6, so the effect cancels. 
    269269        const double spherical_correction = fmax(fabs(cos(M_PI_180*theta)), 1.e-6); 

  • sasmodels/kernelpy.py

    r14a15a3 r0d6e865  
    1010 
    1111import numpy as np  # type: ignore 
    12 from numpy import pi, cos  #type: ignore 
     12from numpy import pi, sin, cos  #type: ignore 
    1313 
    1414from . import details 
    … …  
    220220            partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
    221221            if call_details.theta_par >= 0: 
    222                 cor = cos(pi / 180 * parameters[call_details.theta_par]) 
     222                cor = sin(pi / 180 * parameters[call_details.theta_par]) 
    223223                spherical_correction = max(abs(cor), 1e-6) 
    224224            p0_index = loop_index%p0_length 

  • sasmodels/models/barbell.c

    r2222134 r0d6e865  
    9595    // Compute angle alpha between q and the cylinder axis 
    9696    double sn, cn; // slots to hold sincos function output 
    97     SINCOS(theta*M_PI_180, sn, cn); 
    98     const double q = sqrt(qx*qx+qy*qy); 
    99     const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     97    SINCOS(phi*M_PI_180, sn, cn); 
     98    const double q = sqrt(qx*qx + qy*qy); 
     99    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    100100    const double alpha = acos(cos_val); // rod angle relative to q 
    101101 

  • sasmodels/models/barbell.py

    rb0c4271 r0d6e865  
    7272    Definition of the angles for oriented 2D barbells. 
    7373 
    74 .. figure:: img/cylinder_angle_projection.jpg 
    75     :width: 600px 
    76  
    77     Examples of the angles for oriented pp against the detector plane. 
    7874 
    7975References 

  • sasmodels/models/capped_cylinder.c

    r2222134 r0d6e865  
    116116    // Compute angle alpha between q and the cylinder axis 
    117117    double sn, cn; 
    118     SINCOS(theta*M_PI_180, sn, cn); 
     118    SINCOS(phi*M_PI_180, sn, cn); 
    119119    const double q = sqrt(qx*qx+qy*qy); 
    120     const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     120    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    121121    const double alpha = acos(cos_val); // rod angle relative to q 
    122122 

  • sasmodels/models/capped_cylinder.py

    rb0c4271 r0d6e865  
    7575    Definition of the angles for oriented 2D cylinders. 
    7676 
    77 .. figure:: img/cylinder_angle_projection.jpg 
    78     :width: 600px 
    79  
    80     Examples of the angles for oriented 2D cylinders against the detector plane. 
    8177 
    8278References 
    … …  
    132128              ["radius_cap", "Ang",     20, [0, inf],    "volume", "Cap radius"], 
    133129              ["length",     "Ang",    400, [0, inf],    "volume", "Cylinder length"], 
    134               ["theta",      "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 
    135               ["phi",        "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 
     130              ["theta",      "degrees", 60, [-inf, inf], "orientation", "inclination angle"], 
     131              ["phi",        "degrees", 60, [-inf, inf], "orientation", "deflection angle"], 
    136132             ] 
    137133# pylint: enable=bad-whitespace, line-too-long 

  • sasmodels/models/core_shell_bicelle.c

    r2222134 r0d6e865  
    117117 
    118118    // Cylinder orientation 
    119     const double cyl_x = cos(theta) * cos(phi); 
    120     const double cyl_y = sin(theta); 
     119    const double cyl_x = sin(theta) * cos(phi); 
     120    const double cyl_y = sin(theta) * sin(phi); 
    121121 
    122122    // Compute the angle btw vector q and the axis of the cylinder 

  • sasmodels/models/core_shell_cylinder.c

    r43b7eea r0d6e865  
    6969 
    7070    // Compute angle alpha between q and the cylinder axis 
    71     SINCOS(theta*M_PI_180, sn, cn); 
     71    SINCOS(phi*M_PI_180, sn, cn); 
    7272    // # The following correction factor exists in sasview, but it can't be 
    7373    // # right, so we are leaving it out for now. 
    7474    // const double correction = fabs(cn)*M_PI_2; 
    7575    const double q = sqrt(qx*qx+qy*qy); 
    76     const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     76    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    7777    const double alpha = acos(cos_val); 
    7878 

  • sasmodels/models/core_shell_ellipsoid.c

    r5031ca3 r0d6e865  
    9696 
    9797    // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis. 
    98     const double cyl_x = cos(theta) * cos(phi); 
    99     const double cyl_y = sin(theta); 
     98    const double cyl_x = sin(theta) * cos(phi); 
     99    const double cyl_y = sin(theta) * sin(phi); 
    100100 
    101101    const double sldcs = core_sld - shell_sld; 

  • sasmodels/models/cylinder.c

    re9b1663d r0d6e865  
    11double form_volume(double radius, double length); 
     2double fq(double q, double sn, double cn,double radius, double length); 
     3double orient_avg_1D(double q, double radius, double length); 
    24double Iq(double q, double sld, double solvent_sld, double radius, double length); 
    35double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    … …  
    1113} 
    1214 
     15double fq(double q, double sn, double cn, double radius, double length) 
     16{ 
     17    // precompute qr and qh to save time in the loop 
     18    const double qr = q*radius; 
     19    const double qh = q*0.5*length;  
     20    return  sas_J1c(qr*sn) * sinc(qh*cn) ; 
     21} 
     22 
     23double orient_avg_1D(double q, double radius, double length) 
     24{ 
     25    // translate a point in [-1,1] to a point in [0, pi/2] 
     26    const double zm = M_PI_4; 
     27    const double zb = M_PI_4;  
     28 
     29    double total = 0.0; 
     30    for (int i=0; i<76 ;i++) { 
     31        const double alpha = Gauss76Z[i]*zm + zb; 
     32        double sn, cn; // slots to hold sincos function output 
     33        // alpha(theta,phi) the projection of the cylinder on the detector plane 
     34        SINCOS(alpha, sn, cn); 
     35        total += Gauss76Wt[i] * fq(q, sn, cn, radius, length)* fq(q, sn, cn, radius, length) * sn; 
     36    } 
     37    // translate dx in [-1,1] to dx in [lower,upper] 
     38    return total*zm; 
     39} 
     40 
    1341double Iq(double q, 
    1442    double sld, 
    … …  
    1745    double length) 
    1846{ 
    19     // precompute qr and qh to save time in the loop 
    20     const double qr = q*radius; 
    21     const double qh = q*0.5*length; 
    22  
    23     // translate a point in [-1,1] to a point in [0, pi/2] 
    24     const double zm = M_PI_4; 
    25     const double zb = M_PI_4; 
    26  
    27     double total = 0.0; 
    28     for (int i=0; i<76 ;i++) { 
    29         const double alpha = Gauss76Z[i]*zm + zb; 
    30         double sn, cn; 
    31         SINCOS(alpha, sn, cn); 
    32         const double fq = sinc(qh*cn) * sas_J1c(qr*sn); 
    33         total += Gauss76Wt[i] * fq*fq * sn; 
    34     } 
    35     // translate dx in [-1,1] to dx in [lower,upper] 
    36     const double form = total*zm; 
    3747    const double s = (sld - solvent_sld) * form_volume(radius, length); 
    38     return 1.0e-4 * s * s * form; 
     48    return 1.0e-4 * s * s * orient_avg_1D(q, radius, length); 
    3949} 
    4050 
    … …  
    5161 
    5262    // Compute angle alpha between q and the cylinder axis 
    53     SINCOS(theta*M_PI_180, sn, cn); 
     63    SINCOS(phi*M_PI_180, sn, cn); 
    5464    const double q = sqrt(qx*qx + qy*qy); 
    55     const double cos_val = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); 
     65    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
     66 
    5667    const double alpha = acos(cos_val); 
    5768 
    5869    SINCOS(alpha, sn, cn); 
    59     const double fq = sinc(q*0.5*length*cn) * sas_J1c(q*radius*sn); 
    6070    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    61     return 1.0e-4 * square(s * fq); 
     71    return 1.0e-4 * square(s * fq(q, sn, cn, radius, length)); 
    6272} 

  • sasmodels/models/cylinder.py

    r416f5c7 r0d6e865  
    3535.. math:: 
    3636 
    37     F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
     37    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\theta)\sin(\theta)d\theta} 
    3838 
    3939 
    … …  
    4848    Definition of the angles for oriented cylinders. 
    4949 
    50 .. figure:: img/cylinder_angle_projection.jpg 
    51  
    52     Examples of the angles for oriented cylinders against the detector plane. 
    5350 
    5451NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    … …  
    7774 
    7875    P(q) = \int_0^{\pi/2} d\phi 
    79         \int_0^\pi p(\theta, \phi) P_0(q,\alpha) \sin \theta\ d\theta 
     76        \int_0^\pi p(\alpha) P_0(q,\alpha) \sin \alpha\ d\alpha 
    8077 
    8178 
    82 where $p(\theta,\phi)$ is the probability distribution for the orientation 
     79where $p(\theta,\phi) = 1$ is the probability distribution for the orientation 
    8380and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented 
    8481system, and then comparing to the 1D result. 
    … …  
    8784---------- 
    8885 
    89 None 
    90  
     86J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997). 
     87G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951). 
    9188""" 
    9289 
    … …  
    123120               "Cylinder length"], 
    124121              ["theta", "degrees", 60, [-inf, inf], "orientation", 
    125                "In plane angle"], 
     122               "latitude"], 
    126123              ["phi", "degrees", 60, [-inf, inf], "orientation", 
    127                "Out of plane angle"], 
     124               "longitude"], 
    128125             ] 
    129126 
    130 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
     127source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",  "cylinder.c"] 
    131128 
    132129def ER(radius, length): 
    … …  
    148145 
    149146qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    150 tests = [[{}, 0.2, 0.042761386790780453], 
    151          [{}, [0.2], [0.042761386790780453]], 
    152          [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
    153          [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
    154         ] 
     147# After redefinition of angles, find new tests values  
     148#tests = [[{}, 0.2, 0.042761386790780453], 
     149#         [{}, [0.2], [0.042761386790780453]], 
     150#         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
     151#         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
     152#        ] 
    155153del qx, qy  # not necessary to delete, but cleaner 
    156154# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc 

  • sasmodels/models/ellipsoid.c

    ra807206 r0d6e865  
    5252 
    5353    const double q = sqrt(qx*qx + qy*qy); 
    54     SINCOS(theta*M_PI_180, sn, cn); 
    55     // TODO: check if this is actually sin(alpha), not cos(alpha) 
    56     const double cos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    57     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
     54    SINCOS(phi*M_PI_180, sn, cn); 
     55    const double cos_alpha = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
     56    const double alpha = acos(cos_alpha); 
     57    SINCOS(alpha, sn, cn); 
     58    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sn); 
    5859    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5960 

  • sasmodels/models/ellipsoid.py

    r416f5c7 r0d6e865  
    5757The $\theta$ and $\phi$ parameters are not used for the 1D output. 
    5858 
    59 .. _ellipsoid-geometry: 
    6059 
    61 .. figure:: img/ellipsoid_angle_projection.jpg 
    62  
    63     The angles for oriented ellipsoid, shown here as oblate, $a$ = $R_p$ and $b$ = $R_e$ 
    6460 
    6561Validation 

  • sasmodels/models/hollow_cylinder.c

    raea2e2a r0d6e865  
    5454 
    5555    // Cylinder orientation 
    56     cyl_x = cos(theta) * cos(phi); 
    57     cyl_y = sin(theta); 
     56    cyl_x = sin(theta) * cos(phi); 
     57    cyl_y = sin(theta) * sin(phi); 
    5858    //cyl_z = -cos(theta) * sin(phi); 
    5959 

  • sasmodels/models/stacked_disks.c

    ra807206 r0d6e865  
    142142 
    143143    // parallelepiped orientation 
    144     const double cyl_x = ct * cp; 
    145     const double cyl_y = st; 
     144    const double cyl_x = st * cp; 
     145    const double cyl_y = st * sp; 
    146146 
    147147    // Compute the angle btw vector q and the 

  • sasmodels/models/stacked_disks.py

    r7c57861 r0d6e865  
    7777the axis of the cylinder using two angles $\theta$ and $\varphi$. 
    7878 
    79 .. figure:: img/stacked_disks_angle_definition.jpg 
    80  
    81     Examples of the angles for oriented stacked disks against 
     79.. figure:: img/cylinder_angle_definition.jpg 
     80 
     81    Examples of the angles against 
    8282    the detector plane. 
    8383 
    84 .. figure:: img/stacked_disks_angle_projection.jpg 
    85  
    86     Examples of the angles for oriented pp against the detector plane. 
    8784 
    8885Our model uses the form factor calculations implemented in a c-library provided 
Note: See TracChangeset for help on using the changeset viewer.