Changeset 0d6e865 in sasmodels


Ignore:
Timestamp:
Oct 11, 2016 11:42:00 AM (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:
1fdb555
Parents:
19e7ca7
Message:

Rewriting selected models in spherical coordinates. This fixes ticket #492 for these models.

Location:
sasmodels
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • 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_bicelle.py

    r416f5c7 r0d6e865  
    7171    Definition of the angles for the oriented core shell bicelle tmodel. 
    7272 
    73 .. figure:: img/cylinder_angle_projection.jpg 
    74     :width: 600px 
    75  
    76     Examples of the angles for oriented pp against the detector plane. 
    7773 
    7874References 
    … …  
    161157 
    162158qx, qy = 0.4 * cos(90), 0.5 * sin(0) 
    163 tests = [ 
    164     # Accuracy tests based on content in test/utest_other_models.py 
    165     [{'radius': 20.0, 
    166       'thick_rim': 10.0, 
    167       'thick_face': 10.0, 
    168       'length': 400.0, 
    169       'sld_core': 1.0, 
    170       'sld_face': 4.0, 
    171       'sld_rim': 4.0, 
    172       'sld_solvent': 1.0, 
    173       'background': 0.0, 
    174      }, 0.001, 353.550], 
    175  
    176     [{'radius': 20.0, 
    177       'thick_rim': 10.0, 
    178       'thick_face': 10.0, 
    179       'length': 400.0, 
    180       'sld_core': 1.0, 
    181       'sld_face': 4.0, 
    182       'sld_rim': 4.0, 
    183       'sld_solvent': 1.0, 
    184       'theta': 90.0, 
    185       'phi': 0.0, 
    186       'background': 0.00, 
    187      }, (qx, qy), 24.9167], 
    188  
    189     # Additional tests with larger range of parameters 
    190     [{'radius': 3.0, 
    191       'thick_rim': 100.0, 
    192       'thick_face': 100.0, 
    193       'length': 1200.0, 
    194       'sld_core': 5.0, 
    195       'sld_face': 41.0, 
    196       'sld_rim': 42.0, 
    197       'sld_solvent': 21.0, 
    198      }, 0.05, 1670.1828], 
    199     ] 

  • 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.