Changes in / [3fd0499:b2921d0] in sasmodels


Ignore:
Files:
5 added
7 deleted
28 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

    r650c6d2 r650c6d2  
    8484    -edit starts the parameter explorer 
    8585    -default/-demo* use demo vs default parameters 
    86     -help/-html shows the model docs instead of running the model 
     86    -html shows the model docs instead of running the model 
    8787    -title="note" adds note to the plot title, after the model name 
    8888    -data="path" uses q, dq from the data file 
    … …  
    836836    'linear', 'log', 'q4', 
    837837    'hist', 'nohist', 
    838     'edit', 'html', 'help', 
     838    'edit', 'html', 
    839839    'demo', 'default', 
    840840    ]) 
    … …  
    10051005        elif arg == '-default':    opts['use_demo'] = False 
    10061006        elif arg == '-html':    opts['html'] = True 
    1007         elif arg == '-help':    opts['html'] = True 
    10081007    # pylint: enable=bad-whitespace 
    10091008 

  • 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/barbell.py

    r0b56f38 rfcb33e4  
    8787* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
    8888""" 
    89 from numpy import inf, sin, cos, pi 
     89from numpy import inf 
    9090 
    9191name = "barbell" 
    … …  
    125125            phi_pd=15, phi_pd_n=0, 
    126126           ) 
    127 q = 0.1 
    128 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    129 qx = q*cos(pi/6.0) 
    130 qy = q*sin(pi/6.0) 
    131 tests = [[{}, 0.075, 25.5691260532], 
    132         [{'theta':80., 'phi':10.}, (qx, qy), 3.04233067789], 
    133         ] 
    134 del qx, qy  # not necessary to delete, but cleaner 

  • 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/bcc_paracrystal.py

    re2d6e3b r925ad6e  
    9999""" 
    100100 
    101 from numpy import inf, pi 
     101from numpy import inf 
    102102 
    103103name = "bcc_paracrystal" 
    … …  
    141141    psi_pd=15, psi_pd_n=0, 
    142142    ) 
    143 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    144 # add 2d test later 
    145 q =4.*pi/220. 
    146 tests = [ 
    147     [{ }, 
    148      [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]], 
    149 ] 

  • sasmodels/models/capped_cylinder.py

    r0b56f38 rfcb33e4  
    9191 
    9292""" 
    93 from numpy import inf, sin, cos, pi 
     93from numpy import inf 
    9494 
    9595name = "capped_cylinder" 
    … …  
    145145            theta_pd=15, theta_pd_n=45, 
    146146            phi_pd=15, phi_pd_n=1) 
    147 q = 0.1 
    148 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    149 qx = q*cos(pi/6.0) 
    150 qy = q*sin(pi/6.0) 
    151 tests = [[{}, 0.075, 26.0698570695], 
    152         [{'theta':80., 'phi':10.}, (qx, qy), 0.561811990502], 
    153         ] 
    154 del qx, qy  # not necessary to delete, but cleaner 

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

    r0b56f38 r8afefae  
    8888""" 
    8989 
    90 from numpy import inf, sin, cos, pi 
     90from numpy import inf, sin, cos 
    9191 
    9292name = "core_shell_bicelle" 
    … …  
    155155            theta=90, 
    156156            phi=0) 
    157 q = 0.1 
    158 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    159 qx = q*cos(pi/6.0) 
    160 qy = q*sin(pi/6.0) 
    161 tests = [[{}, 0.05, 7.4883545957], 
    162         [{'theta':80., 'phi':10.}, (qx, qy), 2.81048892474 ], 
    163         ] 
    164 del qx, qy  # not necessary to delete, but cleaner 
    165157 
     158#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 

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

    r15a90c1 r8afefae  
    8080 
    8181 
    82 .. figure:: img/elliptical_cylinder_angle_definition.png 
     82.. figure:: img/elliptical_cylinder_angle_definition.jpg 
    8383 
    84     Definition of the angles for the oriented core_shell_bicelle_elliptical particles.    
    85  
     84    Definition of the angles for the oriented core_shell_bicelle_elliptical model. 
     85    Note that *theta* and *phi* are currently defined differently to those for the core_shell_bicelle model. 
    8686 
    8787 
    … …  
    9999""" 
    100100 
    101 from numpy import inf, sin, cos, pi 
     101from numpy import inf, sin, cos 
    102102 
    103103name = "core_shell_bicelle_elliptical" 
    … …  
    150150            psi=0) 
    151151 
    152 q = 0.1 
    153 # april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 
    154 qx = q*cos(pi/6.0) 
    155 qy = q*sin(pi/6.0) 
     152#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 
    156153 
    157154tests = [ 
    … …  
    162159    'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0}, 
    163160    0.015, 286.540286], 
    164 #    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
    165         ] 
    166  
    167 del qx, qy  # not necessary to delete, but cleaner 
     161] 

  • sasmodels/models/core_shell_cylinder.py

    r0b56f38 r8e68ea0  
    7373""" 
    7474 
    75 from numpy import pi, inf, sin, cos 
     75from numpy import pi, inf 
    7676 
    7777name = "core_shell_cylinder" 
    … …  
    151151            theta_pd=15, theta_pd_n=45, 
    152152            phi_pd=15, phi_pd_n=1) 
    153 q = 0.1 
    154 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    155 qx = q*cos(pi/6.0) 
    156 qy = q*sin(pi/6.0) 
    157 tests = [[{}, 0.075, 10.8552692237], 
    158         [{}, (qx, qy), 0.444618752741 ], 
    159         ] 
    160 del qx, qy  # not necessary to delete, but cleaner 
     153 

  • 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

    r1916c52 rcb0dc22  
    1010 
    1111.. note:: 
    12    This model was originally ported from NIST IGOR macros. However, it is not 
    13    yet fully understood by the SasView developers and is currently under review. 
     12   This model was originally ported from NIST IGOR macros. However,t is not 
     13   yet fully understood by the SasView developers and is currently review. 
    1414 
    1515The form factor is normalized by the particle volume $V$ such that 
    … …  
    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(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta)}{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta} 
    54     - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right] 
    55     \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right] 
    56     \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right] 
    5750 
    5851.. note:: 
    5952 
    60     Why does t_B not appear in the above equation? 
    6153    For the calculation of the form factor to be valid, the sides of the solid 
    62     MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
     54    MUST be chosen such that** $A < B < C$. 
    6355    If this inequality is not satisfied, the model will not report an error, 
    6456    but the calculation will not be correct and thus the result wrong. 
    6557 
    6658FITTING NOTES 
    67 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 
    6860value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    6961However, **no interparticle interference effects are included in this 
    … …  
    8173NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    8274based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    83 and length $(C+2t_C)$ values, after appropriately 
    84 sorting the three dimensions to give an oblate or prolate particle, to give an  
    85 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     75and length $(C+2t_C)$ values, and used as the effective radius 
     76for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    8677 
    8778To provide easy access to the orientation of the parallelepiped, we define the 
    … …  
    9283*x*-axis of the detector. 
    9384 
    94 .. figure:: img/parallelepiped_angle_definition.png 
     85.. figure:: img/parallelepiped_angle_definition.jpg 
    9586 
    9687    Definition of the angles for oriented core-shell parallelepipeds. 
    9788 
    98 .. figure:: img/parallelepiped_angle_projection.png 
     89.. figure:: img/parallelepiped_angle_projection.jpg 
    9990 
    10091    Examples of the angles for oriented core-shell parallelepipeds against the 
    … …  
    121112 
    122113import numpy as np 
    123 from numpy import pi, inf, sqrt, cos, sin 
     114from numpy import pi, inf, sqrt 
    124115 
    125116name = "core_shell_parallelepiped" 
    … …  
    195186            psi_pd=10, psi_pd_n=1) 
    196187 
    197 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
    198 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
     188qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    199189tests = [[{}, 0.2, 0.533149288477], 
    200190         [{}, [0.2], [0.533149288477]], 
    201          [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
    202          [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
     191         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.032102135569], 
     192         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.032102135569]], 
    203193        ] 
    204194del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/models/cylinder.py

    r3330bb4 r3330bb4  
    3838 
    3939 
    40 Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with 
    41 $sin(\alpha)=\sqrt{1-u^2}$. 
     40Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with  
     41$sin(\alpha)=\sqrt{1-u^2}$.  
    4242 
    4343The output of the 1D scattering intensity function for randomly oriented 
    … …  
    6464 
    6565    Definition of the angles for oriented cylinders. 
    66  
    67 .. figure:: img/cylinder_angle_projection.png 
    68  
    69     Examples for oriented cylinders. 
    7066 
    7167The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 
    … …  
    156152tests = [[{}, 0.2, 0.042761386790780453], 
    157153        [{}, [0.2], [0.042761386790780453]], 
    158 #  new coords 
     154#  new coords     
    159155        [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 
    160156        [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 

  • 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

    r0b56f38 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 
    122 from numpy import inf, sin, cos, pi 
     100from numpy import inf 
    123101 
    124102name = "ellipsoid" 
    … …  
    161139def ER(radius_polar, radius_equatorial): 
    162140    import numpy as np 
    163 # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     141 
    164142    ee = np.empty_like(radius_polar) 
    165143    idx = radius_polar > radius_equatorial 
    … …  
    190168            theta_pd=15, theta_pd_n=45, 
    191169            phi_pd=15, phi_pd_n=1) 
    192 q = 0.1 
    193 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    194 qx = q*cos(pi/6.0) 
    195 qy = q*sin(pi/6.0) 
    196 tests = [[{}, 0.05, 54.8525847025], 
    197         [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026 ], 
    198         ] 
    199 del qx, qy  # not necessary to delete, but cleaner 

  • 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/elliptical_cylinder.py

    r15a90c1 rfcb33e4  
    6464oriented system. 
    6565 
    66 .. figure:: img/elliptical_cylinder_angle_definition.png 
     66.. figure:: img/elliptical_cylinder_angle_definition.jpg 
    6767 
    68     Definition of angles for oriented elliptical cylinder, where axis_ratio >1, 
    69     and angle $\Psi$ is a rotation around the axis of the cylinder. 
     68    Definition of angles for 2D 
    7069 
    71 .. figure:: img/elliptical_cylinder_angle_projection.png 
     70.. figure:: img/cylinder_angle_projection.jpg 
    7271 
    7372    Examples of the angles for oriented elliptical cylinders against the 
    74     detector plane, with $\Psi$ = 0. 
     73    detector plane. 
    7574 
    7675NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    … …  
    109108""" 
    110109 
    111 from numpy import pi, inf, sqrt, sin, cos 
     110from numpy import pi, inf, sqrt 
    112111 
    113112name = "elliptical_cylinder" 
    … …  
    150149                           + (length + radius) * (length + pi * radius)) 
    151150    return 0.5 * (ddd) ** (1. / 3.) 
    152 q = 0.1 
    153 # april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 
    154 qx = q*cos(pi/6.0) 
    155 qy = q*sin(pi/6.0) 
    156151 
    157152tests = [ 
    … …  
    163158      'sld_solvent':1.0, 'background':0.0}, 
    164159     0.001, 675.504402], 
    165 #    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
    166160] 

  • 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/fcc_paracrystal.py

    re2d6e3b r925ad6e  
    9090""" 
    9191 
    92 from numpy import inf, pi 
     92from numpy import inf 
    9393 
    9494name = "fcc_paracrystal" 
    … …  
    128128            psi_pd=15, psi_pd_n=0, 
    129129           ) 
    130 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    131 # add 2d test later 
    132 q =4.*pi/220. 
    133 tests = [ 
    134     [{ }, 
    135      [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]], 
    136 ] 

  • sasmodels/models/hollow_cylinder.py

    r0b56f38 raea2e2a  
    6060""" 
    6161 
    62 from numpy import pi, inf, sin, cos 
     62from numpy import pi, inf 
    6363 
    6464name = "hollow_cylinder" 
    … …  
    129129            theta_pd=10, theta_pd_n=5, 
    130130           ) 
    131 q = 0.1 
    132 # april 6 2017, rkh added a 2d unit test, assume correct! 
    133 qx = q*cos(pi/6.0) 
    134 qy = q*sin(pi/6.0) 
     131 
    135132# Parameters for unit tests 
    136133tests = [ 
    137134    [{}, 0.00005, 1764.926], 
    138135    [{}, 'VR', 1.8], 
    139     [{}, 0.001, 1756.76], 
    140     [{}, (qx, qy), 2.36885476192  ], 
    141         ] 
    142 del qx, qy  # not necessary to delete, but cleaner 
     136    [{}, 0.001, 1756.76] 
     137    ] 

  • 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/parallelepiped.py

    r3330bb4 r3330bb4  
    1515.. _parallelepiped-image: 
    1616 
    17  
    1817.. figure:: img/parallelepiped_geometry.jpg 
    1918 
    … …  
    2221.. note:: 
    2322 
    24    The edge of the solid used to have to satisfy the condition that $A < B < C$. 
    25    After some improvements to the effective radius calculation, used with an S(Q), 
    26    it is beleived that this is no longer the case. 
     23   The edge of the solid must satisfy the condition that $A < B < C$. 
     24   This requirement is not enforced in the model, so it is up to the 
     25   user to check this during the analysis. 
    2726 
    2827The 1D scattering intensity $I(q)$ is calculated as: 
    … …  
    7271 
    7372NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    74 the averaged effective radius, after appropriately 
    75 sorting the three dimensions, to give an oblate or prolate particle, $(=\sqrt{A B / \pi})$ and 
     73the averaged effective radius $(=\sqrt{A B / \pi})$ and 
    7674length $(= C)$ values, and used as the effective radius for 
    7775$S(q)$ when $P(q) \cdot S(q)$ is applied. 
    … …  
    104102.. _parallelepiped-orientation: 
    105103 
    106 .. figure:: img/parallelepiped_angle_definition.png 
    107  
    108     Definition of the angles for oriented parallelepiped, shown with $A < B < C$. 
    109  
    110 .. figure:: img/parallelepiped_angle_projection.png 
    111  
    112     Examples of the angles for an oriented parallelepiped against the 
     104.. figure:: img/parallelepiped_angle_definition.jpg 
     105 
     106    Definition of the angles for oriented parallelepipeds. 
     107 
     108.. figure:: img/parallelepiped_angle_projection.jpg 
     109 
     110    Examples of the angles for oriented parallelepipeds against the 
    113111    detector plane. 
    114112 
    … …  
    118116.. math:: 
    119117 
    120     P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2 
    121                   \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2 
    122                   \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2 
     118    P(q_x, q_y) = \left[\frac{\sin(qA\cos\alpha/2)}{(qA\cos\alpha/2)}\right]^2 
     119                  \left[\frac{\sin(qB\cos\beta/2)}{(qB\cos\beta/2)}\right]^2 
     120                  \left[\frac{\sin(qC\cos\gamma/2)}{(qC\cos\gamma/2)}\right]^2 
    123121 
    124122with 
    … …  
    156154angles. 
    157155 
     156This model is based on form factor calculations implemented in a c-library 
     157provided by the NIST Center for Neutron Research (Kline, 2006). 
    158158 
    159159References 
    … …  
    163163 
    164164R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    165  
    166 Authorship and Verification 
    167 ---------------------------- 
    168  
    169 * **Author:** This model is based on form factor calculations implemented in a c-library 
    170 provided by the NIST Center for Neutron Research (Kline, 2006). 
    171 * **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017 
    172 * **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017 
    173  
    174165""" 
    175166 
    176167import numpy as np 
    177 from numpy import pi, inf, sqrt, sin, cos 
     168from numpy import pi, inf, sqrt 
    178169 
    179170name = "parallelepiped" 
    … …  
    189180            mu = q*B 
    190181        V: Volume of the rectangular parallelepiped 
    191         alpha: angle between the long axis of the 
     182        alpha: angle between the long axis of the  
    192183            parallelepiped and the q-vector for 1D 
    193184""" 
    … …  
    217208def ER(length_a, length_b, length_c): 
    218209    """ 
    219     Return effective radius (ER) for P(q)*S(q) 
     210        Return effective radius (ER) for P(q)*S(q) 
    220211    """ 
    221     # now that axes can be in any size order, need to sort a,b,c where a~b and c is either much smaller 
    222     # or much larger 
    223     abc = np.vstack((length_a, length_b, length_c)) 
    224     abc = np.sort(abc, axis=0) 
    225     selector = (abc[1] - abc[0]) > (abc[2] - abc[1]) 
    226     length = np.where(selector, abc[0], abc[2]) 
     212 
    227213    # surface average radius (rough approximation) 
    228     radius = np.sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi) 
    229  
    230     ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius)) 
     214    surf_rad = sqrt(length_a * length_b / pi) 
     215 
     216    ddd = 0.75 * surf_rad * (2 * surf_rad * length_c + (length_c + surf_rad) * (length_c + pi * surf_rad)) 
    231217    return 0.5 * (ddd) ** (1. / 3.) 
    232218 
    … …  
    244230            phi_pd=10, phi_pd_n=1, 
    245231            psi_pd=10, psi_pd_n=10) 
    246 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
    247 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
     232 
     233qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    248234tests = [[{}, 0.2, 0.17758004974], 
    249235         [{}, [0.2], [0.17758004974]], 
    250          [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475], 
    251          [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]], 
     236         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00560296014], 
     237         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00560296014]], 
    252238        ] 
    253239del qx, qy  # not necessary to delete, but cleaner 

  • 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/stacked_disks.py

    r0b56f38 rb7e8b94  
    103103""" 
    104104 
    105 from numpy import inf, sin, cos, pi 
     105from numpy import inf 
    106106 
    107107name = "stacked_disks" 
    … …  
    152152# After redefinition of spherical coordinates - 
    153153# tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 
    154 q = 0.1 
    155 # april 6 2017, rkh added a 2d unit test, assume correct! 
    156 qx = q*cos(pi/6.0) 
    157 qy = q*sin(pi/6.0) 
     154# but should not matter here as so far all the tests are 1D not 2D 
    158155tests = [ 
    159156# Accuracy tests based on content in test/utest_extra_models.py. 
    … …  
    189186    [{'thick_core': 10.0, 
    190187      'thick_layer': 15.0, 
    191       'radius': 100.0, 
    192       'n_stacking': 5, 
    193       'sigma_d': 0.0, 
    194       'sld_core': 4.0, 
    195       'sld_layer': -0.4, 
    196       'solvent_sd': 5.0, 
    197       'theta': 90.0, 
    198       'phi': 20.0, 
    199       'scale': 0.01, 
    200       'background': 0.001}, 
    201       (qx, qy), 0.0491167089952  ], 
    202     [{'thick_core': 10.0, 
    203       'thick_layer': 15.0, 
    204188      'radius': 3000.0, 
    205189      'n_stacking': 5, 
    … …  
    244228      'background': 0.001, 
    245229     }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 
    246     [{'thick_core': 10.0, 
    247       'thick_layer': 15.0, 
    248       'radius': 3000.0, 
    249       'n_stacking': 1.0, 
    250       'sigma_d': 0.0, 
    251       'sld_core': 4.0, 
    252       'sld_layer': -0.4, 
    253       'solvent_sd': 5.0, 
    254       'theta': 90.0, 
    255       'phi': 20.0, 
    256       'scale': 0.01, 
    257       'background': 0.001, 
    258      }, (qx, qy), 0.0341738733124 ], 
    259230 
    260231    [{'thick_core': 10.0, 

  • 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

    r3330bb4 r3330bb4  
    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 for randomly oriented particles is defined by the average over all orientations $\Omega$ of: 
    19  
    20 .. math:: 
    21  
    22     P(q) = \text{scale}(\Delta\rho)^2\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, 
    6636we define the axis of the cylinder using the angles $\theta$, $\phi$ 
    67 and $\psi$. These angles are defined analogously to the elliptical_cylinder below 
    68  
    69 .. figure:: img/elliptical_cylinder_angle_definition.png 
    70  
    71     Definition of angles for oriented triaxial ellipsoid, where radii shown here are $a < b << c$ 
    72     and angle $\Psi$ is a rotation around the axis of the particle. 
    73  
     37and $\psi$. These angles are defined on 
     38:numref:`triaxial-ellipsoid-angles` . 
    7439The angle $\psi$ is the rotational angle around its own $c$ axis 
    7540against the $q$ plane. For example, $\psi = 0$ when the 
    … …  
    7843.. _triaxial-ellipsoid-angles: 
    7944 
    80 .. figure:: img/triaxial_ellipsoid_angle_projection.png 
     45.. figure:: img/triaxial_ellipsoid_angle_projection.jpg 
    8146 
    82     Some example angles for oriented ellipsoid. 
     47    The angles for oriented ellipsoid. 
    8348 
    8449The radius-of-gyration for this system is  $R_g^2 = (R_a R_b R_c)^2/5$. 
    8550 
    86 The contrast $\Delta\rho$ is defined as SLD(ellipsoid) - SLD(solvent).  In the 
     51The contrast is defined as SLD(ellipsoid) - SLD(solvent).  In the 
    8752parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major 
    8853equatorial radius, and $R_c$ is the polar radius of the ellipsoid. 
    … …  
    10469---------- 
    10570 
    106 [1] Finnigan, J.A., Jacobs, D.J., 1971. 
    107 *Light scattering by ellipsoidal particles in solution*, 
    108 J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310 
    109  
    110 Authorship and Verification 
    111 ---------------------------- 
    112  
    113 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    114 * **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017 
    115 * **Last Reviewed by:** Paul Kienzle & Richard Heenan **Date:**  April 4, 2017 
    116  
     71L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 
     72and Neutron Scattering*, Plenum, New York, 1987. 
    11773""" 
    11874 
    119 from numpy import inf, sin, cos, pi 
     75from numpy import inf 
    12076 
    12177name = "triaxial_ellipsoid" 
    … …  
    13591               "Solvent scattering length density"], 
    13692              ["radius_equat_minor", "Ang", 20, [0, inf], "volume", 
    137                "Minor equatorial radius, Ra"], 
     93               "Minor equatorial radius"], 
    13894              ["radius_equat_major", "Ang", 400, [0, inf], "volume", 
    139                "Major equatorial radius, Rb"], 
     95               "Major equatorial radius"], 
    14096              ["radius_polar", "Ang", 10, [0, inf], "volume", 
    141                "Polar radius, Rc"], 
     97               "Polar radius"], 
    14298              ["theta", "degrees", 60, [-inf, inf], "orientation", 
    14399               "In plane angle"], 
    … …  
    152108def ER(radius_equat_minor, radius_equat_major, radius_polar): 
    153109    """ 
    154     Returns the effective radius used in the S*P calculation 
     110        Returns the effective radius used in the S*P calculation 
    155111    """ 
    156112    import numpy as np 
    157113    from .ellipsoid import ER as ellipsoid_ER 
    158  
    159     # now that radii can be in any size order, radii need sorting a,b,c where a~b and c is either much smaller 
    160     # or much larger 
    161     radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar)) 
    162     radii = np.sort(radii, axis=0) 
    163     selector = (radii[1] - radii[0]) > (radii[2] - radii[1]) 
    164     polar = np.where(selector, radii[0], radii[2]) 
    165     equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2])) 
    166     return ellipsoid_ER(polar, equatorial) 
     114    return ellipsoid_ER(radius_polar, np.sqrt(radius_equat_minor * radius_equat_major)) 
    167115 
    168116demo = dict(scale=1, background=0, 
    … …  
    176124            phi_pd=15, phi_pd_n=1, 
    177125            psi_pd=15, psi_pd_n=1) 
    178  
    179 q = 0.1 
    180 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    181 # add 2d test after pull #890 
    182 qx = q*cos(pi/6.0) 
    183 qy = q*sin(pi/6.0) 
    184 tests = [[{}, 0.05, 24.8839548033], 
    185 #        [{'theta':80., 'phi':10.}, (qx, qy), 9999. ], 
    186         ] 
    187 del qx, qy  # not necessary to delete, but cleaner 
Note: See TracChangeset for help on using the changeset viewer.