relax disp.spin lock offset+field figure

Reference to original figure

Comparing to Figure 1 and 10 the reference:

• Palmer, 3rd, A. G. and Massi, F. (2006). Characterization of the dynamics of biomacromolecules using rotating-frame spin relaxation NMR spectroscopy. Chem. Rev., 106(5), 1700-1719. (DOI: 10.1021/cr0404287)

Try to reproduce Figure 1.

Script to produce figure

Python script: Reproduce figure 1 of Palmer and Massi, 2006 using matplotlib.

#############
# 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)
    tX_t, tY_t, _ = proj3d.proj_transform(dataX_t, dataY_t, dataZ_t, ax.get_proj())
    for i in range(len(dataX_t)):
        label = labels_t[i]
        label.xy = tX_t[i],tY_t[i]
        label.update_positions(fig.canvas.renderer)
    fig.canvas.draw()
    return

#########

fig = plt.figure(figsize=(12, 12))
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 = [1,0,0]
Sy = [0,-1,0]
Sz = [0,0,1]
Szp = [1.1,0,1.1]
Sxp = [0.75,0,-0.75]
w1 = [0.7,0,0]
O1 = [0,0,0.3]
we1 = [w1[0],0,O1[2]]
O_av = [0.0,0,Szp[2]/Szp[0]*w1[0]]
we = [w1[0],0,Szp[2]/Szp[0]*w1[0]]

data = np.array([Sx, Sy, Sz , Szp , Sxp, w1, O1, we1, O_av, we])
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)
#######

#######
#Input 3D Data for Ttheta

t1 = [0.5*we1[0],0,0.5*we1[2]]
t = [0.5*we[0],0,0.5*we[2]]

data_t = np.array([t1, t])
textlab_t = [r"$\theta_1$", r"$\theta$"]

#Separate into X, Y, Z for greater clarity
dataX_t = data_t[:,0]
dataY_t = data_t[:,1]
dataZ_t = data_t[:,2]

#3D scatter plot
ax.scatter(dataX_t, dataY_t, dataZ_t, marker = 'o', c='k', s=1)

#Transform co-ordinates to get initial 2D projection
tX_t, tY_t, _ = proj3d.proj_transform(dataX_t, dataY_t, dataZ_t, ax.get_proj())

#Array of labels
labels_t = []
#Loop through data points to initially annotate scatter plot
#and populate labels array
for i in range(len(dataX_t)):
    #text='['+str(int(dataX[i]))+','+str(int(dataY[i]))+','+str(int(dataZ[i]))+']'
    text=textlab_t[i]
    label = ax.annotate(text,
            xycoords='data',
            xy = (tX[i], tY[i]), xytext = (-100, 50),
            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 = 'angle3,angleA=0,angleB=90'))
    labels_t.append(label)
#######



# Make Sx
x, y, z = Sx
Sx = Arrow3D([-1, x],[0, y], [0, z], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(Sx)

# Make Sy
x, y, z = Sy
fSy = Arrow3D([0, x],[1, y], [0, z], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(fSy)

# Make Sz
x, y, z = Sz
fSz = Arrow3D([0, x],[0, y], [-1, z], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(fSz)

# Make Szp
x, y, z = Szp
fSzp = Arrow3D([0, x],[0, y], [0, z], mutation_scale=20, lw=3, arrowstyle="-|>", color="k")
ax.add_artist(fSzp)

# Make Sxp
x, y, z = Sxp
fSxp = Arrow3D([0, x],[0, y], [0, z], mutation_scale=20, lw=3, arrowstyle="-|>", color="k")
ax.add_artist(fSxp)

# Make w1
x, y, z = w1
fw1 = Arrow3D([0, x],[0, y], [0, z], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(fw1)

# Make O1
x, y, z = O1
fO1 = Arrow3D([0, x],[0, y], [0, y], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(fO1)

# Make we1
x, y, z = we1
fwe1 = Arrow3D([0, x],[0, y], [0, z], mutation_scale=20, lw=1, arrowstyle="-|>", color="b")
ax.add_artist(fwe1)

# Make O1 -> we1
fO1_we1 = Arrow3D([O1[0], we1[0]],[O1[1], we1[1]], [O1[2], we1[2]], mutation_scale=20, lw=1, arrowstyle="-", color="b")
ax.add_artist(fO1_we1)

# Make O_av
x, y, z = O_av
fO_av = Arrow3D([0, x],[0, y], [0, z], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(fO_av)

# Make we
x, y, z = we
fwe = Arrow3D([0, x],[0, y], [0, z], mutation_scale=20, lw=1, arrowstyle="-|>", color="r")
ax.add_artist(fwe)

# Make w1 -> we
fw1_we = Arrow3D([w1[0], we[0]],[w1[1], we[1]], [w1[2], we[2]], mutation_scale=20, lw=1, arrowstyle="-", color="r")
ax.add_artist(fw1_we)

# Make O_ave -> we
fO_av_we = Arrow3D([O_av[0], we[0]],[O_av[1], we[1]], [O_av[2], we[2]], mutation_scale=20, lw=1, arrowstyle="-", color="r")
ax.add_artist(fO_av_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()