Changeset 6e5c0b7 in sasmodels


Ignore:
Timestamp:
Apr 4, 2017 3:28:58 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
3a45c2c
Parents:
933af72 (diff), c4e3215 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-890

Files:
1 added
21 edited

Legend:

Unmodified
Added
Removed
  • LICENSE.txt

    rc8de1bd rf68e2a5  
    1 This program is in the public domain. 
     1Copyright (c) 2009-2017, SasView Developers 
     2All rights reserved. 
    23 
    3 Individual files may be copyright different authors, with licences that 
    4 allow modification and redistribution in source or binary form with or 
    5 without modification, and a disclaimer of warranty. 
     4Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 
     5 
     61. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 
     7 
     82. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 
     9 
     103. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. 
     11 
     12THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
  • sasmodels/compare.py

    r01ea374 r6e5c0b7  
    3030 
    3131import sys 
     32import os 
    3233import math 
    3334import datetime 
     
    4041from . import kerneldll 
    4142from . import exception 
    42 from .data import plot_theory, empty_data1D, empty_data2D 
     43from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4344from .direct_model import DirectModel 
    4445from .convert import revert_name, revert_pars, constrain_new_to_old 
     
    8586    -help/-html shows the model docs instead of running the model 
    8687    -title="note" adds note to the plot title, after the model name 
     88    -data="path" uses q, dq from the data file 
    8789 
    8890Any two calculation engines can be selected for comparison: 
     
    747749    comp = opts['engines'][1] if have_comp else None 
    748750    data = opts['data'] 
     751    use_data = have_base ^ have_comp 
    749752 
    750753    # Plot if requested 
     
    763766    if have_base: 
    764767        if have_comp: plt.subplot(131) 
    765         plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
     768        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    766769        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    767770        #cbar_title = "log I" 
     
    769772        if have_base: plt.subplot(132) 
    770773        if not opts['is2d'] and have_base: 
    771             plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
    772         plot_theory(data, comp_value, view=view, use_data=False, limits=limits) 
     774            plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
     775        plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits) 
    773776        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    774777        #cbar_title = "log I" 
     
    784787            err[err>cutoff] = cutoff 
    785788        #err,errstr = base/comp,"ratio" 
    786         plot_theory(data, None, resid=err, view=errview, use_data=False) 
     789        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
    787790        if view == 'linear': 
    788791            plt.xscale('linear') 
     
    834837VALUE_OPTIONS = [ 
    835838    # Note: random is both a name option and a value option 
    836     'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 
     839    'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 
    837840    ] 
    838841 
     
    951954        'html'      : False, 
    952955        'title'     : None, 
     956        'data'      : None, 
    953957    } 
    954958    engines = [] 
     
    971975        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:]) 
    972976        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    973         elif arg.startswith('-title'):     opts['title'] = arg[7:] 
     977        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
     978        elif arg.startswith('-data='):     opts['data'] = arg[6:] 
    974979        elif arg == '-random':  opts['seed'] = np.random.randint(1000000) 
    975980        elif arg == '-preset':  opts['seed'] = -1 
     
    11131118 
    11141119    # Create the computational engines 
    1115     data, _ = make_data(opts) 
     1120    if opts['data'] is not None: 
     1121        data = load_data(os.path.expanduser(opts['data'])) 
     1122    else: 
     1123        data, _ = make_data(opts) 
    11161124    if n1: 
    11171125        base = make_engine(model_info, data, engines[0], opts['cutoff']) 
     
    11431151    show html docs for the model 
    11441152    """ 
    1145     import wx  # type: ignore 
    1146     from .generate import view_html_from_info 
    1147     app = wx.App() if wx.GetApp() is None else None 
    1148     view_html_from_info(opts['def'][0]) 
    1149     if app: app.MainLoop() 
    1150  
     1153    import os 
     1154    from .generate import make_html 
     1155    from . import rst2html 
     1156 
     1157    info = opts['def'][0] 
     1158    html = make_html(info) 
     1159    path = os.path.dirname(info.filename) 
     1160    url = "file://"+path.replace("\\","/")[2:]+"/" 
     1161    rst2html.view_html_qtapp(html, url) 
    11511162 
    11521163def explore(opts): 
  • sasmodels/data.py

    r40a87fa ra769b54  
    5454    if data is None: 
    5555        raise IOError("Data %r could not be loaded" % filename) 
     56    if hasattr(data, 'x'): 
     57        data.qmin, data.qmax = data.x.min(), data.x.max() 
     58        data.mask = (np.isnan(data.y) if data.y is not None 
     59                     else np.zeros_like(data.x, dtype='bool')) 
    5660    return data 
    5761 
     
    348352    # data, but they already handle the masking and graph markup already, so 
    349353    # do not repeat. 
    350     if hasattr(data, 'lam'): 
     354    if hasattr(data, 'isSesans') and data.isSesans: 
    351355        _plot_result_sesans(data, None, None, use_data=True, limits=limits) 
    352356    elif hasattr(data, 'qx_data'): 
     
    376380    *Iq_calc* is the raw theory values without resolution smearing 
    377381    """ 
    378     if hasattr(data, 'lam'): 
     382    if hasattr(data, 'isSesans') and data.isSesans: 
    379383        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 
    380384    elif hasattr(data, 'qx_data'): 
  • sasmodels/direct_model.py

    rb397165 ra769b54  
    192192 
    193193        # interpret data 
    194         if hasattr(data, 'lam'): 
     194        if hasattr(data, 'isSesans') and data.isSesans: 
    195195            self.data_type = 'sesans' 
    196196        elif hasattr(data, 'qx_data'): 
  • sasmodels/generate.py

    r1e7b0db0 rc4e3215  
    928928    from . import rst2html 
    929929    url = "file://"+dirname(info.filename)+"/" 
    930     rst2html.wxview(make_html(info), url=url) 
     930    rst2html.view_html(make_html(info), url=url) 
    931931 
    932932def demo_time(): 
  • sasmodels/rst2html.py

    r0890871 rc4e3215  
    2929from docutils.writers.html4css1 import HTMLTranslator 
    3030from docutils.nodes import SkipNode 
    31  
    32 def wxview(html, url="", size=(850, 540)): 
    33     import wx 
    34     from wx.html2 import WebView 
    35     frame = wx.Frame(None, -1, size=size) 
    36     view = WebView.New(frame) 
    37     view.SetPage(html, url) 
    38     frame.Show() 
    39     return frame 
    40  
    41 def view_rst(filename): 
    42     from os.path import expanduser 
    43     with open(expanduser(filename)) as fid: 
    44         rst = fid.read() 
    45     html = rst2html(rst) 
    46     wxview(html) 
    4731 
    4832def rst2html(rst, part="whole", math_output="mathjax"): 
     
    155139    assert replace_dollar(u"a (again $in parens$) a") == u"a (again :math:`in parens`) a" 
    156140 
    157 def view_rst_app(filename): 
     141def load_rst_as_html(filename): 
     142    from os.path import expanduser 
     143    with open(expanduser(filename)) as fid: 
     144        rst = fid.read() 
     145    html = rst2html(rst) 
     146    return html 
     147 
     148def wxview(html, url="", size=(850, 540)): 
     149    import wx 
     150    from wx.html2 import WebView 
     151    frame = wx.Frame(None, -1, size=size) 
     152    view = WebView.New(frame) 
     153    view.SetPage(html, url) 
     154    frame.Show() 
     155    return frame 
     156 
     157def qtview(html, url=""): 
     158    try: 
     159        from PyQt5.QtWebKitWidgets import QWebView 
     160        from PyQt5.QtCore import QUrl 
     161    except ImportError: 
     162        from PyQt4.QtWebkit import QWebView 
     163        from PyQt4.QtCore import QUrl 
     164    helpView = QWebView() 
     165    helpView.setHtml(html, QUrl(url)) 
     166    helpView.show() 
     167    return helpView 
     168 
     169def view_html_wxapp(html, url=""): 
    158170    import wx  # type: ignore 
    159171    app = wx.App() 
    160     view_rst(filename) 
     172    frame = wxview(html, url) 
    161173    app.MainLoop() 
    162174 
     175def view_html_qtapp(html, url=""): 
     176    import sys 
     177    try: 
     178        from PyQt5.QtWidgets import QApplication 
     179    except ImportError: 
     180        from PyQt4.QtGui import QApplication 
     181    app = QApplication([]) 
     182    frame = qtview(html, url) 
     183    sys.exit(app.exec_()) 
     184 
     185def view_rst_app(filename, qt=False): 
     186    import os 
     187    html = load_rst_as_html(filename) 
     188    url="file://"+os.path.abspath(filename)+"/" 
     189    if qt: 
     190        view_html_qtapp(html, url) 
     191    else: 
     192        view_html_wxapp(html, url) 
    163193 
    164194if __name__ == "__main__": 
    165195    import sys 
    166     view_rst_app(sys.argv[1]) 
     196    view_rst_app(sys.argv[1], qt=True) 
    167197 
  • 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/kernel_header.c

    rdaeef4c 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/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/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_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_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 r933af72  
    4848amplitudes of the core and shell, in the same manner as a core-shell model. 
    4949 
     50.. math:: 
     51 
     52    F_{a}(Q,\alpha,\beta)= 
     53    \left[\frac{\sin(Q(L_A+2t_A)/2\sin\alpha \sin\beta)}{Q(L_A+2t_A)/2\sin\alpha\sin\beta} 
     54    - \frac{\sin(QL_A/2\sin\alpha \sin\beta)}{QL_A/2\sin\alpha \sin\beta} \right] 
     55    \left[\frac{\sin(QL_B/2\sin\alpha \sin\beta)}{QL_B/2\sin\alpha \sin\beta} \right] 
     56    \left[\frac{\sin(QL_C/2\sin\alpha \sin\beta)}{QL_C/2\sin\alpha \sin\beta} \right] 
    5057 
    5158.. note:: 
     
    5764 
    5865FITTING NOTES 
    59 If the scale is set equal to the particle volume fraction, |phi|, the returned 
     66If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    6067value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    6168However, **no interparticle interference effects are included in this 
  • 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 r3b571ae  
    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 
  • 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/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/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/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/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

    r925ad6e r28d3067  
    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 is defined by the average over all orientations $\Omega$, 
    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}\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, 
     
    6999---------- 
    70100 
    71 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 
    72 and Neutron Scattering*, Plenum, New York, 1987. 
     101[1] Finnigan, J.A., Jacobs, D.J., 1971. 
     102*Light scattering by ellipsoidal particles in solution*, 
     103J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310 
     104 
    73105""" 
    74106 
     
    91123               "Solvent scattering length density"], 
    92124              ["radius_equat_minor", "Ang", 20, [0, inf], "volume", 
    93                "Minor equatorial radius"], 
     125               "Minor equatorial radius, Ra"], 
    94126              ["radius_equat_major", "Ang", 400, [0, inf], "volume", 
    95                "Major equatorial radius"], 
     127               "Major equatorial radius, Rb"], 
    96128              ["radius_polar", "Ang", 10, [0, inf], "volume", 
    97                "Polar radius"], 
     129               "Polar radius, Rc"], 
    98130              ["theta", "degrees", 60, [-inf, inf], "orientation", 
    99131               "In plane angle"], 
Note: See TracChangeset for help on using the changeset viewer.