Relax disp.spin lock offset+field figure
Jump to navigation
Jump to search
Script to produce figure
#############
# Made by Troels E. Linnet
# PhD student
# Copenhagen University
# SBiNLab, Structural Biology and NMR Laboratory.
# March 2014
#
# Trying to reproduce Figure 1 in:
#
# Palmer, A.G. & Massi, F. (2006). Characterization of the dynamics of biomacromolecules using rotating-frame spin relaxation NMR spectroscopy.
# Chem. Rev. 106, 1700-1719 http://dx.doi.org/10.1021/cr04042875
#
# This script is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#######################
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Ellipse, PathPatch, FancyArrowPatch
from mpl_toolkits.mplot3d import Axes3D, proj3d
import mpl_toolkits.mplot3d.art3d as art3d
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import pylab
######## Define helper functions
class Arrow3D(FancyArrowPatch):
def __init__(self, xs, ys, zs, *args, **kwargs):
FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
self._verts3d = xs, ys, zs
def draw(self, renderer):
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
FancyArrowPatch.draw(self, renderer)
def rotation_matrix(d):
"""
Calculates a rotation matrix given a vector d. The direction of d
corresponds to the rotation axis. The length of d corresponds to
the sin of the angle of rotation.
Variant of: http://mail.scipy.org/pipermail/numpy-discussion/2009-March/040806.html
"""
sin_angle = np.linalg.norm(d)
if sin_angle == 0:
return np.identity(3)
d /= sin_angle
eye = np.eye(3)
ddt = np.outer(d, d)
skew = np.array([[ 0, d[2], -d[1]],
[-d[2], 0, d[0]],
[d[1], -d[0], 0]], dtype=np.float64)
M = ddt + np.sqrt(1 - sin_angle**2) * (eye - ddt) + sin_angle * skew
return M
def pathpatch_2d_to_3d(pathpatch, z = 0, normal = 'z'):
"""
Transforms a 2D Patch to a 3D patch using the given normal vector.
The patch is projected into they XY plane, rotated about the origin
and finally translated by z.
"""
if type(normal) is str: #Translate strings to normal vectors
index = "xyz".index(normal)
normal = np.roll((1,0,0), index)
normal /= np.linalg.norm(normal) #Make sure the vector is normalised
path = pathpatch.get_path() #Get the path and the associated transform
trans = pathpatch.get_patch_transform()
path = trans.transform_path(path) #Apply the transform
pathpatch.__class__ = art3d.PathPatch3D #Change the class
pathpatch._code3d = path.codes #Copy the codes
pathpatch._facecolor3d = pathpatch.get_facecolor #Get the face color
verts = path.vertices #Get the vertices in 2D
d = np.cross(normal, (0, 0, 1)) #Obtain the rotation vector
M = rotation_matrix(d) #Get the rotation matrix
pathpatch._segment3d = np.array([np.dot(M, (x, y, 0)) + (0, 0, z) for x, y in verts])
def pathpatch_translate(pathpatch, delta):
"""
Translates the 3D pathpatch by the amount delta.
"""
pathpatch._segment3d += delta
def update_position(e):
print "From update position"
#Transform co-ordinates to get new 2D projection
tX, tY, _ = proj3d.proj_transform(dataX, dataY, dataZ, ax.get_proj())
for i in range(len(dataX)):
label = labels[i]
label.xy = tX[i],tY[i]
label.update_positions(fig.canvas.renderer)
fig.canvas.draw()
return
#########
fig = plt.figure()
ax=fig.gca(projection='3d')
# Draw a Center point
ax.scatter([0],[0],[0],c="k",s=50)
# Draw circle
circle = Circle((0, 0), 1)
circle.set_color("k")
circle.set_alpha(80)
#circle.set_fill(False)
ax.add_patch(circle)
#art3d.pathpatch_2d_to_3d(circle, z=0, zdir='y')
pathpatch_2d_to_3d(circle, z=0, normal = (1, 0, 1) )
#Input 3D Data for Text
# Sx Sy Sx Szp Sxp w1 O1 we1 O_av w_e
data = np.array([[1,0,0], [0,-1,0], [0,0,1], [0.975,0,1.17], [0.75,0,-0.75], [0.7,0,0], [0,0,0.3], [0.7,0,0.3], [0.0,0,0.84], [0.7,0,0.84]])
textlab = [r"S$_x$", r"S$_y$=S$_y$'", r"S$_z$", r"S$_z$'", r"S$_x$'", r"$\omega_1$", r"$\Omega_1$", r"$\omega_{e1}$", r"$\bar{\Omega}$", r"$\omega_{e}$"]
#Separate into X, Y, Z for greater clarity
dataX = data[:,0]
dataY = data[:,1]
dataZ = data[:,2]
#3D scatter plot
ax.scatter(dataX, dataY, dataZ, marker = 'o', c='k', s=1)
#Transform co-ordinates to get initial 2D projection
tX, tY, _ = proj3d.proj_transform(dataX, dataY, dataZ, ax.get_proj())
#Array of labels
labels = []
#Loop through data points to initially annotate scatter plot
#and populate labels array
for i in range(len(dataX)):
#text='['+str(int(dataX[i]))+','+str(int(dataY[i]))+','+str(int(dataZ[i]))+']'
text=textlab[i]
label = ax.annotate(text,
xycoords='data',
xy = (tX[i], tY[i]), xytext = (-20, 20),
textcoords = 'offset points', ha = 'right', va = 'top', fontsize=12,
bbox = dict(boxstyle = 'round,pad=0.5', fc = 'grey', alpha = 0.8),
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))
labels.append(label)
# Make Sx
i = 0
Sx = Arrow3D([-1, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(Sx)
# Make Sy
i += 1
Sy = Arrow3D([0, dataX[i]],[1, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(Sy)
# Make Sz
i += 1
Sz = Arrow3D([0, dataX[i]],[0, dataY[i]], [-1, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(Sz)
# Make Szp
i += 1
Szp = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=3, arrowstyle="-|>", color="k")
ax.add_artist(Szp)
# Make Sxp
i += 1
Sxp = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=3, arrowstyle="-|>", color="k")
ax.add_artist(Sxp)
# Make w1
i += 1
w1 = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(w1)
# Make O1
i += 1
O1 = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(O1)
# Make we1
i += 1
we1 = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="b")
ax.add_artist(we1)
# Make O_av
i += 1
O_av = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(O_av)
# Make we
i += 1
we = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="r")
ax.add_artist(we)
ax.set_xlim3d(-1, 1)
ax.set_ylim3d(-1, 1)
ax.set_zlim3d(-1, 1)
fig.canvas.mpl_connect('button_release_event', update_position)
plt.show()