Changeset 5778279 in sasmodels
 Timestamp:
 Sep 25, 2018 7:39:08 PM (4 weeks ago)
 Branches:
 ticket_1156
 Parents:
 2a39ca4
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

explore/realspace.py
r2a39ca4 r5778279 9 9 import numpy as np 10 10 from numpy import pi, radians, sin, cos, sqrt 11 from numpy.random import poisson, uniform, randn, rand 11 from numpy.random import poisson, uniform, randn, rand, randint 12 12 from numpy.polynomial.legendre import leggauss 13 13 from scipy.integrate import simps … … 84 84 center = np.array([0., 0., 0.])[:, None] 85 85 r_max = None 86 lattice_size = np.array((1, 1, 1)) 87 lattice_spacing = np.array((1., 1., 1.)) 88 lattice_distortion = 0.0 89 lattice_rotation = 0.0 90 lattice_type = "" 86 91 87 92 def volume(self): … … 105 110 self.center = self.center + np.array([x, y, z])[:, None] 106 111 return self 112 113 def lattice(self, size=(1, 1, 1), spacing=(2, 2, 2), type="sc", 114 distortion=0.0, rotation=0.0): 115 self.lattice_size = np.asarray(size, 'i') 116 self.lattice_spacing = np.asarray(spacing, 'd') 117 self.lattice_type = type 118 self.lattice_distortion = distortion 119 self.lattice_rotation = rotation 107 120 108 121 def _adjust(self, points): … … 111 124 else: 112 125 points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 126 if self.lattice_type: 127 points = self._apply_lattice(points) 113 128 return points.T 114 129 115 def r_bins(self, q, over_sampling=1, r_step=0.): 116 r_max = min(2 * pi / q[0], self.r_max) 130 def r_bins(self, q, over_sampling=10, r_step=0.): 131 if self.lattice_type: 132 r_max = np.sqrt(np.sum(self.lattice_size*self.lattice_spacing*self.dims)**2)/2 133 else: 134 r_max = self.r_max 135 #r_max = min(2 * pi / q[0], r_max) 117 136 if r_step == 0.: 118 137 r_step = 2 * pi / q[1] / over_sampling 119 138 #r_step = 0.01 120 139 return np.arange(r_step, r_max, r_step) 140 141 def _apply_lattice(self, points): 142 """Spread points to different lattice positions""" 143 size = self.lattice_size 144 spacing = self.lattice_spacing 145 shuffle = self.lattice_distortion 146 rotate = self.lattice_rotation 147 lattice = self.lattice_type 148 149 if rotate != 0: 150 # To vectorize the rotations we will need to unwrap the matrix multiply 151 raise NotImplementedError("don't handle rotations yet") 152 153 # Determine the number of lattice points in the lattice 154 shapes_per_cell = 2 if lattice == "bcc" else 4 if lattice == "fcc" else 1 155 number_of_lattice_points = np.prod(size) * shapes_per_cell 156 157 # For each point in the original shape, figure out which lattice point 158 # to translate it to. This is both cell index (i*ny*nz + j*nz + k) as 159 # well as the point in the cell (corner, body center or face center). 160 nsamples = points.shape[1] 161 lattice_point = randint(number_of_lattice_points, size=nsamples) 162 163 # Translate the cell index into the i,j,k coordinates of the senter 164 cell_index = lattice_point // shapes_per_cell 165 center = np.vstack((cell_index//(size[1]*size[2]), 166 (cell_index%(size[1]*size[2]))//size[2], 167 cell_index%size[2])) 168 center = np.asarray(center, dtype='d') 169 if lattice == "bcc": 170 center[:, lattice_point % shapes_per_cell == 1] += [[0.5], [0.5], [0.5]] 171 elif lattice == "fcc": 172 center[:, lattice_point % shapes_per_cell == 1] += [[0.0], [0.5], [0.5]] 173 center[:, lattice_point % shapes_per_cell == 2] += [[0.5], [0.0], [0.5]] 174 center[:, lattice_point % shapes_per_cell == 3] += [[0.5], [0.5], [0.0]] 175 176 # Each lattice point has its own displacement from the ideal position. 177 # Not checking that shapes do not overlap if displacement is too large. 178 offset = shuffle*(randn(3, number_of_lattice_points) if shuffle < 0.3 179 else rand(3, number_of_lattice_points)) 180 center += offset[:, cell_index] 181 182 # Each lattice point has its own rotation. Rotate the point prior to 183 # applying any displacement. 184 # rotation = rotate*(randn(size=(shapes, 3)) if shuffle < 30 else rand(size=(nsamples, 3))) 185 # for k in shapes: points[k] = rotation[k]*points[k] 186 points += center*(np.array([spacing])*np.array(self.dims)).T 187 return points 121 188 122 189 class Composite(Shape): … … 867 934 import pylab 868 935 if show_points: 869 936 plot_points(rho, points); pylab.figure() 870 937 plot_calc(r, Pr, q, Iq, theory=theory, title=title) 871 938 pylab.gcf().canvas.set_window_title(title) … … 937 1004 nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 938 1005 dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 939 shuffle, rotate= opts.shuffle, opts.rotate1006 distortion, rotation = opts.shuffle, opts.rotate 940 1007 shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 941 1008 view = tuple(float(v) for v in opts.view.split(',')) 942 if nx > 1 or ny > 1 or nz > 1: 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) 1009 # If comparing a sphere in a cubic lattice, compare against the 1010 # corresponding paracrystalline model. 1011 if opts.shape == "sphere" and dx == dy == dz and nx*ny*nz > 1: 1012 radius = pars.get('radius', DEFAULT_SPHERE_RADIUS) 1013 model_name = opts.type + "_paracrystal" 1014 model_pars = { 1015 "scale": 1., 1016 "background": 0., 1017 "lattice_spacing": 2*radius*dx, 1018 "lattice_distortion": distortion, 1019 "radius": radius, 1020 "sld": pars.get('rho', DEFAULT_SPHERE_CONTRAST), 1021 "sld_solvent": 0., 1022 "theta": view[0], 1023 "phi": view[1], 1024 "psi": view[2], 1025 } 1026 fn, fn_xy = wrap_sasmodel(model_name, **model_pars) 1027 if nx*ny*nz > 1: 1028 if rotation != 0: 1029 print("building %s lattice"%opts.type) 1030 build_lattice = LATTICE_FUNCTIONS[opts.type] 1031 shape = build_lattice(shape, nx, ny, nz, dx, dy, dz, 1032 distortion, rotation) 1033 else: 1034 shape.lattice(size=(nx, ny, nz), spacing=(dx, dy, dz), 1035 type=opts.type, 1036 rotation=rotation, distortion=distortion) 1037 964 1038 title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 965 1039 if opts.dim == 1:
Note: See TracChangeset
for help on using the changeset viewer.