This is an ipython notebook that generates an interactive plot of $p_b$ versus $t$ for magnetic resonance (Created by Larry Engelhardt, 2015)

To jump directly to the interactive plot, click here:

Note: It might take your web browser a minute to load the interactive plot, so if the link above might not work immediately upon opening this document.

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.

In [1]:
# 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':50,'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

The following cell defines a new Python function which is described immediately below the cell.

In [2]:
# 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$.

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.

In [3]:
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 as input:

  1. psi0 will contain the initial values: $a_R(0)$, $a_I(0)$, $b_R(0)$, and $b_I(0)$
  2. tArray will contain an array of many uniformly-spaced values of time, $t$.
  3. 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.

The following cell defines the parameters that will be used for the simulations:

In [4]:
tMax = 200; dt = 0.05  # (in nanoseconds)
tArray = arange(0, tMax, dt)    # 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 (-)
initialCond = (aR, aI, bR, bI)  # Initial conditions

The following cell defines a function that will create and return a plot of $p_b(t)$ for a given value of the parameters.

In [5]:
def probability_plot(omega, omega1, omega0):
    param = (omega, omega0, omega1)
    p_b = prob_down(initialCond, tArray, param)
    # The following lines create and return the plot
    fig, ax = subplots(figsize=(15, 7))
    ax.plot(tArray, p_b, linewidth=2)
    title('$\omega_{0} = $ '+str(omega0)+' GHZ,  '\
          +'$\omega = $ '+str(omega)+' GHZ,  and  '\
          + '$\omega_{1} = $ ' +str(omega1)+' GHz',
             fontsize=26)
    xlabel('Time (nanoseconds)', fontsize=26)
    ylabel('Probability, $P_b(t)$', fontsize=26)
    xlim([0, tMax])
    ylim([0, 1.0])
    grid(True)
    return fig

Note that the probability_plot function takes two things as input: $\omega$ and $\omega_1$.

The following cell provides ranges of values for $\omega$ and $\omega_1$ which will be used to create sliders in order to interactively explore the results from probability_plot.

In [6]:
from ipywidgets import StaticInteract, RangeWidget
StaticInteract(probability_plot, omega =RangeWidget(0.9, 0.99, 0.01),
                                 omega1=RangeWidget(0.02, 0.1, 0.02),
                                 omega0=RangeWidget(0.9, 1.0, 0.05))
Out[6]:
omega:
omega0:
omega1:

Interactive Figure: Drag the sliders to explore how $p_b(t)$ dependes on the values of the parameters $\omega$, $\omega_1$, and $\omega_0$.