Changeset 893bea2 in sasmodels
 Timestamp:
 Feb 28, 2018 10:23:19 AM (6 years ago)
 Branches:
 master, core_shell_microgels, magnetic_model, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 4f6f9431
 Parents:
 199bd07
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

explore/realspace.py
r3db96b0 r893bea2 171 171 points = uniform(1, 1, size=(num_points, 3)) 172 172 radius = points[:, 0]**2 + points[:, 1]**2 173 points = self._scale*points[radius <= 1]173 points = points[radius <= 1] 174 174 values = self.value.repeat(points.shape[0]) 175 return values, self._adjust(points) 175 return values, self._adjust(self._scale*points) 176 177 class EllipticalBicelle(Shape): 178 def __init__(self, ra, rb, length, 179 thick_rim, thick_face, 180 value_core, value_rim, value_face, 181 center=(0, 0, 0), orientation=(0, 0, 0)): 182 self.rotate(*orientation) 183 self.shift(*center) 184 self.value = value_core 185 self.ra, self.rb, self.length = ra, rb, length 186 self.thick_rim, self.thick_face = thick_rim, thick_face 187 self.value_rim, self.value_face = value_rim, value_face 188 189 # reset cylinder to outer dimensions for calculating scale, etc. 190 ra = self.ra + self.thick_rim 191 rb = self.rb + self.thick_rim 192 length = self.length + 2*self.thick_face 193 self._scale = np.array([ra, rb, length/2])[None, :] 194 self.r_max = sqrt(4*max(ra, rb)**2 + length**2) 195 self.dims = 2*ra, 2*rb, length 196 self.volume = pi*ra*rb*length 197 198 def sample(self, density): 199 # randomly sample from a box of side length 2*r, excluding anything 200 # not in the cylinder 201 ra = self.ra + self.thick_rim 202 rb = self.rb + self.thick_rim 203 length = self.length + 2*self.thick_face 204 num_points = poisson(density*4*ra*rb*length) 205 points = uniform(1, 1, size=(num_points, 3)) 206 radius = points[:, 0]**2 + points[:, 1]**2 207 points = points[radius <= 1] 208 # set all to core value first 209 values = np.ones_like(points[:, 0])*self.value 210 # then set value to face value if z > face/(length/2)) 211 values[abs(points[:, 2]) > self.length/(self.length + 2*self.thick_face)] = self.value_face 212 # finally set value to rim value if outside the core ellipse 213 radius = (points[:, 0]**2*(1+(self.thick_rim/self.ra)**2) 214 + points[:, 1]**2*(1+(self.thick_rim/self.rb)**2)) 215 values[radius>1] = self.value_rim 216 return values, self._adjust(self._scale*points) 176 217 177 218 class TriaxialEllipsoid(Shape): … … 479 520 #low, high = points.min(axis=0), points.max(axis=0) 480 521 #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]]) 522 ax.set_xlabel("x") 523 ax.set_ylabel("y") 524 ax.set_zlabel("z") 481 525 ax.autoscale(True) 482 526 … … 609 653 return Iq.reshape(qx.shape) 610 654 655 def sasmodels_Iq(kernel, q, pars): 656 from sasmodels.data import empty_data1D 657 from sasmodels.direct_model import DirectModel 658 data = empty_data1D(q) 659 calculator = DirectModel(data, kernel) 660 Iq = calculator(**pars) 661 return Iq 662 663 def sasmodels_Iqxy(kernel, qx, qy, pars, view): 664 from sasmodels.data import Data2D 665 from sasmodels.direct_model import DirectModel 666 Iq = 100 * np.ones_like(qx) 667 data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq)) 668 data.x_bins = qx[0,:] 669 data.y_bins = qy[:,0] 670 data.filename = "fake data" 671 672 calculator = DirectModel(data, kernel) 673 pars_plus_view = pars.copy() 674 pars_plus_view.update(theta=view[0], phi=view[1], psi=view[2]) 675 Iqxy = calculator(**pars_plus_view) 676 return Iqxy.reshape(qx.shape) 677 678 def wrap_sasmodel(name, **pars): 679 from sasmodels.core import load_model 680 kernel = load_model(name) 681 fn = lambda q: sasmodels_Iq(kernel, q, pars) 682 fn_xy = lambda qx, qy, view: sasmodels_Iqxy(kernel, qx, qy, pars, view) 683 return fn, fn_xy 684 685 611 686 #  Test cases  612 687 … … 634 709 fn_xy = lambda qx, qy, view: csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 635 710 slda, sldb, sldc, sld_core, view=view) 711 return shape, fn, fn_xy 712 713 def build_ellcyl(ra=25, rb=50, length=125, rho=2.): 714 shape = EllipticalCylinder(ra, rb, length, rho) 715 fn, fn_xy = wrap_sasmodel( 716 'elliptical_cylinder', 717 scale=1, 718 background=0, 719 radius_minor=ra, 720 axis_ratio=rb/ra, 721 length=length, 722 sld=rho, 723 sld_solvent=0, 724 ) 725 return shape, fn, fn_xy 726 727 def build_cscyl(ra=30, rb=90, length=30, thick_rim=8, thick_face=14, 728 sld_core=4, sld_rim=1, sld_face=7): 729 shape = EllipticalBicelle( 730 ra=ra, rb=rb, length=length, 731 thick_rim=thick_rim, thick_face=thick_face, 732 value_core=sld_core, value_rim=sld_rim, value_face=sld_face, 733 ) 734 fn, fn_xy = wrap_sasmodel( 735 'core_shell_bicelle_elliptical', 736 scale=1, 737 background=0, 738 radius=rb, 739 x_core=ra/rb, 740 length=length, 741 sld_core=sld_core, 742 sld_face=sld_face, 743 sld_rim=sld_rim, 744 sld_solvent=0, 745 ) 636 746 return shape, fn, fn_xy 637 747 … … 652 762 653 763 SHAPE_FUNCTIONS = OrderedDict([ 654 ("cylinder", build_cylinder), 764 ("cyl", build_cylinder), 765 ("ellcyl", build_ellcyl), 655 766 ("sphere", build_sphere), 656 767 ("box", build_box), 657 768 ("csbox", build_csbox), 769 ("cscyl", build_cscyl), 658 770 ]) 659 771 SHAPES = list(SHAPE_FUNCTIONS.keys()) … … 683 795 mesh=100, qmax=1.0, samples=5000): 684 796 rho_solvent = 0 685 qx = np.linspace(0.0, qmax, mesh) 686 qy = np.linspace(0.0, qmax, mesh) 797 #qx = np.linspace(0.0, qmax, mesh) 798 #qy = np.linspace(0.0, qmax, mesh) 799 qx = np.linspace(qmax, qmax, mesh) 800 qy = np.linspace(qmax, qmax, mesh) 687 801 Qx, Qy = np.meshgrid(qx, qy) 688 802 sampling_density = samples / shape.volume … … 693 807 Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 694 808 print("calc Iqxy time", time.time()  t0) 809 t0 = time.time() 695 810 theory = fn(Qx, Qy, view) if fn is not None else None 811 print("calc theory time", time.time()  t0) 696 812 Iqxy += 0.001 * Iqxy.max() 697 813 if theory is not None:
Note: See TracChangeset
for help on using the changeset viewer.