relax disp.spin lock offset+field figure

Reference to original figure

Comparing to Figure 1 and 10 the reference.

Palmer, A.G. & Massi, F. (2006). Characterization of the dynamics of biomacromolecules using rotating-frame spin relaxation NMR spectroscopy. Chem. Rev. 106, 1700-1719 

DOI

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="-