Changeset 2a39ca4 in sasmodels
 Timestamp:
 Sep 25, 2018 12:20:12 PM (5 months ago)
 Branches:
 ticket_1156
 Children:
 5778279
 Parents:
 642046e
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

explore/realspace.py
r362a66f r2a39ca4 78 78 79 79 80 I3 = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]]) 81 80 82 class Shape: 81 rotation = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]])83 rotation = I3 82 84 center = np.array([0., 0., 0.])[:, None] 83 85 r_max = None … … 96 98 97 99 def rotate(self, theta, phi, psi): 98 self.rotation = rotation(theta, phi, psi) * self.rotation 100 if theta != 0. or phi != 0. or psi != 0.: 101 self.rotation = rotation(theta, phi, psi) * self.rotation 99 102 return self 100 103 … … 104 107 105 108 def _adjust(self, points): 106 points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 109 if self.rotation is I3: 110 points = points.T + self.center 111 else: 112 points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 107 113 return points.T 108 114 … … 669 675 Iq = 100 * np.ones_like(qx) 670 676 data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq)) 671 data.x_bins = qx[0, :]672 data.y_bins = qy[:, 0]677 data.x_bins = qx[0, :] 678 data.y_bins = qy[:, 0] 673 679 data.filename = "fake data" 674 680 … … 695 701 return shape, fn, fn_xy 696 702 697 def build_sphere(radius=125, rho=2): 703 DEFAULT_SPHERE_RADIUS = 125 704 DEFAULT_SPHERE_CONTRAST = 2 705 def build_sphere(radius=DEFAULT_SPHERE_RADIUS, rho=DEFAULT_SPHERE_CONTRAST): 698 706 shape = TriaxialEllipsoid(radius, radius, radius, rho) 699 707 fn = lambda q: sphere_Iq(q, radius)*rho**2 … … 751 759 return shape, fn, fn_xy 752 760 753 def build_ cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,754 shuffle=0, rotate=0):761 def build_sc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 762 shuffle=0, rotate=0): 755 763 a, b, c = shape.dims 756 shapes= [copy(shape)764 corners= [copy(shape) 757 765 .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 758 766 (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, … … 762 770 for iy in range(ny) 763 771 for iz in range(nz)] 764 lattice = Composite( shapes)772 lattice = Composite(corners) 765 773 return lattice 766 774 775 def build_bcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 776 shuffle=0, rotate=0): 777 a, b, c = shape.dims 778 corners = [copy(shape) 779 .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 780 (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 781 (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 782 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 783 for ix in range(nx) 784 for iy in range(ny) 785 for iz in range(nz)] 786 centers = [copy(shape) 787 .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 788 (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 789 (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 790 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 791 for ix in range(nx) 792 for iy in range(ny) 793 for iz in range(nz)] 794 lattice = Composite(corners + centers) 795 return lattice 796 797 def build_fcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 798 shuffle=0, rotate=0): 799 a, b, c = shape.dims 800 corners = [copy(shape) 801 .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 802 (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 803 (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 804 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 805 for ix in range(nx) 806 for iy in range(ny) 807 for iz in range(nz)] 808 faces_a = [copy(shape) 809 .shift((ix+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 810 (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 811 (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 812 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 813 for ix in range(nx) 814 for iy in range(ny) 815 for iz in range(nz)] 816 faces_b = [copy(shape) 817 .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 818 (iy+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 819 (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 820 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 821 for ix in range(nx) 822 for iy in range(ny) 823 for iz in range(nz)] 824 faces_c = [copy(shape) 825 .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 826 (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 827 (iz+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 828 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 829 for ix in range(nx) 830 for iy in range(ny) 831 for iz in range(nz)] 832 lattice = Composite(corners + faces_a + faces_b + faces_c) 833 return lattice 767 834 768 835 SHAPE_FUNCTIONS = OrderedDict([ … … 775 842 ]) 776 843 SHAPES = list(SHAPE_FUNCTIONS.keys()) 844 LATTICE_FUNCTIONS = OrderedDict([ 845 ("sc", build_sc_lattice), 846 ("bcc", build_bcc_lattice), 847 ("fcc", build_fcc_lattice), 848 ]) 849 LATTICE_TYPES = list(LATTICE_FUNCTIONS.keys()) 777 850 778 851 def check_shape(title, shape, fn=None, show_points=False, … … 783 856 r = shape.r_bins(q, r_step=r_step) 784 857 sampling_density = samples / shape.volume 858 print("sampling points") 785 859 rho, points = shape.sample(sampling_density) 860 print("calculating Pr") 786 861 t0 = time.time() 787 862 Pr = calc_Pr(r, rhorho_solvent, points) … … 806 881 Qx, Qy = np.meshgrid(qx, qy) 807 882 sampling_density = samples / shape.volume 883 print("sampling points") 808 884 t0 = time.time() 809 885 rho, points = shape.sample(sampling_density) … … 844 920 help='lattice size') 845 921 parser.add_argument('z', 'spacing', type=str, default='2,2,2', 846 help='lattice spacing') 922 help='lattice spacing (relative to shape)') 923 parser.add_argument('t', 'type', choices=LATTICE_TYPES, 924 default=LATTICE_TYPES[0], 925 help='lattice type') 847 926 parser.add_argument('r', 'rotate', type=float, default=0., 848 927 help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") … … 860 939 shuffle, rotate = opts.shuffle, opts.rotate 861 940 shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 941 view = tuple(float(v) for v in opts.view.split(',')) 862 942 if nx > 1 or ny > 1 or nz > 1: 863 shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 943 print("building %s lattice"%opts.type) 944 lattice = LATTICE_FUNCTIONS[opts.type] 945 shape = lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 946 # If comparing a sphere in a cubic lattice, compare against the 947 # corresponding paracrystalline model. 948 if opts.shape == "sphere" and dx == dy == dz: 949 radius = pars.get('radius', DEFAULT_SPHERE_RADIUS) 950 model_name = opts.type + "_paracrystal" 951 model_pars = { 952 "scale": 1., 953 "background": 0., 954 "lattice_spacing": 2*radius*dx, 955 "lattice_distortion": shuffle, 956 "radius": radius, 957 "sld": pars.get('rho', DEFAULT_SPHERE_CONTRAST), 958 "sld_solvent": 0., 959 "theta": view[0], 960 "phi": view[1], 961 "psi": view[2], 962 } 963 fn, fn_xy = wrap_sasmodel(model_name, **model_pars) 864 964 title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 865 965 if opts.dim == 1: … … 867 967 mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 868 968 else: 869 view = tuple(float(v) for v in opts.view.split(','))870 969 check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 871 970 mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples)
Note: See TracChangeset
for help on using the changeset viewer.