This is an ipython notebook that generates interactive plots of the power spectrum for magnetic resonance (Created by Larry Engelhardt, 2015)

To jump directly to the interactive plots, 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.

In the cell below, we define a function, computePowerSpectrum, that will compute and return the power spectrum, $Y(\omega^{\prime})$ (the modulus squared of the Fourier transform), as a function of angular frequency, $\omega^{\prime}$. 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).

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

tMax = 10000; dt = 0.1       # (in nanoseconds)
tArray = arange(0, tMax, dt) # Array of time values
In [6]:
omega0 = 1.0 # in GHz

def spectrum1_plot(omega, omega1, yAxisScale):
    param = (omega, omega0, omega1)
    omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
    # The following lines create and return the plot
    fig, ax = subplots(figsize=(15, 7), dpi=50)
    ax.plot(omegaPrime, Y, linewidth=2)
    if (yAxisScale == 'Linear'):
        ax.set_yscale("linear")
    else:
        ax.set_yscale("log")
        ylim([1e-1, 1e9])
    title( '$\omega_{0} = $ '+str(omega0)+' GHZ,  '\
          +'$\omega_{1} = $ '+str(omega1)+' GHz,  and  '\
          +'$\omega = $ '+str(omega)+' GHZ',\
             fontsize=26)
    xlabel('$\omega^{\prime}$ (GHz)', fontsize=26)
    ylabel('Power Spectrum (arb. units)', fontsize=26)
    xlim([0, 1.5])# ;  ylim([1e-5, 1.1])
    grid(True)
    return fig

Note that the spectrum1_plot function takes three things as input: $\omega$, $\omega_1$, and whether or not to use a logarithmic scale for the plot. The following cell provides ranges of values for $\omega$ and $\omega_1$ which will be used to create a slider in order to interactively explore the results from spectrum1_plot.

In [7]:
from ipywidgets import StaticInteract, RangeWidget, RadioWidget
StaticInteract(spectrum1_plot,
               omega =RangeWidget(0.1, 0.9, 0.1),
               omega1=RangeWidget(0.1, 0.2, 0.1),
               yAxisScale=RadioWidget(['Linear', 'Logarithmic']))
Out[7]:
omega:
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 for $\omega \neq \omega_0$: Drag the sliders to explore how the power spectrum dependes on $\omega$ and $\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. To explore what happens when $\omega=\omega_0$, use the other interactive figure that is provided below.

In [8]:
def spectrum2_plot(omega0, omega1, yAxisScale):
    omega = omega0
    param = (omega, omega0, omega1)
    omegaPrime, Y = computePowerSpectrum(initialCond, tArray, param)
    # The following lines create and return the plot
    fig, ax = subplots(figsize=(15, 7), dpi=50)
    ax.plot(omegaPrime, Y, linewidth=2)
    if (yAxisScale == 'Linear'):
        ax.set_yscale("linear")
    else:
        ax.set_yscale("log")
        ylim([1e-1, 1e9])
    title( '$\omega_{0} = $ '+str(omega0)+' GHZ,  '\
          +'$\omega_{1} = $ '+str(omega1)+' GHz,  and  '\
          +'$\omega = $ '+str(omega)+' GHZ',\
             fontsize=26)
    xlabel('$\omega^{\prime}$ (GHz)', fontsize=26)
    ylabel('Power Spectrum (arb. units)', fontsize=26)
    xlim([0, 2])
    grid(True)
    return fig

Note that the spectrum2_plot function takes three things as input: $\omega_0$, $\omega_1$, and whether or not to use a logarithmic scale for the plot. Within the function, $\omega$ is set equal to $\omega_0$. The following cell provides ranges of values for $\omega_0$ and $\omega_1$ which will be used to create a slider in order to interactively explore the results from spectrum2_plot.

In [9]:
from ipywidgets import StaticInteract, RangeWidget, RadioWidget
StaticInteract(spectrum2_plot,
               omega0=RangeWidget(0.7, 0.9, 0.1),
               omega1=RangeWidget(0.05, 0.5, 0.05),
               yAxisScale=RadioWidget(['Linear', 'Logarithmic']))
Out[9]:
omega0:
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 for $\omega = \omega_0$: Drag the sliders to explore how the power spectrum dependes on $\omega_0$ and $\omega_1$ for $\omega=\omega_0$. You can also use the radio button to switch to a logarithmic scale for the vertical axis. To explore what happens when $\omega \neq \omega_0$, use the other interactive figure that is provided above.