Commit f275165a authored by Sean Fitzgibbon's avatar Sean Fitzgibbon
Browse files

Formatting

parent c216fda1
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"Imports" "# Imports"
] ]
}, },
{ {
...@@ -216,7 +216,7 @@ ...@@ -216,7 +216,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython3", "pygments_lexer": "ipython3",
"version": "3.7.3" "version": "3.7.6"
} }
}, },
"nbformat": 4, "nbformat": 4,
......
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Imports # Imports
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import numpy as np import numpy as np
from scipy.integrate import ode, solve_ivp from scipy.integrate import ode, solve_ivp
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Define the Bloch equation # Define the Bloch equation
$$\frac{d\vec{M}}{dt} = \vec{M}\times \vec{B} - \frac{M_x + M_y}{T2} - \frac{M_z - M_0}{T1}$$ $$\frac{d\vec{M}}{dt} = \vec{M}\times \vec{B} - \frac{M_x + M_y}{T2} - \frac{M_z - M_0}{T1}$$
In this section: In this section:
- define a function - define a function
- numpy functions like np.cross - numpy functions like np.cross
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# define bloch equation # define bloch equation
def bloch(t, M, T1, T2): def bloch(t, M, T1, T2):
# get effective B-field at time t # get effective B-field at time t
B = B_eff(t) B = B_eff(t)
# cross product of M and B, add T1 and T2 relaxation terms # cross product of M and B, add T1 and T2 relaxation terms
return np.array([M[1]*B[2] - M[2]*B[1] - M[0]/T2, return np.array([M[1]*B[2] - M[2]*B[1] - M[0]/T2,
M[2]*B[0] - M[0]*B[2] - M[1]/T2, M[2]*B[0] - M[0]*B[2] - M[1]/T2,
M[0]*B[1] - M[1]*B[0] - (M[2]-1)/T1]) M[0]*B[1] - M[1]*B[0] - (M[2]-1)/T1])
# alternatively # alternatively
#return np.cross(M,B) - np.array([M[0]/T2, M[1]/T2, (M[2]-1)/T1]) #return np.cross(M,B) - np.array([M[0]/T2, M[1]/T2, (M[2]-1)/T1])
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Define the pulse sequence # Define the pulse sequence
We work in the rotating frame, so we only need the amplitude envelope of the RF pulse We work in the rotating frame, so we only need the amplitude envelope of the RF pulse
Typically, B1 excitation fields point in the x- and/or y-directions Typically, B1 excitation fields point in the x- and/or y-directions
Gradient fields point in the z-direction Gradient fields point in the z-direction
In this simple example, a simple sinc-pulse excitation pulse is applied for 1 ms along the x-axis In this simple example, a simple sinc-pulse excitation pulse is applied for 1 ms along the x-axis
then a gradient is turned on for 1.5 ms after that. then a gradient is turned on for 1.5 ms after that.
In this section: In this section:
- constants such as np.pi - constants such as np.pi
- functions like np.sinc - functions like np.sinc
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# define effective B-field # define effective B-field
def B_eff(t): def B_eff(t):
# Do nothing for 0.25 ms # Do nothing for 0.25 ms
if t < 0.25: if t < 0.25:
return np.array([0, 0, 0]) return np.array([0, 0, 0])
# Sinc RF along x-axis and slice-select gradient on for 1.00 ms # Sinc RF along x-axis and slice-select gradient on for 1.00 ms
elif t < 1.25: elif t < 1.25:
return np.array([np.pi*np.sinc((t-0.75)*4), 0, np.pi]) return np.array([np.pi*np.sinc((t-0.75)*4), 0, np.pi])
# Do nothing for 0.25 ms # Do nothing for 0.25 ms
elif t < 1.50: elif t < 1.50:
return np.array([0, 0, 0]) return np.array([0, 0, 0])
# Slice refocusing gradient on for 1.50 ms # Slice refocusing gradient on for 1.50 ms
# Half the area of the slice-select gradient lobe # Half the area of the slice-select gradient lobe
elif t < 3.00: elif t < 3.00:
return np.array([0, 0, -(1/3)*np.pi]) return np.array([0, 0, -(1/3)*np.pi])
# Pulse sequence finished # Pulse sequence finished
else: else:
return np.array([0, 0, 0]) return np.array([0, 0, 0])
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Plot the pulse sequence # Plot the pulse sequence
In this section: In this section:
- unpacking return values - unpacking return values
- unwanted return values - unwanted return values
- list comprehension - list comprehension
- basic plotting - basic plotting
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Create 2 vertical subplots # Create 2 vertical subplots
_, ax = plt.subplots(2, 1, figsize=(12,12)) _, ax = plt.subplots(2, 1, figsize=(12,12))
# Get pulse sequence B-fields from 0 - 5 ms # Get pulse sequence B-fields from 0 - 5 ms
pulseq = [B_eff(t) for t in np.linspace(0,5,1000)] pulseq = [B_eff(t) for t in np.linspace(0,5,1000)]
pulseq = np.array(pulseq) pulseq = np.array(pulseq)
# Plot RF # Plot RF
ax[0].plot(pulseq[:,0]) ax[0].plot(pulseq[:,0])
ax[0].set_ylabel('B1') ax[0].set_ylabel('B1')
# Plot gradient # Plot gradient
ax[1].plot(pulseq[:,2]) ax[1].plot(pulseq[:,2])
ax[1].set_ylabel('Gradient') ax[1].set_ylabel('Gradient')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Integrate ODE # Integrate ODE
This uses a Runge-Kutta variant called the "Dormand-Prince method" by default. Other ode integration methods are available. This uses a Runge-Kutta variant called the "Dormand-Prince method" by default. Other ode integration methods are available.
In this section: In this section:
- ode solvers - ode solvers
- lambdas (anonymous functions) - lambdas (anonymous functions)
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Set the initial conditions # Set the initial conditions
# equilibrium magnetization (M) = (0, 0, 1) # equilibrium magnetization (M) = (0, 0, 1)
M = [0, 0, 1] M = [0, 0, 1]
# Set time interval for integration # Set time interval for integration
t = [0, 5] t = [0, 5]
# Set max step size # Set max step size
dt = 0.005 dt = 0.005
# Set T1 and T2 # Set T1 and T2
T1, T2 = 1500, 50 T1, T2 = 1500, 50
# Integrate ODE # Integrate ODE
# In Scipy 1.2.0, the first argument to solve_ivp must be a function that takes exactly 2 arguments # In Scipy 1.2.0, the first argument to solve_ivp must be a function that takes exactly 2 arguments
sol = solve_ivp(lambda t, M : bloch(t, M, T1, T2), t, M, max_step=dt) sol = solve_ivp(lambda t, M : bloch(t, M, T1, T2), t, M, max_step=dt)
# Grab output # Grab output
t = sol.t t = sol.t
M = sol.y M = sol.y
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Plot Results # Plot Results
In this section: In this section:
- more plotting - more plotting
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Create single axis # Create single axis
_, ax = plt.subplots(figsize=(12,12)) _, ax = plt.subplots(figsize=(12,12))
# Plot x, y and z components of Magnetisation # Plot x, y and z components of Magnetisation
ax.plot(t, M[0,:], label='Mx') ax.plot(t, M[0,:], label='Mx')
ax.plot(t, M[1,:], label='My') ax.plot(t, M[1,:], label='My')
ax.plot(t, M[2,:], label='Mz') ax.plot(t, M[2,:], label='Mz')
# Add legend and grid # Add legend and grid
ax.legend() ax.legend()
ax.grid() ax.grid()
``` ```
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment