Changeset 6e5c0b7 in sasmodels
- Timestamp:
- Apr 4, 2017 5:28:58 PM (8 years ago)
- 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. - Files:
-
- 1 added
- 21 edited
Legend:
- Unmodified
- Added
- Removed
-
LICENSE.txt
rc8de1bd rf68e2a5 1 This program is in the public domain. 1 Copyright (c) 2009-2017, SasView Developers 2 All rights reserved. 2 3 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. 4 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 5 6 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 7 8 2. 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 10 3. 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 12 THIS 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 30 30 31 31 import sys 32 import os 32 33 import math 33 34 import datetime … … 40 41 from . import kerneldll 41 42 from . import exception 42 from .data import plot_theory, empty_data1D, empty_data2D 43 from .data import plot_theory, empty_data1D, empty_data2D, load_data 43 44 from .direct_model import DirectModel 44 45 from .convert import revert_name, revert_pars, constrain_new_to_old … … 85 86 -help/-html shows the model docs instead of running the model 86 87 -title="note" adds note to the plot title, after the model name 88 -data="path" uses q, dq from the data file 87 89 88 90 Any two calculation engines can be selected for comparison: … … 747 749 comp = opts['engines'][1] if have_comp else None 748 750 data = opts['data'] 751 use_data = have_base ^ have_comp 749 752 750 753 # Plot if requested … … 763 766 if have_base: 764 767 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) 766 769 plt.title("%s t=%.2f ms"%(base.engine, base_time)) 767 770 #cbar_title = "log I" … … 769 772 if have_base: plt.subplot(132) 770 773 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) 773 776 plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 774 777 #cbar_title = "log I" … … 784 787 err[err>cutoff] = cutoff 785 788 #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) 787 790 if view == 'linear': 788 791 plt.xscale('linear') … … 834 837 VALUE_OPTIONS = [ 835 838 # 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', 837 840 ] 838 841 … … 951 954 'html' : False, 952 955 'title' : None, 956 'data' : None, 953 957 } 954 958 engines = [] … … 971 975 elif arg.startswith('-cutoff='): opts['cutoff'] = float(arg[8:]) 972 976 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:] 974 979 elif arg == '-random': opts['seed'] = np.random.randint(1000000) 975 980 elif arg == '-preset': opts['seed'] = -1 … … 1113 1118 1114 1119 # 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) 1116 1124 if n1: 1117 1125 base = make_engine(model_info, data, engines[0], opts['cutoff']) … … 1143 1151 show html docs for the model 1144 1152 """ 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) 1151 1162 1152 1163 def explore(opts): -
sasmodels/data.py
r40a87fa ra769b54 54 54 if data is None: 55 55 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')) 56 60 return data 57 61 … … 348 352 # data, but they already handle the masking and graph markup already, so 349 353 # do not repeat. 350 if hasattr(data, ' lam'):354 if hasattr(data, 'isSesans') and data.isSesans: 351 355 _plot_result_sesans(data, None, None, use_data=True, limits=limits) 352 356 elif hasattr(data, 'qx_data'): … … 376 380 *Iq_calc* is the raw theory values without resolution smearing 377 381 """ 378 if hasattr(data, ' lam'):382 if hasattr(data, 'isSesans') and data.isSesans: 379 383 _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 380 384 elif hasattr(data, 'qx_data'): -
sasmodels/direct_model.py
rb397165 ra769b54 192 192 193 193 # interpret data 194 if hasattr(data, ' lam'):194 if hasattr(data, 'isSesans') and data.isSesans: 195 195 self.data_type = 'sesans' 196 196 elif hasattr(data, 'qx_data'): -
sasmodels/generate.py
r1e7b0db0 rc4e3215 928 928 from . import rst2html 929 929 url = "file://"+dirname(info.filename)+"/" 930 rst2html. wxview(make_html(info), url=url)930 rst2html.view_html(make_html(info), url=url) 931 931 932 932 def demo_time(): -
sasmodels/rst2html.py
r0890871 rc4e3215 29 29 from docutils.writers.html4css1 import HTMLTranslator 30 30 from docutils.nodes import SkipNode 31 32 def wxview(html, url="", size=(850, 540)):33 import wx34 from wx.html2 import WebView35 frame = wx.Frame(None, -1, size=size)36 view = WebView.New(frame)37 view.SetPage(html, url)38 frame.Show()39 return frame40 41 def view_rst(filename):42 from os.path import expanduser43 with open(expanduser(filename)) as fid:44 rst = fid.read()45 html = rst2html(rst)46 wxview(html)47 31 48 32 def rst2html(rst, part="whole", math_output="mathjax"): … … 155 139 assert replace_dollar(u"a (again $in parens$) a") == u"a (again :math:`in parens`) a" 156 140 157 def view_rst_app(filename): 141 def 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 148 def 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 157 def 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 169 def view_html_wxapp(html, url=""): 158 170 import wx # type: ignore 159 171 app = wx.App() 160 view_rst(filename)172 frame = wxview(html, url) 161 173 app.MainLoop() 162 174 175 def 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 185 def 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) 163 193 164 194 if __name__ == "__main__": 165 195 import sys 166 view_rst_app(sys.argv[1] )196 view_rst_app(sys.argv[1], qt=True) 167 197 -
explore/angular_pd.py
r12eb36b r8267e0b 47 47 48 48 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 49 theta_center = radians( theta)49 theta_center = radians(90-theta) 50 50 phi_center = radians(phi) 51 51 flow_center = radians(flow) … … 137 137 radius=11., dist=dist) 138 138 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) 140 140 plt.gcf().canvas.draw() 141 141 -
sasmodels/kernel_header.c
rdaeef4c rb00a646 148 148 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 149 149 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 150 161 #if 1 151 162 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 … … 166 177 #endif 167 178 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 168 199 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 169 200 q = sqrt(qx*qx + qy*qy); \ … … 180 211 cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 181 212 } while (0) 213 #endif -
sasmodels/models/bcc_paracrystal.c
r4962519 r50beefe 90 90 double theta, double phi, double psi) 91 91 { 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); 94 94 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; 98 98 99 99 const double qd = 0.5*q*dnn; -
sasmodels/models/core_shell_bicelle.c
r592343f rb260926 30 30 31 31 static double 32 bicelle_kernel(double q q,32 bicelle_kernel(double q, 33 33 double rad, 34 34 double radthick, 35 35 double facthick, 36 double length,36 double halflength, 37 37 double rhoc, 38 38 double rhoh, … … 42 42 double cos_alpha) 43 43 { 44 double si1,si2,be1,be2;45 46 44 const double dr1 = rhoc-rhoh; 47 45 const double dr2 = rhor-rhosolv; 48 46 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); 56 50 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); 61 55 62 56 const double t = vol1*dr1*si1*be1 + … … 64 58 vol3*dr3*si2*be1; 65 59 66 const double retval = t*t *sin_alpha;60 const double retval = t*t; 67 61 68 62 return retval; … … 71 65 72 66 static double 73 bicelle_integration(double q q,67 bicelle_integration(double q, 74 68 double rad, 75 69 double radthick, … … 83 77 // set up the integration end points 84 78 const double uplim = M_PI_4; 85 const double half height= 0.5*length;79 const double halflength = 0.5*length; 86 80 87 81 double summ = 0.0; … … 90 84 double sin_alpha, cos_alpha; // slots to hold sincos function output 91 85 SINCOS(alpha, sin_alpha, cos_alpha); 92 double yyy = Gauss76Wt[i] * bicelle_kernel(q q, rad, radthick, facthick,93 half height, rhoc, rhoh, rhor, rhosolv,86 double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 87 halflength, rhoc, rhoh, rhor, rhosolv, 94 88 sin_alpha, cos_alpha); 95 summ += yyy ;89 summ += yyy*sin_alpha; 96 90 } 97 91 … … 119 113 double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 120 114 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; 126 117 } 127 118 -
sasmodels/models/core_shell_bicelle_elliptical.c
r592343f rf4f85b3 32 32 } 33 33 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) 34 double Iq(double q, 35 double rad, 36 double x_core, 37 double radthick, 38 double facthick, 39 double length, 40 double rhoc, 41 double rhoh, 42 double rhor, 43 double rhosolv) 45 44 { 46 45 double si1,si2,be1,be2; … … 74 73 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 75 74 double inner_sum=0; 76 double sinarg1 = q q*halfheight*cos_alpha;77 double sinarg2 = q q*(halfheight+facthick)*cos_alpha;75 double sinarg1 = q*halfheight*cos_alpha; 76 double sinarg2 = q*(halfheight+facthick)*cos_alpha; 78 77 si1 = sas_sinx_x(sinarg1); 79 78 si2 = sas_sinx_x(sinarg2); … … 83 82 const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 84 83 const double rr = sqrt(rA - rB*cos(beta)); 85 double besarg1 = q q*rr*sin_alpha;86 double besarg2 = q q*(rr+radthick)*sin_alpha;84 double besarg1 = q*rr*sin_alpha; 85 double besarg2 = q*(rr+radthick)*sin_alpha; 87 86 be1 = sas_2J1x_x(besarg1); 88 87 be2 = sas_2J1x_x(besarg2); … … 114 113 { 115 114 // THIS NEEDS TESTING 116 double q q, cos_val, cos_mu, cos_nu;117 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q q, cos_val, cos_mu, cos_nu);115 double q, xhat, yhat, zhat; 116 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 118 117 const double dr1 = rhoc-rhoh; 119 118 const double dr2 = rhor-rhosolv; … … 125 124 const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick); 126 125 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) 128 127 // Given: radius_major = r_ratio * radius_minor 129 128 // 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 ); 135 138 const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 + vol3*dr3*si2*be1); 136 139 //const double vol = form_volume(radius_minor, r_ratio, length); -
sasmodels/models/core_shell_parallelepiped.c
r1e7b0db0 r92dfe0c 134 134 double psi) 135 135 { 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); 138 138 139 139 // cspkernel in csparallelepiped recoded here … … 160 160 double tc = length_a + 2.0*thick_rim_c; 161 161 //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); 168 168 169 169 -
sasmodels/models/core_shell_parallelepiped.py
rcb0dc22 r933af72 48 48 amplitudes of the core and shell, in the same manner as a core-shell model. 49 49 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] 50 57 51 58 .. note:: … … 57 64 58 65 FITTING NOTES 59 If the scale is set equal to the particle volume fraction, |phi|, the returned66 If the scale is set equal to the particle volume fraction, $\phi$, the returned 60 67 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 61 68 However, **no interparticle interference effects are included in this -
sasmodels/models/ellipsoid.c
r130d4c7 r3b571ae 3 3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 4 4 double radius_polar, double radius_equatorial, double theta, double phi); 5 6 static double7 _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 the11 // 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 be18 // slightly better for 2D which has already computed it, but it introduces19 // an extra sqrt and square for 1-D not required by the current form, so20 // leave it as is.21 const double r = radius_equatorial22 * 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 }27 5 28 6 double form_volume(double radius_polar, double radius_equatorial) … … 37 15 double radius_equatorial) 38 16 { 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 39 26 // translate a point in [-1,1] to a point in [0, 1] 27 // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 40 28 const double zm = 0.5; 41 29 const double zb = 0.5; 42 30 double total = 0.0; 43 31 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; 47 36 } 48 37 // translate dx in [-1,1] to dx in [lower,upper] … … 62 51 double q, sin_alpha, cos_alpha; 63 52 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); 65 56 const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 66 57 67 return 1.0e-4 * form * s * s;58 return 1.0e-4 * square(f * s); 68 59 } 69 60 -
sasmodels/models/ellipsoid.py
r925ad6e r3b571ae 18 18 .. math:: 19 19 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} 23 21 24 and 22 for 25 23 26 24 .. math:: 27 25 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} 30 27 31 28 32 29 $\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 31 radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial 32 radius perpendicular to the rotational axis of the ellipsoid and 33 $\Delta \rho$ (contrast) is the scattering length density difference between 34 the scatterer and the solvent. 37 35 38 For randomly oriented particles :36 For randomly oriented particles use the orientational average, 39 37 40 38 .. math:: 41 39 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} 43 41 42 43 computed via substitution of $u=\sin(\alpha)$, $du=\cos(\alpha)\,d\alpha$ as 44 45 .. math:: 46 47 \langle F^2(q) \rangle = \int_0^1{F^2(q, u)\,du} 48 49 with 50 51 .. math:: 52 53 r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2} 44 54 45 55 To provide easy access to the orientation of the ellipsoid, we define … … 48 58 :ref:`cylinder orientation figure <cylinder-angle-definition>`. 49 59 For the ellipsoid, $\theta$ is the angle between the rotational axis 50 and the $z$ -axis. 60 and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$ 61 in the $xy$ plane. 51 62 52 63 NB: The 2nd virial coefficient of the solid ellipsoid is calculated based … … 90 101 than 500. 91 102 103 Model was also tested against the triaxial ellipsoid model with equal major 104 and minor equatorial radii. It is also consistent with the cyclinder model 105 with polar radius equal to length and equatorial radius equal to radius. 106 92 107 References 93 108 ---------- … … 96 111 *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 97 112 Plenum Press, New York, 1987. 113 114 Authorship and Verification 115 ---------------------------- 116 117 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 118 * **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014 119 * **Last Modified by:** Paul Kienzle **Date:** March 22, 2017 98 120 """ 99 121 -
sasmodels/models/elliptical_cylinder.c
r592343f r61104c8 67 67 double theta, double phi, double psi) 68 68 { 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); 71 71 72 72 // Compute: r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 73 73 // 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)); 75 75 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); 77 77 const double Aq = be * si; 78 78 const double delrho = sld - solvent_sld; -
sasmodels/models/fcc_paracrystal.c
r4962519 r50beefe 90 90 double theta, double phi, double psi) 91 91 { 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); 94 94 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; 98 98 const double qd = 0.5*q*dnn; 99 99 const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); -
sasmodels/models/parallelepiped.c
r1e7b0db0 rd605080 67 67 double psi) 68 68 { 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); 71 71 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); 75 75 const double V = form_volume(length_a, length_b, length_c); 76 76 const double drho = (sld - solvent_sld); -
sasmodels/models/sc_paracrystal.c
r4962519 r50beefe 111 111 double psi) 112 112 { 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); 115 115 116 116 const double qd = q*dnn; … … 118 118 const double tanh_qd = tanh(arg); 119 119 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); 123 123 124 124 const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; -
sasmodels/models/triaxial_ellipsoid.c
r925ad6e r68dd6a9 20 20 double radius_polar) 21 21 { 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; 26 27 double outer = 0.0; 27 28 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)); 34 32 35 33 double inner = 0.0; 34 const double um = 0.5; 35 const double ub = 0.5; 36 36 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; 41 42 } 42 outer += Gauss76Wt[i] * 0.5 * inner;43 outer += Gauss76Wt[i] * inner; // correcting for dx later 43 44 } 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); 46 47 const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 47 48 return 1.0e-4 * s * s * fqsq; … … 58 59 double psi) 59 60 { 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); 62 63 63 const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu64 + radius_equat_major*radius_equat_major*cmu*cmu65 + 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); 67 68 const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 68 69 -
sasmodels/models/triaxial_ellipsoid.py
r925ad6e r28d3067 2 2 # Note: model title and parameter table are inserted automatically 3 3 r""" 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**. 4 Definition 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 12 Given an ellipsoid 6 13 7 14 .. math:: 8 15 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 10 17 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 18 the scattering is defined by the average over all orientations $\Omega$, 22 19 23 20 .. math:: 24 21 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} 28 23 29 24 where … … 31 26 .. math:: 32 27 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 32 The $e$, $f$ and $g$ terms are the projections of the orientation vector on $X$, 33 $Y$ and $Z$ respectively. Keeping the orientation fixed at the canonical 34 axes, we can integrate over the incident direction using polar angle 35 $-\pi/2 \le \gamma \le \pi/2$ and equatorial angle $0 \le \phi \le 2\pi$ 36 (as defined in ref [1]), 37 38 .. math:: 39 40 \langle\Phi^2\rangle = \int_0^{2\pi} \int_{-\pi/2}^{\pi/2} \Phi^2(qr) \cos \gamma\,d\gamma d\phi 41 42 with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$. 43 A little algebra yields 44 45 .. math:: 46 47 r^2 = b^2(p_a \sin^2 \phi \cos^2 \gamma + 1 + p_c \sin^2 \gamma) 48 49 for 50 51 .. math:: 52 53 p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 1 54 55 Due to symmetry, the ranges can be restricted to a single quadrant 56 $0 \le \gamma \le \pi/2$ and $0 \le \phi \le \pi/2$, scaling the resulting 57 integral by 8. The computation is done using the substitution $u = \sin\gamma$, 58 $du = \cos\gamma\,d\gamma$, giving 59 60 .. math:: 61 62 \langle\Phi^2\rangle &= 8 \int_0^{\pi/2} \int_0^1 \Phi^2(qr) du d\phi \\ 63 r^2 &= b^2(p_a \sin^2(\phi)(1 - u^2) + 1 + p_c u^2) 34 64 35 65 To provide easy access to the orientation of the triaxial ellipsoid, … … 69 99 ---------- 70 100 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*, 103 J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310 104 73 105 """ 74 106 … … 91 123 "Solvent scattering length density"], 92 124 ["radius_equat_minor", "Ang", 20, [0, inf], "volume", 93 "Minor equatorial radius "],125 "Minor equatorial radius, Ra"], 94 126 ["radius_equat_major", "Ang", 400, [0, inf], "volume", 95 "Major equatorial radius "],127 "Major equatorial radius, Rb"], 96 128 ["radius_polar", "Ang", 10, [0, inf], "volume", 97 "Polar radius "],129 "Polar radius, Rc"], 98 130 ["theta", "degrees", 60, [-inf, inf], "orientation", 99 131 "In plane angle"],
Note: See TracChangeset
for help on using the changeset viewer.