Newer
Older
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.integrate import ode\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Define the Bloch equation\n",
"\n",
"$$\\frac{d\\vec{M}}{dt} = \\vec{M}\\times \\vec{B} - \\frac{M_x + M_y}{T2} - \\frac{M_z - M_0}{T1}$$\n",
"\n",
"In this section:\n",
"- define a function\n",
"- numpy functions like np.cross"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# define bloch equation\n",
"def bloch_ode(t, M, T1, T2):\n",
" # get effective B-field at time t\n",
" B = B_eff(t)\n",
" # cross product of M and B, add T1 and T2 relaxation terms\n",
" return np.array([M[1]*B[2] - M[2]*B[1] - M[0]/T2,\n",
" M[2]*B[0] - M[0]*B[2] - M[1]/T2,\n",
" M[0]*B[1] - M[1]*B[0] - (M[2]-1)/T1])\n",
" # alternatively\n",
" #return np.cross(M,B) - np.array([M[0]/T2, M[1]/T2, (M[2]-1)/T1])"
]
},
{
"cell_type": "markdown",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Define the pulse sequence \n",
"\n",
"We work in the rotating frame, so we only need the amplitude envelope of the RF pulse\n",
"\n",
"Typically, B1 excitation fields point in the x- and/or y-directions \n",
"Gradient fields point in the z-direction\n",
"\n",
"In this simple example, a simple sinc-pulse excitation pulse is applied for 1 ms along the x-axis \n",
"then a gradient is turned on for 1.5 ms after that.\n",
"In this section:\n",
"- constants such as np.pi\n",
"- functions like np.sinc"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# define effective B-field\n",
"def B_eff(t):\n",
" if t < 0.25:\n",
" return np.array([0, 0, 0])\n",
" # Sinc RF along x-axis and slice-select gradient on for 1.00 ms\n",
" elif t < 1.25:\n",
" return np.array([np.pi*np.sinc((t-0.75)*4), 0, np.pi])\n",
" # Do nothing for 0.25 ms\n",
" elif t < 1.50:\n",
" return np.array([0, 0, 0])\n",
" # Slice refocusing gradient on for 1.50 ms\n",
" # Half the area of the slice-select gradient lobe\n",
" elif t < 3.00:\n",
" return np.array([0, 0, -(1/3)*np.pi])\n",
" # Pulse sequence finished\n",
" else:\n",
" return np.array([0, 0, 0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Plot the pulse sequence\n",
"\n",
"In this section:\n",
"- unpacking return values\n",
"- unwanted return values\n",
"- list comprehension\n",
"- basic plotting"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create 2 vertical subplots\n",
"_, ax = plt.subplots(2, 1, figsize=(12,12))\n",
"\n",
"# Get pulse sequence B-fields from 0 - 5 ms\n",
"pulseq = [B_eff(t) for t in np.linspace(0,5,1000)]\n",
"pulseq = np.array(pulseq)\n",
"\n",
"# Plot RF\n",
"ax[0].plot(pulseq[:,0])\n",
"ax[0].set_ylabel('B1')\n",
"\n",
"# Plot gradient\n",
"ax[1].plot(pulseq[:,2])\n",
"ax[1].set_ylabel('Gradient')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Integrate ODE \n",
"\n",
"This uses a Runge-Kutta variant called the \"Dormand-Prince method\"\n",
"\n",
"In this section:\n",
"- list of arrays\n",
"- ode solvers\n",
"- list appending"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
"# Set the initial conditions\n",
"# time (t) = 0\n",
"# equilibrium magnetization (M) = (0, 0, 1)\n",
"t = [0]\n",
"M = [np.array([0, 0, 1])]\n",
"\n",
"# Set integrator time-step\n",
"dt= 0.005\n",
"\n",
"# Set up ODE integrator object\n",
"r = ode(bloch_ode)\n",
"\n",
"# Choose the integrator method\n",
"r.set_integrator('dopri5')\n",
"\n",
"# Pass in initial values\n",
"r.set_initial_value(M[0], t[0])\n",
"\n",
"# Set T1 and T2\n",
"T1, T2 = 1500, 50\n",
"r.set_f_params(T1, T2)\n",
"\n",
"# Integrate by looping over time, moving dt by step size each iteration\n",
"# Append new time point and Magnetisation vector at every step to t and M\n",
"while r.successful() and r.t < 5:\n",
" t.append(r.t + dt)\n",
" M.append(r.integrate(t[-1]))\n",
"\n",
"# Convert M to 2-D numpy array from list of arrays\n",
"M = np.array(M)"
]
},
{
"cell_type": "markdown",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot Results\n",
"\n",
"In this section:\n",
"- more plotting"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create single axis\n",
"_, ax = plt.subplots(figsize=(12,12))\n",
"\n",
"# Plot x, y and z components of Magnetisation\n",
"ax.plot(t, M[:,0], label='Mx')\n",
"ax.plot(t, M[:,1], label='My')\n",
"ax.plot(t, M[:,2], label='Mz')\n",
"\n",
"# Add legend and grid\n",
"ax.legend()\n",
"ax.grid()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
}
},
"nbformat": 4,
"nbformat_minor": 4