2D multiple spin echo imaging

Image matrix = 256 x 256, slice thickness = 5 mm, TR = 4000 ms, Echo spacing = 12 ms, Echo train length = 16, subvoxel = 1 x 1 x 4, number of isochromats = 13,572,096, calculation time = 435.372 s

Python sequence code:

from psdk import *
import numpy as np

gamma = 42.57747892 # [MHz/T]
TR = 4000.0e+3 # [us]
TE = 12.0e+3 # [us]
NR = 256 # Number of readout points
NPE1 = 256 # Number of 1st phase encoding
NEcho = 16 # Number of echoes
fov = [256.0, 256.0, 256.0] # [mm]
dwell_time = 5.0 # [us]
slice_width = 5.0 # [mm]
gx_value = 1e+6 / (dwell_time * gamma * fov[0]) # [mT/m]
gy_value = 2e+6 / (dwell_time * gamma * fov[1]) * NPE1 / NR # [mT/m]
gz_value = 1.25 / (slice_width * 1.0e-3) / gamma # [mT/m]
gx_rise_time = 300.0 # [us]
gy_rise_time = 300.0 # [us]
gz_rise_time = 300.0 # [us]
ex_pulse_width = 1600.0 # [us]
excitation_pulse_flip_angle = 90.0 # [degree]

def sinc_with_hamming(flip_angle, pulse_width, points, *, min=-2.0*np.pi, max=2.0*np.pi):
    x0 = np.arange(min, max, (max - min) / points)
    x1 = x0 + (max - min) / points
    y = (np.sinc(x0 / np.pi) + np.sinc(x1 / np.pi)) * 0.5 * np.hamming(points)
    return flip_angle * y * points / (y.sum() * pulse_width * 360.0e-6 * gamma)

with Sequence("2D multiple SpinEcho"):

    with Block("Excitation", ex_pulse_width + 2.0*gz_rise_time):
        GZ(0.0, gz_value, gz_rise_time)
        RF(gz_rise_time, sinc_with_hamming(excitation_pulse_flip_angle, ex_pulse_width, 160), ex_pulse_width / 160, phase = 0.0)
        GZ(ex_pulse_width + gz_rise_time, 0.0, gz_rise_time)

    with Block("Slice_refocus", ex_pulse_width * 0.5 + gz_rise_time * 2.0) :
        GZ(0.0, -gz_value, gz_rise_time)
        GZ(ex_pulse_width * 0.5 + gz_rise_time, 0.0, gz_rise_time) 

    with Block("ReadPredephasing", NR // 2 * dwell_time * 2.0 + gx_rise_time ):
        GX(0.0, gx_value, gx_rise_time)
        GX(NR // 2 * dwell_time * 2.0, 0.0, gx_rise_time)         

    with Block("Refocus", ex_pulse_width + 2.0*gz_rise_time):
        GZ(0.0, gz_value, gz_rise_time)
        RF(gz_rise_time, 2.0 * sinc_with_hamming(excitation_pulse_flip_angle, ex_pulse_width, 160), ex_pulse_width / 160, phase=np.pi * 0.5)
        GZ(ex_pulse_width + gz_rise_time, 0.0, gz_rise_time)        

    with Block("PhaseEncoding", NR // 2 * dwell_time + gy_rise_time * 2.5):
        GY(0.0, ([gy_value * (i - NPE1 // 2) / NPE1 for i in range(NPE1)], ['PE1']), gy_rise_time)
        GY(NR // 2 * dwell_time, 0.0, gy_rise_time)

    with Block("Readout", NR * dwell_time * 2.0 + gx_rise_time):
        GX(0.0, gx_value, gx_rise_time)
        AD(NR // 2 * dwell_time + gx_rise_time * 0.5, NR, dwell_time)
        GX(NR * dwell_time * 2.0, 0, gx_rise_time)
        
    with Block("Rewinding", NR // 2 * dwell_time + gy_rise_time * 2.5):
        GY(gx_rise_time * 2.5 - gy_rise_time, ([gy_value * (NPE1 // 2 - i) / NPE1 for i in range(NPE1)], ['PE1']), gy_rise_time)
        GY(NR // 2 * dwell_time + gy_rise_time * 2.5 - gy_rise_time, 0.0, gy_rise_time)
        
    with Main():
        with Loop("PE1", NPE1):
            BlockRef("Excitation")
            BlockRef("Slice_refocus")
            BlockRef("ReadPredephasing")
            WaitUntil(TE * 0.5)
            with Loop("Echo", NEcho):
                BlockRef("Refocus")
                BlockRef("PhaseEncoding")
                WaitUntil(TE * 0.5 + gz_rise_time + ex_pulse_width * 0.5 - NR * dwell_time - gx_rise_time) 
                BlockRef("Readout")
                BlockRef("Rewinding")
                WaitUntil(TE * 1.0)                
            WaitUntil(TR)