Magnetic resonance: Computer simulations and visualization to connect quantum theory with classical concepts -- Online supplement (Created by Larry Engelhardt, 2015)¶
The Python code used to compute and plot the data for this article are provided below.
To access interactive/animated versions of the plots from this article, use these links:
This is an IPython notebook, which consists of multiple "cells" much like a Mathematica or Maple notebook. The cells in this notebook contain the following items:
- Source code, written in the Python programming language, that is used to simulate and visualize the phenomenon of magnetic resonance.
- Each of the plots that are included in the article, Magnetic resonance: Computer simulations and visualization to connect quantum theory with classical concepts.
- Text and equations describing what is being done in each cell of this notebook.
- Interactive plots and animations that can be accessed using the links above.
The following cell contains the first lines of Python code. These lines import the necessary libraries and set the preferences that will be used to create plots later in the notebook.
# The following line specifies to display plot "inline" within the ipython notebook
%matplotlib inline
import matplotlib as mpl # Used on the following line to set plot preferences
mpl.rcParams.update({'figure.dpi':200,'figure.figsize':(12,8),'figure.max_open_warning':200,\
'font.size':26,'xtick.labelsize':20,'ytick.labelsize':20})
from matplotlib.pyplot import * # Plotting library
from numpy import * # Numerical library
from scipy.integrate import odeint # Used to numerically solve systems of ODEs
Sec II: Theory¶
The following cell defines a new Python function which is described immediately below the cell.
# Defines the system of ODEs in terms of RATES:
def odeRates(variables, t, omega, omega0, omega1):
"""
To be used in odeint. Returns the values of the DERIVATIVES.
Takes as input: (1) dynamical variables, "variables"
(2) independent variable "t" (3) 3 parameters
"""
# Unpack the variables to make the math clearer below
aR, aI, bR, bI = variables
# DEFINES ODES:
aR_dot = +0.5 * (omega0*aI + omega1*cos(omega*t)*bI)
aI_dot = -0.5 * (omega0*aR + omega1*cos(omega*t)*bR)
bR_dot = -0.5 * (omega0*bI - omega1*cos(omega*t)*aI)
bI_dot = +0.5 * (omega0*bR - omega1*cos(omega*t)*aR)
# Return the list of rates, in the SAME ORDER as the variables:
return [aR_dot, aI_dot, bR_dot, bI_dot]
In the cell above, a Python function named odeRates is defined which computes the rates, $$\dot{a}_R = +\frac{1}{2}\left[ \omega_0 a_I + \omega_1 \cos(\omega t)b_I \right],$$ $$\dot{a}_I = -\frac{1}{2}\left[ \omega_0 a_R + \omega_1 \cos(\omega t) b_R \right],$$ $$\dot{b}_R = -\frac{1}{2}\left[ \omega_0 b_I - \omega_1 \cos(\omega t) a_I \right],$$ $$\dot{b}_I = +\frac{1}{2}\left[ \omega_0 b_R - \omega_1 \cos(\omega t) a_R \right]$$ at a given moment in time, $t$, using the following four lines of code:
aR_dot = +0.5 * (omega0*aI + omega1*cos(omega*t)*bI)
aI_dot = -0.5 * (omega0*aR + omega1*cos(omega*t)*bR)
bR_dot = -0.5 * (omega0*bI - omega1*cos(omega*t)*aI)
bI_dot = +0.5 * (omega0*bR - omega1*cos(omega*t)*aR)
The correspondence between the mathematical equations and these four lines of code should be very clear. Note that all of the variables and parameters that appear on the right-hand side of these equations are provided as input to this function on the line
def odeRates(variables, t, omega, omega0, omega1):
and the computed rates are returned as the output of the function on the line
return [aR_dot, aI_dot, bR_dot, bI_dot]
The following cell defines another function, named prob_down, in which these rates will be used to compute $a_R(t)$, $a_I(t)$, $b_R(t)$, and $b_I(t)$. The function prob_down is described in more detail below.
def prob_down(psi0, tArray, parameters):
# Solves the ODEs (odeint imported from scipy.integrate)
results = odeint(odeRates, psi0, tArray, args=parameters)
return results[:,2]**2 + results[:,3]**2 # p_b (probability of spin down)
The function above takes three variables (also listed in Sec. III of the paper) as input:
- psi0 will contain the initial values: $a_R(0)$, $a_I(0)$, $b_R(0)$, and $b_I(0)$
- tArray will contain an array of many uniformly-spaced values of time, $t$.
- parameters will contain the values of the constant parameters: $\omega$, $\omega_0$, and $\omega_1$.
Using these variables, and function odeRates that is defined above, the following line computes the value of $a_R(t)$, $a_I(t)$, $b_R(t)$, and $b_I(t)$ for every value of $t$ in the array tArray:
results = odeint(odeRates, psi0, tArray, args=parameters)
After this line is executed, results is a two-dimensional array. The first index of the array corresponds to time, and the second index of the array identifies the variable of interest:
- Index = 0 for $a_R$
- Index = 1 for $a_I$
- Index = 2 for $b_R$
- Index = 3 for $b_I$.
Finally, this function returns the following:
return results[:,2]**2 + results[:,3]**2 # p_b (probability of spin down)
This returns an array containing $p_b(t) = b_R^2(t)+b_I^2(t)$, which is the probability of measuring "spin down" at time $t$ for every value of $t$ in the array tArray.
Sec. III: Simulations¶
The following cell defines the parameters that will be used for the simulations:
## PARAMETERS TO BE ADJUSTED ##
# Parameters defining the magnetic field:
omega0 = 1.0; omega1 = 0.1 # Magnitudes (GHz)
omega = 1.0 # Angular frequency (GHz)
# Will be used to define time array:
tMax = 200; dt = 0.1 # (in nanoseconds)
# Initial conditions (real & imaginary parts):
aR = 1.0; aI = 0.0 # a: Up (+)
bR = 0.0; bI = 0.0 # b: Down (-)
################################
# The parameters above generate these:
tArray = arange(0, tMax, dt) # Array of time values
initialCond = (aR, aI, bR, bI) # Initial conditions
param = (omega, omega0, omega1) # Other parameters
p_b = prob_down(initialCond, tArray, param)
plot(tArray, p_b, linewidth=3)
xlabel('Time (nanoseconds)')
ylabel('Probability, $p_b(t)$')
grid(True)
savefig('fig1.eps'); show() # Saves image to disk and shows in notebook
Figure 1: Simulation of $p_b(t)$ using the initial conditions $a_R=1$ and $a_I=b_R=b_I=0$ (spin up at time $t=0$) and parameter values $\omega = \omega_0 = 1.0$ GHz and $\omega_1=0.1$ GHz. An interactive version of this plot is provided here.
Below we repeat the simulation, but now we calculate the flipping probability, $p_{flip}$, which is the maximum value of $p_b(t)$ that is achieved for any value of $t$. This is repeated for many different values of $\omega$ in order to plot $p_{flip}$ versus $\omega$.
The cell below defines the parameter values to be used in order to calculate and plot $p_{flip}(\omega)$.
omegaMin = 0.9; omegaMax = 1.1 # Used to create array below
omega0 = 1.0 # Magnitude in GHz
# Use logarithmic spacing to get more data points in the middle of the range:
omegaRight = logspace(log10(0.001), log10(0.5*(omegaMax-omegaMin)), 20)
omegaLeft = -omegaRight[::-1]
omegaLeft = append(omegaLeft, 0)
omegaArray = append(omegaLeft, omegaRight)
omegaArray = omegaArray + 0.5*(omegaMax+omegaMin)
# Initial conditions (real & imaginary parts):
aR = 1.0; aI = 0.0 # a: Up (+)
bR = 0.0; bI = 0.0 # b: Down (-)
initialCond = (aR, aI, bR, bI) # Initial conditions
N = size(omegaArray) # Number of omega values
pArray = zeros(N) # To store probability for each omega
The following cell calculates and plots $p_{flip}$ versus $\omega$ for two different values of $\omega_1$.
omega1 = 0.10 # First simulate using this value of omega1
tMax = 2.1*pi/omega1
tArray = linspace(0, tMax, 100)
for i in range(N): # Loop through many value of omega
omega = omegaArray[i]
param = (omega, omega0, omega1)
p_b = prob_down(initialCond, tArray, param)
pArray[i] = max(p_b)
plot(omegaArray, pArray, '-', linewidth=3)
omega1 = 0.01 # Repeat using a different value of omega1
tMax = 2.1*pi/omega1
tArray = linspace(0, tMax, 100)
for i in range(N): # Loop through many value of omega
omega = omegaArray[i]
param = (omega, omega0, omega1)
p_b = prob_down(initialCond, tArray, param)
pArray[i] = max(p_b)
plot(omegaArray, pArray, 'r--', linewidth=3)
xlabel('$\omega$ (GHz)')
ylabel('Flipping Probability')
xlim([0.9, 1.1]); ylim([1e-5, 1.1])
grid(True)
savefig('fig2.eps'); show() # Saves image to disk and shows in notebook
Figure 2: Simulation of $p_{flip}(\omega)$ using the initial conditions $a_R=1$ and $a_I=b_R=b_I=0$ (spin up at time $t=0$) with parameter values $\omega_0=1.0$ GHz and two different values of $\omega_1$: $\omega_1=0.1$ GHz (solid curve) and 0.01 GHz (dashed curve). An interactive version of this plot is provided here.
Sec IV: Spectral analysis¶
In the cell below, we define a function, computePowerSpectrum, that will compute and return the power spectrum, $Y(f)$ (the modulus squared of the Fourier transform), as a function of frequency, $f$. In the cell below, the line
powerSpectrum[0] = 0
subtracts off the DC offset that is always present in $p_b(t)$ since $p_b(t)$ can never be negative (since it is a probability).
from numpy.fft import rfft, rfftfreq # Import functions for Fourier transform
def computePowerSpectrum(initialCond, tArray, param):
p_b = prob_down(initialCond, tArray, param)
fourierTransform = rfft(p_b)
powerSpectrum = abs(fourierTransform)**2
powerSpectrum[0] = 0 # Remove the DC (zero frequency) vertical offset
n = size(tArray)
frequencies = rfftfreq(n, d=dt)
omegaPrime = 2*pi*frequencies
return omegaPrime, powerSpectrum
# Initial conditions (real & imaginary parts):
aR = 1.0; aI = 0.0 # a: Up (+)
bR = 0.0; bI = 0.0 # b: Down (-)
initialCond = (aR, aI, bR, bI) # Initial conditions
omega0 = 1.0; omega1 = 0.10 # Magnitudes (GHz)
tMax = 10000; dt = 0.1 # (in nanoseconds)
tArray = arange(0, tMax, dt) # Array of time values
figure() # Creates a single figure that will hold multiple curves
omega = 0.1
param = (omega, omega0, omega1)
omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
plot(omegaPrime, Y, '-', linewidth=2, label='$\omega$ = '+str(omega)+' GHz')
omega = 0.2
param = (omega, omega0, omega1)
omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
plot(omegaPrime, Y, '--', linewidth=3, label='$\omega$ = '+str(omega)+' GHz')
box = legend(fontsize=22, frameon = 1)
box.get_frame().set_facecolor('white') # Solid white background for legend
xlabel('$\omega^{\prime}$ (GHz)')
ylabel('Power Spectrum (arb. units)')
xlim([0, 1.4]); ylim([0, 17000])
grid(True)
Figure 3 (main figure): The power spectrum, $Y(\omega^{\prime})$, versus $\omega^{\prime}$ using parameter values $\omega_0=1.0$ GHz and $\omega_1=0.1$ GHz for two different values of $\omega$ and arbitrary initial conditions. An interactive version of this plot is provided here.
omegaArray = [0.6, 0.7, 0.8, 0.9]
lineStyles = ['-.', '--', ':', '-']
numCurves = len(omegaArray)
omega0 = 1.0; omega1 = 0.10 # Magnitudes (GHz)
tMax = 10000; dt = 0.1 # (in nanoseconds)
tArray = arange(0, tMax, dt) # Array of time values
figure()
for i in range(numCurves):
omega = omegaArray[i]
param = (omega, omega0, omega1)
omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
plot(omegaPrime, Y, lineStyles[i], linewidth=3, label='$\omega$ = '+str(omega)+' GHz')
box = legend(fontsize=22, frameon = 1)
box.get_frame().set_facecolor('white') # Solid white background for legend
yscale('log')
xlabel('$\omega^{\prime}$ (GHz)')
ylabel('Power Spectrum (arb. units)')
xlim([0, 0.5]); ylim([10**-1, 10**9])
grid(True)
show()
Figure 3 (inset): The power spectrum, $Y(\omega^{\prime})$, versus $\omega^{\prime}$ using parameter values $\omega_0=1.0$ GHz and $\omega_1=0.1$ GHz for four different values of $\omega$ and arbitrary initial conditions. An interactive version of this plot is provided here.
In the next cell, we now combine the two previous cell in order to create a single figure with one plot as the main figure and one plot as an inset.
omega0 = 1.0; omega1 = 0.10 # Magnitudes (GHz)
tMax = 10000; dt = 0.1 # (in nanoseconds)
tArray = arange(0, tMax, dt) # Array of time values
figure() # Creates a single figure that will hold multiple curves
###### MAIN FIGURE #######
omega = 0.1
param = (omega, omega0, omega1)
omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
plot(omegaPrime, Y, '-', linewidth=2, label='$\omega$ = '+str(omega)+' GHz')
omega = 0.2
param = (omega, omega0, omega1)
omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
plot(omegaPrime, Y, '--', linewidth=3, label='$\omega$ = '+str(omega)+' GHz')
box = legend(fontsize=22, frameon = 1)
box.get_frame().set_facecolor('white') # Solid white background for legend
xlabel('$\omega^{\prime}$ (GHz)')
ylabel('Power Spectrum (arb. units)')
xlim([0, 1.4]); ylim([0, 17000])
####### INSET ########
a = axes([0.21, 0.35, .3, .5], axisbg='white') # Inset: position & size
omegaArray = [0.6, 0.7, 0.8, 0.9]
lineStyles = ['-.', '--', ':', '-']
numCurves = len(omegaArray)
for i in range(numCurves):
omega = omegaArray[i]
param = (omega, omega0, omega1)
omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
if i == 2: # Make the dotted curve thicker & black to show up better
plot(omegaPrime, Y, lineStyles[i], color='black', linewidth=4, label='$\omega$ = '+str(omega)+' GHz')
else:
plot(omegaPrime, Y, lineStyles[i], linewidth=3, label='$\omega$ = '+str(omega)+' GHz')
box = legend(fontsize=14, frameon = 1)
box.get_frame().set_facecolor('white') # Solid white background for legend
yscale('log')
xlim([0, 0.5]); ylim([10**-1, 10**10])
grid(True)
savefig('fig3.eps'); show() # Saves image to disk and shows in notebook
Figure 3: The power spectrum, $Y(\omega^{\prime})$, versus $\omega^{\prime}$ using parameter values $\omega_0=1.0$ GHz and $\omega_1=0.1$ GHz for several different values of $\omega$ and arbitrary initial conditions. An interactive version of this plot is provided here.
In the following cell, we will again compute the power spectrum, just as we did above, but now for $\omega=\omega_0$.
omega1_Array = [0.1, 0.2, 0.5]
lineStyles = ['-', '--', ':']
numCurves = len(omega1_Array)
omega0 = 1.0; omega = 1.0 # Magnitudes (GHz)
tMax = 10000; dt = 0.1 # (in nanoseconds)
tArray = arange(0, tMax, dt) # Array of time values
figure()
for i in range(numCurves):
omega1 = omega1_Array[i]
param = (omega, omega0, omega1)
omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
if i == 2: # Make the dotted curve thicker & black to show up better
plot(omegaPrime, Y, lineStyles[i], color='black', linewidth=4, label='$\omega_1$ = '+str(omega1)+' GHz')
else:
plot(omegaPrime, Y, lineStyles[i], linewidth=3, label='$\omega_1$ = '+str(omega1)+' GHz')
box = legend(fontsize=22, frameon = 1)
box.get_frame().set_facecolor('white') # Solid white background for legend
yscale('log')
xlabel('$\omega^{\prime}$ (GHz)')
ylabel('Power Spectrum (arb. units)')
xlim([0, 2.5]); ylim([10**1, 10**9])
grid(True)
savefig('fig4.eps'); show() # Saves image to disk and shows in notebook
Figure 4: The power spectrum, $Y(\omega^{\prime})$, versus $\omega^{\prime}$ for $\omega=\omega_0=1.0$ GHz (at resonance) for three different values of $\omega_1$ and arbitrary initial conditions. An interactive version of this plot is provided here.
Section V: Quantum spin vector and classical connection¶
from IPython.display import IFrame
IFrame('MagneticResonance_3D.html', width='100%', height=720)
Interactive version of Figures 5 and 6: The simulation above can be used to create animated versions of Figs. 5 and 6 and to explore how these dynamics depend on each of the parameters. This simulation was written in GlowScript which produces HTML+javascript code that can be run in any browser.
The GlowScript source code is provided below as well as Python versions that use the ivisual library in order be able to modify and execute the code within an IPython notebook. (The Python versions will not work in the HTML version of this notebook; they require IPython.)
GlowScript 1.1 JavaScript
scene.width = 400
scene.height = 250
scene.background = color.white
scene.range = 0.55
scene.background = color.white
scene.up = vec(0,0,1)
scene.forward = vec(0,1,0)
distant_light( {direction:vec(0,-1,1)} )
distant_light( {direction:vec(-1,0,1)} )
var time_label = label( {pixel_pos:true, pos:vec(scene.width-50,10,0), box:false, text:'t = 0 ns'} )
// Parameters
var omega1 = 0.1;
var omega0 = 1.0;
var omega = 0.0;
var aR = 1.0;
var aI = 0.0;
var bR = 0.0;
var bI = 0.0;
var psi = [aR, aI, bR, bI];
var t = 0, dt = 0.1;
function make_title() {
scene.title.html('Simulates the quantum spin vector <<b>s</b>> (black arrow) in\n'+
'a magnetic field (red arrow): <font color="red"><b>H</b></font> =[H<sub>1</sub>cos(ωt), 0, H<sub>0</sub>]\n\n'+
'To vary H<sub>0</sub> and H<sub>1</sub>, we define '+
'ω<sub>0</sub> = gμ<sub>B</sub>H<sub>0</sub>/ℏ and '+
'ω<sub>1</sub> = gμ<sub>B</sub>H<sub>1</sub>/ℏ.\n'+
'- CLICK "BEGIN", THEN VARY ω<sub>1</sub> USING THE FIRST SLIDER.\n'+
'- TO OSCILLATE AUTOMATICALLY, SET ω > 0.')
}
make_title()
// 3D objects
var spin = arrow( {pos:vec(0,0,0), axis_and_length:vec(0,0,0.5), color:color.black, shaftwidth:0.02} )
var field = arrow( {pos:vec(0,0,0), axis_and_length:vec(0.5*omega1,0,0.5*omega0), color:color.red, shaftwidth:0.01} )
var trail = curve({color:color.green, radius: 0.01})
//// CREATING THE CONTROLS ////
$("<div id='controls1' style='display:flex; align-items:center; height:25px' />").appendTo(scene.caption)
// Button to pause (running = false) or play (running = true)
var running = false;
$("<button id='play_button' style='width:60px; margin:10px; display:inline-block;' />").text("BEGIN")
.appendTo(controls1)
.click( function() {
running = !running
if (running) $(this).text("Pause")
else $(this).text("Resume")
})
// Controlling omega1
$("<p id='omega1_slider' style='width:180px; margin:10px; display:inline-block;'></p>")
.appendTo(controls1)
$("<p id='omega1_label' style='display:inline-block; margin:10px;'></p>")
.html("ω<sub>1</sub> = "+omega1+" GHz")
.appendTo(controls1)
$("<div id='controls2' style='display:flex; align-items: center; height:30px;' />").appendTo(scene.caption)
$("<p style='display:inline-block; margin:10px;'></p>")
.html("Enter a number to change the max value of ω<sub>1</sub>:")
.appendTo(controls2)
$("<input id='omega1_max_input' type='text' style='width:20px; margin:10px; display:inline-block;' />")
.appendTo(controls2)
.keyup(function(){
var slider_edges = parseFloat($(this).val());
$( "#omega1_slider" ).slider( "option", "min", -slider_edges);
$( "#omega1_slider" ).slider( "option", "max", +slider_edges);
});
$("<div id='controls3' style='display:flex; align-items: center; height:30px;' />").appendTo(scene.caption)
// Reset button
$("<button style='width:60px; margin:10px;' />").text("RESET")
.appendTo(controls3)
.click( function() {
running = false; $( "#play_button" ).text("BEGIN");
t = 0.0; time_label.text = 't = 0 ns';
omega1 = 0.1; omega0 = 1.0; omega = 0.0;
aR=1.0; aI=0; bR=0; bI=0;
psi = [aR, aI, bR, bI];
$( "#omega1_slider" ).slider( "value", omega1);
$( "#omega0_slider" ).slider( "value", omega0);
$( "#omega_slider" ).slider( "value", omega );
$( "#omega1_label").html("ω<sub>1</sub> = "+omega1+" GHz");
$( "#omega0_label").html("ω<sub>0</sub> = "+omega0+" GHz");
$( "#omega_label").html("ω = "+omega+" GHz");
field.axis_and_length = vec(0.5*omega1,0,0.5*omega0);
spin.axis_and_length = vec(0, 0, 0.5);
trail.clear();
scene.forward = vec(0,1,0);
$('select>option:eq(0)').prop('selected', true);
})
// Drop-down menu to select which axis points into the screen
$("<p id='axis_menu' style='display:inline-block;' > Axis pointing into the screen: </p>").appendTo(controls3)
var s = ""
s += "<option select=selected>+y</option>"
s += "<option>+x</option>"
s += "<option>-y</option>"
s += "<option>-x</option>"
$('<select/>').html(s).css({font:"sans"})
.change(function() { // Picking an axis
switch ($(this).val()) {
case "+y":
scene.forward = vec(0,1,0)
break;
case "+x":
scene.forward = vec(1,0,0)
break;
case "-y":
scene.forward = vec(0,-1,0)
break;
case "-x":
scene.forward = vec(-1,0,0)
break;
}
})
.appendTo(controls3)
$("<div id='controls4' style='display:flex; align-items: center; height:25px;' />").appendTo(scene.caption)
// Controlling omega
$("<p id='omega_slider' style='width:180px; display:inline-block; margin:10px;'></p>")
.appendTo(controls4)
$("<p id='omega_label' style='display:inline-block; margin:10px;'></p>")
.html("ω = "+omega+" GHz")
.appendTo(controls4)
$("<div id='controls5' style='display:flex; align-items: center; height:25px;' />").appendTo(scene.caption)
// Controlling omega0
$("<p id='omega0_slider' style='width:180px; display:inline-block; margin:10px;'></p>")
.appendTo(controls5)
$("<p id='omega0_label' style='display:inline-block; margin:10px;'></p>")
.html("ω<sub>0</sub> = "+omega0+" GHz")
.appendTo(controls5)
var ready = 0 // To delay
// Define sliders properties
$(function(event, ui) { $( '#omega1_slider' ).slider({
orientation: "horizontal",
min: -0.1,
max: 0.1,
value: 0.1,
step: 0.01,
slide: function(event, ui) {
omega1 = ui.value;
$( "#omega1_label").html("ω<sub>1</sub> = "+omega1+" GHz");
field.axis_and_length = vec(0.5*omega1,0,0.5*omega0)
}
})
ready++ // to prevent premature setting of slider value below
})
$(function() { $( '#omega_slider' ).slider({
orientation: "horizontal",
min: 0.0,
max: 1.5,
value: 0.0,
step: 0.01,
slide: function(event, ui) {
omega = ui.value;
$( "#omega_label").html("ω = "+omega+" GHz");
}
})
ready++ // to prevent premature setting of slider value below
})
$(function() { $( '#omega0_slider' ).slider({
orientation: "horizontal",
min: 0.0,
max: 1.5,
value: 1.0,
step: 0.01,
slide: function(event, ui) {
omega0 = ui.value;
$( "#omega0_label").html("ω<sub>0</sub> = "+omega0+" GHz");
field.axis_and_length = vec(0.5*omega1,0,0.5*omega0)
}
})
ready++ // to prevent premature setting of slider value below
})
while (ready < 2) {
rate(30,wait)
}
function display_instructions() {
var s = "In GlowScript programs:\n\n"
s += " To zoom, drag with the left+right mouse buttons,\n
or hold down the Alt key and drag,\n
or use the mouse wheel.\n\n"
s += " Rotate the camera by dragging with the right mouse button,\n
or hold down the Ctrl key and drag.\n\n"
s += "Touch screen: pinch/extend to zoom, swipe or two-finger rotate."
$("<p/>").html(s).appendTo(scene.caption)
}
// Display text below the 3D graphics:
display_instructions()
//// The rest of the code does the MATH ////
// For computations using complex matrices and vectors
get_library("http://www.numericjs.com/lib/numeric-1.2.6.js", wait)
// Computes and returns the derivative for each of the 4 ODEs
function f(psi, t) {
var aR = psi[0], aI = psi[1], bR = psi[2], bI = psi[3];
var aR_dot = +0.5 * (omega0*aI + omega1*Math.cos(omega*t)*bI);
var aI_dot = -0.5 * (omega0*aR + omega1*Math.cos(omega*t)*bR);
var bR_dot = -0.5 * (omega0*bI - omega1*Math.cos(omega*t)*aI);
var bI_dot = +0.5 * (omega0*bR - omega1*Math.cos(omega*t)*aR);
return [aR_dot, aI_dot, bR_dot, bI_dot]
}
// For 4th order Runge-Kutta (RK4)
var k1 = [0, 0, 0, 0];
var k2 = [0, 0, 0, 0];
var k3 = [0, 0, 0, 0];
var k4 = [0, 0, 0, 0];
var psi1 = [0,0,0,0];
var psi2 = [0,0,0,0];
var psi3 = [0,0,0,0];
var derivs = [0, 0, 0, 0];
// Takes one time step forward using RK4
function RK4_step(psi, t, dt) {
var i;
// Compute k1
derivs = f(psi, t);
for (i = 0; i<4; i++) {
k1[i] = dt*derivs[i];
}
// Compute k2
for (i = 0; i<4; i++) {
psi1[i] = psi[i] + 0.5*k1[i];
}
derivs = f(psi1, t+0.5*dt)
for (i = 0; i<4; i++) {
k2[i] = dt*derivs[i];
}
// Compute k3
for (i = 0; i<4; i++) {
psi2[i] = psi[i] + 0.5*k2[i];
}
derivs = f(psi2, t+0.5*dt)
for (i = 0; i<4; i++) {
k3[i] = dt*derivs[i];
}
// Compute k4
for (i = 0; i<4; i++) {
psi3[i] = psi[i] + k3[i];
}
derivs = f(psi3, t+dt)
for (i = 0; i<4; i++) {
k4[i] = dt*derivs[i];
}
for (i = 0; i<4; i++) {
psi[i] = psi[i] + (k1[i] + 2*k2[i] + 2*k3[i] + k4[i])/6.0;
}
return psi
}
// sx, sy, sz operators (in terms of Pauli spin matrices)
var c = 0.5
var sx = new numeric.T( [[0,c], [c,0]], [[0,0], [0,0]] )
var sy = new numeric.T( [[0,0], [0,0]], [[0,-c], [c,0]] )
var sz = new numeric.T( [[c,0], [0,-c]], [[0,0], [0,0]] )
while (true) {
rate(30,wait)
if (running) {
aR = psi[0]; aI = psi[1]; bR = psi[2]; bI = psi[3];
// Computing expectation value
var ket = new numeric.T([[aR],[bR]], [[aI],[bI]])
var bra = ket.transjugate()
var spin_x = bra.dot(sx.dot(ket)).x
var spin_y = bra.dot(sy.dot(ket)).x
var spin_z = bra.dot(sz.dot(ket)).x
var spin_vec = vec(spin_x, spin_y, spin_z)
trail.push(spin_vec)
spin.axis_and_length = spin_vec
field.axis_and_length = vec(0.5*omega1*Math.cos(omega*t),0,0.5*omega0)
psi = RK4_step(psi, t, dt)
time_label.text = 't = '+t.toFixed(1)+' ns'
t = t + dt;
}
}In order to execute GlowScript code directly within an IPython notebook, this needs to be done once to install "glowscriptmagic":
%install_ext https://raw.githubusercontent.com/jcoady/GlowScriptMagic/master/glowscriptmagic.py
Then in each new notebook, this line needs to be executed before using glowscript:
%load_ext glowscriptmagic
Then you can execute VPython code in a cell using cell magic notation %% on the first line. The first line of a GlowScript version of a VPython program in a notebook cell will contain the following line:
%%GlowScript 1.1 VPython
Or a Javascript version will contain the following line:
%%GlowScript 1.1 JavaScript
This cell magic will take the contents of the cell and perform a Python-to-Javascript compilation on it and display the results of the 3D simulation running in javascript.
The remainder of this notebooks provides Python code to compute and animate the results shown in Figs. 5 and 6 using the ivisual library.¶
Note: This requires you to be running an ipython notebook.¶
For more about ivisual, see https://pypi.python.org/pypi/IVisual/¶
Larmor precession (static field, $\omega=0$)¶
The following cell computes $\langle \hat{\vec{s}} \rangle$ at several different times for the case of a static field, $\vec{H} = H_0\mathbf{z} + H_1\mathbf{x}$; and then the next cell will create an animation showing how this vector evolves in three-dimensional space.
%matplotlib inline
from numpy import *
# Parameters defining the magnetic field:
omega = 0.0 # STATIC FIELD
omega1 = 0.1 # Magnitude of H0 (GHz)
omega0 = 1.0 # Magnitude of H0 (GHz)
parameters = (omega, omega0, omega1)
# Will be used to define time array:
tMax = 1.5 * 2*pi/sqrt(omega0**2 + omega1**2)
tArray = linspace(0, tMax, 500) # Array of time values
# Initial conditions (real & imaginary parts):
aR = 1.0; aI = 0.0 # a: Up (+)
bR = 0.0; bI = 0.0 # b: Down (-)
psi0 = (aR, aI, bR, bI) # Initial conditions
# Matrix representations of sx, sy, sz via Pauli spin matrices
sx = 0.5 * array( [[0, 1], [1, 0]] )
sy = 0.5 * array( [[0, -1j], [1j, 0]] )
sz = 0.5 * array( [[1, 0], [0, -1]] )
# Computing the quantum state vs. time
results = odeint(odeRates, psi0, tArray, args=parameters)
# These variables are used in the loop below:
bra = zeros(2, dtype=complex)
ket = zeros(2, dtype=complex)
numSteps = size(tArray)
spinVector = zeros([numSteps,3])
# Computing expectation value of quantum spin vector for each time step
for i in range(numSteps):
ket[0] = results[i,0] + 1j*results[i,1]
ket[1] = results[i,2] + 1j*results[i,3]
bra = conj(ket)
x = real(dot(bra, dot(sx, ket)))
y = real(dot(bra, dot(sy, ket)))
z = real(dot(bra, dot(sz, ket)))
spinVector[i, :] = array([x, y, z])
The next cell creates an animation showing how the $\langle \hat{\vec{s}} \rangle$ vector evolves in three-dimensional space for the case of a static field ($\omega=0$).
NOTE: The following cell creates an interactive, animated version of Figure 5. This cell can be run from the IPython notebook (.IPYNB) version of this notebook, but it cannot be run from the HTML version.
%autosave 86400
# The line above will prevent autosaving during the animation
from ivisual import *
from IPython.html import widgets
from IPython.display import display
larmor_scene = canvas(title = 'Larmor Precession') # Create scene
distant_light(direction=(0,-1,1), color=color.white)
# Create objects to be visualized
tip = sphere(pos=(0,0,0.5), radius=0.0, color=color.green,
opacity=0.8, make_trail=True, retain=numSteps)
spin = arrow(pos=(0,0,0), axis=(0,0,0.5), color=color.black, shaftwidth=0.02)
field = arrow(pos=(0, 0, 0), axis=(0.5*omega1, 0, 0.5*omega0), color=color.red, shaftwidth=0.01)
timeLabel = label(pos=(0.6,0,0))
# Configure scene
larmor_scene.background = color.white
larmor_scene.up = vector(0,0,1)
larmor_scene.forward = vector(0,1,0) # Change to (1,0,0) to view along the x axis
larmor_scene.autoscale = 1
def animation(b): # 'b' will represent a button click
for i in range(numSteps): # Loop through all time steps
rate(20) # Max of 20 frames per second
spin.axis = spinVector[i,:]
tip.pos = spinVector[i,:]
message = "Time: t = " + "%.2f ns" % tArray[i]
timeLabel.text = message
button = widgets.Button(description="Start")
display(button, larmor_scene)
button.on_click(animation)
Magnetic resonance (oscillating field, $\omega\neq0$)¶
The following cell computes $\langle \hat{\vec{s}} \rangle$ at several different times for the case of an oscillating field, $\vec{H} = H_0\mathbf{z} + H_1\cos(\omega t)\mathbf{x}$, and then the next cell will create an animation showing how this vector evolves in three-dimensional space.
%matplotlib inline
from numpy import *
# Parameters defining the magnetic field:
omega = 1.0 # OSCILLATING FIELD
omega0 = 1.0; omega1 = 0.1 # Magnitudes of H0 and H1 (GHz)
parameters = (omega, omega0, omega1)
# Will be used to define time array:
tMax = 4*pi/omega1; dt = 0.01 # (in nanoseconds)
tArray = linspace(0, tMax, 2000) # Array of time values
# Initial conditions (real & imaginary parts):
aR = 1.0; aI = 0.0 # a: Up (+)
bR = 0.0; bI = 0.0 # b: Down (-)
psi0 = (aR, aI, bR, bI) # Initial conditions
# Matrix representations of sx, sy, sz via Pauli spin matrices
sx = 0.5 * array( [[0, 1], [1, 0]] )
sy = 0.5 * array( [[0, -1j], [1j, 0]] )
sz = 0.5 * array( [[1, 0], [0, -1]] )
# Computing the quantum state vs. time
results = odeint(odeRates, psi0, tArray, args=parameters)
# These variables are used in the loop below:
bra = zeros(2, dtype=complex)
ket = zeros(2, dtype=complex)
numSteps = size(tArray)
spinVector = zeros([numSteps,3])
# Computing expectation value of quantum spin vector for each time step
for i in range(numSteps):
ket[0] = results[i,0] + 1j*results[i,1]
ket[1] = results[i,2] + 1j*results[i,3]
bra = conj(ket)
x = real(dot(bra, dot(sx, ket)))
y = real(dot(bra, dot(sy, ket)))
z = real(dot(bra, dot(sz, ket)))
spinVector[i, :] = array([x, y, z])
The next cell creates an animation showing how the $\langle \hat{\vec{s}} \rangle$ vector evolves in three-dimensional space for the case of an oscillating field ($\omega\neq 0$).
NOTE: The following cell creates an interactive, animated version of Figure 6. This cell can be run from the IPython notebook (.IPYNB) version of this notebook, but it cannot be run from the HTML version.
%autosave 86400
# The line above will prevent autosaving during the animation
from ivisual import *
from IPython.html import widgets
from IPython.display import display
titleString = 'Magnetic Resonance using parameters'+\
': omega0 = ' + str(omega0)+\
', omega1 = ' + str(omega1)+' and'+\
', omega = ' + str(omega)
resonance_scene = canvas(title = titleString) # Create scene
distant_light(direction=(0,-1,1), color=color.white)
distant_light(direction=(-1,0,1), color=color.white)
# Create objects to be visualized
tip = sphere(pos=(0,0,0.5), radius=0.0, color=color.green,
opacity=0.8, make_trail=True, retain=numSteps)
spin = arrow(pos=(0,0,0), axis=(0,0,0.5), color=color.black, shaftwidth=0.02)
field = arrow(pos=(0, 0, 0), axis=(0.5*omega1, 0, 0.5*omega0), color=color.red, shaftwidth=0.01)
#timeLabel = label(pos=(0.6,0,0))
# Configure scene
resonance_scene.background = color.white
resonance_scene.up = vector(0,0,1)
resonance_scene.forward = vector(0,1,0) # VIEWING ALONG THE Y AXIS
#resonance_scene.forward = vector(1,0,0) # VIEWING ALONG THE X AXIS
resonance_scene.autoscale = 1
def animation(b): # 'b' will represent a button click
for i in range(numSteps): # Loop through all time steps
rate(20) # Max of 20 frames per second
spin.axis = spinVector[i,:]
tip.pos = spinVector[i,:]
Hx = omega1*cos(omega*tArray[i]); Hz = omega0
field.axis = (0.5*Hx, 0, 0.5*Hz)
message = "Time: t = " + "%.2f ns" % tArray[i]
#timeLabel.text = message
start_button = widgets.Button(description="Start")
axis_button = widgets.ToggleButtons(
description='Select which axis points into the screen:',
options=['+y', '+x', '-y', '-x']
)
# Changes which axis is oriented into the screen
# and is called whenever "axis_button" is clicked:
def change_axis(name, value):
if value == '+x':
ahead = vector(1,0,0)
elif value == '+y':
ahead = vector(0,1,0)
elif value == '-x':
ahead = vector(-1,0,0)
elif value == '-y':
ahead = vector(0,-1,0)
resonance_scene.forward = ahead
display(start_button, axis_button, resonance_scene)
start_button.on_click(animation)
axis_button.on_trait_change(change_axis, 'value')