# Changeset 802c412 in sasmodels

Ignore:
Timestamp:
Mar 26, 2018 4:55:17 PM (6 years ago)
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
b3703f5
Parents:
d86f0fc
Message:

lint reduction

File:
1 edited

### Legend:

Unmodified
 rd86f0fc import matplotlib.pyplot as plt from matplotlib.widgets import Slider, CheckButtons from matplotlib import cm from matplotlib.widgets import Slider import numpy as np from numpy import pi, cos, sin, sqrt, exp, degrees, radians def draw_beam(ax, view=(0, 0)): def draw_beam(axes, view=(0, 0)): """ Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) """ #ax.plot([0,0],[0,0],[1,-1]) #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) #axes.plot([0,0],[0,0],[1,-1]) #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) steps = 25 x, y, z = [v.reshape(shape) for v in points] ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1): """Draw an ellipsoid.""" a, b, c = size x, y, z = transform_xyz(view, jitter, x, y, z) ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) draw_labels(ax, view, jitter, [ axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) draw_labels(axes, view, jitter, [ ('c+', [+0, +0, +c], [+1, +0, +0]), ('c-', [+0, +0, -c], [+0, +0, -1]), ]) def draw_sc(ax, size, view, jitter, steps=None, alpha=1): def draw_sc(axes, size, view, jitter, steps=None, alpha=1): """Draw points for simple cubic paracrystal""" atoms = _build_sc() _draw_crystal(ax, size, view, jitter, atoms=atoms) def draw_fcc(ax, size, view, jitter, steps=None, alpha=1): _draw_crystal(axes, size, view, jitter, atoms=atoms) def draw_fcc(axes, size, view, jitter, steps=None, alpha=1): """Draw points for face-centered cubic paracrystal""" # Build the simple cubic crystal # y and z planes can be generated by substituting x for y and z respectively atoms.extend(zip(x+y+z, y+z+x, z+x+y)) _draw_crystal(ax, size, view, jitter, atoms=atoms) def draw_bcc(ax, size, view, jitter, steps=None, alpha=1): _draw_crystal(axes, size, view, jitter, atoms=atoms) def draw_bcc(axes, size, view, jitter, steps=None, alpha=1): """Draw points for body-centered cubic paracrystal""" # Build the simple cubic crystal ) atoms.extend(zip(x, y, z)) _draw_crystal(ax, size, view, jitter, atoms=atoms) def _draw_crystal(ax, size, view, jitter, steps=None, alpha=1, atoms=None): _draw_crystal(axes, size, view, jitter, atoms=atoms) def _draw_crystal(axes, size, view, jitter, steps=None, alpha=1, atoms=None): atoms, size = np.asarray(atoms, 'd').T, np.asarray(size, 'd') x, y, z = atoms*size[:, None] x, y, z = transform_xyz(view, jitter, x, y, z) ax.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o') ax.scatter(x[1:], y[1:], z[1:], c='r', marker='o') axes.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o') axes.scatter(x[1:], y[1:], z[1:], c='r', marker='o') def _build_sc(): return atoms def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1): """Draw a parallelepiped.""" a, b, c = size x, y, z = transform_xyz(view, jitter, x, y, z) ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) # Draw pink face on box. z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) draw_labels(ax, view, jitter, [ axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) draw_labels(axes, view, jitter, [ ('c+', [+0, +0, +c], [+1, +0, +0]), ('c-', [+0, +0, -c], [+0, +0, -1]), ]) def draw_sphere(ax, radius=10., steps=100): def draw_sphere(axes, radius=10., steps=100): """Draw a sphere""" u = np.linspace(0, 2 * np.pi, steps) y = radius * np.outer(np.sin(u), np.sin(v)) z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w') def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), draw_shape=draw_parallelepiped): """ if dpsi == 0: cloud = [v for v in cloud if v[2] == 0] draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8) scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] for point in cloud: delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] draw_shape(ax, size, view, delta, alpha=0.8) draw_shape(axes, size, view, delta, alpha=0.8) for v in 'xyz': a, b, c = size lim = np.sqrt(a**2 + b**2 + c**2) getattr(ax, 'set_'+v+'lim')([-lim, lim]) getattr(ax, v+'axis').label.set_text(v) getattr(axes, 'set_'+v+'lim')([-lim, lim]) getattr(axes, v+'axis').label.set_text(v) PROJECTIONS = [ 'azimuthal_equal_area', ] def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', projection='equirectangular'): """ Should allow free movement in theta, but phi is distorted. """ # TODO: try Kent distribution instead of a gaussian warped by projection t = np.linspace(-1, 1, n) weights = np.ones_like(t) def _rotate(theta_i, phi_j): return Rx(phi_j)*Ry(theta_i) def _weight(theta_i, phi_j, wi, wj): return wi*wj*abs(cos(radians(theta_i))) def _weight(theta_i, phi_j, w_i, w_j): return w_i*w_j*abs(cos(radians(theta_i))) elif projection == 'sinusoidal':  #define PROJECTION 2 def _rotate(theta_i, phi_j): #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) return Rx(longitude)*Ry(latitude) def _weight(theta_i, phi_j, wi, wj): def _weight(theta_i, phi_j, w_i, w_j): latitude = theta_i scale = cos(radians(latitude)) w = 1 if abs(phi_j) < abs(scale)*180 else 0 return w*wi*wj active = 1 if abs(phi_j) < abs(scale)*180 else 0 return active*w_i*w_j elif projection == 'guyou':  #define PROJECTION 3  (eventually?) def _rotate(theta_i, phi_j): from guyou import guyou_invert from .guyou import guyou_invert #latitude, longitude = guyou_invert([theta_i], [phi_j]) longitude, latitude = guyou_invert([phi_j], [theta_i]) return Rx(longitude[0])*Ry(latitude[0]) def _weight(theta_i, phi_j, wi, wj): return wi*wj def _weight(theta_i, phi_j, w_i, w_j): return w_i*w_j elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry def _rotate(theta_i, phi_j): #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) return Rz(longitude)*Ry(latitude) def _weight(theta_i, phi_j, wi, wj): def _weight(theta_i, phi_j, w_i, w_j): # Weighting for each point comes from the integral: #     \int\int I(q, lat, log) sin(lat) dlat dlog # the entire sphere, and treats theta and phi identically. latitude = sqrt(theta_i**2 + phi_j**2) w = sin(radians(latitude))/latitude if latitude != 0 else 1 return w*wi*wj if latitude < 180 else 0 weight = sin(radians(latitude))/latitude if latitude != 0 else 1 return weight*w_i*w_j if latitude < 180 else 0 elif projection == 'azimuthal_equal_area': def _rotate(theta_i, phi_j): R = min(1, sqrt(theta_i**2 + phi_j**2)/180) latitude = 180-degrees(2*np.arccos(R)) radius = min(1, sqrt(theta_i**2 + phi_j**2)/180) latitude = 180-degrees(2*np.arccos(radius)) longitude = degrees(np.arctan2(phi_j, theta_i)) #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) return Rz(longitude)*Ry(latitude) def _weight(theta_i, phi_j, wi, wj): def _weight(theta_i, phi_j, w_i, w_j): latitude = sqrt(theta_i**2 + phi_j**2) w = sin(radians(latitude))/latitude if latitude != 0 else 1 return w*wi*wj if latitude < 180 else 0 weight = sin(radians(latitude))/latitude if latitude != 0 else 1 return weight*w_i*w_j if latitude < 180 else 0 else: raise ValueError("unknown projection %r"%projection) for theta_i in dtheta*t for phi_j in dphi*t]) # select just the active points (i.e., those with phi < 180 w = np.array([_weight(theta_i, phi_j, wi, wj) for wi, theta_i in zip(weights, dtheta*t) for wj, phi_j in zip(weights, dphi*t)]) w = np.array([_weight(theta_i, phi_j, w_i, w_j) for w_i, theta_i in zip(weights, dtheta*t) for w_j, phi_j in zip(weights, dphi*t)]) #print(max(w), min(w), min(w[w>0])) points = points[:, w > 0] w /= max(w) if 0: # Kent distribution points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t]) xp, yp, zp = [np.array(v).flatten() for v in points] kappa = max(1e6, radians(dtheta)/(2*pi)) beta = 1/max(1e-6, radians(dphi)/(2*pi))/kappa w = exp(kappa*zp) #+ beta*(xp**2 + yp**2) print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w)) #w /= abs(cos(radians( #w /= sum(w) # rotate relative to beam points = orient_relative_to_beam(view, points) x, y, z = [np.array(v).flatten() for v in points] #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51)) ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) def draw_labels(ax, view, jitter, text): axes.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) def draw_labels(axes, view, jitter, text): """ Draw text at a particular location. for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): zdir = np.asarray(zdir).flatten() ax.text(p[0], p[1], p[2], label, zdir=zdir) axes.text(p[0], p[1], p[2], label, zdir=zdir) # Definition of rotation matrices comes from wikipedia: def Rx(angle): """Construct a matrix to rotate points about *x* by *angle* degrees.""" a = radians(angle) angle = radians(angle) R = [[1, 0, 0], [0, +cos(a), -sin(a)], [0, +sin(a), +cos(a)]] [0, +cos(angle), -sin(angle)], [0, +sin(angle), +cos(angle)]] return np.matrix(R) def Ry(angle): """Construct a matrix to rotate points about *y* by *angle* degrees.""" a = radians(angle) R = [[+cos(a), 0, +sin(a)], angle = radians(angle) R = [[+cos(angle), 0, +sin(angle)], [0, 1, 0], [-sin(a), 0, +cos(a)]] [-sin(angle), 0, +cos(angle)]] return np.matrix(R) def Rz(angle): """Construct a matrix to rotate points about *z* by *angle* degrees.""" a = radians(angle) R = [[+cos(a), -sin(a), 0], [+sin(a), +cos(a), 0], angle = radians(angle) R = [[+cos(angle), -sin(angle), 0], [+sin(angle), +cos(angle), 0], [0, 0, 1]] return np.matrix(R) return data[offset], data[-1] def draw_scattering(calculator, ax, view, jitter, dist='gaussian'): def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): """ Plot the scattering for the particular view. *calculator* is returned from :func:build_model.  *ax* are the 3D axes *calculator* is returned from :func:build_model.  *axes* are the 3D axes on which the data will be plotted.  *view* and *jitter* are the current orientation and orientation dispersity.  *dist* is one of the sasmodels level[level < 0] = 0 colors = plt.get_cmap()(level) ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) elif 1: ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, levels=np.linspace(vmin, vmax, 24)) else: ax.pcolormesh(qx, qy, Iqxy) axes.pcolormesh(qx, qy, Iqxy) def build_model(model_name, n=150, qmax=0.5, **pars): plt.gcf().canvas.set_window_title(projection) #gs = gridspec.GridSpec(2,1,height_ratios=[4,1]) #ax = plt.subplot(gs[0], projection='3d') ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') #axes = plt.subplot(gs[0], projection='3d') axes = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection ax.axis('square') axes.axis('square') except Exception: pass ## add control widgets to plot axtheta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) axpsi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) stheta = Slider(axtheta, 'Theta', -90, 90, valinit=theta) sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) axdtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) axdpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) # Note: using ridiculous definition of rectangle distribution, whose width # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep # the maximum width to 90. dlimit = DIST_LIMITS[dist] sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi) sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi) sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) limit = [0, 0, 0.5, 5][dims] jitter = [0 if v < limit else v for v in jitter] ax.cla() draw_beam(ax, (0, 0)) draw_jitter(ax, view, jitter, dist=dist, size=size, draw_shape=draw_shape) #draw_jitter(ax, view, (0,0,0)) draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection) draw_scattering(calculator, ax, view, jitter, dist=dist) axes.cla() draw_beam(axes, (0, 0)) draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape) #draw_jitter(axes, view, (0,0,0)) draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection) draw_scattering(calculator, axes, view, jitter, dist=dist) plt.gcf().canvas.draw()