Changes in / [933af72:cb0dc22] in sasmodels


Ignore:
Files:
1 deleted
16 edited

Legend:

Unmodified
Added
Removed

  • explore/angular_pd.py

    r8267e0b r12eb36b  
    4747 
    4848def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    49     theta_center = radians(90-theta) 
     49    theta_center = radians(theta) 
    5050    phi_center = radians(phi) 
    5151    flow_center = radians(flow) 
    … …  
    137137                              radius=11., dist=dist) 
    138138        if not axis.startswith('d'): 
    139             ax.view_init(elev=90-theta if use_new else theta, azim=phi) 
     139            ax.view_init(elev=theta, azim=phi) 
    140140        plt.gcf().canvas.draw() 
    141141 

  • sasmodels/compare.py

    r01ea374 rf72d70a  
    8383    -edit starts the parameter explorer 
    8484    -default/-demo* use demo vs default parameters 
    85     -help/-html shows the model docs instead of running the model 
     85    -html shows the model docs instead of running the model 
    8686    -title="note" adds note to the plot title, after the model name 
    8787 
    … …  
    829829    'linear', 'log', 'q4', 
    830830    'hist', 'nohist', 
    831     'edit', 'html', 'help', 
     831    'edit', 'html', 
    832832    'demo', 'default', 
    833833    ]) 
    … …  
    996996        elif arg == '-default':    opts['use_demo'] = False 
    997997        elif arg == '-html':    opts['html'] = True 
    998         elif arg == '-help':    opts['html'] = True 
    999998    # pylint: enable=bad-whitespace 
    1000999 

  • sasmodels/kernel_header.c

    rb00a646 r1e7b0db0  
    148148inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149149 
    150 // To rotate from the canonical position to theta, phi, psi, first rotate by 
    151 // psi about the major axis, oriented along z, which is a rotation in the 
    152 // detector plane xy. Next rotate by theta about the y axis, aligning the major 
    153 // axis in the xz plane. Finally, rotate by phi in the detector plane xy. 
    154 // To compute the scattering, undo these rotations in reverse order: 
    155 //     rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi 
    156 // The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit 
    157 // vector in the q direction. 
    158 // To change between counterclockwise and clockwise rotation, change the 
    159 // sign of phi and psi. 
    160  
    161150#if 1 
    162151//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
    … …  
    177166#endif 
    178167 
    179 #if 1 
    180 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \ 
    181     q = sqrt(qx*qx + qy*qy); \ 
    182     const double qxhat = qx/q; \ 
    183     const double qyhat = qy/q; \ 
    184     double sin_theta, cos_theta; \ 
    185     double sin_phi, cos_phi; \ 
    186     double sin_psi, cos_psi; \ 
    187     SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
    188     SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
    189     SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
    190     xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \ 
    191          + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \ 
    192     yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \ 
    193          + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \ 
    194     zhat = qxhat*(-sin_theta*cos_phi) \ 
    195          + qyhat*(-sin_theta*sin_phi); \ 
    196     } while (0) 
    197 #else 
    198 // SasView 3.x definition of orientation 
    199168#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
    200169    q = sqrt(qx*qx + qy*qy); \ 
    … …  
    211180    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
    212181    } while (0) 
    213 #endif 

  • sasmodels/models/bcc_paracrystal.c

    r50beefe r4962519  
    9090    double theta, double phi, double psi) 
    9191{ 
    92     double q, zhat, yhat, xhat; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     92    double q, cos_a1, cos_a2, cos_a3; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    9494 
    95     const double a1 = +xhat - zhat + yhat; 
    96     const double a2 = +xhat + zhat - yhat; 
    97     const double a3 = -xhat + zhat + yhat; 
     95    const double a1 = +cos_a3 - cos_a1 + cos_a2; 
     96    const double a2 = +cos_a3 + cos_a1 - cos_a2; 
     97    const double a3 = -cos_a3 + cos_a1 + cos_a2; 
    9898 
    9999    const double qd = 0.5*q*dnn; 

  • sasmodels/models/core_shell_bicelle.c

    rb260926 r592343f  
    3030 
    3131static double 
    32 bicelle_kernel(double q, 
     32bicelle_kernel(double qq, 
    3333              double rad, 
    3434              double radthick, 
    3535              double facthick, 
    36               double halflength, 
     36              double length, 
    3737              double rhoc, 
    3838              double rhoh, 
    … …  
    4242              double cos_alpha) 
    4343{ 
     44    double si1,si2,be1,be2; 
     45 
    4446    const double dr1 = rhoc-rhoh; 
    4547    const double dr2 = rhor-rhosolv; 
    4648    const double dr3 = rhoh-rhor; 
    47     const double vol1 = M_PI*square(rad)*2.0*(halflength); 
    48     const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick); 
    49     const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick); 
     49    const double vol1 = M_PI*rad*rad*(2.0*length); 
     50    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 
     51    const double vol3 = M_PI*rad*rad*2.0*(length+facthick); 
     52    double besarg1 = qq*rad*sin_alpha; 
     53    double besarg2 = qq*(rad+radthick)*sin_alpha; 
     54    double sinarg1 = qq*length*cos_alpha; 
     55    double sinarg2 = qq*(length+facthick)*cos_alpha; 
    5056 
    51     const double be1 = sas_2J1x_x(q*(rad)*sin_alpha); 
    52     const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha); 
    53     const double si1 = sas_sinx_x(q*(halflength)*cos_alpha); 
    54     const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha); 
     57    be1 = sas_2J1x_x(besarg1); 
     58    be2 = sas_2J1x_x(besarg2); 
     59    si1 = sas_sinx_x(sinarg1); 
     60    si2 = sas_sinx_x(sinarg2); 
    5561 
    5662    const double t = vol1*dr1*si1*be1 + 
    … …  
    5864                     vol3*dr3*si2*be1; 
    5965 
    60     const double retval = t*t; 
     66    const double retval = t*t*sin_alpha; 
    6167 
    6268    return retval; 
    … …  
    6571 
    6672static double 
    67 bicelle_integration(double q, 
     73bicelle_integration(double qq, 
    6874                   double rad, 
    6975                   double radthick, 
    … …  
    7783    // set up the integration end points 
    7884    const double uplim = M_PI_4; 
    79     const double halflength = 0.5*length; 
     85    const double halfheight = 0.5*length; 
    8086 
    8187    double summ = 0.0; 
    … …  
    8490        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8591        SINCOS(alpha, sin_alpha, cos_alpha); 
    86         double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 
    87                              halflength, rhoc, rhoh, rhor, rhosolv, 
     92        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 
     93                             halfheight, rhoc, rhoh, rhor, rhosolv, 
    8894                             sin_alpha, cos_alpha); 
    89         summ += yyy*sin_alpha; 
     95        summ += yyy; 
    9096    } 
    9197 
    … …  
    113119    double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
    114120                           0.5*length, core_sld, face_sld, rim_sld, 
    115                            solvent_sld, sin_alpha, cos_alpha); 
    116     return 1.0e-4*answer; 
     121                           solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha); 
     122 
     123    answer *= 1.0e-4; 
     124 
     125    return answer; 
    117126} 
    118127 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rf4f85b3 r592343f  
    3232} 
    3333 
    34 double Iq(double q, 
    35           double rad, 
    36           double x_core, 
    37           double radthick, 
    38           double facthick, 
    39           double length, 
    40           double rhoc, 
    41           double rhoh, 
    42           double rhor, 
    43           double rhosolv) 
     34double  
     35                Iq(double qq, 
     36                   double rad, 
     37                   double x_core, 
     38                   double radthick, 
     39                   double facthick, 
     40                   double length, 
     41                   double rhoc, 
     42                   double rhoh, 
     43                   double rhor, 
     44                   double rhosolv) 
    4445{ 
    4546    double si1,si2,be1,be2; 
    … …  
    7374        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    7475        double inner_sum=0; 
    75         double sinarg1 = q*halfheight*cos_alpha; 
    76         double sinarg2 = q*(halfheight+facthick)*cos_alpha; 
     76        double sinarg1 = qq*halfheight*cos_alpha; 
     77        double sinarg2 = qq*(halfheight+facthick)*cos_alpha; 
    7778        si1 = sas_sinx_x(sinarg1); 
    7879        si2 = sas_sinx_x(sinarg2); 
    … …  
    8283            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    8384            const double rr = sqrt(rA - rB*cos(beta)); 
    84             double besarg1 = q*rr*sin_alpha; 
    85             double besarg2 = q*(rr+radthick)*sin_alpha; 
     85            double besarg1 = qq*rr*sin_alpha; 
     86            double besarg2 = qq*(rr+radthick)*sin_alpha; 
    8687            be1 = sas_2J1x_x(besarg1); 
    8788            be2 = sas_2J1x_x(besarg2); 
    … …  
    113114{ 
    114115       // THIS NEEDS TESTING 
    115     double q, xhat, yhat, zhat; 
    116     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     116    double qq, cos_val, cos_mu, cos_nu; 
     117    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu); 
    117118    const double dr1 = rhoc-rhoh; 
    118119    const double dr2 = rhor-rhosolv; 
    … …  
    124125    const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick); 
    125126 
    126     // Compute:  r = sqrt((radius_major*zhat)^2 + (radius_minor*yhat)^2) 
     127    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    127128    // Given:    radius_major = r_ratio * radius_minor   
    128129    // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle 
    129     const double rad_minor = rad; 
    130     const double rad_major = rad*x_core; 
    131     const double r_hat = sqrt(square(rad_major*xhat) + square(rad_minor*yhat)); 
    132     const double rshell_hat = sqrt(square((rad_major+radthick)*xhat) 
    133                                    + square((rad_minor+radthick)*yhat)); 
    134     const double be1 = sas_2J1x_x( q*r_hat ); 
    135     const double be2 = sas_2J1x_x( q*rshell_hat ); 
    136     const double si1 = sas_sinx_x( q*halfheight*zhat ); 
    137     const double si2 = sas_sinx_x( q*(halfheight + facthick)*zhat ); 
     130    const double r = rad*sqrt(square(x_core*cos_nu) + cos_mu*cos_mu); 
     131    const double be1 = sas_2J1x_x( qq*r ); 
     132    const double be2 = sas_2J1x_x( qq*(r + radthick ) ); 
     133    const double si1 = sas_sinx_x( qq*halfheight*cos_val ); 
     134    const double si2 = sas_sinx_x( qq*(halfheight + facthick)*cos_val ); 
    138135    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1); 
    139136    //const double vol = form_volume(radius_minor, r_ratio, length); 

  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c r1e7b0db0  
    134134    double psi) 
    135135{ 
    136     double q, zhat, yhat, xhat; 
    137     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     136    double q, cos_val_a, cos_val_b, cos_val_c; 
     137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
    138138 
    139139    // cspkernel in csparallelepiped recoded here 
    … …  
    160160    double tc = length_a + 2.0*thick_rim_c; 
    161161    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162     double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    163     double siB = sas_sinx_x(0.5*q*length_b*yhat); 
    164     double siC = sas_sinx_x(0.5*q*length_c*zhat); 
    165     double siAt = sas_sinx_x(0.5*q*ta*xhat); 
    166     double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167     double siCt = sas_sinx_x(0.5*q*tc*zhat); 
     162    double siA = sas_sinx_x(0.5*q*length_a*cos_val_a); 
     163    double siB = sas_sinx_x(0.5*q*length_b*cos_val_b); 
     164    double siC = sas_sinx_x(0.5*q*length_c*cos_val_c); 
     165    double siAt = sas_sinx_x(0.5*q*ta*cos_val_a); 
     166    double siBt = sas_sinx_x(0.5*q*tb*cos_val_b); 
     167    double siCt = sas_sinx_x(0.5*q*tc*cos_val_c); 
    168168     
    169169 

  • sasmodels/models/core_shell_parallelepiped.py

    rcb0dc22 rcb0dc22  
    4848amplitudes of the core and shell, in the same manner as a core-shell model. 
    4949 
    50 .. math:: 
    51  
    52     F_{a}(Q,\alpha,\beta)= 
    53     \left[\frac{\sin(Q(L_A+2t_A)/2\sin\alpha \sin\beta)}{Q(L_A+2t_A)/2\sin\alpha\sin\beta} 
    54     - \frac{\sin(QL_A/2\sin\alpha \sin\beta)}{QL_A/2\sin\alpha \sin\beta} \right] 
    55     \left[\frac{\sin(QL_B/2\sin\alpha \sin\beta)}{QL_B/2\sin\alpha \sin\beta} \right] 
    56     \left[\frac{\sin(QL_C/2\sin\alpha \sin\beta)}{QL_C/2\sin\alpha \sin\beta} \right] 
    5750 
    5851.. note:: 
    … …  
    6457 
    6558FITTING NOTES 
    66 If the scale is set equal to the particle volume fraction, $\phi$, the returned 
     59If the scale is set equal to the particle volume fraction, |phi|, the returned 
    6760value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    6861However, **no interparticle interference effects are included in this 

  • sasmodels/models/ellipsoid.c

    r3b571ae r130d4c7  
    33double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
     5 
     6static double 
     7_ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha) 
     8{ 
     9    double ratio = radius_polar/radius_equatorial; 
     10    // Using ratio v = Rp/Re, we can expand the following to match the 
     11    // form given in Guinier (1955) 
     12    //     r = Re * sqrt(1 + cos^2(T) (v^2 - 1)) 
     13    //       = Re * sqrt( (1 - cos^2(T)) + v^2 cos^2(T) ) 
     14    //       = Re * sqrt( sin^2(T) + v^2 cos^2(T) ) 
     15    //       = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) ) 
     16    // 
     17    // Instead of using pythagoras we could pass in sin and cos; this may be 
     18    // slightly better for 2D which has already computed it, but it introduces 
     19    // an extra sqrt and square for 1-D not required by the current form, so 
     20    // leave it as is. 
     21    const double r = radius_equatorial 
     22                     * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 
     23    const double f = sas_3j1x_x(q*r); 
     24 
     25    return f*f; 
     26} 
    527 
    628double form_volume(double radius_polar, double radius_equatorial) 
    … …  
    1537    double radius_equatorial) 
    1638{ 
    17     // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955) 
    18     //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT 
    19     //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT 
    20     //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT 
    21     // u-substitution of 
    22     //     u = sin, du = cos dT 
    23     //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
    24     const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
    25  
    2639    // translate a point in [-1,1] to a point in [0, 1] 
    27     // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    2840    const double zm = 0.5; 
    2941    const double zb = 0.5; 
    3042    double total = 0.0; 
    3143    for (int i=0;i<76;i++) { 
    32         const double u = Gauss76Z[i]*zm + zb; 
    33         const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
    34         const double f = sas_3j1x_x(q*r); 
    35         total += Gauss76Wt[i] * f * f; 
     44        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
     45        const double cos_alpha = Gauss76Z[i]*zm + zb; 
     46        total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    3647    } 
    3748    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    5162    double q, sin_alpha, cos_alpha; 
    5263    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    53     const double r = sqrt(square(radius_equatorial*sin_alpha) 
    54                           + square(radius_polar*cos_alpha)); 
    55     const double f = sas_3j1x_x(q*r); 
     64    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    5665    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5766 
    58     return 1.0e-4 * square(f * s); 
     67    return 1.0e-4 * form * s * s; 
    5968} 
    6069 

  • sasmodels/models/ellipsoid.py

    r3b571ae r925ad6e  
    1818.. math:: 
    1919 
    20     F(q,\alpha) = \Delta \rho V \frac{3(\sin qr  - qr \cos qr)}{(qr)^3} 
     20    F(q,\alpha) = \frac{3 \Delta \rho V (\sin[qr(R_p,R_e,\alpha)] 
     21                - \cos[qr(R_p,R_e,\alpha)])} 
     22                {[qr(R_p,R_e,\alpha)]^3} 
    2123 
    22 for 
     24and 
    2325 
    2426.. math:: 
    2527 
    26     r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2} 
     28    r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha 
     29        + R_p^2 \cos^2 \alpha \right]^{1/2} 
    2730 
    2831 
    2932$\alpha$ is the angle between the axis of the ellipsoid and $\vec q$, 
    30 $V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar 
    31 radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial 
    32 radius perpendicular to the rotational axis of the ellipsoid and 
    33 $\Delta \rho$ (contrast) is the scattering length density difference between 
    34 the scatterer and the solvent. 
     33$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the 
     34rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular 
     35to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the 
     36scattering length density difference between the scatterer and the solvent. 
    3537 
    36 For randomly oriented particles use the orientational average, 
     38For randomly oriented particles: 
    3739 
    3840.. math:: 
    3941 
    40    \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha} 
     42   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
    4143 
    42  
    43 computed via substitution of $u=\sin(\alpha)$, $du=\cos(\alpha)\,d\alpha$ as 
    44  
    45 .. math:: 
    46  
    47     \langle F^2(q) \rangle = \int_0^1{F^2(q, u)\,du} 
    48  
    49 with 
    50  
    51 .. math:: 
    52  
    53     r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2} 
    5444 
    5545To provide easy access to the orientation of the ellipsoid, we define 
    … …  
    5848:ref:`cylinder orientation figure <cylinder-angle-definition>`. 
    5949For the ellipsoid, $\theta$ is the angle between the rotational axis 
    60 and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$ 
    61 in the $xy$ plane. 
     50and the $z$ -axis. 
    6251 
    6352NB: The 2nd virial coefficient of the solid ellipsoid is calculated based 
    … …  
    10190than 500. 
    10291 
    103 Model was also tested against the triaxial ellipsoid model with equal major 
    104 and minor equatorial radii.  It is also consistent with the cyclinder model 
    105 with polar radius equal to length and equatorial radius equal to radius. 
    106  
    10792References 
    10893---------- 
    … …  
    11196*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 
    11297Plenum Press, New York, 1987. 
    113  
    114 Authorship and Verification 
    115 ---------------------------- 
    116  
    117 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    118 * **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014 
    119 * **Last Modified by:** Paul Kienzle **Date:** March 22, 2017 
    12098""" 
    12199 

  • sasmodels/models/elliptical_cylinder.c

    r61104c8 r592343f  
    6767     double theta, double phi, double psi) 
    6868{ 
    69     double q, xhat, yhat, zhat; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     69    double q, cos_val, cos_mu, cos_nu; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val, cos_mu, cos_nu); 
    7171 
    7272    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    7373    // Given:    radius_major = r_ratio * radius_minor 
    74     const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat)); 
     74    const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu); 
    7575    const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*zhat*0.5*length); 
     76    const double si = sas_sinx_x(q*0.5*length*cos_val); 
    7777    const double Aq = be * si; 
    7878    const double delrho = sld - solvent_sld; 

  • sasmodels/models/fcc_paracrystal.c

    r50beefe r4962519  
    9090    double theta, double phi, double psi) 
    9191{ 
    92     double q, zhat, yhat, xhat; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     92    double q, cos_a1, cos_a2, cos_a3; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    9494 
    95     const double a1 = yhat + xhat; 
    96     const double a2 = xhat + zhat; 
    97     const double a3 = yhat + zhat; 
     95    const double a1 = cos_a2 + cos_a3; 
     96    const double a2 = cos_a3 + cos_a1; 
     97    const double a3 = cos_a2 + cos_a1; 
    9898    const double qd = 0.5*q*dnn; 
    9999    const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 

  • sasmodels/models/parallelepiped.c

    rd605080 r1e7b0db0  
    6767    double psi) 
    6868{ 
    69     double q, xhat, yhat, zhat; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     69    double q, cos_val_a, cos_val_b, cos_val_c; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
    7171 
    72     const double siA = sas_sinx_x(0.5*length_a*q*xhat); 
    73     const double siB = sas_sinx_x(0.5*length_b*q*yhat); 
    74     const double siC = sas_sinx_x(0.5*length_c*q*zhat); 
     72    const double siA = sas_sinx_x(0.5*q*length_a*cos_val_a); 
     73    const double siB = sas_sinx_x(0.5*q*length_b*cos_val_b); 
     74    const double siC = sas_sinx_x(0.5*q*length_c*cos_val_c); 
    7575    const double V = form_volume(length_a, length_b, length_c); 
    7676    const double drho = (sld - solvent_sld); 

  • sasmodels/models/sc_paracrystal.c

    r50beefe r4962519  
    111111          double psi) 
    112112{ 
    113     double q, zhat, yhat, xhat; 
    114     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     113    double q, cos_a1, cos_a2, cos_a3; 
     114    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    115115 
    116116    const double qd = q*dnn; 
    … …  
    118118    const double tanh_qd = tanh(arg); 
    119119    const double cosh_qd = cosh(arg); 
    120     const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd) 
    121                     * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 
    122                     * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 
     120    const double Zq = tanh_qd/(1. - cos(qd*cos_a1)/cosh_qd) 
     121                    * tanh_qd/(1. - cos(qd*cos_a2)/cosh_qd) 
     122                    * tanh_qd/(1. - cos(qd*cos_a3)/cosh_qd); 
    123123 
    124124    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 

  • sasmodels/models/triaxial_ellipsoid.c

    r68dd6a9 r925ad6e  
    2020    double radius_polar) 
    2121{ 
    22     const double pa = square(radius_equat_minor/radius_equat_major) - 1.0; 
    23     const double pc = square(radius_polar/radius_equat_major) - 1.0; 
    24     // translate a point in [-1,1] to a point in [0, pi/2] 
    25     const double zm = M_PI_4; 
    26     const double zb = M_PI_4; 
     22    double sn, cn; 
     23    // translate a point in [-1,1] to a point in [0, 1] 
     24    const double zm = 0.5; 
     25    const double zb = 0.5; 
    2726    double outer = 0.0; 
    2827    for (int i=0;i<76;i++) { 
    29         //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2; 
    30         const double phi = Gauss76Z[i]*zm + zb; 
    31         const double pa_sinsq_phi = pa*square(sin(phi)); 
     28        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
     29        const double x = 0.5*(Gauss76Z[i] + 1.0); 
     30        SINCOS(M_PI_2*x, sn, cn); 
     31        const double acosx2 = radius_equat_minor*radius_equat_minor*cn*cn; 
     32        const double bsinx2 = radius_equat_major*radius_equat_major*sn*sn; 
     33        const double c2 = radius_polar*radius_polar; 
    3234 
    3335        double inner = 0.0; 
    34         const double um = 0.5; 
    35         const double ub = 0.5; 
    3636        for (int j=0;j<76;j++) { 
    37             // translate a point in [-1,1] to a point in [0, 1] 
    38             const double usq = square(Gauss76Z[j]*um + ub); 
    39             const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 
    40             const double fq = sas_3j1x_x(q*r); 
    41             inner += Gauss76Wt[j] * fq * fq; 
     37            const double ysq = square(Gauss76Z[j]*zm + zb); 
     38            const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq); 
     39            const double fq = sas_3j1x_x(t); 
     40            inner += Gauss76Wt[j] * fq * fq ; 
    4241        } 
    43         outer += Gauss76Wt[i] * inner;  // correcting for dx later 
     42        outer += Gauss76Wt[i] * 0.5 * inner; 
    4443    } 
    45     // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    46     const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
     44    // translate dx in [-1,1] to dx in [lower,upper] 
     45    const double fqsq = outer*zm; 
    4746    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    4847    return 1.0e-4 * s * s * fqsq; 
    … …  
    5958    double psi) 
    6059{ 
    61     double q, xhat, yhat, zhat; 
    62     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     60    double q, calpha, cmu, cnu; 
     61    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu); 
    6362 
    64     const double r = sqrt(square(radius_equat_minor*xhat) 
    65                           + square(radius_equat_major*yhat) 
    66                           + square(radius_polar*zhat)); 
    67     const double fq = sas_3j1x_x(q*r); 
     63    const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu 
     64                          + radius_equat_major*radius_equat_major*cmu*cmu 
     65                          + radius_polar*radius_polar*calpha*calpha); 
     66    const double fq = sas_3j1x_x(t); 
    6867    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    6968 

  • sasmodels/models/triaxial_ellipsoid.py

    r28d3067 r925ad6e  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
     4All three axes are of different lengths with $R_a \leq R_b \leq R_c$ 
     5**Users should maintain this inequality for all calculations**. 
     6 
     7.. math:: 
     8 
     9    P(q) = \text{scale} V \left< F^2(q) \right> + \text{background} 
     10 
     11where the volume $V = 4/3 \pi R_a R_b R_c$, and the averaging 
     12$\left<\ldots\right>$ is applied over all orientations for 1D. 
     13 
     14.. figure:: img/triaxial_ellipsoid_geometry.jpg 
     15 
     16    Ellipsoid schematic. 
     17 
    418Definition 
    519---------- 
    620 
    7 .. figure:: img/triaxial_ellipsoid_geometry.jpg 
    8  
    9     Ellipsoid with $R_a$ as *radius_equat_minor*, $R_b$ as *radius_equat_major* 
    10     and $R_c$ as *radius_polar*. 
    11  
    12 Given an ellipsoid 
     21The form factor calculated is 
    1322 
    1423.. math:: 
    1524 
    16     \frac{X^2}{R_a^2} + \frac{Y^2}{R_b^2} + \frac{Z^2}{R_c^2} = 1 
    17  
    18 the scattering is defined by the average over all orientations $\Omega$, 
    19  
    20 .. math:: 
    21  
    22     P(q) = \text{scale}\frac{V}{4 \pi}\int_\Omega \Phi^2(qr) d\Omega + \text{background} 
     25    P(q) = \frac{\text{scale}}{V}\int_0^1\int_0^1 
     26        \Phi^2(qR_a^2\cos^2( \pi x/2) + qR_b^2\sin^2(\pi y/2)(1-y^2) + R_c^2y^2) 
     27        dx dy 
    2328 
    2429where 
    … …  
    2631.. math:: 
    2732 
    28     \Phi(qr) &= 3 j_1(qr)/qr = 3 (\sin qr - qr \cos qr)/(qr)^3 \\ 
    29     r^2 &= R_a^2e^2 + R_b^2f^2 + R_c^2g^2 \\ 
    30     V &= \tfrac{4}{3} \pi R_a R_b R_c 
    31  
    32 The $e$, $f$ and $g$ terms are the projections of the orientation vector on $X$, 
    33 $Y$ and $Z$ respectively.  Keeping the orientation fixed at the canonical 
    34 axes, we can integrate over the incident direction using polar angle 
    35 $-\pi/2 \le \gamma \le \pi/2$ and equatorial angle $0 \le \phi \le 2\pi$ 
    36 (as defined in ref [1]), 
    37  
    38  .. math:: 
    39  
    40      \langle\Phi^2\rangle = \int_0^{2\pi} \int_{-\pi/2}^{\pi/2} \Phi^2(qr) \cos \gamma\,d\gamma d\phi 
    41  
    42 with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$. 
    43 A little algebra yields 
    44  
    45 .. math:: 
    46  
    47     r^2 = b^2(p_a \sin^2 \phi \cos^2 \gamma + 1 + p_c \sin^2 \gamma) 
    48  
    49 for 
    50  
    51 .. math:: 
    52  
    53     p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 1 
    54  
    55 Due to symmetry, the ranges can be restricted to a single quadrant 
    56 $0 \le \gamma \le \pi/2$ and $0 \le \phi \le \pi/2$, scaling the resulting 
    57 integral by 8. The computation is done using the substitution $u = \sin\gamma$, 
    58 $du = \cos\gamma\,d\gamma$, giving 
    59  
    60 .. math:: 
    61  
    62     \langle\Phi^2\rangle &= 8 \int_0^{\pi/2} \int_0^1 \Phi^2(qr) du d\phi \\ 
    63     r^2 &= b^2(p_a \sin^2(\phi)(1 - u^2) + 1 + p_c u^2) 
     33    \Phi(u) = 3 u^{-3} (\sin u - u \cos u) 
    6434 
    6535To provide easy access to the orientation of the triaxial ellipsoid, 
    … …  
    9969---------- 
    10070 
    101 [1] Finnigan, J.A., Jacobs, D.J., 1971. 
    102 *Light scattering by ellipsoidal particles in solution*, 
    103 J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310 
    104  
     71L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 
     72and Neutron Scattering*, Plenum, New York, 1987. 
    10573""" 
    10674 
    … …  
    12391               "Solvent scattering length density"], 
    12492              ["radius_equat_minor", "Ang", 20, [0, inf], "volume", 
    125                "Minor equatorial radius, Ra"], 
     93               "Minor equatorial radius"], 
    12694              ["radius_equat_major", "Ang", 400, [0, inf], "volume", 
    127                "Major equatorial radius, Rb"], 
     95               "Major equatorial radius"], 
    12896              ["radius_polar", "Ang", 10, [0, inf], "volume", 
    129                "Polar radius, Rc"], 
     97               "Polar radius"], 
    13098              ["theta", "degrees", 60, [-inf, inf], "orientation", 
    13199               "In plane angle"], 
Note: See TracChangeset for help on using the changeset viewer.