Changes in / [f6fd413:13bd583] in sasmodels
- Files:
-
- 3 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/guide/pd/polydispersity.rst
reda8b30 r92d330fd 42 42 calculations are generally more robust with more data points or more angles. 43 43 44 The following fivedistribution functions are provided:44 The following distribution functions are provided: 45 45 46 46 * *Rectangular Distribution* 47 * *Uniform Distribution* 47 48 * *Gaussian Distribution* 48 49 * *Lognormal Distribution* 49 50 * *Schulz Distribution* 50 51 * *Array Distribution* 52 * *Boltzmann Distribution* 51 53 52 54 These are all implemented as *number-average* distributions. … … 85 87 Rectangular distribution. 86 88 89 90 91 Uniform Distribution 92 ^^^^^^^^^^^^^^^^^^^^^^^^ 93 94 The Uniform Distribution is defined as 95 96 .. math:: 97 98 f(x) = \frac{1}{\text{Norm}} 99 \begin{cases} 100 1 & \text{for } |x - \bar x| \leq \sigma \\ 101 0 & \text{for } |x - \bar x| > \sigma 102 \end{cases} 103 104 where $\bar x$ is the mean of the distribution, $\sigma$ is the half-width, and 105 *Norm* is a normalization factor which is determined during the numerical 106 calculation. 107 108 Note that the polydispersity is given by 109 110 .. math:: \text{PD} = \sigma / \bar x 111 112 .. figure:: pd_uniform.jpg 113 114 Uniform distribution. 115 116 The value $N_\sigma$ is ignored for this distribution. 117 87 118 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 88 119 … … 183 214 ^^^^^^^^^^^^^^^^^^ 184 215 185 This user-definable distribution should be given as a s asimple ASCII text216 This user-definable distribution should be given as a simple ASCII text 186 217 file where the array is defined by two columns of numbers: $x$ and $f(x)$. 187 218 The $f(x)$ will be normalized to 1 during the computation. … … 202 233 given for the model will have no affect, and will be ignored when computing 203 234 the average. This means that any parameter with an array distribution will 204 not be fittable. 235 not be fitable. 236 237 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 238 239 Boltzmann Distribution 240 ^^^^^^^^^^^^^^^^^^^^^^ 241 242 The Boltzmann Distribution is defined as 243 244 .. math:: 245 246 f(x) = \frac{1}{\text{Norm}} 247 \exp\left(-\frac{ | x - \bar x | }{\sigma}\right) 248 249 where $\bar x$ is the mean of the distribution and *Norm* is a normalization 250 factor which is determined during the numerical calculation. 251 The width is defined as 252 253 .. math:: \sigma=\frac{k T}{E} 254 255 which is the inverse Boltzmann factor, 256 where $k$ is the Boltzmann constant, $T$ the temperature in Kelvin and $E$ a 257 characteristic energy per particle. 258 259 .. figure:: pd_boltzmann.jpg 260 261 Boltzmann distribution. 205 262 206 263 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ -
explore/realspace.py
r8fb2a94 r032646d 3 3 import time 4 4 from copy import copy 5 import os 6 import argparse 7 from collections import OrderedDict 5 8 6 9 import numpy as np 7 10 from numpy import pi, radians, sin, cos, sqrt 8 from numpy.random import poisson, uniform 11 from numpy.random import poisson, uniform, randn, rand 9 12 from numpy.polynomial.legendre import leggauss 10 13 from scipy.integrate import simps 11 14 from scipy.special import j1 as J1 15 16 try: 17 import numba 18 USE_NUMBA = True 19 except ImportError: 20 USE_NUMBA = False 12 21 13 22 # Definition of rotation matrices comes from wikipedia: … … 82 91 raise NotImplementedError() 83 92 93 def dims(self): 94 # type: () -> float, float, float 95 raise NotImplementedError() 96 84 97 def rotate(self, theta, phi, psi): 85 98 self.rotation = rotation(theta, phi, psi) * self.rotation … … 116 129 for s2 in shapes] 117 130 self.r_max = max(distances + [s.r_max for s in shapes]) 118 119 def volume(self): 120 return sum(shape.volume() for shape in self.shapes) 131 self.volume = sum(shape.volume for shape in self.shapes) 121 132 122 133 def sample(self, density): … … 133 144 self._scale = np.array([a/2, b/2, c/2])[None, :] 134 145 self.r_max = sqrt(a**2 + b**2 + c**2) 135 136 def volume(self): 137 return self.a*self.b*self.c 146 self.dims = a, b, c 147 self.volume = a*b*c 138 148 139 149 def sample(self, density): 140 num_points = poisson(density*self. a*self.b*self.c)150 num_points = poisson(density*self.volume) 141 151 points = self._scale*uniform(-1, 1, size=(num_points, 3)) 142 152 values = self.value.repeat(points.shape[0]) … … 152 162 self._scale = np.array([ra, rb, length/2])[None, :] 153 163 self.r_max = sqrt(4*max(ra, rb)**2 + length**2) 154 155 def volume(self): 156 return pi*self.ra*self.rb*self.length 164 self.dims = 2*ra, 2*rb, length 165 self.volume = pi*ra*rb*length 157 166 158 167 def sample(self, density): 159 # density of the bounding box 168 # randomly sample from a box of side length 2*r, excluding anything 169 # not in the cylinder 160 170 num_points = poisson(density*4*self.ra*self.rb*self.length) 161 171 points = uniform(-1, 1, size=(num_points, 3)) … … 174 184 self._scale = np.array([ra, rb, rc])[None, :] 175 185 self.r_max = 2*max(ra, rb, rc) 176 177 def volume(self): 178 return 4*pi/3 * self.ra * self.rb * self.rc 186 self.dims = 2*ra, 2*rb, 2*rc 187 self.volume = 4*pi/3 * ra * rb * rc 179 188 180 189 def sample(self, density): … … 194 203 self.rotate(*orientation) 195 204 self.shift(*center) 205 helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 206 total_radius = self.helix_radius + self.tube_radius 196 207 self.helix_radius, self.helix_pitch = helix_radius, helix_pitch 197 208 self.tube_radius, self.tube_length = tube_radius, tube_length 198 helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 199 self.r_max = sqrt((2*helix_radius + 2*tube_radius)*2 200 + (helix_length + 2*tube_radius)**2) 201 202 def volume(self): 209 self.r_max = sqrt(4*total_radius + (helix_length + 2*tube_radius)**2) 210 self.dims = 2*total_radius, 2*total_radius, helix_length 203 211 # small tube radius approximation; for larger tubes need to account 204 212 # for the fact that the inner length is much shorter than the outer 205 213 # length 206 returnpi*self.tube_radius**2*self.tube_length214 self.volume = pi*self.tube_radius**2*self.tube_length 207 215 208 216 def points(self, density): … … 244 252 return values, self._adjust(points) 245 253 246 NUMBA = False 247 if NUMBA: 254 def csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 255 core = Box(a, b, c, sld_core) 256 side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0)) 257 side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0)) 258 side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2)) 259 side_a2 = copy(side_a).shift(-a-da, 0, 0) 260 side_b2 = copy(side_b).shift(0, -b-db, 0) 261 side_c2 = copy(side_c).shift(0, 0, -c-dc) 262 shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 263 shape.dims = 2*da+a, 2*db+b, 2*dc+c 264 return shape 265 266 def _Iqxy(values, x, y, z, qa, qb, qc): 267 """I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r)""" 268 Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 269 for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 270 return Iq 271 272 if USE_NUMBA: 273 # Override simple numpy solution with numba if available 248 274 from numba import njit 249 275 @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])") … … 252 278 for j in range(len(Iq)): 253 279 total = 0. + 0j 254 for k in range(len( Iq)):280 for k in range(len(values)): 255 281 total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k])) 256 282 Iq[j] = abs(total)**2 … … 263 289 values = rho*volume 264 290 x, y, z = points.T 291 values, x, y, z, qa, qb, qc = [np.asarray(v, 'd') 292 for v in (values, x, y, z, qa, qb, qc)] 265 293 266 294 # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r) 267 if NUMBA: 268 Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 269 else: 270 Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 271 for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 295 Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 272 296 return np.asarray(Iq).reshape(qx.shape) / np.sum(volume) 273 297 … … 338 362 """ 339 363 340 if NUMBA: 364 if USE_NUMBA: 365 # Override simple numpy solution with numba if available 341 366 @njit("f8[:](f8[:], f8[:], f8[:,:])") 342 def _calc_Pr_uniform _jit(r, rho, points):367 def _calc_Pr_uniform(r, rho, points): 343 368 dr = r[0] 344 369 n_max = len(r) … … 362 387 # P(r) with uniform steps in r is 3x faster; check if we are uniform 363 388 # before continuing 389 r, rho, points = [np.asarray(v, 'd') for v in (r, rho, points)] 364 390 if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 365 391 Pr = _calc_Pr_nonuniform(r, rho, points) 366 392 else: 367 if NUMBA: 368 Pr = _calc_Pr_uniform_jit(r, rho, points) 369 else: 370 Pr = _calc_Pr_uniform(r, rho, points) 393 Pr = _calc_Pr_uniform(r, rho, points) 371 394 return Pr / Pr.max() 372 395 … … 399 422 return edges 400 423 424 # -------------- plotters ---------------- 401 425 def plot_calc(r, Pr, q, Iq, theory=None): 402 426 import matplotlib.pyplot as plt … … 410 434 plt.ylabel('Iq') 411 435 if theory is not None: 412 plt.loglog(theory[0], theory[1] , '-', label='analytic')436 plt.loglog(theory[0], theory[1]/theory[1][0], '-', label='analytic') 413 437 plt.legend() 414 438 … … 444 468 ax.autoscale(True) 445 469 446 def check_shape(shape, fn=None): 447 rho_solvent = 0 448 q = np.logspace(-3, 0, 200) 449 r = shape.r_bins(q, r_step=0.01) 450 sampling_density = 6*5000 / shape.volume() 451 rho, points = shape.sample(sampling_density) 452 t0 = time.time() 453 Pr = calc_Pr(r, rho-rho_solvent, points) 454 print("calc Pr time", time.time() - t0) 455 Iq = calc_Iq(q, r, Pr) 456 theory = (q, fn(q)) if fn is not None else None 457 458 import pylab 459 #plot_points(rho, points); pylab.figure() 460 plot_calc(r, Pr, q, Iq, theory=theory) 461 pylab.show() 462 463 def check_shape_2d(shape, fn=None, view=(0, 0, 0)): 464 rho_solvent = 0 465 nq, qmax = 100, 1.0 466 qx = np.linspace(0.0, qmax, nq) 467 qy = np.linspace(0.0, qmax, nq) 468 Qx, Qy = np.meshgrid(qx, qy) 469 sampling_density = 50000 / shape.volume() 470 #t0 = time.time() 471 rho, points = shape.sample(sampling_density) 472 #print("sample time", time.time() - t0) 473 t0 = time.time() 474 Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 475 print("calc time", time.time() - t0) 476 theory = fn(Qx, Qy) if fn is not None else None 477 Iqxy += 0.001 * Iqxy.max() 478 if theory is not None: 479 theory += 0.001 * theory.max() 480 481 import pylab 482 #plot_points(rho, points); pylab.figure() 483 plot_calc_2d(qx, qy, Iqxy, theory=theory) 484 pylab.show() 485 470 # ----------- Analytic models -------------- 486 471 def sas_sinx_x(x): 487 472 with np.errstate(all='ignore'): … … 510 495 for k, qk in enumerate(q): 511 496 qab, qc = qk*sin_alpha, qk*cos_alpha 512 Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2)497 Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 513 498 Iq[k] = np.sum(w*Fq**2) 514 Iq = Iq /Iq[0]499 Iq = Iq 515 500 return Iq 516 501 517 502 def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)): 518 503 qa, qb, qc = invert_view(qx, qy, view) 519 qab = np.sqrt(qa**2 + qb**2)520 Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2)504 qab = sqrt(qa**2 + qb**2) 505 Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 521 506 Iq = Fq**2 522 507 return Iq.reshape(qx.shape) … … 524 509 def sphere_Iq(q, radius): 525 510 Iq = sas_3j1x_x(q*radius)**2 526 return Iq/Iq[0] 511 return Iq 512 513 def box_Iq(q, a, b, c): 514 z, w = leggauss(76) 515 outer_sum = np.zeros_like(q) 516 for cos_alpha, outer_w in zip((z+1)/2, w): 517 sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 518 qc = q*cos_alpha 519 siC = c*sas_sinx_x(c*qc/2) 520 inner_sum = np.zeros_like(q) 521 for beta, inner_w in zip((z + 1)*pi/4, w): 522 qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta) 523 siA = a*sas_sinx_x(a*qa/2) 524 siB = b*sas_sinx_x(b*qb/2) 525 Fq = siA*siB*siC 526 inner_sum += inner_w * Fq**2 527 outer_sum += outer_w * inner_sum 528 Iq = outer_sum / 4 # = outer*um*zm*8.0/(4.0*M_PI) 529 return Iq 530 531 def box_Iqxy(qx, qy, a, b, c, view=(0, 0, 0)): 532 qa, qb, qc = invert_view(qx, qy, view) 533 sia = sas_sinx_x(qa*a/2) 534 sib = sas_sinx_x(qb*b/2) 535 sic = sas_sinx_x(qc*c/2) 536 Fq = sia*sib*sic 537 Iq = Fq**2 538 return Iq.reshape(qx.shape) 527 539 528 540 def csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core): … … 539 551 sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 540 552 qc = q*cos_alpha 541 siC = c* j0(c*qc/2)542 siCt = tC* j0(tC*qc/2)553 siC = c*sas_sinx_x(c*qc/2) 554 siCt = tC*sas_sinx_x(tC*qc/2) 543 555 inner_sum = np.zeros_like(q) 544 556 for beta, inner_w in zip((z + 1)*pi/4, w): 545 557 qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta) 546 siA = a* j0(a*qa/2)547 siB = b* j0(b*qb/2)548 siAt = tA* j0(tA*qa/2)549 siBt = tB* j0(tB*qb/2)558 siA = a*sas_sinx_x(a*qa/2) 559 siB = b*sas_sinx_x(b*qb/2) 560 siAt = tA*sas_sinx_x(tA*qa/2) 561 siBt = tB*sas_sinx_x(tB*qb/2) 550 562 if overlapping: 551 563 Fq = (dr0*siA*siB*siC … … 584 596 return Iq.reshape(qx.shape) 585 597 586 def check_cylinder(radius=25, length=125, rho=2.): 598 # --------- Test cases ----------- 599 600 def build_cylinder(radius=25, length=125, rho=2.): 587 601 shape = EllipticalCylinder(radius, radius, length, rho) 588 fn = lambda q: cylinder_Iq(q, radius, length) 589 check_shape(shape, fn) 590 591 def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)): 592 shape = EllipticalCylinder(radius, radius, length, rho) 593 fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 594 check_shape_2d(shape, fn, view=view) 595 596 def check_cylinder_2d_lattice(radius=25, length=125, rho=2., 597 view=(0, 0, 0)): 598 nx, dx = 1, 2*radius 599 ny, dy = 30, 2*radius 600 nz, dz = 30, length 601 dx, dy, dz = 2*dx, 2*dy, 2*dz 602 def center(*args): 603 sigma = 0.333 604 space = 2 605 return [(space*n+np.random.randn()*sigma)*x for n, x in args] 606 shapes = [EllipticalCylinder(radius, radius, length, rho, 607 #center=(ix*dx, iy*dy, iz*dz) 608 orientation=np.random.randn(3)*0, 609 center=center((ix, dx), (iy, dy), (iz, dz)) 610 ) 602 fn = lambda q: cylinder_Iq(q, radius, length)*rho**2 603 fn_xy = lambda qx, qy, view: cylinder_Iqxy(qx, qy, radius, length, view=view)*rho**2 604 return shape, fn, fn_xy 605 606 def build_sphere(radius=125, rho=2): 607 shape = TriaxialEllipsoid(radius, radius, radius, rho) 608 fn = lambda q: sphere_Iq(q, radius)*rho**2 609 fn_xy = lambda qx, qy, view: sphere_Iq(np.sqrt(qx**2+qy**2), radius)*rho**2 610 return shape, fn, fn_xy 611 612 def build_box(a=10, b=20, c=30, rho=2.): 613 shape = Box(a, b, c, rho) 614 fn = lambda q: box_Iq(q, a, b, c)*rho**2 615 fn_xy = lambda qx, qy, view: box_Iqxy(qx, qy, a, b, c, view=view)*rho**2 616 return shape, fn, fn_xy 617 618 def build_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 619 shape = csbox(a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 620 fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 621 fn_xy = lambda qx, qy, view: csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 622 slda, sldb, sldc, sld_core, view=view) 623 return shape, fn, fn_xy 624 625 def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 626 shuffle=0, rotate=0): 627 a, b, c = shape.dims 628 shapes = [copy(shape) 629 .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 630 (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 631 (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 632 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 611 633 for ix in range(nx) 612 634 for iy in range(ny) 613 635 for iz in range(nz)] 614 shape = Composite(shapes) 615 fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 616 check_shape_2d(shape, fn, view=view) 617 618 def check_sphere(radius=125, rho=2): 619 shape = TriaxialEllipsoid(radius, radius, radius, rho) 620 fn = lambda q: sphere_Iq(q, radius) 621 check_shape(shape, fn) 622 623 def check_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 624 core = Box(a, b, c, sld_core) 625 side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0)) 626 side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0)) 627 side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2)) 628 side_a2 = copy(side_a).shift(-a-da, 0, 0) 629 side_b2 = copy(side_b).shift(0, -b-db, 0) 630 side_c2 = copy(side_c).shift(0, 0, -c-dc) 631 shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 632 def fn(q): 633 return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 634 #check_shape(shape, fn) 635 636 view = (20, 30, 40) 637 def fn_xy(qx, qy): 638 return csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 639 slda, sldb, sldc, sld_core, view=view) 640 check_shape_2d(shape, fn_xy, view=view) 636 lattice = Composite(shapes) 637 return lattice 638 639 640 SHAPE_FUNCTIONS = OrderedDict([ 641 ("cylinder", build_cylinder), 642 ("sphere", build_sphere), 643 ("box", build_box), 644 ("csbox", build_csbox), 645 ]) 646 SHAPES = list(SHAPE_FUNCTIONS.keys()) 647 648 def check_shape(title, shape, fn=None, show_points=False, 649 mesh=100, qmax=1.0, r_step=0.01, samples=5000): 650 rho_solvent = 0 651 qmin = qmax/100. 652 q = np.logspace(np.log10(qmin), np.log10(qmax), mesh) 653 r = shape.r_bins(q, r_step=r_step) 654 sampling_density = samples / shape.volume 655 rho, points = shape.sample(sampling_density) 656 t0 = time.time() 657 Pr = calc_Pr(r, rho-rho_solvent, points) 658 print("calc Pr time", time.time() - t0) 659 Iq = calc_Iq(q, r, Pr) 660 theory = (q, fn(q)) if fn is not None else None 661 662 import pylab 663 if show_points: 664 plot_points(rho, points); pylab.figure() 665 plot_calc(r, Pr, q, Iq, theory=theory) 666 pylab.gcf().canvas.set_window_title(title) 667 pylab.show() 668 669 def check_shape_2d(title, shape, fn=None, view=(0, 0, 0), show_points=False, 670 mesh=100, qmax=1.0, samples=5000): 671 rho_solvent = 0 672 qx = np.linspace(0.0, qmax, mesh) 673 qy = np.linspace(0.0, qmax, mesh) 674 Qx, Qy = np.meshgrid(qx, qy) 675 sampling_density = samples / shape.volume 676 t0 = time.time() 677 rho, points = shape.sample(sampling_density) 678 print("point generation time", time.time() - t0) 679 t0 = time.time() 680 Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 681 print("calc Iqxy time", time.time() - t0) 682 theory = fn(Qx, Qy, view) if fn is not None else None 683 Iqxy += 0.001 * Iqxy.max() 684 if theory is not None: 685 theory += 0.001 * theory.max() 686 687 import pylab 688 if show_points: 689 plot_points(rho, points); pylab.figure() 690 plot_calc_2d(qx, qy, Iqxy, theory=theory) 691 pylab.gcf().canvas.set_window_title(title) 692 pylab.show() 693 694 def main(): 695 parser = argparse.ArgumentParser( 696 description="Compute scattering from realspace sampling", 697 formatter_class=argparse.ArgumentDefaultsHelpFormatter, 698 ) 699 parser.add_argument('-d', '--dim', type=int, default=1, 700 help='dimension 1 or 2') 701 parser.add_argument('-m', '--mesh', type=int, default=100, 702 help='number of mesh points') 703 parser.add_argument('-s', '--samples', type=int, default=5000, 704 help="number of sample points") 705 parser.add_argument('-q', '--qmax', type=float, default=0.5, 706 help='max q') 707 parser.add_argument('-v', '--view', type=str, default='0,0,0', 708 help='theta,phi,psi angles') 709 parser.add_argument('-n', '--lattice', type=str, default='1,1,1', 710 help='lattice size') 711 parser.add_argument('-z', '--spacing', type=str, default='2,2,2', 712 help='lattice spacing') 713 parser.add_argument('-r', '--rotate', type=float, default=0., 714 help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") 715 parser.add_argument('-w', '--shuffle', type=float, default=0., 716 help="position relative to lattice, gaussian < 0.3, uniform otherwise") 717 parser.add_argument('-p', '--plot', action='store_true', 718 help='plot points') 719 parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], 720 help='oriented shape') 721 parser.add_argument('pars', type=str, nargs='*', help='shape parameters') 722 opts = parser.parse_args() 723 pars = {key: float(value) for p in opts.pars for key, value in [p.split('=')]} 724 nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 725 dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 726 shuffle, rotate = opts.shuffle, opts.rotate 727 shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 728 if nx > 1 or ny > 1 or nz > 1: 729 shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 730 title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 731 if opts.dim == 1: 732 check_shape(title, shape, fn, show_points=opts.plot, 733 mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 734 else: 735 view = tuple(float(v) for v in opts.view.split(',')) 736 check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 737 mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 738 641 739 642 740 if __name__ == "__main__": 643 check_cylinder(radius=10, length=40) 644 #check_cylinder_2d(radius=10, length=40, view=(90,30,0)) 645 #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0)) 646 #check_sphere() 647 #check_csbox() 648 #check_csbox(da=100, db=200, dc=300) 741 main() -
sasmodels/kernel_iq.c
rec8d4ac raadec17 93 93 // Compute the magnetic sld 94 94 static double mag_sld( 95 const int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag95 const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 96 96 const double qx, const double qy, 97 97 const double px, const double py, … … 103 103 const double perp = qy*mx - qx*my; 104 104 switch (xs) { 105 default: // keep compiler happy; condition ensures xs in [0,1,2,3] 105 106 case 0: // uu => sld - D M_perpx 106 107 return sld - px*perp; … … 659 660 660 661 // loop over uu, ud real, du real, dd, ud imag, du imag 661 for ( int xs=0; xs<6; xs++) {662 for (unsigned int xs=0; xs<6; xs++) { 662 663 const double xs_weight = spins[xs]; 663 664 if (xs_weight > 1.e-8) { -
sasmodels/weights.py
r2d81cfe r3d58247 97 97 return x, px 98 98 99 class UniformDispersion(Dispersion): 100 r""" 101 Uniform dispersion, with width $\sigma$. 102 103 .. math:: 104 105 w = 1 106 """ 107 type = "uniform" 108 default = dict(npts=35, width=0, nsigmas=None) 109 def _weights(self, center, sigma, lb, ub): 110 x = np.linspace(center-sigma, center+sigma, self.npts) 111 x = x[(x >= lb) & (x <= ub)] 112 return x, np.ones_like(x) 99 113 100 114 class RectangleDispersion(Dispersion): … … 107 121 """ 108 122 type = "rectangle" 109 default = dict(npts=35, width=0, nsigmas=1.70325) 110 def _weights(self, center, sigma, lb, ub): 111 x = self._linspace(center, sigma, lb, ub) 112 x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] 113 return x, np.ones_like(x) 114 123 default = dict(npts=35, width=0, nsigmas=1.73205) 124 def _weights(self, center, sigma, lb, ub): 125 x = self._linspace(center, sigma, lb, ub) 126 x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] 127 return x, np.ones_like(x) 115 128 116 129 class LogNormalDispersion(Dispersion): … … 190 203 return x, px 191 204 205 class BoltzmannDispersion(Dispersion): 206 r""" 207 Boltzmann dispersion, with $\sigma=k T/E$. 208 209 .. math:: 210 211 w = \exp\left( -|x - c|/\sigma\right) 212 """ 213 type = "boltzmann" 214 default = dict(npts=35, width=0, nsigmas=3) 215 def _weights(self, center, sigma, lb, ub): 216 x = self._linspace(center, sigma, lb, ub) 217 px = np.exp(-np.abs(x-center) / np.abs(sigma)) 218 return x, px 192 219 193 220 # dispersion name -> disperser lookup table. … … 196 223 MODELS = OrderedDict((d.type, d) for d in ( 197 224 RectangleDispersion, 225 UniformDispersion, 198 226 ArrayDispersion, 199 227 LogNormalDispersion, 200 228 GaussianDispersion, 201 229 SchulzDispersion, 230 BoltzmannDispersion 202 231 )) 203 232
Note: See TracChangeset
for help on using the changeset viewer.