From 102ea6fce9108b378a64395c0d7d3945490c7d7b Mon Sep 17 00:00:00 2001
From: Christoph Kolbitsch <christoph.kolbitsch@ptb.de>
Date: Tue, 24 Feb 2026 11:45:02 +0000
Subject: [PATCH 1/6] sat pulse added

---
 src/mrseq/preparations/adiabatic_sat_prep.py  | 161 ++++++++++++++++++
 .../{b0_b1_wasabi.py => b0_b1_wasabiti.py}    |  30 +++-
 src/mrseq/utils/__init__.py                   |   1 +
 src/mrseq/utils/hyperbolic_secant_pulse.py    | 116 +++++++++++++
 ...b0_b1_wasabi.py => test_b0_b1_wasabiti.py} |   2 +-
 5 files changed, 300 insertions(+), 10 deletions(-)
 create mode 100644 src/mrseq/preparations/adiabatic_sat_prep.py
 rename src/mrseq/scripts/{b0_b1_wasabi.py => b0_b1_wasabiti.py} (93%)
 create mode 100644 src/mrseq/utils/hyperbolic_secant_pulse.py
 rename tests/scripts/{test_b0_b1_wasabi.py => test_b0_b1_wasabiti.py} (85%)

diff --git a/src/mrseq/preparations/adiabatic_sat_prep.py b/src/mrseq/preparations/adiabatic_sat_prep.py
new file mode 100644
index 0000000..8683b14
--- /dev/null
+++ b/src/mrseq/preparations/adiabatic_sat_prep.py
@@ -0,0 +1,161 @@
+"""Generate composite or adiabatic saturation pulse trains."""
+
+import numpy as np
+from pypulseq import Opts
+from pypulseq import Sequence
+from pypulseq import make_sinc_pulse
+from pypulseq import make_trapezoid
+
+from mrseq.utils.hyperbolic_secant_pulse import make_hypsec_90
+
+
+def add_adia_sat_block(
+    seq: Sequence | None = None,
+    n_pulses: int = 3,
+    max_b1: float = 20,
+    sys: Opts | None = None,
+) -> tuple[Sequence, float | None]:
+    """Add adiabatic saturation pulse train to a PyPulseq Sequence object.
+
+    Parameters
+    ----------
+    seq
+        PyPulseq Sequence object
+    n_pulses
+        number of rf pulses in pulse train
+    max_b1
+        maximum B1 field amplitude in µT
+    sys
+        system limits
+
+    Returns
+    -------
+    seq
+        PyPulseq Sequence object
+    last_spoil_dur
+        duration of the last spoiler gradient in s
+    """
+    if not seq:
+        seq = Sequence()
+
+    if not sys:
+        sys = Opts(max_grad=24, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s')
+
+    # spoilers
+    spoil_amp0 = 0.8 * sys.max_grad  # Hz/m
+    spoil_amp1 = -0.7 * sys.max_grad  # Hz/m
+    spoil_amp2 = 0.6 * sys.max_grad  # Hz/m
+
+    rise_time = 1.0e-3  # spoiler rise time in seconds
+    spoil_dur = 5.5e-3  # complete spoiler duration in seconds
+
+    gx_spoil0, gy_spoil0, gz_spoil0 = (
+        make_trapezoid(channel=c, system=sys, amplitude=spoil_amp0, duration=spoil_dur, rise_time=rise_time)
+        for c in ['x', 'y', 'z']
+    )
+    gx_spoil1, gy_spoil1, gz_spoil1 = (
+        make_trapezoid(channel=c, system=sys, amplitude=spoil_amp1, duration=spoil_dur, rise_time=rise_time)
+        for c in ['x', 'y', 'z']
+    )
+    gx_spoil2, gy_spoil2, gz_spoil2 = (
+        make_trapezoid(channel=c, system=sys, amplitude=spoil_amp2, duration=spoil_dur, rise_time=rise_time)
+        for c in ['x', 'y', 'z']
+    )
+
+    hypsec_90 = make_hypsec_90(amp=max_b1, system=sys)
+
+    for i in range(n_pulses):
+        seq.add_block(hypsec_90)
+        if i % 3 == 0:
+            seq.add_block(gx_spoil0, gy_spoil1, gz_spoil2)
+        elif i % 2 == 0:
+            seq.add_block(gx_spoil2, gy_spoil0, gz_spoil1)
+        else:
+            seq.add_block(gx_spoil1, gy_spoil2, gz_spoil0)
+
+    return seq, spoil_dur
+
+
+def add_sat_pulse_train(
+    seq: Sequence = None,
+    n_pulses: int = 6,
+    max_b1: float = 14,
+    sys: Opts = None,
+    grad_strength_factor: float = 1.0,
+) -> tuple[Sequence, float | None]:
+    """Create a saturation pulse train according to Chow, K. et al.
+
+    For more information see:
+    Saturation pulse design for quantitative myocardial T1 mapping. J. Cardiovasc. Magn. Reason. 17, 1-15 (2015).
+
+    Parameters
+    ----------
+    seq
+        PyPulseq Sequence object
+    n_pulses
+        number of rf pulses in pulse train
+    max_b1
+        maximum B1 field amplitude in µT
+    sys
+        system limits
+    grad_strength_factor
+        increase the spoiler gradient strength by this factor
+
+    Returns
+    -------
+    seq
+        PyPulseq Sequence object
+    last_spoil_dur
+        duration of the last spoiler gradient in s
+    """
+    if not seq:
+        seq = Sequence()
+
+    if not sys:
+        sys = Opts(max_grad=24, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s')
+
+    flip_angles = [115, 90, 125, 85, 176, 223]
+    grad_dur = [5, 5.25, 4.5, 3.75, 3.5, 3, 4]
+    grad_area_values = [153, 198, 136, 111, 72, 60, 119]
+    grad_area = np.array(grad_area_values) * grad_strength_factor
+    grad_ro = [-1, 1, -1, 0, 0, 0, 1]
+    grad_pe = [0, 1, -1, 1, -1, 0, 0]
+    grad_ss = [-1, 1, 0, -1, 0, 1, -1]
+
+    dur_block = np.array(flip_angles) / (360 * sys.gamma * 10e-7 * max_b1)
+    dur_sinc = np.ceil(dur_block / 0.22570566672775255 * 1e5) * 1e-5  # factor equals mean of sinc shape
+    dur_block = np.ceil(dur_block * 1e5) * 1e-5
+
+    for i in range(n_pulses + 1):
+        if grad_ro[i] != 0:
+            gx = make_trapezoid(
+                channel='x', area=grad_ro[i] * grad_area[i], duration=grad_dur[i] / 1000, rise_time=5e-4
+            )
+        if grad_pe[i] != 0:
+            gy = make_trapezoid(
+                channel='y', area=grad_pe[i] * grad_area[i], duration=grad_dur[i] / 1000, rise_time=5e-4
+            )
+        if grad_ss[i] != 0:
+            gz = make_trapezoid(
+                channel='z', area=grad_ss[i] * grad_area[i], duration=grad_dur[i] / 1000, rise_time=5e-4
+            )
+
+        if grad_ro[i] != 0 and grad_pe[i] != 0 and grad_ss[i] != 0:
+            seq.add_block(gx, gy, gz)
+        elif grad_ro[i] != 0 and grad_pe[i] != 0:
+            seq.add_block(gx, gy)
+        elif grad_ro[i] != 0 and grad_ss[i] != 0:
+            seq.add_block(gx, gz)
+        elif grad_ro[i] != 0:
+            seq.add_block(gx)
+        elif grad_pe[i] != 0 and grad_ss[i] != 0:
+            seq.add_block(gy, gz)
+        elif grad_pe[i] != 0:
+            seq.add_block(gy)
+        elif grad_ss[i] != 0:
+            seq.add_block(gz)
+
+        if i < n_pulses:
+            seq.add_block(make_sinc_pulse(flip_angle=flip_angles[i] * np.pi / 180, duration=dur_sinc[i], system=sys))
+
+    return seq, grad_dur[-1] * 1e-3
diff --git a/src/mrseq/scripts/b0_b1_wasabi.py b/src/mrseq/scripts/b0_b1_wasabiti.py
similarity index 93%
rename from src/mrseq/scripts/b0_b1_wasabi.py
rename to src/mrseq/scripts/b0_b1_wasabiti.py
index 6a8fd0f..5100730 100644
--- a/src/mrseq/scripts/b0_b1_wasabi.py
+++ b/src/mrseq/scripts/b0_b1_wasabiti.py
@@ -5,6 +5,7 @@
 import numpy as np
 import pypulseq as pp
 
+from mrseq.preparations.adiabatic_sat_prep import add_adia_sat_block
 from mrseq.utils import find_gx_flat_time_on_adc_raster
 from mrseq.utils import round_to_raster
 from mrseq.utils import sys_defaults
@@ -13,11 +14,11 @@
 from mrseq.utils.trajectory import cartesian_phase_encoding
 
 
-def wasabi_gre_centric_kernel(
+def wasabiti_gre_centric_kernel(
     system: pp.Opts,
     frequency_offsets: np.ndarray,
     norm_offset: float | None,
-    t_recovery: float,
+    t_recovery: np.ndarray,
     t_recovery_norm: float,
     fov_xy: float,
     n_readout: int,
@@ -85,9 +86,15 @@ def wasabi_gre_centric_kernel(
     if readout_oversampling < 1:
         raise ValueError('Readout oversampling factor must be >= 1.')
 
+    if len(t_recovery) != len(frequency_offsets):
+        raise ValueError(
+            f'Number of recovery times ({len(t_recovery)}) and offsets ({len(frequency_offsets)}) have to match.'
+        )
+
     # add normalization offset to beginning of frequency offsets if specified
     if norm_offset is not None:
         frequency_offsets = np.concatenate(([norm_offset], frequency_offsets))
+        t_recovery = np.concatenate(([t_recovery_norm], t_recovery))
 
     # create WASABI block pulse
     rf_prep_flipangle_rad = rf_prep_amplitude * 1e-6 * rf_prep_duration * GYROMAGNETIC_RATIO_PROTON * 2 * np.pi
@@ -167,11 +174,13 @@ def wasabi_gre_centric_kernel(
         # set repetition ('REP') label for current frequency offset
         rep_label = pp.make_label(type='SET', label='REP', value=int(rep_idx))
 
-        # add delay for normalization offset
-        if rep_idx == 0 and norm_offset is not None:
-            seq.add_block(pp.make_delay(t_recovery_norm))
-        else:
-            seq.add_block(pp.make_delay(t_recovery))
+        # add adiabatic saturation pulse train and recovery time
+        seq, last_spoil_dur = add_adia_sat_block(seq=seq, sys=system)
+        seq.add_block(
+            pp.make_delay(
+                round_to_raster(t_recovery[rep_idx] - last_spoil_dur, raster_time=system.block_duration_raster)
+            )
+        )
 
         # update frequency offset of WASABI block pulse and add it to sequence
         rf_prep.freq_offset = freq_offset_hz
@@ -233,7 +242,7 @@ def main(
     system: pp.Opts | None = None,
     frequency_offsets: np.ndarray | None = None,
     norm_offset: float | None = -35e3,
-    t_recovery: float = 3.0,
+    t_recovery: float | np.ndarray = 3.0,
     t_recovery_norm: float = 12.0,
     fov_xy: float = 200e-3,
     n_readout: int = 128,
@@ -288,6 +297,9 @@ def main(
     if frequency_offsets is None:
         frequency_offsets = np.linspace(-240, 240, 31)
 
+    if not isinstance(t_recovery, np.ndarray):
+        t_recovery = np.ones_like(frequency_offsets) * t_recovery
+
     # define settings of rf excitation pulse
     rf_duration = 1.28e-3  # duration of the rf excitation pulse [s]
     rf_flip_angle = 10.0  # flip angle of rf excitation pulse [°]
@@ -309,7 +321,7 @@ def main(
     rf_prep_duration: float = 5e-3
     rf_prep_amplitude: float = 3.75
 
-    seq = wasabi_gre_centric_kernel(
+    seq = wasabiti_gre_centric_kernel(
         system=system,
         frequency_offsets=frequency_offsets,
         norm_offset=norm_offset,
diff --git a/src/mrseq/utils/__init__.py b/src/mrseq/utils/__init__.py
index 00d413e..2f2b4fd 100644
--- a/src/mrseq/utils/__init__.py
+++ b/src/mrseq/utils/__init__.py
@@ -3,3 +3,4 @@
 from mrseq.utils.vds import variable_density_spiral_trajectory
 from mrseq.utils.trajectory import spiral_acquisition
 from mrseq.utils.ismrmrd import create_header, combine_ismrmrd_files
+from mrseq.utils.hyperbolic_secant_pulse import make_hypsec_180, make_hypsec_90
diff --git a/src/mrseq/utils/hyperbolic_secant_pulse.py b/src/mrseq/utils/hyperbolic_secant_pulse.py
new file mode 100644
index 0000000..14d39fb
--- /dev/null
+++ b/src/mrseq/utils/hyperbolic_secant_pulse.py
@@ -0,0 +1,116 @@
+"""Functions to create an adiabatic hyperbolic secant full passage pulse."""
+
+from types import SimpleNamespace
+
+import numpy as np
+from pypulseq import Opts
+
+
+# Todo: import from BMCTool instead ?!
+def calculate_phase(
+    frequency: np.ndarray, duration: float, samples: int, shift_idx: int = -1, pos_offsets: bool = False
+) -> np.ndarray:
+    """Calculate phase modulation for a given frequency modulation."""
+    phase = np.zeros_like(frequency)
+    for i in range(1, samples):
+        phase[i] = phase[i - 1] + (frequency[i] * duration / samples)
+    phase_shift = phase[shift_idx]
+    for i in range(samples):
+        phase[i] = np.fmod(phase[i] + 1e-12 - phase_shift, 2 * np.pi)
+    if not pos_offsets:
+        phase += 2 * np.pi
+    return phase
+
+
+# Todo: import from BMCTool instead ?!
+def create_arbitrary_pulse_with_phase(
+    signal: np.ndarray,
+    flip_angle: float,
+    freq_offset: float = 0,
+    phase_offset: float = 0,
+    system: Opts | None = None,
+) -> SimpleNamespace:
+    """Create a RF pulse with arbitrary shape and phase."""
+    if not system:
+        system = Opts()
+
+    signal *= flip_angle / (2 * np.pi)
+    t = np.arange(1, len(signal) + 1) * system.rf_raster_time
+
+    rf = SimpleNamespace()
+    rf.type = 'rf'
+    rf.signal = signal
+    rf.t = t
+    rf.shape_dur = t[-1]
+    rf.freq_offset = freq_offset
+    rf.phase_offset = phase_offset
+    rf.dead_time = system.rf_dead_time
+    rf.ringdown_time = system.rf_ringdown_time
+    rf.delay = system.rf_dead_time
+    """If rf.ringdown_time > 0:
+
+    t_fill = np.arange(1, round(rf.ringdown_time /
+    system.rf_raster_time) + 1) * system.rf_raster_time rf.t =
+    np.concatenate((rf.t, rf.t[-1] + t_fill)) rf.signal =
+    np.concatenate((rf.signal, np.zeros(len(t_fill))))
+    """
+
+    return rf
+
+
+def calculate_amplitude(t: np.ndarray, amp: float, mu: float, bandwidth: float) -> np.ndarray:
+    """Calculate amplitude modulation for a hypsec inversion pulse."""
+    return np.divide(amp, np.cosh((bandwidth * np.pi / mu) * t))
+
+
+def calculate_frequency(t: np.ndarray, mu: float, bandwidth: float) -> np.ndarray:
+    """Calculate frequency modulation for a hypsec inversion pulse."""
+    beta = bandwidth * np.pi / mu
+    return bandwidth * np.pi * np.tanh(beta * t)
+
+
+def make_hypsec_180(
+    amp: float,
+    pulse_duration: float = 8e-3,
+    mu: float = 6,
+    bandwidth: float = 1200,
+    system: Opts | None = None,
+) -> SimpleNamespace:
+    """Create a hypsec full passage pulse."""
+    if not system:
+        system = Opts()
+
+    samples = int(pulse_duration / system.rf_raster_time)
+    t_pulse = np.arange(-samples // 2, samples // 2) / samples * pulse_duration
+    w1 = calculate_amplitude(t=t_pulse, amp=1, mu=mu, bandwidth=bandwidth)
+    freq = calculate_frequency(t=t_pulse, mu=mu, bandwidth=bandwidth)
+    phase = calculate_phase(frequency=freq, duration=pulse_duration, samples=samples)
+    signal = np.multiply(w1, np.exp(1j * phase))
+    flip_angle = amp * 1e-6 * system.gamma * 2 * np.pi  # factor 1e-6 converts from µT to T
+    hypsec_180 = create_arbitrary_pulse_with_phase(signal=signal, flip_angle=flip_angle, system=system)
+    return hypsec_180
+
+
+def make_hypsec_90(
+    amp: float,
+    pulse_duration: float = 8e-3,
+    mu: float = 6,
+    bandwidth: float = 1200,
+    system: Opts | None = None,
+) -> SimpleNamespace:
+    """Create a hypsec half passage pulse."""
+    if not system:
+        system = Opts()
+
+    samples = int(pulse_duration / system.rf_raster_time)
+    t_pulse = np.divide(np.arange(1, samples + 1), samples) * pulse_duration
+    t_pulse -= t_pulse[-1]
+    w1 = calculate_amplitude(t=t_pulse, amp=1, mu=mu, bandwidth=bandwidth)
+    freq = calculate_frequency(t=t_pulse, mu=mu, bandwidth=bandwidth)
+    freq = freq - freq[-1]  # ensure phase ends with 0 for tip-down pulse
+    phase = calculate_phase(frequency=freq, duration=pulse_duration, samples=samples)
+    signal = np.multiply(w1, np.exp(1j * phase))
+    flip_angle = amp * 1e-6 * system.gamma * 2 * np.pi  # factor 1e-6 converts from µT to T
+    hypsec_90 = create_arbitrary_pulse_with_phase(signal=signal, flip_angle=flip_angle, system=system)
+    # hypsec_90 = make_arbitrary_rf(signal=signal, flip_angle=flip_angle, system=system)
+    return hypsec_90
diff --git a/tests/scripts/test_b0_b1_wasabi.py b/tests/scripts/test_b0_b1_wasabiti.py
similarity index 85%
rename from tests/scripts/test_b0_b1_wasabi.py
rename to tests/scripts/test_b0_b1_wasabiti.py
index 78efe09..d6047fb 100644
--- a/tests/scripts/test_b0_b1_wasabi.py
+++ b/tests/scripts/test_b0_b1_wasabiti.py
@@ -1,7 +1,7 @@
 """Tests for WASABI sequence."""
 
 import pytest
-from mrseq.scripts.b0_b1_wasabi import main as create_seq
+from mrseq.scripts.b0_b1_wasabiti import main as create_seq
 
 EXPECTED_DUR = 122.1104  # defined 2026-02-10
 

From cf2220aebad947385d9e8f3684b6ced02d41910d Mon Sep 17 00:00:00 2001
From: Christoph Kolbitsch <christoph.kolbitsch@ptb.de>
Date: Tue, 24 Feb 2026 11:47:21 +0000
Subject: [PATCH 2/6] test fixed

---
 tests/scripts/test_b0_b1_wasabiti.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/tests/scripts/test_b0_b1_wasabiti.py b/tests/scripts/test_b0_b1_wasabiti.py
index d6047fb..5cba3ac 100644
--- a/tests/scripts/test_b0_b1_wasabiti.py
+++ b/tests/scripts/test_b0_b1_wasabiti.py
@@ -1,9 +1,10 @@
 """Tests for WASABI sequence."""
 
+import numpy as np
 import pytest
 from mrseq.scripts.b0_b1_wasabiti import main as create_seq
 
-EXPECTED_DUR = 122.1104  # defined 2026-02-10
+EXPECTED_DUR = 123.2429  # defined 2026-02-24
 
 
 def test_default_seq_duration(system_defaults):
@@ -11,3 +12,9 @@ def test_default_seq_duration(system_defaults):
     seq, _ = create_seq(system=system_defaults, show_plots=False)
     duration = seq.duration()[0]
     assert duration == pytest.approx(EXPECTED_DUR)
+
+
+def test_seq_error_mismatch_t_offset(system_defaults):
+    """Test if error is raised on mismatch between recovery times and offsets."""
+    with pytest.raises(ValueError):
+        create_seq(system=system_defaults, frequency_offsets=np.ones(3), t_recovery=np.ones(2), show_plots=False)

From cfd09a2266c8e312c1da92ee4f65e23340b5c53a Mon Sep 17 00:00:00 2001
From: Christoph Kolbitsch <christoph.kolbitsch@ptb.de>
Date: Tue, 24 Feb 2026 11:59:30 +0000
Subject: [PATCH 3/6] tidy

---
 src/mrseq/preparations/adiabatic_sat_prep.py | 75 +++++++++++---------
 src/mrseq/scripts/b0_b1_wasabiti.py          | 19 +++--
 src/mrseq/utils/hyperbolic_secant_pulse.py   |  4 +-
 3 files changed, 54 insertions(+), 44 deletions(-)

diff --git a/src/mrseq/preparations/adiabatic_sat_prep.py b/src/mrseq/preparations/adiabatic_sat_prep.py
index 8683b14..5bf1146 100644
--- a/src/mrseq/preparations/adiabatic_sat_prep.py
+++ b/src/mrseq/preparations/adiabatic_sat_prep.py
@@ -1,32 +1,31 @@
-"""Generate composite or adiabatic saturation pulse trains."""
+"""Cmposite or adiabatic saturation pulse trains."""
 
 import numpy as np
-from pypulseq import Opts
-from pypulseq import Sequence
-from pypulseq import make_sinc_pulse
-from pypulseq import make_trapezoid
+import pypulseq as pp
 
+from mrseq.utils import sys_defaults
 from mrseq.utils.hyperbolic_secant_pulse import make_hypsec_90
 
 
 def add_adia_sat_block(
-    seq: Sequence | None = None,
+    seq: pp.Sequence | None = None,
+    system: pp.Opts | None = None,
     n_pulses: int = 3,
     max_b1: float = 20,
-    sys: Opts | None = None,
-) -> tuple[Sequence, float | None]:
+) -> tuple[pp.Sequence, float]:
     """Add adiabatic saturation pulse train to a PyPulseq Sequence object.
 
     Parameters
     ----------
     seq
         PyPulseq Sequence object
+    system
+        system limits
     n_pulses
         number of rf pulses in pulse train
     max_b1
         maximum B1 field amplitude in µT
-    sys
-        system limits
+
 
     Returns
     -------
@@ -35,34 +34,36 @@ def add_adia_sat_block(
     last_spoil_dur
         duration of the last spoiler gradient in s
     """
-    if not seq:
-        seq = Sequence()
+    # set system to default if not provided
+    if system is None:
+        system = sys_defaults
 
-    if not sys:
-        sys = Opts(max_grad=24, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s')
+    # create new sequence if not provided
+    if seq is None:
+        seq = pp.Sequence(system=system)
 
     # spoilers
-    spoil_amp0 = 0.8 * sys.max_grad  # Hz/m
-    spoil_amp1 = -0.7 * sys.max_grad  # Hz/m
-    spoil_amp2 = 0.6 * sys.max_grad  # Hz/m
+    spoil_amp0 = 0.8 * system.max_grad  # Hz/m
+    spoil_amp1 = -0.7 * system.max_grad  # Hz/m
+    spoil_amp2 = 0.6 * system.max_grad  # Hz/m
 
     rise_time = 1.0e-3  # spoiler rise time in seconds
     spoil_dur = 5.5e-3  # complete spoiler duration in seconds
 
     gx_spoil0, gy_spoil0, gz_spoil0 = (
-        make_trapezoid(channel=c, system=sys, amplitude=spoil_amp0, duration=spoil_dur, rise_time=rise_time)
+        pp.make_trapezoid(channel=c, system=system, amplitude=spoil_amp0, duration=spoil_dur, rise_time=rise_time)
         for c in ['x', 'y', 'z']
     )
     gx_spoil1, gy_spoil1, gz_spoil1 = (
-        make_trapezoid(channel=c, system=sys, amplitude=spoil_amp1, duration=spoil_dur, rise_time=rise_time)
+        pp.make_trapezoid(channel=c, system=system, amplitude=spoil_amp1, duration=spoil_dur, rise_time=rise_time)
         for c in ['x', 'y', 'z']
     )
     gx_spoil2, gy_spoil2, gz_spoil2 = (
-        make_trapezoid(channel=c, system=sys, amplitude=spoil_amp2, duration=spoil_dur, rise_time=rise_time)
+        pp.make_trapezoid(channel=c, system=system, amplitude=spoil_amp2, duration=spoil_dur, rise_time=rise_time)
         for c in ['x', 'y', 'z']
     )
 
-    hypsec_90 = make_hypsec_90(amp=max_b1, system=sys)
+    hypsec_90 = make_hypsec_90(amp=max_b1, system=system)
 
     for i in range(n_pulses):
         seq.add_block(hypsec_90)
@@ -77,12 +78,12 @@ def add_adia_sat_block(
 
 
 def add_sat_pulse_train(
-    seq: Sequence = None,
+    seq: pp.Sequence | None = None,
+    system: pp.Opts | None = None,
     n_pulses: int = 6,
     max_b1: float = 14,
-    sys: Opts = None,
     grad_strength_factor: float = 1.0,
-) -> tuple[Sequence, float | None]:
+) -> tuple[pp.Sequence, float]:
     """Create a saturation pulse train according to Chow, K. et al.
 
     For more information see:
@@ -92,12 +93,12 @@ def add_sat_pulse_train(
     ----------
     seq
         PyPulseq Sequence object
+    system
+        system limits
     n_pulses
         number of rf pulses in pulse train
     max_b1
         maximum B1 field amplitude in µT
-    sys
-        system limits
     grad_strength_factor
         increase the spoiler gradient strength by this factor
 
@@ -108,11 +109,13 @@ def add_sat_pulse_train(
     last_spoil_dur
         duration of the last spoiler gradient in s
     """
-    if not seq:
-        seq = Sequence()
+    # set system to default if not provided
+    if system is None:
+        system = sys_defaults
 
-    if not sys:
-        sys = Opts(max_grad=24, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s')
+    # create new sequence if not provided
+    if seq is None:
+        seq = pp.Sequence(system=system)
 
     flip_angles = [115, 90, 125, 85, 176, 223]
     grad_dur = [5, 5.25, 4.5, 3.75, 3.5, 3, 4]
@@ -122,21 +125,21 @@ def add_sat_pulse_train(
     grad_pe = [0, 1, -1, 1, -1, 0, 0]
     grad_ss = [-1, 1, 0, -1, 0, 1, -1]
 
-    dur_block = np.array(flip_angles) / (360 * sys.gamma * 10e-7 * max_b1)
+    dur_block = np.array(flip_angles) / (360 * system.gamma * 10e-7 * max_b1)
     dur_sinc = np.ceil(dur_block / 0.22570566672775255 * 1e5) * 1e-5  # factor equals mean of sinc shape
     dur_block = np.ceil(dur_block * 1e5) * 1e-5
 
     for i in range(n_pulses + 1):
         if grad_ro[i] != 0:
-            gx = make_trapezoid(
+            gx = pp.make_trapezoid(
                 channel='x', area=grad_ro[i] * grad_area[i], duration=grad_dur[i] / 1000, rise_time=5e-4
             )
         if grad_pe[i] != 0:
-            gy = make_trapezoid(
+            gy = pp.make_trapezoid(
                 channel='y', area=grad_pe[i] * grad_area[i], duration=grad_dur[i] / 1000, rise_time=5e-4
             )
         if grad_ss[i] != 0:
-            gz = make_trapezoid(
+            gz = pp.make_trapezoid(
                 channel='z', area=grad_ss[i] * grad_area[i], duration=grad_dur[i] / 1000, rise_time=5e-4
             )
 
@@ -156,6 +159,8 @@ def add_sat_pulse_train(
             seq.add_block(gz)
 
         if i < n_pulses:
-            seq.add_block(make_sinc_pulse(flip_angle=flip_angles[i] * np.pi / 180, duration=dur_sinc[i], system=sys))
+            seq.add_block(
+                pp.make_sinc_pulse(flip_angle=flip_angles[i] * np.pi / 180, duration=dur_sinc[i], system=system)
+            )
 
     return seq, grad_dur[-1] * 1e-3
diff --git a/src/mrseq/scripts/b0_b1_wasabiti.py b/src/mrseq/scripts/b0_b1_wasabiti.py
index 5100730..a185069 100644
--- a/src/mrseq/scripts/b0_b1_wasabiti.py
+++ b/src/mrseq/scripts/b0_b1_wasabiti.py
@@ -33,6 +33,7 @@ def wasabiti_gre_centric_kernel(
     rf_apodization: float,
     rf_spoiling_inc: float,
     adc_dwell_time: float,
+    sat_pulse_max_b1: float,
 ) -> pp.Sequence:
     """Generate a WASABI sequence for simultaneous B0 and B1 mapping using a centric-out cartesian GRE readout.
 
@@ -71,9 +72,11 @@ def wasabiti_gre_centric_kernel(
     rf_apodization
         Apodization factor of rf excitation pulse
     rf_spoiling_inc
-        Phase increment used for RF spoiling. Set to 0 to disable RF spoiling.
+        Phase increment used for RF spoiling. Set to 0 to disable RF spoiling
     adc_dwell_time
         Dwell time of ADC.
+    sat_pulse_max_b1
+        Maximum B1 field amplitude of adiabatic saturation pulse (in µT)
 
     Returns
     -------
@@ -151,7 +154,7 @@ def wasabiti_gre_centric_kernel(
     # create post preparation spoiler gradient
     prep_spoil = pp.make_trapezoid(channel='z', area=5 * gz_spoil.area, system=system)
 
-    # calculate minimum echo time
+    # calculate minimum echo time - only for sequence definitions
     min_te = (
         rf.shape_dur / 2  # time from center to end of RF pulse
         + max(rf.ringdown_time, gz.fall_time)  # RF ringdown time or gradient fall time
@@ -161,7 +164,7 @@ def wasabiti_gre_centric_kernel(
         + (k0_center_id + 0.5) * adc.dwell  # time from beginning of ADC to time point of k-space center sample
     ).item()
 
-    # calculate minimum repetition time
+    # calculate minimum repetition time - only for sequence definitions
     min_tr = (
         pp.calc_duration(gz)  # rf pulse
         + pp.calc_duration(gzr, gx_pre)  # slice selection re-phasing gradient and readout pre-winder
@@ -175,7 +178,7 @@ def wasabiti_gre_centric_kernel(
         rep_label = pp.make_label(type='SET', label='REP', value=int(rep_idx))
 
         # add adiabatic saturation pulse train and recovery time
-        seq, last_spoil_dur = add_adia_sat_block(seq=seq, sys=system)
+        seq, last_spoil_dur = add_adia_sat_block(seq=seq, system=system, max_b1=sat_pulse_max_b1)
         seq.add_block(
             pp.make_delay(
                 round_to_raster(t_recovery[rep_idx] - last_spoil_dur, raster_time=system.block_duration_raster)
@@ -318,8 +321,11 @@ def main(
     )
 
     # WASABI block pulse
-    rf_prep_duration: float = 5e-3
-    rf_prep_amplitude: float = 3.75
+    rf_prep_duration = 5e-3
+    rf_prep_amplitude = 3.75
+
+    # Saturation pulse
+    sat_pulse_max_b1 = 20  # (µT)
 
     seq = wasabiti_gre_centric_kernel(
         system=system,
@@ -340,6 +346,7 @@ def main(
         rf_apodization=rf_apodization,
         rf_spoiling_inc=rf_spoiling_inc,
         adc_dwell_time=adc_dwell_time,
+        sat_pulse_max_b1=sat_pulse_max_b1,
     )
 
     # check timing of the sequence
diff --git a/src/mrseq/utils/hyperbolic_secant_pulse.py b/src/mrseq/utils/hyperbolic_secant_pulse.py
index 14d39fb..038f5ad 100644
--- a/src/mrseq/utils/hyperbolic_secant_pulse.py
+++ b/src/mrseq/utils/hyperbolic_secant_pulse.py
@@ -1,4 +1,4 @@
-"""Functions to create an adiabatic hyperbolic secant full passage pulse."""
+"""Adiabatic hyperbolic secant full passage pulse."""
 
 from types import SimpleNamespace
 
@@ -6,7 +6,6 @@
 from pypulseq import Opts
 
 
-# Todo: import from BMCTool instead ?!
 def calculate_phase(
     frequency: np.ndarray, duration: float, samples: int, shift_idx: int = -1, pos_offsets: bool = False
 ) -> np.ndarray:
@@ -22,7 +21,6 @@ def calculate_phase(
     return phase
 
 
-# Todo: import from BMCTool instead ?!
 def create_arbitrary_pulse_with_phase(
     signal: np.ndarray,
     flip_angle: float,

From eec35d2e85e549ef20205896b2255001ae44e3a7 Mon Sep 17 00:00:00 2001
From: Christoph Kolbitsch <christoph.kolbitsch@ptb.de>
Date: Tue, 24 Feb 2026 13:40:56 +0000
Subject: [PATCH 4/6] noise added and gy scaling

---
 src/mrseq/scripts/b0_b1_wasabiti.py | 43 ++++++++++++++++++-----------
 1 file changed, 27 insertions(+), 16 deletions(-)

diff --git a/src/mrseq/scripts/b0_b1_wasabiti.py b/src/mrseq/scripts/b0_b1_wasabiti.py
index a185069..3028d05 100644
--- a/src/mrseq/scripts/b0_b1_wasabiti.py
+++ b/src/mrseq/scripts/b0_b1_wasabiti.py
@@ -23,6 +23,7 @@ def wasabiti_gre_centric_kernel(
     fov_xy: float,
     n_readout: int,
     readout_oversampling: float,
+    gx_pre_duration: float,
     n_phase_encoding: int,
     slice_thickness: float,
     rf_prep_duration: float,
@@ -55,6 +56,8 @@ def wasabiti_gre_centric_kernel(
         Number of frequency encoding steps.
     readout_oversampling
         Readout oversampling factor, commonly 2. This reduces aliasing artifacts.
+    gx_pre_duration
+        Duration of readout pre-winder gradient (in seconds)
     n_phase_encoding
         Number of phase encoding steps.
     slice_thickness
@@ -122,13 +125,6 @@ def wasabiti_gre_centric_kernel(
         use='excitation',
     )
 
-    # define centric-out phase encoding steps
-    kpe, _ = cartesian_phase_encoding(
-        n_phase_encoding=n_phase_encoding,
-        acceleration=1,
-        sampling_order='low_high',
-    )
-
     # create readout gradient and ADC
     delta_k = 1 / (fov_xy * readout_oversampling)
     n_readout_with_oversampling = int(n_readout * readout_oversampling)
@@ -142,12 +138,24 @@ def wasabiti_gre_centric_kernel(
     adc = pp.make_adc(num_samples=n_readout_with_oversampling, duration=gx.flat_time, delay=gx.rise_time, system=system)
 
     # create frequency encoding pre- and re-winder gradient
-    gx_pre = pp.make_trapezoid(channel='x', area=-gx.area / 2 - delta_k / 2, system=system)
-    gx_post = pp.make_trapezoid(channel='x', area=-gx.area / 2 + delta_k / 2, system=system)
+    gx_pre = pp.make_trapezoid(channel='x', area=-gx.area / 2 - delta_k / 2, duration=gx_pre_duration, system=system)
+    gx_post = pp.make_trapezoid(channel='x', area=-gx.area / 2 + delta_k / 2, duration=gx_pre_duration, system=system)
     k0_center_id = np.where((np.arange(n_readout_with_oversampling) - n_readout_with_oversampling / 2) * delta_k == 0)[
         0
     ][0]
 
+    # define centric-out phase encoding steps
+    kpe, _ = cartesian_phase_encoding(
+        n_phase_encoding=n_phase_encoding,
+        acceleration=1,
+        sampling_order='low_high',
+    )
+
+    # phase encoding gradient
+    gy_pre_max = pp.make_trapezoid(
+        channel='y', area=delta_k * n_phase_encoding / 2, duration=gx_pre_duration, system=system
+    )
+
     # create readout spoiler gradient
     gz_spoil = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system)
 
@@ -210,9 +218,7 @@ def wasabiti_gre_centric_kernel(
                 rf_phase = divmod(rf_phase + rf_inc, 360.0)[1]
 
             # calculate phase encoding gradient for current phase encoding step
-            gy_pre = pp.make_trapezoid(
-                channel='y', area=pe_idx * delta_k, duration=pp.calc_duration(gzr, gx_pre), system=system
-            )
+            gy_pre = pp.scale_grad(gy_pre_max, (pe_idx - n_phase_encoding / 2) / (n_phase_encoding / 2))
 
             # add slice-selective rf pulse
             seq.add_block(rf, gz)
@@ -224,10 +230,13 @@ def wasabiti_gre_centric_kernel(
             seq.add_block(gx, adc, pe_label, rep_label)
 
             # add x and y re-winder and spoiler gradient in z-direction
-            gy_post = pp.make_trapezoid(
-                channel='y', amplitude=-gy_pre.amplitude, rise_time=gy_pre.rise_time, flat_time=gy_pre.flat_time
-            )
-            seq.add_block(gx_post, gy_post, gz_spoil)
+            seq.add_block(gx_post, pp.scale_grad(gy_pre, -1), gz_spoil)
+
+    # obtain noise samples
+    seq.add_block(pp.make_label(label='LIN', type='SET', value=0), pp.make_label(label='REP', type='SET', value=0))
+    seq.add_block(adc, pp.make_label(label='NOISE', type='SET', value=True))
+    seq.add_block(pp.make_label(label='NOISE', type='SET', value=False))
+    seq.add_block(pp.make_delay(system.rf_dead_time))
 
     # write all required parameters in the seq-file header/definitions
     seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness])
@@ -313,6 +322,7 @@ def main(
 
     # define ADC and gradient timing
     n_readout_with_oversampling = int(n_readout * readout_oversampling)
+    gx_pre_duration = 1.0e-3  # duration of readout pre-winder gradient [s]
     adc_dwell_time = round_to_raster(
         1.0 / (receiver_bandwidth_per_pixel * n_readout_with_oversampling), system.adc_raster_time
     )
@@ -336,6 +346,7 @@ def main(
         fov_xy=fov_xy,
         n_readout=n_readout,
         readout_oversampling=readout_oversampling,
+        gx_pre_duration=gx_pre_duration,
         n_phase_encoding=n_phase_encoding,
         slice_thickness=slice_thickness,
         rf_prep_duration=rf_prep_duration,

From 037fa2673ab75793530250bf15dfdf9e92d91734 Mon Sep 17 00:00:00 2001
From: Christoph Kolbitsch <christoph.kolbitsch@ptb.de>
Date: Tue, 24 Feb 2026 13:50:05 +0000
Subject: [PATCH 5/6] kpe fixes

---
 src/mrseq/scripts/b0_b1_wasabiti.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mrseq/scripts/b0_b1_wasabiti.py b/src/mrseq/scripts/b0_b1_wasabiti.py
index 3028d05..482633c 100644
--- a/src/mrseq/scripts/b0_b1_wasabiti.py
+++ b/src/mrseq/scripts/b0_b1_wasabiti.py
@@ -153,7 +153,7 @@ def wasabiti_gre_centric_kernel(
 
     # phase encoding gradient
     gy_pre_max = pp.make_trapezoid(
-        channel='y', area=delta_k * n_phase_encoding / 2, duration=gx_pre_duration, system=system
+        channel='y', area=delta_k * readout_oversampling * n_phase_encoding / 2, duration=gx_pre_duration, system=system
     )
 
     # create readout spoiler gradient
@@ -218,7 +218,7 @@ def wasabiti_gre_centric_kernel(
                 rf_phase = divmod(rf_phase + rf_inc, 360.0)[1]
 
             # calculate phase encoding gradient for current phase encoding step
-            gy_pre = pp.scale_grad(gy_pre_max, (pe_idx - n_phase_encoding / 2) / (n_phase_encoding / 2))
+            gy_pre = pp.scale_grad(gy_pre_max, pe_idx / (n_phase_encoding / 2))
 
             # add slice-selective rf pulse
             seq.add_block(rf, gz)

From c7c17a40f459ecd8243812d90882389ad2790a2f Mon Sep 17 00:00:00 2001
From: Christoph Kolbitsch <christoph.kolbitsch@ptb.de>
Date: Tue, 24 Feb 2026 14:15:46 +0000
Subject: [PATCH 6/6] scale spoiling

---
 src/mrseq/scripts/b0_b1_wasabiti.py | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/mrseq/scripts/b0_b1_wasabiti.py b/src/mrseq/scripts/b0_b1_wasabiti.py
index 482633c..56d41f9 100644
--- a/src/mrseq/scripts/b0_b1_wasabiti.py
+++ b/src/mrseq/scripts/b0_b1_wasabiti.py
@@ -35,6 +35,8 @@ def wasabiti_gre_centric_kernel(
     rf_spoiling_inc: float,
     adc_dwell_time: float,
     sat_pulse_max_b1: float,
+    gz_spoil_duration: float,
+    gz_spoil_area: float,
 ) -> pp.Sequence:
     """Generate a WASABI sequence for simultaneous B0 and B1 mapping using a centric-out cartesian GRE readout.
 
@@ -80,6 +82,10 @@ def wasabiti_gre_centric_kernel(
         Dwell time of ADC.
     sat_pulse_max_b1
         Maximum B1 field amplitude of adiabatic saturation pulse (in µT)
+    gz_spoil_duration
+        Duration of spoiler gradient (in seconds).
+    gz_spoil_area
+        Area of spoiler gradient (in 1/meters = Hz/m * s).
 
     Returns
     -------
@@ -157,10 +163,10 @@ def wasabiti_gre_centric_kernel(
     )
 
     # create readout spoiler gradient
-    gz_spoil = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system)
+    gz_spoil = pp.make_trapezoid(channel='z', system=system, area=gz_spoil_area, duration=gz_spoil_duration)
 
     # create post preparation spoiler gradient
-    prep_spoil = pp.make_trapezoid(channel='z', area=5 * gz_spoil.area, system=system)
+    prep_spoil = pp.make_trapezoid(channel='z', area=5 * gz_spoil.area, system=system, duration=5 * gz_spoil_duration)
 
     # calculate minimum echo time - only for sequence definitions
     min_te = (
@@ -330,6 +336,10 @@ def main(
         n_readout_with_oversampling, adc_dwell_time, system.grad_raster_time, system.adc_raster_time
     )
 
+    # define spoiling
+    gz_spoil_duration = 0.8e-3  # duration of spoiler gradient [s]
+    gz_spoil_area = 4 / slice_thickness  # area / zeroth gradient moment of spoiler gradient
+
     # WASABI block pulse
     rf_prep_duration = 5e-3
     rf_prep_amplitude = 3.75
@@ -358,6 +368,8 @@ def main(
         rf_spoiling_inc=rf_spoiling_inc,
         adc_dwell_time=adc_dwell_time,
         sat_pulse_max_b1=sat_pulse_max_b1,
+        gz_spoil_area=gz_spoil_area,
+        gz_spoil_duration=gz_spoil_duration,
     )
 
     # check timing of the sequence