This is an ipython notebook that generates an interactive plot of flipping probability versus $\omega$ 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.

Below 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)$.

In [4]:
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

Next, we make an interactive version of Figure 2 by first defining a function that will return a figure and then providing a range of values for $\omega_1$.

In [5]:
def flipping_plot(log10_omega1, yAxisScale):
    omega1 = 10**log10_omega1 # Converts from log-scale of slider
    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)
    fig, ax = subplots(figsize=(15, 7))
    ax.plot(omegaArray, pArray, '-o', linewidth=2)
    if (yAxisScale == 'Linear'):
        ax.set_yscale("linear")
    else:
        ax.set_yscale("log")
    omega1text = '%.3f' % omega1 # SHOW 3 DIGITS AFTER THE DECIMAL
    title('$\omega_{1} = $ ' + omega1text +' GHz', fontsize=26)
    xlabel('$\omega$ (GHz)', fontsize=26)
    ylabel('Flipping Probability', fontsize=26)
    xlim([0.9, 1.1]);  ylim([1e-5, 1.1])
    grid(True)
    return fig

Note that the flipping_plot function takes two things as input: $\log_{10}(\omega_1)$ and whether or not to use a logarithmic scale for the plot. The following cell provides a range of values for $\log_{10}(\omega_1)$ which will be used to create a slider in order to interactively explore the results from flipping_plot.

In [6]:
from ipywidgets import StaticInteract, RangeWidget, RadioWidget
StaticInteract(flipping_plot,
               log10_omega1=RangeWidget(-3, -1, 0.5),
               yAxisScale=RadioWidget(['Linear', 'Logarithmic']))
Out[6]:
log10_omega1:
yAxisScale: Linear: Logarithmic:

NOTE: Make sure that one of the options (either Linear or Logarithmic) is selected from the radio buttons above.

Interactive Figure: Drag the sliders to explore how the flipping probability dependes on $\omega_1$ for a fixed value of $\omega_0=1.0$ GHz. You can also use the radio button to switch to a logarithmic scale for the vertical axis.