Changeset 110f69c in sasmodels
- Timestamp:
- Nov 28, 2017 11:17:57 AM (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:
- e65c3ba
- Parents:
- 167d0f1
- Location:
- sasmodels
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/compare.py
r3bfd924 r110f69c 150 150 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True 151 151 152 # list of math functions for use in evaluating parameters 153 MATH = dict((k,getattr(math, k)) for k in dir(math) if not k.startswith('_')) 152 def build_math_context(): 153 # type: () -> Dict[str, Callable] 154 """build dictionary of functions from math module""" 155 return dict((k, getattr(math, k)) 156 for k in dir(math) if not k.startswith('_')) 157 158 #: list of math functions for use in evaluating parameters 159 MATH = build_math_context() 154 160 155 161 # CRUFT python 2.6 … … 231 237 pass 232 238 233 def __exit__(self, exc_type, exc_value, t raceback):239 def __exit__(self, exc_type, exc_value, tb): 234 240 # type: (Any, BaseException, Any) -> None 235 241 # TODO: better typing for __exit__ method … … 374 380 375 381 def _random_pd(model_info, pars): 382 # type: (ModelInfo, Dict[str, float]) -> None 383 """ 384 Generate a random dispersity distribution for the model. 385 386 1% no shape dispersity 387 85% single shape parameter 388 13% two shape parameters 389 1% three shape parameters 390 391 If oriented, then put dispersity in theta, add phi and psi dispersity 392 with 10% probability for each. 393 """ 376 394 pd = [p for p in model_info.parameters.kernel_parameters if p.polydisperse] 377 395 pd_volume = [] … … 444 462 value = pars[p.name] 445 463 if p.units == 'Ang' and value > maxdim: 446 pars[p.name] = maxdim*10**np.random.uniform(-3, 0)464 pars[p.name] = maxdim*10**np.random.uniform(-3, 0) 447 465 448 466 def constrain_pars(model_info, pars): … … 490 508 if pars['radius'] < pars['thick_string']: 491 509 pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius'] 492 pass493 510 494 511 elif name == 'rpa': … … 608 625 return pars 609 626 627 # TODO: remove support for sasview 3.x models 610 628 def eval_sasview(model_info, data): 611 629 # type: (Modelinfo, Data) -> Calculator … … 621 639 from sas.models.dispersion_models import models as dispersers 622 640 623 def get_model_class(name):641 def _get_model_class(name): 624 642 # type: (str) -> "sas.models.BaseComponent" 625 643 #print("new",sorted(_pars.items())) … … 641 659 composition_type, parts = model_info.composition 642 660 if composition_type == 'product': 643 P, S = [ get_model_class(revert_name(p))() for p in parts]661 P, S = [_get_model_class(revert_name(p))() for p in parts] 644 662 model = [MultiplicationModel(P, S)] 645 663 else: … … 649 667 if old_name is None: 650 668 raise ValueError("model %r does not exist in old sasview" 651 % model_info.id)652 ModelClass = get_model_class(old_name)669 % model_info.id) 670 ModelClass = _get_model_class(old_name) 653 671 model = [ModelClass()] 654 672 model[0].disperser_handles = {} … … 847 865 # print a separate seed for each dataset for better reproducibility 848 866 new_seed = np.random.randint(1000000) 849 print("Set %d uses -random=%i"%(k+1, new_seed))867 print("Set %d uses -random=%i"%(k+1, new_seed)) 850 868 np.random.seed(new_seed) 851 869 opts['pars'] = parse_pars(opts, maxdim=maxdim) … … 868 886 def run_models(opts, verbose=False): 869 887 # type: (Dict[str, Any]) -> Dict[str, Any] 888 """ 889 Process a parameter set, return calculation results and times. 890 """ 870 891 871 892 base, comp = opts['engines'] … … 941 962 def plot_models(opts, result, limits=None, setnum=0): 942 963 # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 964 """ 965 Plot the results from :func:`run_model`. 966 """ 943 967 import matplotlib.pyplot as plt 944 968 … … 987 1011 errview = 'linear' 988 1012 if 0: # 95% cutoff 989 sorted = np.sort(err.flatten())990 cutoff = sorted [int(sorted.size*0.95)]1013 sorted_err = np.sort(err.flatten()) 1014 cutoff = sorted_err[int(sorted_err.size*0.95)] 991 1015 err[err > cutoff] = cutoff 992 1016 #err,errstr = base/comp,"ratio" … … 1106 1130 1107 1131 INTEGER_RE = re.compile("^[+-]?[1-9][0-9]*$") 1108 def isnumber(str): 1109 match = FLOAT_RE.match(str) 1110 isfloat = (match and not str[match.end():]) 1111 return isfloat or INTEGER_RE.match(str) 1132 def isnumber(s): 1133 # type: (str) -> bool 1134 """Return True if string contains an int or float""" 1135 match = FLOAT_RE.match(s) 1136 isfloat = (match and not s[match.end():]) 1137 return isfloat or INTEGER_RE.match(s) 1112 1138 1113 1139 # For distinguishing pairs of models for comparison … … 1314 1340 1315 1341 def set_spherical_integration_parameters(opts, steps): 1342 # type: (Dict[str, Any], int) -> None 1316 1343 """ 1317 1344 Set integration parameters for spherical integration over the entire … … 1337 1364 'psi_pd_type=rectangle', 1338 1365 ]) 1339 pass1340 1366 1341 1367 def parse_pars(opts, maxdim=np.inf): 1368 # type: (Dict[str, Any], float) -> Tuple[Dict[str, float], Dict[str, float]] 1369 """ 1370 Generate a parameter set. 1371 1372 The default values come from the model, or a randomized model if a seed 1373 value is given. Next, evaluate any parameter expressions, constraining 1374 the value of the parameter within and between models. If *maxdim* is 1375 given, limit parameters with units of Angstrom to this value. 1376 1377 Returns a pair of parameter dictionaries for base and comparison models. 1378 """ 1342 1379 model_info, model_info2 = opts['info'] 1343 1380 … … 1378 1415 print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 1379 1416 return None 1380 v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v, v)1417 v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v, v) 1381 1418 if v1 and k in pars: 1382 1419 presets[k] = float(v1) if isnumber(v1) else v1 … … 1434 1471 html = make_html(info) 1435 1472 path = os.path.dirname(info.filename) 1436 url = "file://" +path.replace("\\","/")[2:]+"/"1473 url = "file://" + path.replace("\\", "/")[2:] + "/" 1437 1474 rst2html.view_html_qtapp(html, url) 1438 1475 … … 1458 1495 frame.panel.Layout() 1459 1496 frame.panel.aui.Split(0, wx.TOP) 1460 def reset_parameters(event):1497 def _reset_parameters(event): 1461 1498 model.revert_values() 1462 1499 signal.update_parameters(problem) 1463 frame.Bind(wx.EVT_TOOL, reset_parameters, frame.ToolBar.GetToolByPos(1))1500 frame.Bind(wx.EVT_TOOL, _reset_parameters, frame.ToolBar.GetToolByPos(1)) 1464 1501 if is_mac: frame.Show() 1465 1502 # If running withing an app, start the main loop … … 1504 1541 1505 1542 def revert_values(self): 1543 # type: () -> None 1544 """ 1545 Restore starting values of the parameters. 1546 """ 1506 1547 for k, v in self.starting_values.items(): 1507 1548 self.pars[k].value = v 1508 1549 1509 1550 def model_update(self): 1551 # type: () -> None 1552 """ 1553 Respond to signal that model parameters have been changed. 1554 """ 1510 1555 pass 1511 1556 -
sasmodels/details.py
rce99754 r110f69c 23 23 # CRUFT: np.meshgrid requires multiple vectors 24 24 def meshgrid(*args): 25 """See docs from a recent version of numpy""" 25 26 if len(args) > 1: 26 27 return np.meshgrid(*args) … … 231 232 npars = kernel.info.parameters.npars 232 233 nvalues = kernel.info.parameters.nvalues 233 scalars = [value for value, dispersity,weight in mesh]234 scalars = [value for value, _dispersity, _weight in mesh] 234 235 # skipping scale and background when building values and weights 235 values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ())236 _values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 236 237 #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights) 237 238 length = np.array([len(w) for w in weights]) … … 287 288 mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 288 289 mag = mag.reshape(-1, 3) 289 scale = mag[:, 0]290 scale = mag[:, 0] 290 291 if np.any(scale): 291 292 theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py
r129bdc4 r110f69c 6 6 core-shell scattering length density profile. Thus this is a variation 7 7 of the core-shell bicelle model, but with an elliptical cylinder for the core. 8 In this version the "rim" or "belt" does NOT extend the full length of the particle,9 but has the same length as the core. 10 Outer shells on the rims and flat ends may be of different thicknesses and 11 scattering lengthdensities. The form factor is normalized by the total particle volume.8 In this version the "rim" or "belt" does NOT extend the full length of 9 the particle, but has the same length as the core. Outer shells on the 10 rims and flat ends may be of different thicknesses and scattering length 11 densities. The form factor is normalized by the total particle volume. 12 12 This version includes an approximate "interfacial roughness". 13 13 … … 15 15 .. figure:: img/core_shell_bicelle_belt_geometry.png 16 16 17 Schematic cross-section of bicelle with belt. Note however that the model here18 calculates for rectangular, not curved, rims as shown below.17 Schematic cross-section of bicelle with belt. Note however that the model 18 here calculates for rectangular, not curved, rims as shown below. 19 19 20 20 .. figure:: img/core_shell_bicelle_belt_parameters.png 21 21 22 Cross section of model used here. Users will have 23 to decide how to distribute "heads" and "tails" between the rim, face 22 Cross section of model used here. Users will have 23 to decide how to distribute "heads" and "tails" between the rim, face 24 24 and core regions in order to estimate appropriate starting parameters. 25 25 … … 30 30 .. math:: 31 31 32 \rho(r) = 33 \begin{cases} 32 \rho(r) = 33 \begin{cases} 34 34 &\rho_c \text{ for } 0 \lt r \lt R; -L/2 \lt z\lt L/2 \\[1.5ex] 35 &\rho_f \text{ for } 0 \lt r \lt R; -(L/2 +t_ {face}) \lt z\lt -L/2;36 L/2 \lt z\lt (L/2+t_ {face}) \\[1.5ex]37 &\rho_r\text{ for } R \lt r \lt R+t_ {rim}; -L/2 \lt z\lt L/235 &\rho_f \text{ for } 0 \lt r \lt R; -(L/2 +t_\text{face}) \lt z\lt -L/2; 36 L/2 \lt z\lt (L/2+t_\text{face}) \\[1.5ex] 37 &\rho_r\text{ for } R \lt r \lt R+t_\text{rim}; -L/2 \lt z\lt L/2 38 38 \end{cases} 39 39 40 40 The form factor for the bicelle is calculated in cylindrical coordinates, where 41 $\alpha$ is the angle between the $Q$ vector and the cylinder axis, and $\psi$ is the angle42 for the ellipsoidal cross section core, to give:41 $\alpha$ is the angle between the $Q$ vector and the cylinder axis, and $\psi$ 42 is the angle for the ellipsoidal cross section core, to give: 43 43 44 44 .. math:: 45 45 46 I(Q,\alpha,\psi) = \frac{\text{scale}}{V_t} \cdot 47 F(Q,\alpha, \psi)^2.sin(\alpha).exp\left \{ -\frac{1}{2}Q^2\sigma^2 \right \} + \text{background} 46 I(Q,\alpha,\psi) = \frac{\text{scale}}{V_t} 47 \cdot F(Q,\alpha, \psi)^2 \cdot \sin(\alpha) 48 \cdot\exp\left\{ -\frac{1}{2}Q^2\sigma^2 \right\} + \text{background} 48 49 49 where a numerical integration of $F(Q,\alpha, \psi)^2.sin(\alpha)$ is carried out over \alpha and \psi for: 50 where a numerical integration of $F(Q,\alpha, \psi)^2\sin(\alpha)$ is 51 carried out over $\alpha$ and $\psi$ for: 50 52 51 53 .. math:: 52 54 53 \begin{align} 54 F(Q,\alpha,\psi) = &\bigg[ 55 (\rho_c -\rho_r - \rho_f + \rho_s) V_c \frac{2J_1(QR'sin \alpha)}{QR'sin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\ 56 &+(\rho_f - \rho_s) V_{c+f} \frac{2J_1(QR'sin\alpha)}{QR'sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha} \\ 57 &+(\rho_r - \rho_s) V_{c+r} \frac{2J_1(Q(R'+t_r)sin\alpha)}{Q(R'+t_r)sin\alpha}\frac{sin(Q(L/2)cos\alpha)}{Q(L/2)cos\alpha} 55 F(Q,\alpha,\psi) = &\bigg[ 56 (\rho_c -\rho_r - \rho_f + \rho_s) V_c 57 \frac{2J_1(QR'\sin \alpha)}{QR'\sin\alpha} 58 \frac{\sin(QL\cos\alpha/2)}{Q(L/2)\cos\alpha} \\ 59 &+(\rho_f - \rho_s) V_{c+f} 60 \frac{2J_1(QR'\sin\alpha)}{QR'\sin\alpha} 61 \frac{\sin(Q(L/2+t_f)\cos\alpha)}{Q(L/2+t_f)\cos\alpha} \\ 62 &+(\rho_r - \rho_s) V_{c+r} 63 \frac{2J_1(Q(R'+t_r)\sin\alpha)}{Q(R'+t_r)\sin\alpha} 64 \frac{\sin(Q(L/2)\cos\alpha)}{Q(L/2)\cos\alpha} 58 65 \bigg] 59 \end{align}60 66 61 67 where … … 63 69 .. math:: 64 70 65 R'=\frac{R}{\sqrt{2}}\sqrt{(1+X_{core}^{2}) + (1-X_{core}^{2})cos(\psi)} 66 67 68 and $V_t = \pi.(R+t_r)(Xcore.R+t_r).L + 2.\pi.Xcore.R^2.t_f$ is the total volume of the bicelle, 69 $V_c = \pi.Xcore.R^2.L$ the volume of the core, $V_{c+f} = \pi.Xcore.R^2.(L+2.t_f)$ 70 the volume of the core plus the volume of the faces, $V_{c+r} = \pi.(R+t_r)(Xcore.R+t_r),.L$ 71 the volume of the core plus the rim, $R$ is the radius 72 of the core, $Xcore$ is the axial ratio of the core, $L$ the length of the core, 73 $t_f$ the thickness of the face, $t_r$ the thickness of the rim and $J_1$ the usual 74 first order bessel function. The core has radii $R$ and $Xcore.R$ so is circular, 75 as for the core_shell_bicelle model, for $Xcore$ =1. Note that you may need to 76 limit the range of $Xcore$, especially if using the Monte-Carlo algorithm, as 77 setting radius to $R/Xcore$ and axial ratio to $1/Xcore$ gives an equivalent solution! 71 R' = \frac{R}{\sqrt{2}} 72 \sqrt{(1+X_\text{core}^{2}) + (1-X_\text{core}^{2})\cos(\psi)} 78 73 79 An approximation for the effects of "Gaussian interfacial roughness" $\sigma$ is included, 80 by multiplying $I(Q)$ by $exp\left \{ -\frac{1}{2}Q^2\sigma^2 \right \}$ . This 81 applies, in some way, to all interfaces in the model not just the external ones. 82 (Note that for a one dimensional system convolution of the scattering length density profile 83 with a Gaussian of standard deviation $\sigma$ does exactly this multiplication.) Leave 84 $\sigma$ set to zero for the usual sharp interfaces. 74 75 and $V_t = \pi (R+t_r)(X_\text{core} R+t_r) L + 2 \pi X_\text{core} R^2 t_f$ is 76 the total volume of the bicelle, $V_c = \pi X_\text{core} R^2 L$ the volume of 77 the core, $V_{c+f} = \pi X_\text{core} R^2 (L+2 t_f)$ the volume of the core 78 plus the volume of the faces, $V_{c+r} = \pi (R+t_r)(X_\text{core} R+t_r) L$ 79 the volume of the core plus the rim, $R$ is the radius of the core, 80 $X_\text{core}$ is the axial ratio of the core, $L$ the length of the core, 81 $t_f$ the thickness of the face, $t_r$ the thickness of the rim and $J_1$ the 82 usual first order bessel function. The core has radii $R$ and $X_\text{core} R$ 83 so is circular, as for the core_shell_bicelle model, for $X_\text{core}=1$. 84 Note that you may need to limit the range of $X_\text{core}$, especially if 85 using the Monte-Carlo algorithm, as setting radius to $R/X_\text{core}$ and 86 axial ratio to $1/X_\text{core}$ gives an equivalent solution! 87 88 An approximation for the effects of "Gaussian interfacial roughness" $\sigma$ 89 is included, by multiplying $I(Q)$ by 90 $\exp\left \{ -\frac{1}{2}Q^2\sigma^2 \right \}$. This applies, in some way, to 91 all interfaces in the model not just the external ones. (Note that for a one 92 dimensional system convolution of the scattering length density profile with 93 a Gaussian of standard deviation $\sigma$ does exactly this multiplication.) 94 Leave $\sigma$ set to zero for the usual sharp interfaces. 85 95 86 96 The output of the 1D scattering intensity function for randomly oriented 87 97 bicelles is then given by integrating over all possible $\alpha$ and $\psi$. 88 98 89 For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 90 for further details of the calculation and angular dispersions see :ref:`orientation` . 99 For oriented bicelles the *theta*, *phi* and *psi* orientation parameters 100 will appear when fitting 2D data, for further details of the calculation 101 and angular dispersions see :ref:`orientation` . 91 102 92 103 .. figure:: img/elliptical_cylinder_angle_definition.png 93 104 94 Definition of the angles for the oriented core_shell_bicelle_elliptical particles. 105 Definition of the angles for the oriented core_shell_bicelle_elliptical 106 particles. 95 107 96 108 … … 115 127 description = """ 116 128 core_shell_bicelle_elliptical_belt_rough 117 Elliptical cylinder core, optional shell on the two flat faces, and "belt" shell of 129 Elliptical cylinder core, optional shell on the two flat faces, and "belt" shell of 118 130 uniform thickness on its rim (in this case NOT extending around the end faces). 119 131 with approximate interfacial roughness. … … 173 185 174 186 [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0, 175 'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0},176 0.015, 189.328],177 #[{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ],178 187 'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0}, 188 0.015, 189.328], 189 #[{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 190 ] 179 191 180 192 del qx, qy # not necessary to delete, but cleaner -
sasmodels/models/ellipsoid.py
reda8b30 r110f69c 53 53 r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2} 54 54 55 For 2d data from oriented ellipsoids the direction of the rotation axis of 56 the ellipsoid is defined using two angles $\theta$ and $\phi$ as for the 55 For 2d data from oriented ellipsoids the direction of the rotation axis of 56 the ellipsoid is defined using two angles $\theta$ and $\phi$ as for the 57 57 :ref:`cylinder orientation figure <cylinder-angle-definition>`. 58 58 For the ellipsoid, $\theta$ is the angle between the rotational axis 59 59 and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$ 60 in the $xy$ plane, for further details of the calculation and angular 60 in the $xy$ plane, for further details of the calculation and angular 61 61 dispersions see :ref:`orientation` . 62 62 … … 209 209 qy = q*sin(pi/6.0) 210 210 tests = [[{}, 0.05, 54.8525847025], 211 [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026 211 [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026], 212 212 ] 213 213 del qx, qy # not necessary to delete, but cleaner -
sasmodels/models/stacked_disks.py
reda8b30 r110f69c 74 74 the layers. 75 75 76 2d scattering from oriented stacks is calculated in the same way as for cylinders, 77 for further details of the calculation and angular dispersions see :ref:`orientation` . 76 2d scattering from oriented stacks is calculated in the same way as for 77 cylinders, for further details of the calculation and angular dispersions 78 see :ref:`orientation`. 78 79 79 80 .. figure:: img/cylinder_angle_definition.png 80 81 81 82 Angles $\theta$ and $\phi$ orient the stack of discs relative 82 to the beam line coordinates, where the beam is along the $z$ axis. Rotation $\theta$, initially 83 in the $xz$ plane, is carried out first, then rotation $\phi$ about the $z$ axis. Orientation distributions 84 are described as rotations about two perpendicular axes $\delta_1$ and $\delta_2$ 85 in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes. 83 to the beam line coordinates, where the beam is along the $z$ axis. 84 Rotation $\theta$, initially in the $xz$ plane, is carried out first, 85 then rotation $\phi$ about the $z$ axis. Orientation distributions are 86 described as rotations about two perpendicular axes $\delta_1$ and 87 $\delta_2$ in the frame of the cylinder itself, which when 88 $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes. 86 89 87 90 -
sasmodels/special.py
r706f466 r110f69c 105 105 106 106 p1evl(x, c, n): 107 Evaluat ion ofnormalized polynomial $p(x) = x^n + \sum_{i=0}^{n-1} c_i x^i$107 Evaluate normalized polynomial $p(x) = x^n + \sum_{i=0}^{n-1} c_i x^i$ 108 108 using Horner's method so it is faster and more accurate. 109 109 … … 155 155 .. math:: 156 156 157 \text{Si}(x) \sim \frac{\pi}{2} 158 - \frac{\cos(x)}{x}\left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} - \frac{6!}{x^6} \right) 159 - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 157 \text{Si}(x) \sim \frac{\pi}{2} 158 - \frac{\cos(x)}{x} 159 \left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} - \frac{6!}{x^6} \right) 160 - \frac{\sin(x)}{x} 161 \left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 160 162 161 163 For small arguments, … … 207 209 # erf, erfc, tgamma, lgamma **do not use** 208 210 209 # Rotations of q210 def ORIENT_SYMMETRIC(qx, qy, theta, phi):211 q = sqrt(qx*qx + qy*qy)212 q = sqrt(qx*qx + qy*qy)213 sin_phi, cos_phi = sin(radians(phi)), cos(radians(phi))214 cn = (cn*qx + sn*qy)/q * sin(radians(theta))215 cn[q==0.] = 1.216 sn = sqrt(1 - cn**2)217 return q, cn, sn218 219 def ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi):220 q = sqrt(qx*qx + qy*qy)221 qxhat = qx/q222 qyhat = qy/q223 sin_theta, cos_theta = sin(radians(theta)), cos(radians(theta))224 sin_phi, cos_phi = sin(radians(phi)), cos(radians(phi))225 sin_psi, cos_psi = sin(radians(psi)), cos(radians(psi))226 227 xhat = (qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi)228 + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi))229 yhat = (qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi)230 + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi))231 zhat = (qxhat*(-sin_theta*cos_phi)232 + qyhat*(-sin_theta*sin_phi))233 return q, xhat, yhat, zhat234 235 236 211 # non-standard constants and functions 237 238 212 M_PI_180, M_4PI_3 = M_PI/180, 4*M_PI/3 239 213 … … 262 236 from scipy.special import jn as sas_JN 263 237 264 # missing sas_Si265 238 def sas_Si(x): 266 239 return scipy.special.sici(x)[0] … … 289 262 with np.errstate(all='ignore'): 290 263 retvalue = 2*sas_J1(x)/x 291 retvalue[x ==0] = 1.264 retvalue[x == 0] = 1. 292 265 return retvalue 293 266 … … 500 473 501 474 Gauss150Z = np.array([ 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 475 -0.9998723404457334, 476 -0.9993274305065947, 477 -0.9983473449340834, 478 -0.9969322929775997, 479 -0.9950828645255290, 480 -0.9927998590434373, 481 -0.9900842691660192, 482 -0.9869372772712794, 483 -0.9833602541697529, 484 -0.9793547582425894, 485 -0.9749225346595943, 486 -0.9700655145738374, 487 -0.9647858142586956, 488 -0.9590857341746905, 489 -0.9529677579610971, 490 -0.9464345513503147, 491 -0.9394889610042837, 492 -0.9321340132728527, 493 -0.9243729128743136, 494 -0.9162090414984952, 495 -0.9076459563329236, 496 -0.8986873885126239, 497 -0.8893372414942055, 498 -0.8795995893549102, 499 -0.8694786750173527, 500 -0.8589789084007133, 501 -0.8481048644991847, 502 -0.8368612813885015, 503 -0.8252530581614230, 504 -0.8132852527930605, 505 -0.8009630799369827, 506 -0.7882919086530552, 507 -0.7752772600680049, 508 -0.7619248049697269, 509 -0.7482403613363824, 510 -0.7342298918013638, 511 -0.7198995010552305, 512 -0.7052554331857488, 513 -0.6903040689571928, 514 -0.6750519230300931, 515 -0.6595056411226444, 516 -0.6436719971150083, 517 -0.6275578900977726, 518 -0.6111703413658551, 519 -0.5945164913591590, 520 -0.5776035965513142, 521 -0.5604390262878617, 522 -0.5430302595752546, 523 -0.5253848818220803, 524 -0.5075105815339176, 525 -0.4894151469632753, 526 -0.4711064627160663, 527 -0.4525925063160997, 528 -0.4338813447290861, 529 -0.4149811308476706, 530 -0.3959000999390257, 531 -0.3766465660565522, 532 -0.3572289184172501, 533 -0.3376556177463400, 534 -0.3179351925907259, 535 -0.2980762356029071, 536 -0.2780873997969574, 537 -0.2579773947782034, 538 -0.2377549829482451, 539 -0.2174289756869712, 540 -0.1970082295132342, 541 -0.1765016422258567, 542 -0.1559181490266516, 543 -0.1352667186271445, 544 -0.1145563493406956, 545 -0.0937960651617229, 546 -0.0729949118337358, 547 -0.0521619529078925, 548 -0.0313062657937972, 549 -0.0104369378042598, 550 0.0104369378042598, 551 0.0313062657937972, 552 0.0521619529078925, 553 0.0729949118337358, 554 0.0937960651617229, 555 0.1145563493406956, 556 0.1352667186271445, 557 0.1559181490266516, 558 0.1765016422258567, 559 0.1970082295132342, 560 0.2174289756869712, 561 0.2377549829482451, 562 0.2579773947782034, 563 0.2780873997969574, 564 0.2980762356029071, 565 0.3179351925907259, 566 0.3376556177463400, 567 0.3572289184172501, 568 0.3766465660565522, 569 0.3959000999390257, 570 0.4149811308476706, 571 0.4338813447290861, 572 0.4525925063160997, 573 0.4711064627160663, 574 0.4894151469632753, 575 0.5075105815339176, 576 0.5253848818220803, 577 0.5430302595752546, 578 0.5604390262878617, 579 0.5776035965513142, 580 0.5945164913591590, 581 0.6111703413658551, 582 0.6275578900977726, 583 0.6436719971150083, 584 0.6595056411226444, 585 0.6750519230300931, 586 0.6903040689571928, 587 0.7052554331857488, 588 0.7198995010552305, 589 0.7342298918013638, 590 0.7482403613363824, 591 0.7619248049697269, 592 0.7752772600680049, 593 0.7882919086530552, 594 0.8009630799369827, 595 0.8132852527930605, 596 0.8252530581614230, 597 0.8368612813885015, 598 0.8481048644991847, 599 0.8589789084007133, 600 0.8694786750173527, 601 0.8795995893549102, 602 0.8893372414942055, 603 0.8986873885126239, 604 0.9076459563329236, 605 0.9162090414984952, 606 0.9243729128743136, 607 0.9321340132728527, 608 0.9394889610042837, 609 0.9464345513503147, 610 0.9529677579610971, 611 0.9590857341746905, 612 0.9647858142586956, 613 0.9700655145738374, 614 0.9749225346595943, 615 0.9793547582425894, 616 0.9833602541697529, 617 0.9869372772712794, 618 0.9900842691660192, 619 0.9927998590434373, 620 0.9950828645255290, 621 0.9969322929775997, 622 0.9983473449340834, 623 0.9993274305065947, 624 0.9998723404457334 652 625 ]) 653 626 654 627 Gauss150Wt = np.array([ 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 628 0.0003276086705538, 629 0.0007624720924706, 630 0.0011976474864367, 631 0.0016323569986067, 632 0.0020663664924131, 633 0.0024994789888943, 634 0.0029315036836558, 635 0.0033622516236779, 636 0.0037915348363451, 637 0.0042191661429919, 638 0.0046449591497966, 639 0.0050687282939456, 640 0.0054902889094487, 641 0.0059094573005900, 642 0.0063260508184704, 643 0.0067398879387430, 644 0.0071507883396855, 645 0.0075585729801782, 646 0.0079630641773633, 647 0.0083640856838475, 648 0.0087614627643580, 649 0.0091550222717888, 650 0.0095445927225849, 651 0.0099300043714212, 652 0.0103110892851360, 653 0.0106876814158841, 654 0.0110596166734735, 655 0.0114267329968529, 656 0.0117888704247183, 657 0.0121458711652067, 658 0.0124975796646449, 659 0.0128438426753249, 660 0.0131845093222756, 661 0.0135194311690004, 662 0.0138484622795371, 663 0.0141714592928592, 664 0.0144882814685445, 665 0.0147987907597169, 666 0.0151028518701744, 667 0.0154003323133401, 668 0.0156911024699895, 669 0.0159750356447283, 670 0.0162520081211971, 671 0.0165218992159766, 672 0.0167845913311726, 673 0.0170399700056559, 674 0.0172879239649355, 675 0.0175283451696437, 676 0.0177611288626114, 677 0.0179861736145128, 678 0.0182033813680609, 679 0.0184126574807331, 680 0.0186139107660094, 681 0.0188070535331042, 682 0.0189920016251754, 683 0.0191686744559934, 684 0.0193369950450545, 685 0.0194968900511231, 686 0.0196482898041878, 687 0.0197911283358190, 688 0.0199253434079123, 689 0.0200508765398072, 690 0.0201676730337687, 691 0.0202756819988200, 692 0.0203748563729175, 693 0.0204651529434560, 694 0.0205465323660984, 695 0.0206189591819181, 696 0.0206824018328499, 697 0.0207368326754401, 698 0.0207822279928917, 699 0.0208185680053983, 700 0.0208458368787627, 701 0.0208640227312962, 702 0.0208731176389954, 703 0.0208731176389954, 704 0.0208640227312962, 705 0.0208458368787627, 706 0.0208185680053983, 707 0.0207822279928917, 708 0.0207368326754401, 709 0.0206824018328499, 710 0.0206189591819181, 711 0.0205465323660984, 712 0.0204651529434560, 713 0.0203748563729175, 714 0.0202756819988200, 715 0.0201676730337687, 716 0.0200508765398072, 717 0.0199253434079123, 718 0.0197911283358190, 719 0.0196482898041878, 720 0.0194968900511231, 721 0.0193369950450545, 722 0.0191686744559934, 723 0.0189920016251754, 724 0.0188070535331042, 725 0.0186139107660094, 726 0.0184126574807331, 727 0.0182033813680609, 728 0.0179861736145128, 729 0.0177611288626114, 730 0.0175283451696437, 731 0.0172879239649355, 732 0.0170399700056559, 733 0.0167845913311726, 734 0.0165218992159766, 735 0.0162520081211971, 736 0.0159750356447283, 737 0.0156911024699895, 738 0.0154003323133401, 739 0.0151028518701744, 740 0.0147987907597169, 741 0.0144882814685445, 742 0.0141714592928592, 743 0.0138484622795371, 744 0.0135194311690004, 745 0.0131845093222756, 746 0.0128438426753249, 747 0.0124975796646449, 748 0.0121458711652067, 749 0.0117888704247183, 750 0.0114267329968529, 751 0.0110596166734735, 752 0.0106876814158841, 753 0.0103110892851360, 754 0.0099300043714212, 755 0.0095445927225849, 756 0.0091550222717888, 757 0.0087614627643580, 758 0.0083640856838475, 759 0.0079630641773633, 760 0.0075585729801782, 761 0.0071507883396855, 762 0.0067398879387430, 763 0.0063260508184704, 764 0.0059094573005900, 765 0.0054902889094487, 766 0.0050687282939456, 767 0.0046449591497966, 768 0.0042191661429919, 769 0.0037915348363451, 770 0.0033622516236779, 771 0.0029315036836558, 772 0.0024994789888943, 773 0.0020663664924131, 774 0.0016323569986067, 775 0.0011976474864367, 776 0.0007624720924706, 777 0.0003276086705538 805 778 ])
Note: See TracChangeset
for help on using the changeset viewer.