Pulse sequences for BlochSolver can be developed using the Python pulse sequence development kit (Python PSDK).

The Python PSDK manual is here.

A typical 2D RF spoiled gradient echo sequence is written by the Python PSDK as follows.

from psdk import *
import numpy as np

gamma = 42.57747892 # [MHz/T]
TR = 20.0e+3 # [us]
TE = 6.0e+3 # [us]
NR = 256 # Number of readout points
NPE1 = 256 # Number of 1st phase encoding
fov = [256.0, 256.0, 256.0] # [mm]
dwell_time = 10.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_rt = 300.0 # [us] gx rise time
gy_rt = 300.0 # [us] gy rise time
gz_rt = 300.0 # [us] gz rise time
ex_pulse_width = 3200.0 # [us]
ex_pulse_flip_angle = 30.0 # [degree]

def phase_shift_angle(i):
    phi = i * (i + 1) / 2 * 117 / 360
    phi -= round(phi)
    return 2.0 * np.pi * phi

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 Spoiled GRE'):

    with Block('Excitation', ex_pulse_width + 2.0*gz_rt):
        GZ(0.0, gz_value, gz_rt)
        RF(gz_rt, sinc_with_hamming(ex_pulse_flip_angle, ex_pulse_width, 160), ex_pulse_width / 160, 
        phase=([phase_shift_angle(i) for i in range(NPE1)], ['PE1']))
        GZ(ex_pulse_width + gz_rt, 0.0, gz_rt)
        
    with Block('PhaseEncoding', NR // 2 * dwell_time + gx_rt * 2.5):
        GX(0.0, -gx_value, gx_rt)
        GY(0.0, ([gy_value * (i - NPE1 // 2) / NPE1 for i in range(NPE1)], ['PE1']), gy_rt)
        GY(NR // 2 * dwell_time, 0.0, gy_rt)
        GX(NR // 2 * dwell_time + gx_rt * 0.5, gx_value, gx_rt * 2.0)
        GZ(0.0, -gz_value, gz_rt)
        GZ(ex_pulse_width * 0.5, 0.0, gz_rt) 

    with Block('Readout', NR * dwell_time):
        AD(0.0, NR, dwell_time, phase=([phase_shift_angle(i) for i in range(NPE1)], ['PE1']))
        
    with Block('Rewinding', NR // 2 * dwell_time + gx_rt * 2.5):
        GY(gx_rt * 2.5 - gy_rt, ([gy_value * (NPE1 // 2 - i) / NPE1 for i in range(NPE1)], ['PE1']), gy_rt)
        GX(NR // 2 * dwell_time + gx_rt * 2.5 - gx_rt, 0.0, gx_rt)
        GY(NR // 2 * dwell_time + gx_rt * 2.5 - gy_rt, 0.0, gy_rt)
        
    with Main():
        with Loop('PE1', NPE1):
            BlockRef('Excitation')
            WaitUntil(TE + ex_pulse_width * 0.5 + gz_rt - NR // 2 * 2 * dwell_time - gx_rt * 2.5)
            BlockRef('PhaseEncoding')
            BlockRef('Readout')
            BlockRef('Rewinding')
            WaitUntil(TR)

1. The first part of the source code is definitions and initialization of the parameters and variables.

2. The second part is definitions of functions such as phase_shift_angle() for RF spoiling and sinc_with_hamming() for selective excitation of the slicing plane.

3. Blocks are defined within the with syntax of the Sequence class. The Blocks are defined for RF excitation, Phase encoding, Signal readout, and Gradient rewinding.

4. The loop structure is described after the Main() entry.