Changeset 6ab64c9 in sasmodels
- Timestamp:
- Oct 12, 2017 2:53:26 PM (7 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 0d19f42
- Parents:
- 706f466 (diff), 09141ff (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:
-
- 8 added
- 2 deleted
- 37 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/guide/magnetism/magnetism.rst
r59485a4 r1f058ea 77 77 =========== ================================================================ 78 78 M0_sld = $D_M M_0$ 79 Up_theta = $\theta_ {up}$79 Up_theta = $\theta_\mathrm{up}$ 80 80 M_theta = $\theta_M$ 81 81 M_phi = $\phi_M$ -
doc/guide/pd/polydispersity.rst
rf8a2baa r1f058ea 95 95 \exp\left(-\frac{(x - \bar x)^2}{2\sigma^2}\right) 96 96 97 where $\bar x$ is the mean of the distribution and *Norm* is a normalization factor98 which is determined during the numerical calculation.97 where $\bar x$ is the mean of the distribution and *Norm* is a normalization 98 factor which is determined during the numerical calculation. 99 99 100 100 The polydispersity is … … 122 122 during the numerical calculation. 123 123 124 The median value for the distribution will be the value given for the respective125 size parameter, for example, *radius=60*.124 The median value for the distribution will be the value given for the 125 respective size parameter, for example, *radius=60*. 126 126 127 127 The polydispersity is given by $\sigma$ … … 208 208 209 209 Many commercial Dynamic Light Scattering (DLS) instruments produce a size 210 polydispersity parameter, sometimes even given the symbol $p$ This210 polydispersity parameter, sometimes even given the symbol $p$\ ! This 211 211 parameter is defined as the relative standard deviation coefficient of 212 212 variation of the size distribution and is NOT the same as the polydispersity -
doc/guide/resolution.rst
r30b60d2 r1f058ea 17 17 resolution contribution into a model calculation/simulation (which by definition 18 18 will be exact) to make it more representative of what has been measured 19 experimentally - a process called *smearing*. sasmodels does the latter.19 experimentally - a process called *smearing*. Sasmodels does the latter. 20 20 21 21 Both smearing and desmearing rely on functions to describe the resolution 22 effect. sasmodels provides three smearing algorithms:22 effect. Sasmodels provides three smearing algorithms: 23 23 24 24 * *Slit Smearing* … … 99 99 100 100 For discrete $q$ values, at the $q$ values of the data points and at the $q$ 101 values extended up to $q_N = q_i + \Delta q_ v$ the smeared101 values extended up to $q_N = q_i + \Delta q_u$ the smeared 102 102 intensity can be approximately calculated as 103 103 -
doc/guide/scripting.rst
r2e66ef5 r4aa5dce 69 69 $ bumps example/model.py --preview 70 70 71 Note that bumps and sasmodels are included as part of the SasView 72 distribution. On windows, bumps can be called from the cmd prompt 73 as follows:: 74 75 SasViewCom bumps.cli example/model.py --preview 76 71 77 Using sasmodels directly 72 78 ======================== … … 105 111 plt.loglog(q, Iq) 106 112 plt.show() 113 114 On windows, this can be called from the cmd prompt using sasview as:: 115 116 SasViewCom example/cylinder_eval.py -
example/oriented_usans.py
r1cd24b4 r74b0495 6 6 7 7 # Spherical particle data, not ellipsoids 8 sans, usans = load_data(' ../../sasview/sasview/test/1d_data/latex_smeared.xml')8 sans, usans = load_data('latex_smeared.xml', index='all') 9 9 usans.qmin, usans.qmax = np.min(usans.x), np.max(usans.x) 10 10 usans.mask = (usans.x < 0.0) -
example/simul_fit.py
r1a4d4c0 r74b0495 1 #!/usr/bin/env python2 # -*- coding: utf-8 -*-3 4 # To Sasview/documents/scripts5 6 1 from bumps.names import * 7 2 from sasmodels.core import load_model … … 9 4 from sasmodels.data import load_data, plot_data 10 5 6 # latex data, same sample usans and sans 7 # particles radius ~2300, uniform dispersity 8 datasets = load_data('latex_smeared.xml', index='all') 9 #[print(data) for data in datasets] 11 10 12 """ IMPORT THE DATA USED """ 13 datafiles = ['latex_smeared_out_0.txt', 'latex_smeared_out_1.txt'] 14 datasets = [load_data(el) for el in datafiles] 11 # A single sphere model to share between the datasets. We will use 12 # FreeVariables below to set the parameters that are independent between 13 # the datasets. 14 kernel = load_model('sphere') 15 pars = dict(scale=0.01, background=0.0, sld=5.0, sld_solvent=0.0, radius=1500., 16 #radius_pd=0.1, radius_pd_n=35, 17 ) 18 model = Model(kernel, **pars) 15 19 16 for data in datasets: 17 data.qmin = 0.0 18 data.qmax = 10.0 20 # radius and polydispersity (if any) are shared 21 model.radius.range(0, inf) 22 #model.radius_pd.range(0, 1) 19 23 20 #sphere model 21 kernel = load_model('sphere', dtype="single") 22 pars = dict(scale=0.01, background=0.0, sld=1.0, sld_solvent=6.0, radius=1500.) 23 model = Model(kernel, **pars) 24 model.radius.range(0, inf) 25 #model.background.range(-inf, inf) 26 #model.scale.range(0, inf) 27 model.sld.range(-inf, inf) 28 model.sld_solvent.range(-inf, inf) 24 # Contrast and dilution are the same for both measurements, but are not 25 # separable with a single measurement (i.e., I(q) ~ F(q) contrast^2 Vf), 26 # so fit one of scale, sld or solvent sld. With absolute scaling from 27 # data reduction, can use the same parameter for both datasets. 28 model.scale.range(0, inf) 29 #model.sld.range(-inf, inf) 30 #model.sld_solvent.range(-inf, inf) 29 31 32 # Background is different for sans and usans so set it as a free variable 33 # in the model. 30 34 free = FreeVariables( 31 names=[data. filenamefor data in datasets],35 names=[data.run[0] for data in datasets], 32 36 background=model.background, 33 scale=model.scale,34 37 ) 35 38 free.background.range(-inf, inf) 36 free.scale.range(0, inf)37 39 38 M = [Experiment(data=data, model=model) for data in datasets] 40 # Note: can access the parameters for the individual models using 41 # free.background[0] and free.background[1], setting constraints or 42 # ranges as appropriate. 43 44 # For more complex systems where different datasets require independent models, 45 # separate models can be defined, with parameters tied together using 46 # constraint expressions. For example, the following could be used to fit 47 # data set 1 to spheres and data set 2 to cylinders of the same volume: 48 # model1 = Model(load_model('sphere')) 49 # model2 = Model(load_model('cylinder')) 50 # model1.sld = model2.sld 51 # model1.sld_solvent = model2.sld_solvent 52 # model1.scale = model2.scale 53 # # set cylinders and spheres to the same volume 54 # model1.radius = (3/4*model2.radius**2*model2.length)**(1/3) 55 # model1.background.range(0, 2) 56 # model2.background.range(0, 2) 57 58 # Setup the experiments, sharing the same model across all datasets. 59 M = [Experiment(data=data, model=model, name=data.run[0]) for data in datasets] 39 60 40 61 problem = FitProblem(M, freevars=free) 41 42 print(problem._parameters) -
sasmodels/bumps_model.py
r3330bb4 r74b0495 133 133 """ 134 134 _cache = None # type: Dict[str, np.ndarray] 135 def __init__(self, data, model, cutoff=1e-5 ):135 def __init__(self, data, model, cutoff=1e-5, name=None): 136 136 # type: (Data, Model, float) -> None 137 137 # remember inputs so we can inspect from outside 138 self.name = data.filename if name is None else name 138 139 self.model = model 139 140 self.cutoff = cutoff … … 204 205 """ 205 206 data, theory, resid = self._data, self.theory(), self.residuals() 206 plot_theory(data, theory, resid, view, Iq_calc=self.Iq_calc) 207 # TODO: hack to display oriented usans 2-D pattern 208 Iq_calc = self.Iq_calc if isinstance(self.Iq_calc, tuple) else None 209 plot_theory(data, theory, resid, view, Iq_calc=Iq_calc) 207 210 208 211 def simulate_data(self, noise=None): -
sasmodels/compare.py
r765eb0e rb76191e 533 533 % (pd, n, nsigma, nsigma, pdtype) 534 534 if M0 != 0.: 535 line += " M0:%.3f m phi:%.1f mtheta:%.1f" % (M0, mphi, mtheta)535 line += " M0:%.3f mtheta:%.1f mphi:%.1f" % (M0, mtheta, mphi) 536 536 return line 537 537 -
sasmodels/data.py
r630156b r09141ff 44 44 Data = Union["Data1D", "Data2D", "SesansData"] 45 45 46 def load_data(filename ):46 def load_data(filename, index=0): 47 47 # type: (str) -> Data 48 48 """ … … 55 55 filename, indexstr = filename[:-1].split('[') 56 56 index = int(indexstr) 57 else:58 index = None59 57 datasets = loader.load(filename) 60 if datasets is None:58 if not datasets: # None or [] 61 59 raise IOError("Data %r could not be loaded" % filename) 62 60 if not isinstance(datasets, list): 63 61 datasets = [datasets] 64 if index is None and len(datasets) > 1: 65 raise ValueError("Need to specify filename[index] for multipart data") 66 data = datasets[index if index is not None else 0] 67 if hasattr(data, 'x'): 68 data.qmin, data.qmax = data.x.min(), data.x.max() 69 data.mask = (np.isnan(data.y) if data.y is not None 70 else np.zeros_like(data.x, dtype='bool')) 71 elif hasattr(data, 'qx_data'): 72 data.mask = ~data.mask 73 return data 62 for data in datasets: 63 if hasattr(data, 'x'): 64 data.qmin, data.qmax = data.x.min(), data.x.max() 65 data.mask = (np.isnan(data.y) if data.y is not None 66 else np.zeros_like(data.x, dtype='bool')) 67 elif hasattr(data, 'qx_data'): 68 data.mask = ~data.mask 69 return datasets[index] if index != 'all' else datasets 74 70 75 71 -
doc/genapi.py
r2e66ef5 r706f466 59 59 #('alignment', 'GPU data alignment [unused]'), 60 60 ('bumps_model', 'Bumps interface'), 61 ('compare', 'Compare models on different compute engines'), 61 62 ('compare_many', 'Batch compare models on different compute engines'), 62 ('co mpare', 'Compare models on different compute engines'),63 ('conversion_table', 'Model conversion table'), 63 64 ('convert', 'Sasview to sasmodel converter'), 64 65 ('core', 'Model access'), … … 82 83 ('sasview_model', 'Sasview interface'), 83 84 ('sesans', 'SESANS calculation routines'), 85 ('special', 'Special functions library'), 84 86 ('weights', 'Distribution functions'), 85 87 ] -
explore/jitter.py
r85190c2 rd4c33d6 3 3 dispersity and possible replacement algorithms. 4 4 """ 5 import sys 6 5 7 import mpl_toolkits.mplot3d # Adds projection='3d' option to subplot 6 8 import matplotlib.pyplot as plt 7 9 from matplotlib.widgets import Slider, CheckButtons 8 10 from matplotlib import cm 9 10 11 import numpy as np 11 12 from numpy import pi, cos, sin, sqrt, exp, degrees, radians 12 13 13 def draw_beam(ax ):14 def draw_beam(ax, view=(0, 0)): 14 15 #ax.plot([0,0],[0,0],[1,-1]) 15 16 #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) … … 22 23 x = r*np.outer(np.cos(u), np.ones_like(v)) 23 24 y = r*np.outer(np.sin(u), np.ones_like(v)) 24 z = np.outer(np.ones_like(u), v) 25 z = 1.3*np.outer(np.ones_like(u), v) 26 27 theta, phi = view 28 shape = x.shape 29 points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 30 points = Rz(phi)*Ry(theta)*points 31 x, y, z = [v.reshape(shape) for v in points] 25 32 26 33 ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 27 34 28 def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi): 29 size=[0.1, 0.4, 1.0] 30 view=[theta, phi, psi] 31 shimmy=[0,0,0] 32 #draw_shape = draw_parallelepiped 33 draw_shape = draw_ellipsoid 35 def draw_jitter(ax, view, jitter): 36 size = [0.1, 0.4, 1.0] 37 draw_shape = draw_parallelepiped 38 #draw_shape = draw_ellipsoid 34 39 35 40 #np.random.seed(10) … … 64 69 [ 1, 1, 1], 65 70 ] 71 dtheta, dphi, dpsi = jitter 66 72 if dtheta == 0: 67 73 cloud = [v for v in cloud if v[0] == 0] … … 70 76 if dpsi == 0: 71 77 cloud = [v for v in cloud if v[2] == 0] 72 draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8)78 draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 73 79 for point in cloud: 74 shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]]75 draw_shape(ax, size, view, shimmy, alpha=0.8)80 delta = [dtheta*point[0], dphi*point[1], dpsi*point[2]] 81 draw_shape(ax, size, view, delta, alpha=0.8) 76 82 for v in 'xyz': 77 83 a, b, c = size … … 80 86 getattr(ax, v+'axis').label.set_text(v) 81 87 82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1):88 def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 83 89 a,b,c = size 84 theta, phi, psi = view85 dtheta, dphi, dpsi = shimmy86 87 90 u = np.linspace(0, 2 * np.pi, steps) 88 91 v = np.linspace(0, np.pi, steps) … … 90 93 y = b*np.outer(np.sin(u), np.sin(v)) 91 94 z = c*np.outer(np.ones_like(u), np.cos(v)) 92 93 shape = x.shape 94 points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 95 points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 96 points = Rz(phi)*Ry(theta)*Rz(psi)*points 97 x,y,z = [v.reshape(shape) for v in points] 95 x, y, z = transform_xyz(view, jitter, x, y, z) 98 96 99 97 ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 100 98 101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 99 draw_labels(ax, view, jitter, [ 100 ('c+', [ 0, 0, c], [ 1, 0, 0]), 101 ('c-', [ 0, 0,-c], [ 0, 0,-1]), 102 ('a+', [ a, 0, 0], [ 0, 0, 1]), 103 ('a-', [-a, 0, 0], [ 0, 0,-1]), 104 ('b+', [ 0, b, 0], [-1, 0, 0]), 105 ('b-', [ 0,-b, 0], [-1, 0, 0]), 106 ]) 107 108 def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 102 109 a,b,c = size 103 theta, phi, psi = view104 dtheta, dphi, dpsi = shimmy105 106 110 x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 107 111 y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) … … 118 122 ]) 119 123 120 points = np.matrix([x,y,z]) 121 points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 122 points = Rz(phi)*Ry(theta)*Rz(psi)*points 123 124 x,y,z = [np.array(v).flatten() for v in points] 124 x, y, z = transform_xyz(view, jitter, x, y, z) 125 125 ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 126 126 127 def draw_sphere(ax, radius=10., steps=100): 128 u = np.linspace(0, 2 * np.pi, steps) 129 v = np.linspace(0, np.pi, steps) 130 131 x = radius * np.outer(np.cos(u), np.sin(v)) 132 y = radius * np.outer(np.sin(u), np.sin(v)) 133 z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 134 ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 135 136 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 137 theta_center = radians(theta) 138 phi_center = radians(phi) 139 flow_center = radians(flow) 140 dtheta = radians(dtheta) 141 dphi = radians(dphi) 142 143 # 10 point 3-sigma gaussian weights 144 t = np.linspace(-3., 3., 11) 127 draw_labels(ax, view, jitter, [ 128 ('c+', [ 0, 0, c], [ 1, 0, 0]), 129 ('c-', [ 0, 0,-c], [ 0, 0,-1]), 130 ('a+', [ a, 0, 0], [ 0, 0, 1]), 131 ('a-', [-a, 0, 0], [ 0, 0,-1]), 132 ('b+', [ 0, b, 0], [-1, 0, 0]), 133 ('b-', [ 0,-b, 0], [-1, 0, 0]), 134 ]) 135 136 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gauss'): 137 theta, phi, psi = view 138 dtheta, dphi, dpsi = jitter 145 139 if dist == 'gauss': 140 t = np.linspace(-3, 3, n) 146 141 weights = exp(-0.5*t**2) 147 142 elif dist == 'rect': 143 t = np.linspace(0, 1, n) 148 144 weights = np.ones_like(t) 149 145 else: 150 146 raise ValueError("expected dist to be 'gauss' or 'rect'") 151 theta = dtheta*t 152 phi = dphi*t 153 154 x = radius * np.outer(cos(phi), cos(theta)) 155 y = radius * np.outer(sin(phi), cos(theta)) 156 z = radius * np.outer(np.ones_like(phi), sin(theta)) 157 #w = np.outer(weights, weights*abs(cos(dtheta*t))) 158 w = np.outer(weights, weights*abs(cos(theta))) 159 160 x, y, z, w = [v.flatten() for v in (x,y,z,w)] 161 x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center) 162 163 ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 164 165 def rotate(x, y, z, phi, theta, psi): 166 R = Rz(phi)*Ry(theta)*Rz(psi) 167 p = np.vstack([x,y,z]) 168 return R*p 147 148 # mesh in theta, phi formed by rotating z 149 z = np.matrix([[0], [0], [radius]]) 150 points = np.hstack([Rx(phi_i)*Ry(theta_i)*z 151 for theta_i in dtheta*t 152 for phi_i in dphi*t]) 153 # rotate relative to beam 154 points = orient_relative_to_beam(view, points) 155 156 w = np.outer(weights, weights) 157 158 x, y, z = [np.array(v).flatten() for v in points] 159 ax.scatter(x, y, z, c=w.flatten(), marker='o', vmin=0., vmax=1.) 169 160 170 161 def Rx(angle): … … 188 179 [0., 0., 1.]] 189 180 return np.matrix(R) 181 182 def transform_xyz(view, jitter, x, y, z): 183 x, y, z = [np.asarray(v) for v in (x, y, z)] 184 shape = x.shape 185 points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 186 points = apply_jitter(jitter, points) 187 points = orient_relative_to_beam(view, points) 188 x, y, z = [np.array(v).reshape(shape) for v in points] 189 return x, y, z 190 191 def apply_jitter(jitter, points): 192 dtheta, dphi, dpsi = jitter 193 points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 194 return points 195 196 def orient_relative_to_beam(view, points): 197 theta, phi, psi = view 198 points = Rz(phi)*Ry(theta)*Rz(psi)*points 199 return points 200 201 def draw_labels(ax, view, jitter, text): 202 labels, locations, orientations = zip(*text) 203 px, py, pz = zip(*locations) 204 dx, dy, dz = zip(*orientations) 205 206 px, py, pz = transform_xyz(view, jitter, px, py, pz) 207 dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 208 209 for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 210 zdir = np.asarray(zdir).flatten() 211 ax.text(p[0], p[1], p[2], label, zdir=zdir) 212 213 def draw_sphere(ax, radius=10., steps=100): 214 u = np.linspace(0, 2 * np.pi, steps) 215 v = np.linspace(0, np.pi, steps) 216 217 x = radius * np.outer(np.cos(u), np.sin(v)) 218 y = radius * np.outer(np.sin(u), np.sin(v)) 219 z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 220 ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 190 221 191 222 def main(): … … 206 237 207 238 axcolor = 'lightgoldenrodyellow' 239 208 240 axtheta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 209 241 axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) … … 212 244 sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 213 245 spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 246 214 247 axdtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 215 248 axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) … … 217 250 sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 218 251 sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 219 sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dp hi)252 sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dpsi) 220 253 221 254 def update(val, axis=None): 222 theta, phi, psi= stheta.val, sphi.val, spsi.val223 dtheta, dphi, dpsi= sdtheta.val, sdphi.val, sdpsi.val255 view = stheta.val, sphi.val, spsi.val 256 jitter = sdtheta.val, sdphi.val, sdpsi.val 224 257 ax.cla() 225 draw_beam(ax )226 draw_ shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi)227 # if not axis.startswith('d'):228 # ax.view_init(elev=theta, azim=phi)258 draw_beam(ax, (0, 0)) 259 draw_jitter(ax, view, jitter) 260 #draw_jitter(ax, view, (0,0,0)) 261 draw_mesh(ax, view, jitter) 229 262 plt.gcf().canvas.draw() 230 263 -
sasmodels/details.py
rccd5f01 rf39759c 15 15 16 16 import numpy as np # type: ignore 17 from numpy import pi, cos, sin 17 from numpy import pi, cos, sin, radians 18 18 19 19 try: … … 29 29 30 30 try: 31 from typing import List 31 from typing import List, Tuple, Sequence 32 32 except ImportError: 33 33 pass 34 34 else: 35 35 from .modelinfo import ModelInfo 36 from .modelinfo import ParameterTable 36 37 37 38 … … 53 54 coordinates, the total circumference decreases as latitude varies from 54 55 pi r^2 at the equator to 0 at the pole, and the weight associated 55 with a range of phivalues needs to be scaled by this circumference.56 with a range of latitude values needs to be scaled by this circumference. 56 57 This scale factor needs to be updated each time the theta value 57 58 changes. *theta_par* indicates which of the values in the parameter … … 231 232 nvalues = kernel.info.parameters.nvalues 232 233 scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs] 233 values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) 234 # skipping scale and background when building values and weights 235 values, weights = zip(*pairs[2:npars+2]) if npars else ((), ()) 236 weights = correct_theta_weights(kernel.info.parameters, values, weights) 234 237 length = np.array([len(w) for w in weights]) 235 238 offset = np.cumsum(np.hstack((0, length))) … … 244 247 return call_details, data, is_magnetic 245 248 249 def correct_theta_weights(parameters, values, weights): 250 # type: (ParameterTable, Sequence[np.ndarray], Sequence[np.ndarray]) -> Sequence[np.ndarray] 251 """ 252 If there is a theta parameter, update the weights of that parameter so that 253 the cosine weighting required for polar integration is preserved. Avoid 254 evaluation strictly at the pole, which would otherwise send the weight to 255 zero. 256 257 Note: values and weights do not include scale and background 258 """ 259 # TODO: document code, explaining why scale and background are skipped 260 # given that we don't have scale and background in the list, we 261 # should be calling the variables something other than values and weights 262 # Apparently the parameters.theta_offset similarly skips scale and 263 # and background, so the indexing works out. 264 if parameters.theta_offset >= 0: 265 index = parameters.theta_offset 266 theta = values[index] 267 theta_weight = np.minimum(abs(cos(radians(theta))), 1e-6) 268 # copy the weights list so we can update it 269 weights = list(weights) 270 weights[index] = theta_weight*np.asarray(weights[index]) 271 weights = tuple(weights) 272 return weights 273 246 274 247 275 def convert_magnetism(parameters, values): 276 # type: (ParameterTable, Sequence[np.ndarray]) -> bool 248 277 """ 249 278 Convert magnetism values from polar to rectangular coordinates. … … 255 284 scale = mag[:,0] 256 285 if np.any(scale): 257 theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180.286 theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 258 287 cos_theta = cos(theta) 259 288 mag[:, 0] = scale*cos_theta*cos(phi) # mx … … 269 298 """ 270 299 Create a mesh grid of dispersion parameters and weights. 300 301 pars is a list of pairs (values, weights), where the values are the 302 individual parameter values at which to evaluate the polydispersity 303 distribution and weights are the weights associated with each value. 304 305 Only the volume parameters should be included in this list. Orientation 306 parameters do not affect the calculation of effective radius or volume 307 ratio. 271 308 272 309 Returns [p1,p2,...],w where pj is a vector of values for parameter j -
sasmodels/kernel_header.c
rbb4b509 r73cbc5b 164 164 SINCOS(phi*M_PI_180, sn, cn); \ 165 165 q = sqrt(qx*qx + qy*qy); \ 166 cn 166 cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180)); \ 167 167 sn = sqrt(1 - cn*cn); \ 168 168 } while (0) -
sasmodels/kernel_iq.c
rbde38b5 rd4c33d6 25 25 int32_t num_weights; // total length of the weights vector 26 26 int32_t num_active; // number of non-trivial pd loops 27 int32_t theta_par; // id of spherical correction variable 27 int32_t theta_par; // id of spherical correction variable (not used) 28 28 } ProblemDetails; 29 29 … … 173 173 174 174 175 #if MAX_PD>0176 const int theta_par = details->theta_par;177 const int fast_theta = (theta_par == p0);178 const int slow_theta = (theta_par >= 0 && !fast_theta);179 double spherical_correction = 1.0;180 #else181 // Note: if not polydisperse the weights cancel and we don't need the182 // spherical correction.183 const double spherical_correction = 1.0;184 #endif185 186 175 int step = pd_start; 187 176 … … 220 209 #endif 221 210 #if MAX_PD>0 222 if (slow_theta) { // Theta is not in inner loop223 spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);224 }225 211 while(i0 < n0) { 226 212 local_values.vector[p0] = v0[i0]; 227 213 double weight0 = w0[i0] * weight1; 228 214 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); 229 if (fast_theta) { // Theta is in inner loop230 spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);231 }232 215 #else 233 216 const double weight0 = 1.0; … … 244 227 // Note: weight==0 must always be excluded 245 228 if (weight0 > cutoff) { 246 // spherical correction is set at a minimum of 1e-6, otherwise there 247 // would be problems looking at models with theta=90. 248 const double weight = weight0 * spherical_correction; 249 pd_norm += weight * CALL_VOLUME(local_values.table); 229 pd_norm += weight0 * CALL_VOLUME(local_values.table); 250 230 251 231 #ifdef USE_OPENMP … … 304 284 #endif // !MAGNETIC 305 285 //printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 306 result[q_index] += weight * scattering;286 result[q_index] += weight0 * scattering; 307 287 } 308 288 } -
sasmodels/kernel_iq.cl
rbde38b5 rd4c33d6 25 25 int32_t num_weights; // total length of the weights vector 26 26 int32_t num_active; // number of non-trivial pd loops 27 int32_t theta_par; // id of spherical correction variable 27 int32_t theta_par; // id of spherical correction variable (not used) 28 28 } ProblemDetails; 29 29 … … 169 169 170 170 171 #if MAX_PD>0172 const int theta_par = details->theta_par;173 const bool fast_theta = (theta_par == p0);174 const bool slow_theta = (theta_par >= 0 && !fast_theta);175 double spherical_correction = 1.0;176 #else177 // Note: if not polydisperse the weights cancel and we don't need the178 // spherical correction.179 const double spherical_correction = 1.0;180 #endif181 182 171 int step = pd_start; 183 172 … … 217 206 #endif 218 207 #if MAX_PD>0 219 if (slow_theta) { // Theta is not in inner loop220 spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);221 }222 208 while(i0 < n0) { 223 209 local_values.vector[p0] = v0[i0]; 224 210 double weight0 = w0[i0] * weight1; 225 211 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); 226 if (fast_theta) { // Theta is in inner loop227 spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);228 }229 212 #else 230 213 const double weight0 = 1.0; … … 232 215 233 216 //if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); } 234 //if (q_index == 0) printf("sphcor: %g\n", spherical_correction);235 217 236 218 #ifdef INVALID … … 241 223 // Note: weight==0 must always be excluded 242 224 if (weight0 > cutoff) { 243 // spherical correction is set at a minimum of 1e-6, otherwise there 244 // would be problems looking at models with theta=90. 245 const double weight = weight0 * spherical_correction; 246 pd_norm += weight * CALL_VOLUME(local_values.table); 225 pd_norm += weight0 * CALL_VOLUME(local_values.table); 247 226 248 227 #if defined(MAGNETIC) && NUM_MAGNETIC > 0 … … 296 275 const double scattering = CALL_IQ(q, q_index, local_values.table); 297 276 #endif // !MAGNETIC 298 this_result += weight * scattering;277 this_result += weight0 * scattering; 299 278 } 300 279 } -
sasmodels/kernelpy.py
r3b32f8d r62d7601 201 201 202 202 pd_norm = 0.0 203 spherical_correction = 1.0204 203 partial_weight = np.NaN 205 204 weight = np.NaN 206 205 207 206 p0_par = call_details.pd_par[0] 208 p0_is_theta = (p0_par == call_details.theta_par)209 207 p0_length = call_details.pd_length[0] 210 208 p0_index = p0_length … … 223 221 parameters[pd_par] = pd_value[pd_offset+pd_index] 224 222 partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 225 if call_details.theta_par >= 0:226 cor = sin(pi / 180 * parameters[call_details.theta_par])227 spherical_correction = max(abs(cor), 1e-6)228 223 p0_index = loop_index%p0_length 229 224 230 225 weight = partial_weight * pd_weight[p0_offset + p0_index] 231 226 parameters[p0_par] = pd_value[p0_offset + p0_index] 232 if p0_is_theta:233 cor = cos(pi/180 * parameters[p0_par])234 spherical_correction = max(abs(cor), 1e-6)235 227 p0_index += 1 236 228 if weight > cutoff: … … 244 236 245 237 # update value and norm 246 weight *= spherical_correction247 238 total += weight * Iq 248 239 pd_norm += weight * form_volume() -
sasmodels/model_test.py
r65314f7 r8db25bf 86 86 suite = unittest.TestSuite() 87 87 88 if models[0] == 'all':88 if models[0] in core.KINDS: 89 89 skip = models[1:] 90 models = list_models( )90 models = list_models(models[0]) 91 91 else: 92 92 skip = [] -
sasmodels/models/barbell.c
r592343f r2a0b2b1 10 10 //barbell kernel - same as dumbell 11 11 static double 12 _bell_kernel(double q , double h, double radius_bell,13 double half_length , double sin_alpha, double cos_alpha)12 _bell_kernel(double qab, double qc, double h, double radius_bell, 13 double half_length) 14 14 { 15 15 // translate a point in [-1,1] to a point in [lower,upper] … … 26 26 // m = q R cos(alpha) 27 27 // b = q(L/2-h) cos(alpha) 28 const double m = q*radius_bell*cos_alpha; // cos argument slope29 const double b = q*(half_length-h)*cos_alpha; // cos argument intercept30 const double q rst = q*radius_bell*sin_alpha; // Q*R*sin(theta)28 const double m = radius_bell*qc; // cos argument slope 29 const double b = (half_length-h)*qc; // cos argument intercept 30 const double qab_r = radius_bell*qab; // Q*R*sin(theta) 31 31 double total = 0.0; 32 32 for (int i = 0; i < 76; i++){ 33 33 const double t = Gauss76Z[i]*zm + zb; 34 34 const double radical = 1.0 - t*t; 35 const double bj = sas_2J1x_x(q rst*sqrt(radical));35 const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 36 36 const double Fq = cos(m*t + b) * radical * bj; 37 37 total += Gauss76Wt[i] * Fq; … … 44 44 45 45 static double 46 _fq(double q, double h, 47 double radius_bell, double radius, double half_length, 48 double sin_alpha, double cos_alpha) 46 _fq(double qab, double qc, double h, 47 double radius_bell, double radius, double half_length) 49 48 { 50 const double bell_fq = _bell_kernel(q , h, radius_bell, half_length, sin_alpha, cos_alpha);51 const double bj = sas_2J1x_x( q*radius*sin_alpha);52 const double si = sas_sinx_x( q*half_length*cos_alpha);49 const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length); 50 const double bj = sas_2J1x_x(radius*qab); 51 const double si = sas_sinx_x(half_length*qc); 53 52 const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 54 53 const double Aq = bell_fq + cyl_fq; … … 84 83 double sin_alpha, cos_alpha; // slots to hold sincos function output 85 84 SINCOS(alpha, sin_alpha, cos_alpha); 86 const double Aq = _fq(q , h, radius_bell, radius, half_length, sin_alpha, cos_alpha);85 const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 87 86 total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 88 87 } … … 103 102 double q, sin_alpha, cos_alpha; 104 103 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 104 const double qab = q*sin_alpha; 105 const double qc = q*cos_alpha; 105 106 106 107 const double h = -sqrt(square(radius_bell) - square(radius)); 107 const double Aq = _fq(q , h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha);108 const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 108 109 109 110 // Multiply by contrast^2 and convert to cm-1 -
sasmodels/models/bcc_paracrystal.c
r50beefe r7e0b281 1 double form_volume(double radius); 2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 3 double Iqxy(double qx, double qy, double dnn, 4 double d_factor, double radius,double sld, double solvent_sld, 5 double theta, double phi, double psi); 1 static double 2 bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 3 { 4 #if 0 // Equations as written in Matsuoka 5 const double a1 = (+qa + qb + qc)/2.0; 6 const double a2 = (-qa - qb + qc)/2.0; 7 const double a3 = (-qa + qb - qc)/2.0; 8 #else 9 const double a1 = (+qa + qb - qc)/2.0; 10 const double a2 = (+qa - qb + qc)/2.0; 11 const double a3 = (-qa + qb + qc)/2.0; 12 #endif 6 13 7 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 8 double _BCCeval(double Theta, double Phi, double temp1, double temp3); 9 double _sphereform(double q, double radius, double sld, double solvent_sld); 14 #if 1 15 // Numerator: (1 - exp(a)^2)^3 16 // => (-(exp(2a) - 1))^3 17 // => -expm1(2a)^3 18 // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 19 // => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 20 // => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 21 const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 22 const double exp_arg = exp(arg); 23 const double Zq = -cube(expm1(2.0*arg)) 24 / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 25 * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 26 * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 27 #else 28 // Alternate form, which perhaps is more approachable 29 const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 30 const double sinh_qd = sinh(arg); 31 const double cosh_qd = cosh(arg); 32 const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) 33 * sinh_qd/(cosh_qd - cos(dnn*a2)) 34 * sinh_qd/(cosh_qd - cos(dnn*a3)); 35 #endif 36 37 return Zq; 38 } 10 39 11 40 12 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 13 14 const double Da = d_factor*dnn; 15 const double temp1 = q*q*Da*Da; 16 const double temp3 = q*dnn; 17 18 double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 19 return(retVal); 41 // occupied volume fraction calculated from lattice symmetry and sphere radius 42 static double 43 bcc_volume_fraction(double radius, double dnn) 44 { 45 return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 20 46 } 21 47 22 double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 23 24 double result; 25 double sin_theta,cos_theta,sin_phi,cos_phi; 26 SINCOS(Theta, sin_theta, cos_theta); 27 SINCOS(Phi, sin_phi, cos_phi); 28 29 const double temp6 = sin_theta; 30 const double temp7 = sin_theta*cos_phi + sin_theta*sin_phi + cos_theta; 31 const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 32 const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 33 34 const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 35 result = cube(1.0 - (temp10*temp10))*temp6 36 / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 37 * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 38 * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 39 40 return (result); 41 } 42 43 double form_volume(double radius){ 48 static double 49 form_volume(double radius) 50 { 44 51 return sphere_volume(radius); 45 52 } 46 53 47 54 48 double Iq(double q, double dnn,55 static double Iq(double q, double dnn, 49 56 double d_factor, double radius, 50 double sld, double solvent_sld){ 57 double sld, double solvent_sld) 58 { 59 // translate a point in [-1,1] to a point in [0, 2 pi] 60 const double phi_m = M_PI; 61 const double phi_b = M_PI; 62 // translate a point in [-1,1] to a point in [0, pi] 63 const double theta_m = M_PI_2; 64 const double theta_b = M_PI_2; 51 65 52 //Volume fraction calculated from lattice symmetry and sphere radius 53 const double s1 = dnn/sqrt(0.75); 54 const double latticescale = 2.0*sphere_volume(radius/s1); 55 56 const double va = 0.0; 57 const double vb = 2.0*M_PI; 58 const double vaj = 0.0; 59 const double vbj = M_PI; 60 61 double summ = 0.0; 62 double answer = 0.0; 63 for(int i=0; i<150; i++) { 64 //setup inner integral over the ellipsoidal cross-section 65 double summj=0.0; 66 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy is phi 67 for(int j=0;j<150;j++) { 68 //20 gauss points for the inner integral 69 double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy is theta 70 double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi); 71 summj += yyy; 72 } 73 //now calculate the value of the inner integral 74 double answer = (vbj-vaj)/2.0*summj; 75 76 //now calculate outer integral 77 summ = summ+(Gauss150Wt[i] * answer); 78 } //final scaling is done at the end of the function, after the NT_FP64 case 79 80 answer = (vb-va)/2.0*summ; 81 answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 82 83 return answer; 66 double outer_sum = 0.0; 67 for(int i=0; i<150; i++) { 68 double inner_sum = 0.0; 69 const double theta = Gauss150Z[i]*theta_m + theta_b; 70 double sin_theta, cos_theta; 71 SINCOS(theta, sin_theta, cos_theta); 72 const double qc = q*cos_theta; 73 const double qab = q*sin_theta; 74 for(int j=0;j<150;j++) { 75 const double phi = Gauss150Z[j]*phi_m + phi_b; 76 double sin_phi, cos_phi; 77 SINCOS(phi, sin_phi, cos_phi); 78 const double qa = qab*cos_phi; 79 const double qb = qab*sin_phi; 80 const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 81 inner_sum += Gauss150Wt[j] * form; 82 } 83 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 84 outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 85 } 86 outer_sum *= theta_m; 87 const double Zq = outer_sum/(4.0*M_PI); 88 const double Pq = sphere_form(q, radius, sld, solvent_sld); 89 return bcc_volume_fraction(radius, dnn) * Pq * Zq; 84 90 } 85 91 86 92 87 double Iqxy(double qx, double qy,93 static double Iqxy(double qx, double qy, 88 94 double dnn, double d_factor, double radius, 89 95 double sld, double solvent_sld, … … 92 98 double q, zhat, yhat, xhat; 93 99 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 100 const double qa = q*xhat; 101 const double qb = q*yhat; 102 const double qc = q*zhat; 94 103 95 const double a1 = +xhat - zhat + yhat; 96 const double a2 = +xhat + zhat - yhat; 97 const double a3 = -xhat + zhat + yhat; 98 99 const double qd = 0.5*q*dnn; 100 const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 101 const double tanh_qd = tanh(arg); 102 const double cosh_qd = cosh(arg); 103 const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 104 * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 105 * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 106 107 const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 108 //the occupied volume of the lattice 109 const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 110 return lattice_scale * Fq; 104 q = sqrt(qa*qa + qb*qb + qc*qc); 105 const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); 106 const double Pq = sphere_form(q, radius, sld, solvent_sld); 107 return bcc_volume_fraction(radius, dnn) * Pq * Zq; 111 108 } -
sasmodels/models/capped_cylinder.c
r592343f r2a0b2b1 14 14 // radius_cap is the radius of the lens 15 15 // length is the cylinder length, or the separation between the lens halves 16 // alpha is the angle of the cylinder wrt q.16 // theta is the angle of the cylinder wrt q. 17 17 static double 18 _cap_kernel(double q , double h, double radius_cap,19 double half_length, double sin_alpha, double cos_alpha)18 _cap_kernel(double qab, double qc, double h, double radius_cap, 19 double half_length) 20 20 { 21 21 // translate a point in [-1,1] to a point in [lower,upper] … … 26 26 27 27 // cos term in integral is: 28 // cos (q (R t - h + L/2) cos( alpha))28 // cos (q (R t - h + L/2) cos(theta)) 29 29 // so turn it into: 30 30 // cos (m t + b) 31 31 // where: 32 // m = q R cos( alpha)33 // b = q(L/2-h) cos( alpha)34 const double m = q*radius_cap*cos_alpha; // cos argument slope35 const double b = q*(half_length-h)*cos_alpha; // cos argument intercept36 const double q rst = q*radius_cap*sin_alpha; // Q*R*sin(theta)32 // m = q R cos(theta) 33 // b = q(L/2-h) cos(theta) 34 const double m = radius_cap*qc; // cos argument slope 35 const double b = (half_length-h)*qc; // cos argument intercept 36 const double qab_r = radius_cap*qab; // Q*R*sin(theta) 37 37 double total = 0.0; 38 38 for (int i=0; i<76 ;i++) { 39 39 const double t = Gauss76Z[i]*zm + zb; 40 40 const double radical = 1.0 - t*t; 41 const double bj = sas_2J1x_x(q rst*sqrt(radical));41 const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 42 42 const double Fq = cos(m*t + b) * radical * bj; 43 43 total += Gauss76Wt[i] * Fq; … … 50 50 51 51 static double 52 _fq(double q, double h, double radius_cap, double radius, double half_length, 53 double sin_alpha, double cos_alpha) 52 _fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 54 53 { 55 const double cap_Fq = _cap_kernel(q , h, radius_cap, half_length, sin_alpha, cos_alpha);56 const double bj = sas_2J1x_x( q*radius*sin_alpha);57 const double si = sas_sinx_x( q*half_length*cos_alpha);54 const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length); 55 const double bj = sas_2J1x_x(radius*qab); 56 const double si = sas_sinx_x(half_length*qc); 58 57 const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 59 58 const double Aq = cap_Fq + cyl_Fq; … … 101 100 double total = 0.0; 102 101 for (int i=0; i<76 ;i++) { 103 const double alpha= Gauss76Z[i]*zm + zb; 104 double sin_alpha, cos_alpha; // slots to hold sincos function output 105 SINCOS(alpha, sin_alpha, cos_alpha); 106 107 const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha); 108 // sin_alpha for spherical coord integration 109 total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 102 const double theta = Gauss76Z[i]*zm + zb; 103 double sin_theta, cos_theta; // slots to hold sincos function output 104 SINCOS(theta, sin_theta, cos_theta); 105 const double qab = q*sin_theta; 106 const double qc = q*cos_theta; 107 const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 108 // scale by sin_theta for spherical coord integration 109 total += Gauss76Wt[i] * Aq * Aq * sin_theta; 110 110 } 111 111 // translate dx in [-1,1] to dx in [lower,upper] … … 125 125 double q, sin_alpha, cos_alpha; 126 126 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 127 const double qab = q*sin_alpha; 128 const double qc = q*cos_alpha; 127 129 128 130 const double h = sqrt(radius_cap*radius_cap - radius*radius); 129 const double Aq = _fq(q , h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha);131 const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 130 132 131 133 // Multiply by contrast^2 and convert to cm-1 -
sasmodels/models/core_shell_bicelle.c
rb260926 r2a0b2b1 1 double form_volume(double radius, double thick_rim, double thick_face, double length); 2 double Iq(double q, 3 double radius, 4 double thick_rim, 5 double thick_face, 6 double length, 7 double core_sld, 8 double face_sld, 9 double rim_sld, 10 double solvent_sld); 11 12 13 double Iqxy(double qx, double qy, 14 double radius, 15 double thick_rim, 16 double thick_face, 17 double length, 18 double core_sld, 19 double face_sld, 20 double rim_sld, 21 double solvent_sld, 22 double theta, 23 double phi); 24 25 26 double form_volume(double radius, double thick_rim, double thick_face, double length) 1 static double 2 form_volume(double radius, double thick_rim, double thick_face, double length) 27 3 { 28 return M_PI* (radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face);4 return M_PI*square(radius+thick_rim)*(length+2.0*thick_face); 29 5 } 30 6 31 7 static double 32 bicelle_kernel(double q, 33 double rad, 34 double radthick, 35 double facthick, 36 double halflength, 37 double rhoc, 38 double rhoh, 39 double rhor, 40 double rhosolv, 41 double sin_alpha, 42 double cos_alpha) 8 bicelle_kernel(double qab, 9 double qc, 10 double radius, 11 double thick_radius, 12 double thick_face, 13 double halflength, 14 double sld_core, 15 double sld_face, 16 double sld_rim, 17 double sld_solvent) 43 18 { 44 const double dr1 = rhoc-rhoh;45 const double dr2 = rhor-rhosolv;46 const double dr3 = rhoh-rhor;47 const double vol1 = M_PI*square(rad )*2.0*(halflength);48 const double vol2 = M_PI*square(rad +radthick)*2.0*(halflength+facthick);49 const double vol3 = M_PI*square(rad )*2.0*(halflength+facthick);19 const double dr1 = sld_core-sld_face; 20 const double dr2 = sld_rim-sld_solvent; 21 const double dr3 = sld_face-sld_rim; 22 const double vol1 = M_PI*square(radius)*2.0*(halflength); 23 const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face); 24 const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face); 50 25 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);26 const double be1 = sas_2J1x_x((radius)*qab); 27 const double be2 = sas_2J1x_x((radius+thick_radius)*qab); 28 const double si1 = sas_sinx_x((halflength)*qc); 29 const double si2 = sas_sinx_x((halflength+thick_face)*qc); 55 30 56 31 const double t = vol1*dr1*si1*be1 + … … 58 33 vol3*dr3*si2*be1; 59 34 60 const double retval = t*t; 61 62 return retval; 63 35 return t; 64 36 } 65 37 66 38 static double 67 bicelle_integration(double q,68 double rad,69 double radthick,70 double facthick,71 72 double rhoc,73 double rhoh,74 double rhor,75 double rhosolv)39 Iq(double q, 40 double radius, 41 double thick_radius, 42 double thick_face, 43 double length, 44 double sld_core, 45 double sld_face, 46 double sld_rim, 47 double sld_solvent) 76 48 { 77 49 // set up the integration end points … … 79 51 const double halflength = 0.5*length; 80 52 81 double summ= 0.0;53 double total = 0.0; 82 54 for(int i=0;i<N_POINTS_76;i++) { 83 double alpha = (Gauss76Z[i] + 1.0)*uplim; 84 double sin_alpha, cos_alpha; // slots to hold sincos function output 85 SINCOS(alpha, sin_alpha, cos_alpha); 86 double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 87 halflength, rhoc, rhoh, rhor, rhosolv, 88 sin_alpha, cos_alpha); 89 summ += yyy*sin_alpha; 55 double theta = (Gauss76Z[i] + 1.0)*uplim; 56 double sin_theta, cos_theta; // slots to hold sincos function output 57 SINCOS(theta, sin_theta, cos_theta); 58 double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 59 halflength, sld_core, sld_face, sld_rim, sld_solvent); 60 total += Gauss76Wt[i]*fq*fq*sin_theta; 90 61 } 91 62 92 63 // calculate value of integral to return 93 double answer = uplim*summ;94 return answer;64 double answer = total*uplim; 65 return 1.0e-4*answer; 95 66 } 96 67 97 68 static double 98 bicelle_kernel_2d(double qx, double qy,99 100 101 102 103 104 105 106 107 108 69 Iqxy(double qx, double qy, 70 double radius, 71 double thick_rim, 72 double thick_face, 73 double length, 74 double core_sld, 75 double face_sld, 76 double rim_sld, 77 double solvent_sld, 78 double theta, 79 double phi) 109 80 { 110 81 double q, sin_alpha, cos_alpha; 111 82 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 83 const double qab = q*sin_alpha; 84 const double qc = q*cos_alpha; 112 85 113 double answer = bicelle_kernel(q, radius, thick_rim, thick_face,86 double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face, 114 87 0.5*length, core_sld, face_sld, rim_sld, 115 solvent_sld , sin_alpha, cos_alpha);116 return 1.0e-4* answer;88 solvent_sld); 89 return 1.0e-4*fq*fq; 117 90 } 118 119 double Iq(double q,120 double radius,121 double thick_rim,122 double thick_face,123 double length,124 double core_sld,125 double face_sld,126 double rim_sld,127 double solvent_sld)128 {129 double intensity = bicelle_integration(q, radius, thick_rim, thick_face,130 length, core_sld, face_sld, rim_sld, solvent_sld);131 return intensity*1.0e-4;132 }133 134 135 double Iqxy(double qx, double qy,136 double radius,137 double thick_rim,138 double thick_face,139 double length,140 double core_sld,141 double face_sld,142 double rim_sld,143 double solvent_sld,144 double theta,145 double phi)146 {147 double intensity = bicelle_kernel_2d(qx, qy,148 radius,149 thick_rim,150 thick_face,151 length,152 core_sld,153 face_sld,154 rim_sld,155 solvent_sld,156 theta,157 phi);158 159 return intensity;160 } -
sasmodels/models/core_shell_bicelle_elliptical.c
rdedcf34 r2a0b2b1 17 17 double thick_face, 18 18 double length, 19 double rhoc,20 double rhoh,21 double rhor,22 double rhosolv)19 double sld_core, 20 double sld_face, 21 double sld_rim, 22 double sld_solvent) 23 23 { 24 double si1,si2,be1,be2;25 24 // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 26 25 // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 27 // const double uplim = M_PI_4;28 26 const double halfheight = 0.5*length; 29 //const double va = 0.0;30 //const double vb = 1.0;31 // inner integral limits32 //const double vaj=0.0;33 //const double vbj=M_PI;34 35 27 const double r_major = r_minor * x_core; 36 28 const double r2A = 0.5*(square(r_major) + square(r_minor)); 37 29 const double r2B = 0.5*(square(r_major) - square(r_minor)); 38 const double dr1 = (rhoc-rhoh) *M_PI*r_minor*r_major*(2.0*halfheight);;39 const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);40 const double dr3 = (rhoh-rhor) *M_PI*r_minor*r_major*2.0*(halfheight+thick_face);41 //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);42 //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);43 //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);30 const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 31 const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 32 const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 33 const double dr1 = vol1*(sld_core-sld_face); 34 const double dr2 = vol2*(sld_rim-sld_solvent); 35 const double dr3 = vol3*(sld_face-sld_rim); 44 36 45 37 //initialize integral … … 47 39 for(int i=0;i<76;i++) { 48 40 //setup inner integral over the ellipsoidal cross-section 49 // since we generate these lots of times, why not store them somewhere? 50 //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 51 const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 52 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 53 double inner_sum=0; 54 double sinarg1 = q*halfheight*cos_alpha; 55 double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 56 si1 = sas_sinx_x(sinarg1); 57 si2 = sas_sinx_x(sinarg2); 41 //const double va = 0.0; 42 //const double vb = 1.0; 43 //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 44 const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 45 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 46 const double qab = q*sin_theta; 47 const double qc = q*cos_theta; 48 const double si1 = sas_sinx_x(halfheight*qc); 49 const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 50 double inner_sum=0.0; 58 51 for(int j=0;j<76;j++) { 59 52 //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 60 //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 61 const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 62 const double rr = sqrt(r2A - r2B*cos(beta)); 63 double besarg1 = q*rr*sin_alpha; 64 double besarg2 = q*(rr+thick_rim)*sin_alpha; 65 be1 = sas_2J1x_x(besarg1); 66 be2 = sas_2J1x_x(besarg2); 67 inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 68 dr2*si2*be2 + 69 dr3*si2*be1); 53 // inner integral limits 54 //const double vaj=0.0; 55 //const double vbj=M_PI; 56 //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 57 const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 58 const double rr = sqrt(r2A - r2B*cos(phi)); 59 const double be1 = sas_2J1x_x(rr*qab); 60 const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 61 const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 62 63 inner_sum += Gauss76Wt[j] * fq * fq; 70 64 } 71 65 //now calculate outer integral … … 83 77 double thick_face, 84 78 double length, 85 double rhoc,86 double rhoh,87 double rhor,88 double rhosolv,79 double sld_core, 80 double sld_face, 81 double sld_rim, 82 double sld_solvent, 89 83 double theta, 90 84 double phi, 91 85 double psi) 92 86 { 93 // THIS NEEDS TESTING94 87 double q, xhat, yhat, zhat; 95 88 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 96 const double dr1 = rhoc-rhoh; 97 const double dr2 = rhor-rhosolv; 98 const double dr3 = rhoh-rhor; 89 const double qa = q*xhat; 90 const double qb = q*yhat; 91 const double qc = q*zhat; 92 93 const double dr1 = sld_core-sld_face; 94 const double dr2 = sld_rim-sld_solvent; 95 const double dr3 = sld_face-sld_rim; 99 96 const double r_major = r_minor*x_core; 100 97 const double halfheight = 0.5*length; … … 104 101 105 102 // Compute effective radius in rotated coordinates 106 const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat));107 const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat)108 + square((r_minor+thick_rim)* yhat));109 const double be1 = sas_2J1x_x( q *r_hat );110 const double be2 = sas_2J1x_x( q *rshell_hat );111 const double si1 = sas_sinx_x( q*halfheight*zhat);112 const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat);113 const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 + vol3*dr3*si2*be1);114 return 1.0e-4 * Aq;103 const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 104 const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 105 + square((r_minor+thick_rim)*qb)); 106 const double be1 = sas_2J1x_x( qr_hat ); 107 const double be2 = sas_2J1x_x( qrshell_hat ); 108 const double si1 = sas_sinx_x( halfheight*qc ); 109 const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 110 const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 + vol3*dr3*si2*be1; 111 return 1.0e-4 * fq*fq; 115 112 } 116 -
sasmodels/models/core_shell_cylinder.c
r592343f r2a0b2b1 1 double form_volume(double radius, double thickness, double length);2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld,3 double radius, double thickness, double length);4 double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld,5 double radius, double thickness, double length, double theta, double phi);6 7 1 // vd = volume * delta_rho 8 // besarg = q * R * sin(alpha) 9 // siarg = q * L/2 * cos(alpha) 10 double _cyl(double vd, double besarg, double siarg); 11 double _cyl(double vd, double besarg, double siarg) 2 // besarg = q * R * sin(theta) 3 // siarg = q * L/2 * cos(theta) 4 static double _cyl(double vd, double besarg, double siarg) 12 5 { 13 6 return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 14 7 } 15 8 16 double form_volume(double radius, double thickness, double length) 9 static double 10 form_volume(double radius, double thickness, double length) 17 11 { 18 return M_PI* (radius+thickness)*(radius+thickness)*(length+2.0*thickness);12 return M_PI*square(radius+thickness)*(length+2.0*thickness); 19 13 } 20 14 21 double Iq(double q, 15 static double 16 Iq(double q, 22 17 double core_sld, 23 18 double shell_sld, … … 28 23 { 29 24 // precalculate constants 30 const double core_ qr = q*radius;31 const double core_ qh = q*0.5*length;25 const double core_r = radius; 26 const double core_h = 0.5*length; 32 27 const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 33 const double shell_ qr = q*(radius + thickness);34 const double shell_ qh = q*(0.5*length + thickness);28 const double shell_r = (radius + thickness); 29 const double shell_h = (0.5*length + thickness); 35 30 const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 36 31 double total = 0.0; 37 // double lower=0, upper=M_PI_2;38 32 for (int i=0; i<76 ;i++) { 39 // translate a point in [-1,1] to a point in [lower,upper] 40 //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 41 double sn, cn; 42 const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 43 SINCOS(alpha, sn, cn); 44 const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 45 + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 46 total += Gauss76Wt[i] * fq * fq * sn; 33 // translate a point in [-1,1] to a point in [0, pi/2] 34 //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 35 double sin_theta, cos_theta; 36 const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 37 SINCOS(theta, sin_theta, cos_theta); 38 const double qab = q*sin_theta; 39 const double qc = q*cos_theta; 40 const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 41 + _cyl(shell_vd, shell_r*qab, shell_h*qc); 42 total += Gauss76Wt[i] * fq * fq * sin_theta; 47 43 } 48 44 // translate dx in [-1,1] to dx in [lower,upper] … … 64 60 double q, sin_alpha, cos_alpha; 65 61 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 62 const double qab = q*sin_alpha; 63 const double qc = q*cos_alpha; 66 64 67 const double core_ qr = q*radius;68 const double core_ qh = q*0.5*length;65 const double core_r = radius; 66 const double core_h = 0.5*length; 69 67 const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 70 const double shell_ qr = q*(radius + thickness);71 const double shell_ qh = q*(0.5*length + thickness);68 const double shell_r = (radius + thickness); 69 const double shell_h = (0.5*length + thickness); 72 70 const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 73 71 74 const double fq = _cyl(core_vd, core_ qr*sin_alpha, core_qh*cos_alpha)75 + _cyl(shell_vd, shell_ qr*sin_alpha, shell_qh*cos_alpha);72 const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 73 + _cyl(shell_vd, shell_r*qab, shell_h*qc); 76 74 return 1.0e-4 * fq * fq; 77 75 } -
sasmodels/models/core_shell_ellipsoid.c
r0a3d9b2 r2a0b2b1 1 double form_volume(double radius_equat_core,2 double polar_core,3 double equat_shell,4 double polar_shell);5 double Iq(double q,6 double radius_equat_core,7 double x_core,8 double thick_shell,9 double x_polar_shell,10 double core_sld,11 double shell_sld,12 double solvent_sld);13 1 2 // Converted from Igor function gfn4, using the same pattern as ellipsoid 3 // for evaluating the parts of the integral. 4 // FUNCTION gfn4: CONTAINS F(Q,A,B,MU)**2 AS GIVEN 5 // BY (53) & (58-59) IN CHEN AND 6 // KOTLARCHYK REFERENCE 7 // 8 // <OBLATE ELLIPSOID> 9 static double 10 _cs_ellipsoid_kernel(double qab, double qc, 11 double equat_core, double polar_core, 12 double equat_shell, double polar_shell, 13 double sld_core_shell, double sld_shell_solvent) 14 { 15 const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc)); 16 const double si_core = sas_3j1x_x(qr_core); 17 const double volume_core = M_4PI_3*equat_core*equat_core*polar_core; 18 const double fq_core = si_core*volume_core*sld_core_shell; 14 19 15 double Iqxy(double qx, double qy, 16 double radius_equat_core, 17 double x_core, 18 double thick_shell, 19 double x_polar_shell, 20 double core_sld, 21 double shell_sld, 22 double solvent_sld, 23 double theta, 24 double phi); 20 const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc)); 21 const double si_shell = sas_3j1x_x(qr_shell); 22 const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell; 23 const double fq_shell = si_shell*volume_shell*sld_shell_solvent; 25 24 25 return fq_core + fq_shell; 26 } 26 27 27 double form_volume(double radius_equat_core, 28 double x_core, 29 double thick_shell, 30 double x_polar_shell) 28 static double 29 form_volume(double radius_equat_core, 30 double x_core, 31 double thick_shell, 32 double x_polar_shell) 31 33 { 32 34 const double equat_shell = radius_equat_core + thick_shell; … … 37 39 38 40 static double 39 core_shell_ellipsoid_xt_kernel(double q,40 41 42 43 44 45 46 41 Iq(double q, 42 double radius_equat_core, 43 double x_core, 44 double thick_shell, 45 double x_polar_shell, 46 double core_sld, 47 double shell_sld, 48 double solvent_sld) 47 49 { 48 const double lolim = 0.0; 49 const double uplim = 1.0; 50 51 52 const double delpc = core_sld - shell_sld; //core - shell 53 const double delps = shell_sld - solvent_sld; //shell - solvent 54 50 const double sld_core_shell = core_sld - shell_sld; 51 const double sld_shell_solvent = shell_sld - solvent_sld; 55 52 56 53 const double polar_core = radius_equat_core*x_core; … … 58 55 const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 59 56 60 double summ = 0.0; //initialize intergral 57 // translate from [-1, 1] => [0, 1] 58 const double m = 0.5; 59 const double b = 0.5; 60 double total = 0.0; //initialize intergral 61 61 for(int i=0;i<76;i++) { 62 double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim ); 63 double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell, 64 polar_shell, delpc, delps, q); 65 summ += Gauss76Wt[i] * yyy; 62 const double cos_theta = Gauss76Z[i]*m + b; 63 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 64 double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 65 radius_equat_core, polar_core, 66 equat_shell, polar_shell, 67 sld_core_shell, sld_shell_solvent); 68 total += Gauss76Wt[i] * fq * fq; 66 69 } 67 summ *= 0.5*(uplim-lolim);70 total *= m; 68 71 69 72 // convert to [cm-1] 70 return 1.0e-4 * summ;73 return 1.0e-4 * total; 71 74 } 72 75 73 76 static double 74 core_shell_ellipsoid_xt_kernel_2d(double qx, double qy,75 76 77 78 79 80 81 82 83 77 Iqxy(double qx, double qy, 78 double radius_equat_core, 79 double x_core, 80 double thick_shell, 81 double x_polar_shell, 82 double core_sld, 83 double shell_sld, 84 double solvent_sld, 85 double theta, 86 double phi) 84 87 { 85 88 double q, sin_alpha, cos_alpha; 86 89 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 90 const double qab = q*sin_alpha; 91 const double qc = q*cos_alpha; 87 92 88 const double sld cs= core_sld - shell_sld;89 const double sld ss = shell_sld- solvent_sld;93 const double sld_core_shell = core_sld - shell_sld; 94 const double sld_shell_solvent = shell_sld - solvent_sld; 90 95 91 96 const double polar_core = radius_equat_core*x_core; … … 93 98 const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 94 99 95 // Call the IGOR library function to get the kernel: 96 // MUST use gfn4 not gf2 because of the def of params. 97 double answer = gfn4(cos_alpha, 98 radius_equat_core, 99 polar_core, 100 equat_shell, 101 polar_shell, 102 sldcs, 103 sldss, 104 q); 100 double fq = _cs_ellipsoid_kernel(qab, qc, 101 radius_equat_core, polar_core, 102 equat_shell, polar_shell, 103 sld_core_shell, sld_shell_solvent); 105 104 106 105 //convert to [cm-1] 107 answer *= 1.0e-4; 108 109 return answer; 106 return 1.0e-4 * fq * fq; 110 107 } 111 112 double Iq(double q,113 double radius_equat_core,114 double x_core,115 double thick_shell,116 double x_polar_shell,117 double core_sld,118 double shell_sld,119 double solvent_sld)120 {121 double intensity = core_shell_ellipsoid_xt_kernel(q,122 radius_equat_core,123 x_core,124 thick_shell,125 x_polar_shell,126 core_sld,127 shell_sld,128 solvent_sld);129 130 return intensity;131 }132 133 134 double Iqxy(double qx, double qy,135 double radius_equat_core,136 double x_core,137 double thick_shell,138 double x_polar_shell,139 double core_sld,140 double shell_sld,141 double solvent_sld,142 double theta,143 double phi)144 {145 double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy,146 radius_equat_core,147 x_core,148 thick_shell,149 x_polar_shell,150 core_sld,151 shell_sld,152 solvent_sld,153 theta,154 phi);155 156 return intensity;157 } -
sasmodels/models/core_shell_ellipsoid.py
r30b60d2 r8db25bf 141 141 # pylint: enable=bad-whitespace, line-too-long 142 142 143 source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 144 "core_shell_ellipsoid.c"] 143 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 145 144 146 145 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): -
sasmodels/models/core_shell_parallelepiped.c
r92dfe0c r2a0b2b1 1 double form_volume(double length_a, double length_b, double length_c, 1 double form_volume(double length_a, double length_b, double length_c, 2 2 double thick_rim_a, double thick_rim_b, double thick_rim_c); 3 3 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, … … 9 9 double thick_rim_c, double theta, double phi, double psi); 10 10 11 double form_volume(double length_a, double length_b, double length_c, 11 double form_volume(double length_a, double length_b, double length_c, 12 12 double thick_rim_a, double thick_rim_b, double thick_rim_c) 13 13 { 14 14 //return length_a * length_b * length_c; 15 return length_a * length_b * length_c + 16 2.0 * thick_rim_a * length_b * length_c + 15 return length_a * length_b * length_c + 16 2.0 * thick_rim_a * length_b * length_c + 17 17 2.0 * thick_rim_b * length_a * length_c + 18 18 2.0 * thick_rim_c * length_a * length_b; … … 34 34 // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 35 35 // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 36 36 37 37 const double mu = 0.5 * q * length_b; 38 38 39 39 //calculate volume before rescaling (in original code, but not used) 40 //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 41 //double vol = length_a * length_b * length_c + 42 // 2.0 * thick_rim_a * length_b * length_c + 40 //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 41 //double vol = length_a * length_b * length_c + 42 // 2.0 * thick_rim_a * length_b * length_c + 43 43 // 2.0 * thick_rim_b * length_a * length_c + 44 44 // 2.0 * thick_rim_c * length_a * length_b; 45 45 46 46 // Scale sides by B 47 47 const double a_scaled = length_a / length_b; … … 101 101 // ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3); // correct FF : square of sum of phase factors 102 102 // This is the function called by csparallelepiped_analytical_2D_scaled, 103 // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 104 103 // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 104 105 105 // correct FF : sum of square of phase factors 106 106 inner_total += Gauss76Wt[j] * form * form; … … 136 136 double q, zhat, yhat, xhat; 137 137 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 138 const double qa = q*xhat; 139 const double qb = q*yhat; 140 const double qc = q*zhat; 138 141 139 142 // cspkernel in csparallelepiped recoded here … … 160 163 double tc = length_a + 2.0*thick_rim_c; 161 164 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 162 double siA = sas_sinx_x(0.5* q*length_a*xhat);163 double siB = sas_sinx_x(0.5* q*length_b*yhat);164 double siC = sas_sinx_x(0.5* q*length_c*zhat);165 double siAt = sas_sinx_x(0.5* q*ta*xhat);166 double siBt = sas_sinx_x(0.5* q*tb*yhat);167 double siCt = sas_sinx_x(0.5* q*tc*zhat);168 165 double siA = sas_sinx_x(0.5*length_a*qa); 166 double siB = sas_sinx_x(0.5*length_b*qb); 167 double siC = sas_sinx_x(0.5*length_c*qc); 168 double siAt = sas_sinx_x(0.5*ta*qa); 169 double siBt = sas_sinx_x(0.5*tb*qb); 170 double siCt = sas_sinx_x(0.5*tc*qc); 171 169 172 170 173 // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed … … 174 177 + drB*siA*(siBt-siB)*siC*V2 175 178 + drC*siA*siB*(siCt*siCt-siC)*V3); 176 179 177 180 return 1.0e-4 * f * f; 178 181 } -
sasmodels/models/cylinder.c
r592343f rb34fc77 1 double form_volume(double radius, double length);2 double fq(double q, double sn, double cn,double radius, double length);3 double orient_avg_1D(double q, double radius, double length);4 double Iq(double q, double sld, double solvent_sld, double radius, double length);5 double Iqxy(double qx, double qy, double sld, double solvent_sld,6 double radius, double length, double theta, double phi);7 8 1 #define INVALID(v) (v.radius<0 || v.length<0) 9 2 10 double form_volume(double radius, double length) 3 static double 4 form_volume(double radius, double length) 11 5 { 12 6 return M_PI*radius*radius*length; 13 7 } 14 8 15 double fq(double q, double sn, double cn, double radius, double length) 9 static double 10 fq(double qab, double qc, double radius, double length) 16 11 { 17 // precompute qr and qh to save time in the loop 18 const double qr = q*radius; 19 const double qh = q*0.5*length; 20 return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn); 12 return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 21 13 } 22 14 23 double orient_avg_1D(double q, double radius, double length) 15 static double 16 orient_avg_1D(double q, double radius, double length) 24 17 { 25 18 // translate a point in [-1,1] to a point in [0, pi/2] 26 19 const double zm = M_PI_4; 27 const double zb = M_PI_4; 20 const double zb = M_PI_4; 28 21 29 22 double total = 0.0; 30 23 for (int i=0; i<76 ;i++) { 31 const double alpha = Gauss76Z[i]*zm + zb; 32 double sn, cn; // slots to hold sincos function output 33 // alpha(theta,phi) the projection of the cylinder on the detector plane 34 SINCOS(alpha, sn, cn); 35 total += Gauss76Wt[i] * square( fq(q, sn, cn, radius, length) ) * sn; 24 const double theta = Gauss76Z[i]*zm + zb; 25 double sin_theta, cos_theta; // slots to hold sincos function output 26 // theta (theta,phi) the projection of the cylinder on the detector plane 27 SINCOS(theta , sin_theta, cos_theta); 28 const double form = fq(q*sin_theta, q*cos_theta, radius, length); 29 total += Gauss76Wt[i] * form * form * sin_theta; 36 30 } 37 31 // translate dx in [-1,1] to dx in [lower,upper] … … 39 33 } 40 34 41 double Iq(double q, 35 static double 36 Iq(double q, 42 37 double sld, 43 38 double solvent_sld, … … 49 44 } 50 45 51 52 doubleIqxy(double qx, double qy,46 static double 47 Iqxy(double qx, double qy, 53 48 double sld, 54 49 double solvent_sld, … … 60 55 double q, sin_alpha, cos_alpha; 61 56 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 62 //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha); 57 const double qab = q*sin_alpha; 58 const double qc = q*cos_alpha; 59 63 60 const double s = (sld-solvent_sld) * form_volume(radius, length); 64 const double form = fq(q , sin_alpha, cos_alpha, radius, length);61 const double form = fq(qab, qc, radius, length); 65 62 return 1.0e-4 * square(s * form); 66 63 } -
sasmodels/models/ellipsoid.c
r3b571ae r2a0b2b1 1 double form_volume(double radius_polar, double radius_equatorial); 2 double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 4 double radius_polar, double radius_equatorial, double theta, double phi); 5 6 double form_volume(double radius_polar, double radius_equatorial) 1 static double 2 form_volume(double radius_polar, double radius_equatorial) 7 3 { 8 4 return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 9 5 } 10 6 11 double Iq(double q, 7 static double 8 Iq(double q, 12 9 double sld, 13 10 double sld_solvent, … … 41 38 } 42 39 43 double Iqxy(double qx, double qy, 40 static double 41 Iqxy(double qx, double qy, 44 42 double sld, 45 43 double sld_solvent, … … 51 49 double q, sin_alpha, cos_alpha; 52 50 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 53 const double r = sqrt(square(radius_equatorial*sin_alpha) 54 + square(radius_polar*cos_alpha)); 55 const double f = sas_3j1x_x(q*r); 51 const double qab = q*sin_alpha; 52 const double qc = q*cos_alpha; 53 54 const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc)); 55 const double f = sas_3j1x_x(qr); 56 56 const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 57 57 58 58 return 1.0e-4 * square(f * s); 59 59 } 60 -
sasmodels/models/elliptical_cylinder.c
r61104c8 r2a0b2b1 1 double form_volume(double radius_minor, double r_ratio, double length); 2 double Iq(double q, double radius_minor, double r_ratio, double length, 3 double sld, double solvent_sld); 4 double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length, 5 double sld, double solvent_sld, double theta, double phi, double psi); 6 7 8 double 1 static double 9 2 form_volume(double radius_minor, double r_ratio, double length) 10 3 { … … 12 5 } 13 6 14 double7 static double 15 8 Iq(double q, double radius_minor, double r_ratio, double length, 16 9 double sld, double solvent_sld) … … 61 54 62 55 63 double56 static double 64 57 Iqxy(double qx, double qy, 65 58 double radius_minor, double r_ratio, double length, … … 69 62 double q, xhat, yhat, zhat; 70 63 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 64 const double qa = q*xhat; 65 const double qb = q*yhat; 66 const double qc = q*zhat; 71 67 72 68 // Compute: r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 73 69 // Given: radius_major = r_ratio * radius_minor 74 const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat));75 const double be = sas_2J1x_x(q *r);76 const double si = sas_sinx_x(q *zhat*0.5*length);70 const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb)); 71 const double be = sas_2J1x_x(qr); 72 const double si = sas_sinx_x(qc*0.5*length); 77 73 const double Aq = be * si; 78 74 const double delrho = sld - solvent_sld; -
sasmodels/models/fcc_paracrystal.c
r50beefe r7e0b281 1 double form_volume(double radius); 2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 3 double Iqxy(double qx, double qy, double dnn, 4 double d_factor, double radius,double sld, double solvent_sld, 5 double theta, double phi, double psi); 1 static double 2 fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 3 { 4 #if 0 // Equations as written in Matsuoka 5 const double a1 = ( qa + qb)/2.0; 6 const double a2 = (-qa + qc)/2.0; 7 const double a3 = (-qa + qb)/2.0; 8 #else 9 const double a1 = ( qa + qb)/2.0; 10 const double a2 = ( qa + qc)/2.0; 11 const double a3 = ( qb + qc)/2.0; 12 #endif 6 13 7 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 8 double _FCCeval(double Theta, double Phi, double temp1, double temp3); 14 // Numerator: (1 - exp(a)^2)^3 15 // => (-(exp(2a) - 1))^3 16 // => -expm1(2a)^3 17 // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 18 // => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 19 // => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 20 const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 21 const double exp_arg = exp(arg); 22 const double Zq = -cube(expm1(2.0*arg)) 23 / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 24 * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 25 * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 26 27 return Zq; 28 } 9 29 10 30 11 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 12 13 const double Da = d_factor*dnn; 14 const double temp1 = q*q*Da*Da; 15 const double temp3 = q*dnn; 16 17 double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 18 return(retVal); 31 // occupied volume fraction calculated from lattice symmetry and sphere radius 32 static double 33 fcc_volume_fraction(double radius, double dnn) 34 { 35 return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 19 36 } 20 37 21 double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 22 23 double result; 24 double sin_theta,cos_theta,sin_phi,cos_phi; 25 SINCOS(Theta, sin_theta, cos_theta); 26 SINCOS(Phi, sin_phi, cos_phi); 27 28 const double temp6 = sin_theta; 29 const double temp7 = sin_theta*sin_phi + cos_theta; 30 const double temp8 = -sin_theta*cos_phi + cos_theta; 31 const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi; 32 33 const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 34 result = cube(1.0-(temp10*temp10))*temp6 35 / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 36 * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 37 * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 38 39 return (result); 40 } 41 42 double form_volume(double radius){ 38 static double 39 form_volume(double radius) 40 { 43 41 return sphere_volume(radius); 44 42 } 45 43 46 44 47 double Iq(double q, double dnn,45 static double Iq(double q, double dnn, 48 46 double d_factor, double radius, 49 double sld, double solvent_sld){ 47 double sld, double solvent_sld) 48 { 49 // translate a point in [-1,1] to a point in [0, 2 pi] 50 const double phi_m = M_PI; 51 const double phi_b = M_PI; 52 // translate a point in [-1,1] to a point in [0, pi] 53 const double theta_m = M_PI_2; 54 const double theta_b = M_PI_2; 50 55 51 //Volume fraction calculated from lattice symmetry and sphere radius 52 const double s1 = dnn*sqrt(2.0); 53 const double latticescale = 4.0*sphere_volume(radius/s1); 56 double outer_sum = 0.0; 57 for(int i=0; i<150; i++) { 58 double inner_sum = 0.0; 59 const double theta = Gauss150Z[i]*theta_m + theta_b; 60 double sin_theta, cos_theta; 61 SINCOS(theta, sin_theta, cos_theta); 62 const double qc = q*cos_theta; 63 const double qab = q*sin_theta; 64 for(int j=0;j<150;j++) { 65 const double phi = Gauss150Z[j]*phi_m + phi_b; 66 double sin_phi, cos_phi; 67 SINCOS(phi, sin_phi, cos_phi); 68 const double qa = qab*cos_phi; 69 const double qb = qab*sin_phi; 70 const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 71 inner_sum += Gauss150Wt[j] * form; 72 } 73 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 74 outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 75 } 76 outer_sum *= theta_m; 77 const double Zq = outer_sum/(4.0*M_PI); 78 const double Pq = sphere_form(q, radius, sld, solvent_sld); 54 79 55 const double va = 0.0; 56 const double vb = 2.0*M_PI; 57 const double vaj = 0.0; 58 const double vbj = M_PI; 59 60 double summ = 0.0; 61 double answer = 0.0; 62 for(int i=0; i<150; i++) { 63 //setup inner integral over the ellipsoidal cross-section 64 double summj=0.0; 65 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy is phi 66 for(int j=0;j<150;j++) { 67 //20 gauss points for the inner integral 68 double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy is theta 69 double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi); 70 summj += yyy; 71 } 72 //now calculate the value of the inner integral 73 double answer = (vbj-vaj)/2.0*summj; 74 75 //now calculate outer integral 76 summ = summ+(Gauss150Wt[i] * answer); 77 } //final scaling is done at the end of the function, after the NT_FP64 case 78 79 answer = (vb-va)/2.0*summ; 80 answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 81 82 return answer; 80 return fcc_volume_fraction(radius, dnn) * Pq * Zq; 81 } 83 82 84 83 85 } 86 87 double Iqxy(double qx, double qy, 84 static double Iqxy(double qx, double qy, 88 85 double dnn, double d_factor, double radius, 89 86 double sld, double solvent_sld, … … 92 89 double q, zhat, yhat, xhat; 93 90 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 91 const double qa = q*xhat; 92 const double qb = q*yhat; 93 const double qc = q*zhat; 94 94 95 const double a1 = yhat + xhat; 96 const double a2 = xhat + zhat; 97 const double a3 = yhat + zhat; 98 const double qd = 0.5*q*dnn; 99 const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 100 const double tanh_qd = tanh(arg); 101 const double cosh_qd = cosh(arg); 102 const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 103 * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 104 * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 105 106 //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg); 107 108 const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 109 //the occupied volume of the lattice 110 const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 111 return lattice_scale * Fq; 95 q = sqrt(qa*qa + qb*qb + qc*qc); 96 const double Pq = sphere_form(q, radius, sld, solvent_sld); 97 const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); 98 return fcc_volume_fraction(radius, dnn) * Pq * Zq; 112 99 } -
sasmodels/models/hollow_cylinder.c
r592343f r2a0b2b1 1 1 double form_volume(double radius, double thickness, double length); 2 2 double Iq(double q, double radius, double thickness, double length, double sld, 3 3 double solvent_sld); 4 4 double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld, 5 5 double solvent_sld, double theta, double phi); 6 6 7 7 //#define INVALID(v) (v.radius_core >= v.radius) … … 14 14 } 15 15 16 17 16 static double 18 _ hollow_cylinder_kernel(double q,19 double radius, double thickness, double length , double sin_val, double cos_val)17 _fq(double qab, double qc, 18 double radius, double thickness, double length) 20 19 { 21 const double qs = q*sin_val; 22 const double lam1 = sas_2J1x_x((radius+thickness)*qs); 23 const double lam2 = sas_2J1x_x(radius*qs); 20 const double lam1 = sas_2J1x_x((radius+thickness)*qab); 21 const double lam2 = sas_2J1x_x(radius*qab); 24 22 const double gamma_sq = square(radius/(radius+thickness)); 25 //Note: lim_{thickness -> 0} psi = sas_J0(radius*q s)26 //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*q s)27 const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); 28 const double t2 = sas_sinx_x(0.5* q*length*cos_val);23 //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 24 //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 25 const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 26 const double t2 = sas_sinx_x(0.5*length*qc); 29 27 return psi*t2; 30 28 } … … 43 41 { 44 42 const double lower = 0.0; 45 const double upper = 1.0; 43 const double upper = 1.0; //limits of numerical integral 46 44 47 double summ = 0.0; 45 double summ = 0.0; //initialize intergral 48 46 for (int i=0;i<76;i++) { 49 const double cos_ val= 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );50 const double sin_ val = sqrt(1.0 - cos_val*cos_val);51 const double inter = _hollow_cylinder_kernel(q, radius, thickness, length,52 sin_val, cos_val);53 summ += Gauss76Wt[i] * inter * inter;47 const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 48 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 49 const double form = _fq(q*sin_theta, q*cos_theta, 50 radius, thickness, length); 51 summ += Gauss76Wt[i] * form * form; 54 52 } 55 53 … … 66 64 double q, sin_alpha, cos_alpha; 67 65 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 68 const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 69 sin_alpha, cos_alpha); 66 const double qab = q*sin_alpha; 67 const double qc = q*cos_alpha; 68 69 const double form = _fq(qab, qc, radius, thickness, length); 70 70 71 71 const double vol = form_volume(radius, thickness, length); 72 return _hollow_cylinder_scaling( Aq*Aq, solvent_sld-sld, vol);72 return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 73 73 } 74 -
sasmodels/models/parallelepiped.c
rd605080 r2a0b2b1 20 20 { 21 21 const double mu = 0.5 * q * length_b; 22 22 23 23 // Scale sides by B 24 24 const double a_scaled = length_a / length_b; 25 25 const double c_scaled = length_c / length_b; 26 26 27 27 // outer integral (with gauss points), integration limits = 0, 1 28 28 double outer_total = 0; //initialize integral … … 69 69 double q, xhat, yhat, zhat; 70 70 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 71 const double qa = q*xhat; 72 const double qb = q*yhat; 73 const double qc = q*zhat; 71 74 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 const double siA = sas_sinx_x(0.5*length_a*qa); 76 const double siB = sas_sinx_x(0.5*length_b*qb); 77 const double siC = sas_sinx_x(0.5*length_c*qc); 75 78 const double V = form_volume(length_a, length_b, length_c); 76 79 const double drho = (sld - solvent_sld); -
sasmodels/models/sc_paracrystal.c
r50beefe r7e0b281 1 double form_volume(double radius); 1 static double 2 sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 3 { 4 const double a1 = qa; 5 const double a2 = qb; 6 const double a3 = qc; 2 7 3 double Iq(double q, 4 double dnn, 5 double d_factor, 6 double radius, 7 double sphere_sld, 8 double solvent_sld); 8 // Numerator: (1 - exp(a)^2)^3 9 // => (-(exp(2a) - 1))^3 10 // => -expm1(2a)^3 11 // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 12 // => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 13 // => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 14 const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 15 const double exp_arg = exp(arg); 16 const double Zq = -cube(expm1(2.0*arg)) 17 / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 18 * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 19 * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 9 20 10 double Iqxy(double qx, double qy, 11 double dnn, 12 double d_factor, 13 double radius, 14 double sphere_sld, 15 double solvent_sld, 16 double theta, 17 double phi, 18 double psi); 21 return Zq; 22 } 19 23 20 double form_volume(double radius) 24 // occupied volume fraction calculated from lattice symmetry and sphere radius 25 static double 26 sc_volume_fraction(double radius, double dnn) 27 { 28 return sphere_volume(radius/dnn); 29 } 30 31 static double 32 form_volume(double radius) 21 33 { 22 34 return sphere_volume(radius); 23 35 } 24 36 25 static double 26 sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 37 38 static double Iq(double q, double dnn, 39 double d_factor, double radius, 40 double sld, double solvent_sld) 27 41 { 28 double cnt, snt; 29 SINCOS(theta, cnt, snt); 42 // translate a point in [-1,1] to a point in [0, 2 pi] 43 const double phi_m = M_PI_4; 44 const double phi_b = M_PI_4; 45 // translate a point in [-1,1] to a point in [0, pi] 46 const double theta_m = M_PI_4; 47 const double theta_b = M_PI_4; 30 48 31 double cnp, snp;32 SINCOS(phi, cnp, snp);33 49 34 double temp6 = snt; 35 double temp7 = -1.0*temp3*snt*cnp; 36 double temp8 = temp3*snt*snp; 37 double temp9 = temp3*cnt; 38 double result = temp6/((1.0-temp4*cos((temp7))+temp5)* 39 (1.0-temp4*cos((temp8))+temp5)* 40 (1.0-temp4*cos((temp9))+temp5)); 41 return (result); 50 double outer_sum = 0.0; 51 for(int i=0; i<150; i++) { 52 double inner_sum = 0.0; 53 const double theta = Gauss150Z[i]*theta_m + theta_b; 54 double sin_theta, cos_theta; 55 SINCOS(theta, sin_theta, cos_theta); 56 const double qc = q*cos_theta; 57 const double qab = q*sin_theta; 58 for(int j=0;j<150;j++) { 59 const double phi = Gauss150Z[j]*phi_m + phi_b; 60 double sin_phi, cos_phi; 61 SINCOS(phi, sin_phi, cos_phi); 62 const double qa = qab*cos_phi; 63 const double qb = qab*sin_phi; 64 const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 65 inner_sum += Gauss150Wt[j] * form; 66 } 67 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 68 outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 69 } 70 outer_sum *= theta_m; 71 const double Zq = outer_sum/M_PI_2; 72 const double Pq = sphere_form(q, radius, sld, solvent_sld); 73 74 return sc_volume_fraction(radius, dnn) * Pq * Zq; 42 75 } 43 76 44 static double45 sc_integrand(double dnn, double d_factor, double qq, double xx, double yy)46 {47 //Function to calculate integrand values for simple cubic structure48 77 49 double da = d_factor*dnn; 50 double temp1 = qq*qq*da*da; 51 double temp2 = cube(-expm1(-temp1)); 52 double temp3 = qq*dnn; 53 double temp4 = 2.0*exp(-0.5*temp1); 54 double temp5 = exp(-1.0*temp1); 55 56 double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 57 58 return(integrand); 59 } 60 61 double Iq(double q, 62 double dnn, 63 double d_factor, 64 double radius, 65 double sphere_sld, 66 double solvent_sld) 67 { 68 const double va = 0.0; 69 const double vb = M_PI_2; //orientation average, outer integral 70 71 double summ=0.0; 72 double answer=0.0; 73 74 for(int i=0;i<150;i++) { 75 //setup inner integral over the ellipsoidal cross-section 76 double summj=0.0; 77 double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; 78 for(int j=0;j<150;j++) { 79 //150 gauss points for the inner integral 80 double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0; 81 double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij); 82 summj += tmp; 83 } 84 //now calculate the value of the inner integral 85 answer = (vb-va)/2.0*summj; 86 87 //now calculate outer integral 88 double tmp = Gauss150Wt[i] * answer; 89 summ += tmp; 90 } //final scaling is done at the end of the function, after the NT_FP64 case 91 92 answer = (vb-va)/2.0*summ; 93 94 //Volume fraction calculated from lattice symmetry and sphere radius 95 // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3 96 const double latticeScale = sphere_volume(radius/dnn); 97 98 answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale; 99 100 return answer; 101 } 102 103 double Iqxy(double qx, double qy, 104 double dnn, 105 double d_factor, 106 double radius, 107 double sphere_sld, 108 double solvent_sld, 109 double theta, 110 double phi, 111 double psi) 78 static double Iqxy(double qx, double qy, 79 double dnn, double d_factor, double radius, 80 double sld, double solvent_sld, 81 double theta, double phi, double psi) 112 82 { 113 83 double q, zhat, yhat, xhat; 114 84 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 85 const double qa = q*xhat; 86 const double qb = q*yhat; 87 const double qc = q*zhat; 115 88 116 const double qd = q*dnn; 117 const double arg = 0.5*square(qd*d_factor); 118 const double tanh_qd = tanh(arg); 119 const double cosh_qd = cosh(arg); 120 const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd) 121 * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 122 * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 123 124 const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 125 //the occupied volume of the lattice 126 const double lattice_scale = sphere_volume(radius/dnn); 127 return lattice_scale * Fq; 89 q = sqrt(qa*qa + qb*qb + qc*qc); 90 const double Pq = sphere_form(q, radius, sld, solvent_sld); 91 const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); 92 return sc_volume_fraction(radius, dnn) * Pq * Zq; 128 93 } -
sasmodels/models/sc_paracrystal.py
r8f04da4 r9bc4882 161 161 [{'theta': 10.0, 'phi': 20, 'psi': 30.0}, (0.023, 0.045), 0.0177333171285], 162 162 ] 163 164 -
sasmodels/models/stacked_disks.c
r19f996b rb34fc77 1 1 static double stacked_disks_kernel( 2 double q, 2 double qab, 3 double qc, 3 4 double halfheight, 4 5 double thick_layer, … … 9 10 double layer_sld, 10 11 double solvent_sld, 11 double sin_alpha,12 double cos_alpha,13 12 double d) 14 13 … … 20 19 // zi is the dummy variable for the integration (x in Feigin's notation) 21 20 22 const double besarg1 = q*radius*sin_alpha;23 //const double besarg2 = q*radius*sin_alpha;21 const double besarg1 = radius*qab; 22 //const double besarg2 = radius*qab; 24 23 25 const double sinarg1 = q*halfheight*cos_alpha;26 const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha;24 const double sinarg1 = halfheight*qc; 25 const double sinarg2 = (halfheight+thick_layer)*qc; 27 26 28 27 const double be1 = sas_2J1x_x(besarg1); … … 43 42 44 43 // loop for the structure factor S(q) 45 double qd_cos_alpha = q*d*cos_alpha;44 double qd_cos_alpha = d*qc; 46 45 //d*cos_alpha is the projection of d onto q (in other words the component 47 46 //of d that is parallel to q. … … 84 83 double sin_alpha, cos_alpha; // slots to hold sincos function output 85 84 SINCOS(zi, sin_alpha, cos_alpha); 86 double yyy = stacked_disks_kernel(q ,85 double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 87 86 halfheight, 88 87 thick_layer, … … 93 92 layer_sld, 94 93 solvent_sld, 95 sin_alpha,96 cos_alpha,97 94 d); 98 95 summ += Gauss76Wt[i] * yyy * sin_alpha; … … 152 149 double phi) 153 150 { 154 int n_stacking = (int)(fp_n_stacking + 0.5);155 151 double q, sin_alpha, cos_alpha; 156 152 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 153 const double qab = q*sin_alpha; 154 const double qc = q*cos_alpha; 157 155 156 int n_stacking = (int)(fp_n_stacking + 0.5); 158 157 double d = 2.0 * thick_layer + thick_core; 159 158 double halfheight = 0.5*thick_core; 160 double answer = stacked_disks_kernel(q ,159 double answer = stacked_disks_kernel(qab, qc, 161 160 halfheight, 162 161 thick_layer, … … 167 166 layer_sld, 168 167 solvent_sld, 169 sin_alpha,170 cos_alpha,171 168 d); 172 169 -
sasmodels/models/triaxial_ellipsoid.c
r68dd6a9 r2a0b2b1 1 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar);2 double Iq(double q, double sld, double sld_solvent,3 double radius_equat_minor, double radius_equat_major, double radius_polar);4 double Iqxy(double qx, double qy, double sld, double sld_solvent,5 double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi);6 7 1 //#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 8 2 9 10 doubleform_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)3 static double 4 form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 11 5 { 12 6 return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 13 7 } 14 8 15 double Iq(double q, 9 static double 10 Iq(double q, 16 11 double sld, 17 12 double sld_solvent, … … 45 40 // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 46 41 const double fqsq = outer/4.0; // = outer*um*zm*8.0/(4.0*M_PI); 47 const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 48 return 1.0e-4 * s * s * fqsq; 42 const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 43 const double drho = (sld - sld_solvent); 44 return 1.0e-4 * square(vol*drho) * fqsq; 49 45 } 50 46 51 double Iqxy(double qx, double qy, 47 static double 48 Iqxy(double qx, double qy, 52 49 double sld, 53 50 double sld_solvent, … … 61 58 double q, xhat, yhat, zhat; 62 59 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 60 const double qa = q*xhat; 61 const double qb = q*yhat; 62 const double qc = q*zhat; 63 63 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); 68 const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 64 const double qr = sqrt(square(radius_equat_minor*qa) 65 + square(radius_equat_major*qb) 66 + square(radius_polar*qc)); 67 const double fq = sas_3j1x_x(qr); 68 const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 69 const double drho = (sld - sld_solvent); 69 70 70 return 1.0e-4 * square( s* fq);71 return 1.0e-4 * square(vol * drho * fq); 71 72 } 72 -
sasmodels/weights.py
r41e7f2e rd4c33d6 55 55 """ 56 56 sigma = self.width * center if relative else self.width 57 if not relative: 58 # For orientation, the jitter is relative to 0 not the angle 59 #center = 0 60 pass 57 61 if sigma == 0 or self.npts < 2: 58 62 if lb <= center <= ub:
Note: See TracChangeset
for help on using the changeset viewer.