Changes in / [b2921d0:3fd0499] in sasmodels


Ignore:
Files:
7 added
5 deleted
28 edited

Legend:

Unmodified
Added
Removed
  • explore/angular_pd.py

    r12eb36b r8267e0b  
    4747 
    4848def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    49     theta_center = radians(theta) 
     49    theta_center = radians(90-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=theta, azim=phi) 
     139            ax.view_init(elev=90-theta if use_new else 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     -html shows the model docs instead of running the model 
     86    -help/-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', 
     838    'edit', 'html', 'help', 
    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 
    10071008    # pylint: enable=bad-whitespace 
    10081009 
  • sasmodels/kernel_header.c

    r1e7b0db0 rb00a646  
    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 
    150161#if 1 
    151162//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
     
    166177#endif 
    167178 
     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 
    168199#define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
    169200    q = sqrt(qx*qx + qy*qy); \ 
     
    180211    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
    181212    } while (0) 
     213#endif 
  • sasmodels/models/barbell.py

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

    r4962519 r50beefe  
    9090    double theta, double phi, double psi) 
    9191{ 
    92     double q, cos_a1, cos_a2, cos_a3; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
     92    double q, zhat, yhat, xhat; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    9494 
    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; 
     95    const double a1 = +xhat - zhat + yhat; 
     96    const double a2 = +xhat + zhat - yhat; 
     97    const double a3 = -xhat + zhat + yhat; 
    9898 
    9999    const double qd = 0.5*q*dnn; 
  • sasmodels/models/bcc_paracrystal.py

    r925ad6e re2d6e3b  
    9999""" 
    100100 
    101 from numpy import inf 
     101from numpy import inf, pi 
    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 
     145q =4.*pi/220. 
     146tests = [ 
     147    [{ }, 
     148     [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]], 
     149] 
  • sasmodels/models/capped_cylinder.py

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

    r592343f rb260926  
    3030 
    3131static double 
    32 bicelle_kernel(double qq, 
     32bicelle_kernel(double q, 
    3333              double rad, 
    3434              double radthick, 
    3535              double facthick, 
    36               double length, 
     36              double halflength, 
    3737              double rhoc, 
    3838              double rhoh, 
     
    4242              double cos_alpha) 
    4343{ 
    44     double si1,si2,be1,be2; 
    45  
    4644    const double dr1 = rhoc-rhoh; 
    4745    const double dr2 = rhor-rhosolv; 
    4846    const double dr3 = rhoh-rhor; 
    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; 
     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); 
    5650 
    57     be1 = sas_2J1x_x(besarg1); 
    58     be2 = sas_2J1x_x(besarg2); 
    59     si1 = sas_sinx_x(sinarg1); 
    60     si2 = sas_sinx_x(sinarg2); 
     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); 
    6155 
    6256    const double t = vol1*dr1*si1*be1 + 
     
    6458                     vol3*dr3*si2*be1; 
    6559 
    66     const double retval = t*t*sin_alpha; 
     60    const double retval = t*t; 
    6761 
    6862    return retval; 
     
    7165 
    7266static double 
    73 bicelle_integration(double qq, 
     67bicelle_integration(double q, 
    7468                   double rad, 
    7569                   double radthick, 
     
    8377    // set up the integration end points 
    8478    const double uplim = M_PI_4; 
    85     const double halfheight = 0.5*length; 
     79    const double halflength = 0.5*length; 
    8680 
    8781    double summ = 0.0; 
     
    9084        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    9185        SINCOS(alpha, sin_alpha, cos_alpha); 
    92         double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 
    93                              halfheight, rhoc, rhoh, rhor, rhosolv, 
     86        double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 
     87                             halflength, rhoc, rhoh, rhor, rhosolv, 
    9488                             sin_alpha, cos_alpha); 
    95         summ += yyy; 
     89        summ += yyy*sin_alpha; 
    9690    } 
    9791 
     
    119113    double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
    120114                           0.5*length, core_sld, face_sld, rim_sld, 
    121                            solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha); 
    122  
    123     answer *= 1.0e-4; 
    124  
    125     return answer; 
     115                           solvent_sld, sin_alpha, cos_alpha); 
     116    return 1.0e-4*answer; 
    126117} 
    127118 
  • sasmodels/models/core_shell_bicelle.py

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

    r592343f rf4f85b3  
    3232} 
    3333 
    34 double  
    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) 
     34double 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) 
    4544{ 
    4645    double si1,si2,be1,be2; 
     
    7473        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    7574        double inner_sum=0; 
    76         double sinarg1 = qq*halfheight*cos_alpha; 
    77         double sinarg2 = qq*(halfheight+facthick)*cos_alpha; 
     75        double sinarg1 = q*halfheight*cos_alpha; 
     76        double sinarg2 = q*(halfheight+facthick)*cos_alpha; 
    7877        si1 = sas_sinx_x(sinarg1); 
    7978        si2 = sas_sinx_x(sinarg2); 
     
    8382            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    8483            const double rr = sqrt(rA - rB*cos(beta)); 
    85             double besarg1 = qq*rr*sin_alpha; 
    86             double besarg2 = qq*(rr+radthick)*sin_alpha; 
     84            double besarg1 = q*rr*sin_alpha; 
     85            double besarg2 = q*(rr+radthick)*sin_alpha; 
    8786            be1 = sas_2J1x_x(besarg1); 
    8887            be2 = sas_2J1x_x(besarg2); 
     
    114113{ 
    115114       // THIS NEEDS TESTING 
    116     double qq, cos_val, cos_mu, cos_nu; 
    117     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu); 
     115    double q, xhat, yhat, zhat; 
     116    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    118117    const double dr1 = rhoc-rhoh; 
    119118    const double dr2 = rhor-rhosolv; 
     
    125124    const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick); 
    126125 
    127     // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
     126    // Compute:  r = sqrt((radius_major*zhat)^2 + (radius_minor*yhat)^2) 
    128127    // Given:    radius_major = r_ratio * radius_minor   
    129128    // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle 
    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 ); 
     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 ); 
    135138    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1); 
    136139    //const double vol = form_volume(radius_minor, r_ratio, length); 
  • sasmodels/models/core_shell_bicelle_elliptical.py

    r8afefae r15a90c1  
    8080 
    8181 
    82 .. figure:: img/elliptical_cylinder_angle_definition.jpg 
     82.. figure:: img/elliptical_cylinder_angle_definition.png 
    8383 
    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. 
     84    Definition of the angles for the oriented core_shell_bicelle_elliptical particles.    
     85 
    8686 
    8787 
     
    9999""" 
    100100 
    101 from numpy import inf, sin, cos 
     101from numpy import inf, sin, cos, pi 
    102102 
    103103name = "core_shell_bicelle_elliptical" 
     
    150150            psi=0) 
    151151 
    152 #qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 
     152q = 0.1 
     153# april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 
     154qx = q*cos(pi/6.0) 
     155qy = q*sin(pi/6.0) 
    153156 
    154157tests = [ 
     
    159162    'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0}, 
    160163    0.015, 286.540286], 
    161 ] 
     164#    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
     165        ] 
     166 
     167del qx, qy  # not necessary to delete, but cleaner 
  • sasmodels/models/core_shell_cylinder.py

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

    r1e7b0db0 r92dfe0c  
    134134    double psi) 
    135135{ 
    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); 
     136    double q, zhat, yhat, xhat; 
     137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    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*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); 
     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); 
    168168     
    169169 
  • sasmodels/models/core_shell_parallelepiped.py

    rcb0dc22 r1916c52  
    1010 
    1111.. note:: 
    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. 
     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. 
    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] 
    5057 
    5158.. note:: 
    5259 
     60    Why does t_B not appear in the above equation? 
    5361    For the calculation of the form factor to be valid, the sides of the solid 
    54     MUST be chosen such that** $A < B < C$. 
     62    MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
    5563    If this inequality is not satisfied, the model will not report an error, 
    5664    but the calculation will not be correct and thus the result wrong. 
    5765 
    5866FITTING NOTES 
    59 If the scale is set equal to the particle volume fraction, |phi|, the returned 
     67If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    6068value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    6169However, **no interparticle interference effects are included in this 
     
    7381NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    7482based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    75 and length $(C+2t_C)$ values, and used as the effective radius 
    76 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     83and length $(C+2t_C)$ values, after appropriately 
     84sorting the three dimensions to give an oblate or prolate particle, to give an  
     85effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    7786 
    7887To provide easy access to the orientation of the parallelepiped, we define the 
     
    8392*x*-axis of the detector. 
    8493 
    85 .. figure:: img/parallelepiped_angle_definition.jpg 
     94.. figure:: img/parallelepiped_angle_definition.png 
    8695 
    8796    Definition of the angles for oriented core-shell parallelepipeds. 
    8897 
    89 .. figure:: img/parallelepiped_angle_projection.jpg 
     98.. figure:: img/parallelepiped_angle_projection.png 
    9099 
    91100    Examples of the angles for oriented core-shell parallelepipeds against the 
     
    112121 
    113122import numpy as np 
    114 from numpy import pi, inf, sqrt 
     123from numpy import pi, inf, sqrt, cos, sin 
    115124 
    116125name = "core_shell_parallelepiped" 
     
    186195            psi_pd=10, psi_pd_n=1) 
    187196 
    188 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
     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 
     198qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
    189199tests = [[{}, 0.2, 0.533149288477], 
    190200         [{}, [0.2], [0.533149288477]], 
    191          [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.032102135569], 
    192          [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.032102135569]], 
     201         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
     202         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
    193203        ] 
    194204del 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. 
    6670 
    6771The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 
     
    152156tests = [[{}, 0.2, 0.042761386790780453], 
    153157        [{}, [0.2], [0.042761386790780453]], 
    154 #  new coords     
     158#  new coords 
    155159        [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 
    156160        [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 
  • sasmodels/models/ellipsoid.c

    r130d4c7 r3b571ae  
    33double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
    5  
    6 static 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 } 
    275 
    286double form_volume(double radius_polar, double radius_equatorial) 
     
    3715    double radius_equatorial) 
    3816{ 
     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 
    3926    // translate a point in [-1,1] to a point in [0, 1] 
     27    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    4028    const double zm = 0.5; 
    4129    const double zb = 0.5; 
    4230    double total = 0.0; 
    4331    for (int i=0;i<76;i++) { 
    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); 
     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; 
    4736    } 
    4837    // translate dx in [-1,1] to dx in [lower,upper] 
     
    6251    double q, sin_alpha, cos_alpha; 
    6352    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    64     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, 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); 
    6556    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6657 
    67     return 1.0e-4 * form * s * s; 
     58    return 1.0e-4 * square(f * s); 
    6859} 
    6960 
  • sasmodels/models/ellipsoid.py

    r925ad6e r0b56f38  
    1818.. math:: 
    1919 
    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} 
     20    F(q,\alpha) = \Delta \rho V \frac{3(\sin qr  - qr \cos qr)}{(qr)^3} 
    2321 
    24 and 
     22for 
    2523 
    2624.. math:: 
    2725 
    28     r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha 
    29         + R_p^2 \cos^2 \alpha \right]^{1/2} 
     26    r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2} 
    3027 
    3128 
    3229$\alpha$ is the angle between the axis of the ellipsoid and $\vec q$, 
    33 $V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the 
    34 rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular 
    35 to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the 
    36 scattering length density difference between the scatterer and the solvent. 
     30$V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar 
     31radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial 
     32radius perpendicular to the rotational axis of the ellipsoid and 
     33$\Delta \rho$ (contrast) is the scattering length density difference between 
     34the scatterer and the solvent. 
    3735 
    38 For randomly oriented particles: 
     36For randomly oriented particles use the orientational average, 
    3937 
    4038.. math:: 
    4139 
    42    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
     40   \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha} 
    4341 
     42 
     43computed 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 
     49with 
     50 
     51.. math:: 
     52 
     53    r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2} 
    4454 
    4555To provide easy access to the orientation of the ellipsoid, we define 
     
    4858:ref:`cylinder orientation figure <cylinder-angle-definition>`. 
    4959For the ellipsoid, $\theta$ is the angle between the rotational axis 
    50 and the $z$ -axis. 
     60and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$ 
     61in the $xy$ plane. 
    5162 
    5263NB: The 2nd virial coefficient of the solid ellipsoid is calculated based 
     
    90101than 500. 
    91102 
     103Model was also tested against the triaxial ellipsoid model with equal major 
     104and minor equatorial radii.  It is also consistent with the cyclinder model 
     105with polar radius equal to length and equatorial radius equal to radius. 
     106 
    92107References 
    93108---------- 
     
    96111*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 
    97112Plenum Press, New York, 1987. 
     113 
     114Authorship 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 
    98120""" 
    99121 
    100 from numpy import inf 
     122from numpy import inf, sin, cos, pi 
    101123 
    102124name = "ellipsoid" 
     
    139161def ER(radius_polar, radius_equatorial): 
    140162    import numpy as np 
    141  
     163# see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
    142164    ee = np.empty_like(radius_polar) 
    143165    idx = radius_polar > radius_equatorial 
     
    168190            theta_pd=15, theta_pd_n=45, 
    169191            phi_pd=15, phi_pd_n=1) 
     192q = 0.1 
     193# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
     194qx = q*cos(pi/6.0) 
     195qy = q*sin(pi/6.0) 
     196tests = [[{}, 0.05, 54.8525847025], 
     197        [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026 ], 
     198        ] 
     199del qx, qy  # not necessary to delete, but cleaner 
  • sasmodels/models/elliptical_cylinder.c

    r592343f r61104c8  
    6767     double theta, double phi, double psi) 
    6868{ 
    69     double q, cos_val, cos_mu, cos_nu; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val, cos_mu, cos_nu); 
     69    double q, xhat, yhat, zhat; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    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*cos_nu) + cos_mu*cos_mu); 
     74    const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat)); 
    7575    const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*0.5*length*cos_val); 
     76    const double si = sas_sinx_x(q*zhat*0.5*length); 
    7777    const double Aq = be * si; 
    7878    const double delrho = sld - solvent_sld; 
  • sasmodels/models/elliptical_cylinder.py

    rfcb33e4 r15a90c1  
    6464oriented system. 
    6565 
    66 .. figure:: img/elliptical_cylinder_angle_definition.jpg 
     66.. figure:: img/elliptical_cylinder_angle_definition.png 
    6767 
    68     Definition of angles for 2D 
     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. 
    6970 
    70 .. figure:: img/cylinder_angle_projection.jpg 
     71.. figure:: img/elliptical_cylinder_angle_projection.png 
    7172 
    7273    Examples of the angles for oriented elliptical cylinders against the 
    73     detector plane. 
     74    detector plane, with $\Psi$ = 0. 
    7475 
    7576NB: The 2nd virial coefficient of the cylinder is calculated based on the 
     
    108109""" 
    109110 
    110 from numpy import pi, inf, sqrt 
     111from numpy import pi, inf, sqrt, sin, cos 
    111112 
    112113name = "elliptical_cylinder" 
     
    149150                           + (length + radius) * (length + pi * radius)) 
    150151    return 0.5 * (ddd) ** (1. / 3.) 
     152q = 0.1 
     153# april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 
     154qx = q*cos(pi/6.0) 
     155qy = q*sin(pi/6.0) 
    151156 
    152157tests = [ 
     
    158163      'sld_solvent':1.0, 'background':0.0}, 
    159164     0.001, 675.504402], 
     165#    [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 
    160166] 
  • sasmodels/models/fcc_paracrystal.c

    r4962519 r50beefe  
    9090    double theta, double phi, double psi) 
    9191{ 
    92     double q, cos_a1, cos_a2, cos_a3; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
     92    double q, zhat, yhat, xhat; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    9494 
    95     const double a1 = cos_a2 + cos_a3; 
    96     const double a2 = cos_a3 + cos_a1; 
    97     const double a3 = cos_a2 + cos_a1; 
     95    const double a1 = yhat + xhat; 
     96    const double a2 = xhat + zhat; 
     97    const double a3 = yhat + zhat; 
    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

    r925ad6e re2d6e3b  
    9090""" 
    9191 
    92 from numpy import inf 
     92from numpy import inf, pi 
    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 
     132q =4.*pi/220. 
     133tests = [ 
     134    [{ }, 
     135     [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]], 
     136] 
  • sasmodels/models/hollow_cylinder.py

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

    r1e7b0db0 rd605080  
    6767    double psi) 
    6868{ 
    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); 
     69    double q, xhat, yhat, zhat; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    7171 
    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); 
     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); 
    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 
    1718.. figure:: img/parallelepiped_geometry.jpg 
    1819 
     
    2122.. note:: 
    2223 
    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. 
     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. 
    2627 
    2728The 1D scattering intensity $I(q)$ is calculated as: 
     
    7172 
    7273NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    73 the averaged effective radius $(=\sqrt{A B / \pi})$ and 
     74the averaged effective radius, after appropriately 
     75sorting the three dimensions, to give an oblate or prolate particle, $(=\sqrt{A B / \pi})$ and 
    7476length $(= C)$ values, and used as the effective radius for 
    7577$S(q)$ when $P(q) \cdot S(q)$ is applied. 
     
    102104.. _parallelepiped-orientation: 
    103105 
    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 
     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 
    111113    detector plane. 
    112114 
     
    116118.. math:: 
    117119 
    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 
     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 
    121123 
    122124with 
     
    154156angles. 
    155157 
    156 This model is based on form factor calculations implemented in a c-library 
    157 provided 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 
     166Authorship and Verification 
     167---------------------------- 
     168 
     169* **Author:** This model is based on form factor calculations implemented in a c-library 
     170provided 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 
    165174""" 
    166175 
    167176import numpy as np 
    168 from numpy import pi, inf, sqrt 
     177from numpy import pi, inf, sqrt, sin, cos 
    169178 
    170179name = "parallelepiped" 
     
    180189            mu = q*B 
    181190        V: Volume of the rectangular parallelepiped 
    182         alpha: angle between the long axis of the  
     191        alpha: angle between the long axis of the 
    183192            parallelepiped and the q-vector for 1D 
    184193""" 
     
    208217def ER(length_a, length_b, length_c): 
    209218    """ 
    210         Return effective radius (ER) for P(q)*S(q) 
     219    Return effective radius (ER) for P(q)*S(q) 
    211220    """ 
    212  
     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]) 
    213227    # surface average radius (rough approximation) 
    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)) 
     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)) 
    217231    return 0.5 * (ddd) ** (1. / 3.) 
    218232 
     
    230244            phi_pd=10, phi_pd_n=1, 
    231245            psi_pd=10, psi_pd_n=10) 
    232  
    233 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
     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 
     247qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
    234248tests = [[{}, 0.2, 0.17758004974], 
    235249         [{}, [0.2], [0.17758004974]], 
    236          [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00560296014], 
    237          [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00560296014]], 
     250         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475], 
     251         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]], 
    238252        ] 
    239253del qx, qy  # not necessary to delete, but cleaner 
  • sasmodels/models/sc_paracrystal.c

    r4962519 r50beefe  
    111111          double psi) 
    112112{ 
    113     double q, cos_a1, cos_a2, cos_a3; 
    114     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
     113    double q, zhat, yhat, xhat; 
     114    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    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*cos_a1)/cosh_qd) 
    121                     * tanh_qd/(1. - cos(qd*cos_a2)/cosh_qd) 
    122                     * tanh_qd/(1. - cos(qd*cos_a3)/cosh_qd); 
     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); 
    123123 
    124124    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 
  • sasmodels/models/stacked_disks.py

    rb7e8b94 r0b56f38  
    103103""" 
    104104 
    105 from numpy import inf 
     105from numpy import inf, sin, cos, pi 
    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 # but should not matter here as so far all the tests are 1D not 2D 
     154q = 0.1 
     155# april 6 2017, rkh added a 2d unit test, assume correct! 
     156qx = q*cos(pi/6.0) 
     157qy = q*sin(pi/6.0) 
    155158tests = [ 
    156159# Accuracy tests based on content in test/utest_extra_models.py. 
     
    186189    [{'thick_core': 10.0, 
    187190      '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, 
    188204      'radius': 3000.0, 
    189205      'n_stacking': 5, 
     
    228244      'background': 0.001, 
    229245     }, ([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 ], 
    230259 
    231260    [{'thick_core': 10.0, 
  • sasmodels/models/triaxial_ellipsoid.c

    r925ad6e r68dd6a9  
    2020    double radius_polar) 
    2121{ 
    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; 
     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; 
    2627    double outer = 0.0; 
    2728    for (int i=0;i<76;i++) { 
    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; 
     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)); 
    3432 
    3533        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             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 ; 
     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; 
    4142        } 
    42         outer += Gauss76Wt[i] * 0.5 * inner; 
     43        outer += Gauss76Wt[i] * inner;  // correcting for dx later 
    4344    } 
    44     // translate dx in [-1,1] to dx in [lower,upper] 
    45     const double fqsq = outer*zm; 
     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); 
    4647    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    4748    return 1.0e-4 * s * s * fqsq; 
     
    5859    double psi) 
    5960{ 
    60     double q, calpha, cmu, cnu; 
    61     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu); 
     61    double q, xhat, yhat, zhat; 
     62    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    6263 
    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); 
     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); 
    6768    const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    6869 
  • sasmodels/models/triaxial_ellipsoid.py

    r3330bb4 r3330bb4  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4 All three axes are of different lengths with $R_a \leq R_b \leq R_c$ 
    5 **Users should maintain this inequality for all calculations**. 
     4Definition 
     5---------- 
     6 
     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 
     12Given an ellipsoid 
    613 
    714.. math:: 
    815 
    9     P(q) = \text{scale} V \left< F^2(q) \right> + \text{background} 
     16    \frac{X^2}{R_a^2} + \frac{Y^2}{R_b^2} + \frac{Z^2}{R_c^2} = 1 
    1017 
    11 where 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  
    18 Definition 
    19 ---------- 
    20  
    21 The form factor calculated is 
     18the scattering for randomly oriented particles is defined by the average over all orientations $\Omega$ of: 
    2219 
    2320.. math:: 
    2421 
    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 
     22    P(q) = \text{scale}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega \Phi^2(qr) d\Omega + \text{background} 
    2823 
    2924where 
     
    3126.. math:: 
    3227 
    33     \Phi(u) = 3 u^{-3} (\sin u - u \cos u) 
     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 
     32The $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 
     34axes, 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 
     42with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$. 
     43A 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 
     49for 
     50 
     51.. math:: 
     52 
     53    p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 1 
     54 
     55Due 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 
     57integral 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) 
    3464 
    3565To provide easy access to the orientation of the triaxial ellipsoid, 
    3666we define the axis of the cylinder using the angles $\theta$, $\phi$ 
    37 and $\psi$. These angles are defined on 
    38 :numref:`triaxial-ellipsoid-angles` . 
     67and $\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 
    3974The angle $\psi$ is the rotational angle around its own $c$ axis 
    4075against the $q$ plane. For example, $\psi = 0$ when the 
     
    4378.. _triaxial-ellipsoid-angles: 
    4479 
    45 .. figure:: img/triaxial_ellipsoid_angle_projection.jpg 
     80.. figure:: img/triaxial_ellipsoid_angle_projection.png 
    4681 
    47     The angles for oriented ellipsoid. 
     82    Some example angles for oriented ellipsoid. 
    4883 
    4984The radius-of-gyration for this system is  $R_g^2 = (R_a R_b R_c)^2/5$. 
    5085 
    51 The contrast is defined as SLD(ellipsoid) - SLD(solvent).  In the 
     86The contrast $\Delta\rho$ is defined as SLD(ellipsoid) - SLD(solvent).  In the 
    5287parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major 
    5388equatorial radius, and $R_c$ is the polar radius of the ellipsoid. 
     
    69104---------- 
    70105 
    71 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 
    72 and Neutron Scattering*, Plenum, New York, 1987. 
     106[1] Finnigan, J.A., Jacobs, D.J., 1971. 
     107*Light scattering by ellipsoidal particles in solution*, 
     108J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310 
     109 
     110Authorship 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 
    73117""" 
    74118 
    75 from numpy import inf 
     119from numpy import inf, sin, cos, pi 
    76120 
    77121name = "triaxial_ellipsoid" 
     
    91135               "Solvent scattering length density"], 
    92136              ["radius_equat_minor", "Ang", 20, [0, inf], "volume", 
    93                "Minor equatorial radius"], 
     137               "Minor equatorial radius, Ra"], 
    94138              ["radius_equat_major", "Ang", 400, [0, inf], "volume", 
    95                "Major equatorial radius"], 
     139               "Major equatorial radius, Rb"], 
    96140              ["radius_polar", "Ang", 10, [0, inf], "volume", 
    97                "Polar radius"], 
     141               "Polar radius, Rc"], 
    98142              ["theta", "degrees", 60, [-inf, inf], "orientation", 
    99143               "In plane angle"], 
     
    108152def ER(radius_equat_minor, radius_equat_major, radius_polar): 
    109153    """ 
    110         Returns the effective radius used in the S*P calculation 
     154    Returns the effective radius used in the S*P calculation 
    111155    """ 
    112156    import numpy as np 
    113157    from .ellipsoid import ER as ellipsoid_ER 
    114     return ellipsoid_ER(radius_polar, np.sqrt(radius_equat_minor * radius_equat_major)) 
     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) 
    115167 
    116168demo = dict(scale=1, background=0, 
     
    124176            phi_pd=15, phi_pd_n=1, 
    125177            psi_pd=15, psi_pd_n=1) 
     178 
     179q = 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 
     182qx = q*cos(pi/6.0) 
     183qy = q*sin(pi/6.0) 
     184tests = [[{}, 0.05, 24.8839548033], 
     185#        [{'theta':80., 'phi':10.}, (qx, qy), 9999. ], 
     186        ] 
     187del qx, qy  # not necessary to delete, but cleaner 
Note: See TracChangeset for help on using the changeset viewer.