diff --git a/gre_pypulseq_sigpy_demo1.py b/gre_pypulseq_sigpy_demo1.py
new file mode 100644
index 0000000..2662787
--- /dev/null
+++ b/gre_pypulseq_sigpy_demo1.py
@@ -0,0 +1,176 @@
+import math
+
+import numpy as np
+import pypulseq as pp
+import sigpy.mri.rf as rf_ext
+
+import sigpy2pulseq as sp
+
+# ======
+# SETUP
+# ======
+# Create a new sequence object and library options
+seq = pp.Sequence()
+
+# Define FOV and resolution
+fov = 256e-3
+Nx = 256
+Ny = 256
+alpha = 5  # flip angle
+slice_thickness = 3e-3  # slice
+TE = np.array([4.3e-3])
+TR = 10e-3
+rf_spoiling_inc = 117  # RF spoiling increment
+
+ext_pulse_library = 1  # 0 - pypulseq, 1 - sigpy
+ext_pulse_type = 'hypsec'
+disp_pulse = 1
+
+system = pp.Opts(max_grad=28, grad_unit='mT/m', max_slew=150, slew_unit='T/m/s', rf_ringdown_time=20e-6,
+                 rf_dead_time=100e-6, adc_dead_time=10e-6, rf_raster_time=1e-6)
+
+# Common pulse design parameters
+duration = 3e-3  # seconds
+tb = 4
+
+# ======
+# CREATE EVENTS
+# ======
+if (ext_pulse_library):
+    # Declare pulse configuration inputs here
+    pulse_freq_offset = 0.0  # if designing rf pulse without a FM
+    if ext_pulse_type == 'slr':
+        pulse = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls', d1=0.01,
+                                d2=0.01,
+                                cancel_alpha_phs=True)
+
+    if ext_pulse_type == 'sms':
+        pulse_in = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls',
+                                   d1=0.01, d2=0.01,
+                                   cancel_alpha_phs=False)
+        pulse = rf_ext.multiband.mb_rf(pulse_in, n_bands=3, band_sep=25, phs_0_pt='None')
+
+    if ext_pulse_type == 'hypsec':
+        # Let us design two pulses here
+        # - STA for excitation
+        pulse = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls', d1=0.01,
+                                d2=0.01,
+                                cancel_alpha_phs=True)
+
+        # - Hyperbolic secant pulse for inversion as part of magnetization prep
+        duration_hypsec = 2 * duration  # adiabatic pulses need more time
+        pulse_hypsec, pulse_freq_offset = rf_ext.adiabatic.hypsec(n=2002,
+                                                                  beta=500,
+                                                                  mu=60, dur=duration_hypsec)
+
+        rf_inv = sp.sig_2_seq(pulse=pulse_hypsec, flip_angle=alpha * math.pi / 180, system=system,
+                              duration=duration_hypsec,
+                              slice_thickness=slice_thickness,
+                              return_gz=False, time_bw_product=tb, rf_freq=pulse_freq_offset)
+
+    # Display if needed
+    if (disp_pulse):
+        if (ext_pulse_type == 'hypsec'):
+            sp.disp_pulse(pulse_hypsec, pulse_freq_offset, tb, duration_hypsec, system)
+        else:
+            sp.disp_pulse(pulse, 0, tb, duration, system)
+
+    # Convert and integrate with pypulseq
+    rf, gz, gzr, _ = sp.sig_2_seq(pulse=pulse, flip_angle=alpha * math.pi / 180, system=system, duration=duration,
+                                  slice_thickness=slice_thickness,
+                                  return_gz=True, time_bw_product=tb, rf_freq=0)
+else:
+    rf, gz, gzr = pp.make_sinc_pulse(flip_angle=alpha * math.pi / 180, duration=duration,
+                                     slice_thickness=slice_thickness,
+                                     apodization=0.5, time_bw_product=tb, system=system, return_gz=True)
+
+# Define other gradients and ADC events
+delta_k = 1 / fov
+gx = pp.make_trapezoid(channel='x', flat_area=Nx * delta_k, flat_time=3.2e-3, system=system)
+adc = pp.make_adc(num_samples=Nx, duration=gx.flat_time, delay=gx.rise_time, system=system)
+gx_pre = pp.make_trapezoid(channel='x', area=-gx.area / 2, duration=1e-3, system=system)
+gz_reph = pp.make_trapezoid(channel='z', area=-gz.area / 2, duration=1e-3, system=system)
+phase_areas = (np.arange(Ny) - Ny / 2) * delta_k
+
+# gradient spoiling
+gx_spoil = pp.make_trapezoid(channel='x', area=2 * Nx * delta_k, system=system)
+gz_spoil = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system)
+
+# Calculate timing
+delay_TE = np.ceil((TE - pp.calc_duration(gx_pre) - gz.fall_time - gz.flat_time / 2 - pp.calc_duration(
+    gx) / 2) / seq.grad_raster_time) * seq.grad_raster_time
+delay_TR = np.ceil((TR - pp.calc_duration(gz) - pp.calc_duration(gx_pre) - pp.calc_duration(
+    gx) - delay_TE) / seq.grad_raster_time) * seq.grad_raster_time
+
+assert np.all(delay_TE >= 0)
+assert np.all(delay_TR >= pp.calc_duration(gx_spoil, gz_spoil))
+
+rf_phase = 0
+rf_inc = 0
+
+# ======
+# CONSTRUCT SEQUENCE
+# ======
+# Loop over phase encodes and define sequence blocks
+# Insert adiabatic inversion here if chosen - 12 ms - for demo purposes only
+
+for i in range(Ny):
+
+    for j in range(len(TE)):
+        rf.phase_offset = rf_phase / 180 * np.pi
+        adc.phase_offset = rf_phase / 180 * np.pi
+        rf_inc = divmod(rf_inc + rf_spoiling_inc, 360.0)[1]
+        rf_phase = divmod(rf_phase + rf_inc, 360.0)[1]
+
+        if ext_pulse_type == 'hypsec':
+            # if np.mod(i, 1) == 0:
+            if i==0:
+                print('Mag prep now')
+                for n in range(len(rf_inv)):
+                    seq.add_block(rf_inv[n])
+                seq.add_block(pp.make_delay(300e-3))
+
+        # seq.plot()
+        seq.add_block(rf, gz)
+
+        gy_pre = pp.make_trapezoid(channel='y', area=phase_areas[i], duration=pp.calc_duration(gx_pre), system=system)
+        seq.add_block(gx_pre, gy_pre, gz_reph)
+        seq.add_block(pp.make_delay(delay_TE[j]))
+        seq.add_block(gx, adc)
+        gy_pre.amplitude = -gy_pre.amplitude
+        seq.add_block(pp.make_delay(delay_TR[j]), gx_spoil, gy_pre, gz_spoil)
+
+ok, error_report = seq.check_timing()  # Check whether the timing of the sequence is correct
+if ok:
+    print('Timing check passed successfully')
+else:
+    print('Timing check failed. Error listing follows:')
+    [print(e) for e in error_report]
+
+# ======
+# VISUALIZATION
+# ======
+# seq.plot()
+
+# Trajectory calculation and plotting
+# ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc = seq.calculate_kspace()
+# time_axis = np.arange(1, ktraj.shape[1] + 1) * system.grad_raster_time
+# plt.figure()
+# plt.plot(time_axis, ktraj.T)  # Plot the entire k-space trajectory
+# plt.plot(t_adc, ktraj_adc[0], '.')  # Plot sampling points on the kx-axis
+# plt.figure()
+# plt.plot(ktraj[0], ktraj[1], 'b')  # 2D plot
+# plt.axis('equal')  # Enforce aspect ratio for the correct trajectory display
+# plt.plot(ktraj_adc[0], ktraj_adc[1], 'r.')  # Plot  sampling points
+# plt.show()
+
+# Prepare the sequence output for the scanner
+seq.set_definition('FOV', [fov, fov, slice_thickness])
+seq.set_definition('Name', 'gre')
+
+seq.write('gre_pypulseq_ext_' + str(ext_pulse_type) + '.seq')
+
+# Very optional slow step, but useful for testing during development e.g. for the real TE, TR or for staying within
+# slew-rate limits
+# rep = seq.test_report()
+# print(rep)
diff --git a/gre_pypulseq_sigpy_demo2.py b/gre_pypulseq_sigpy_demo2.py
new file mode 100644
index 0000000..eaeecad
--- /dev/null
+++ b/gre_pypulseq_sigpy_demo2.py
@@ -0,0 +1,176 @@
+import math
+
+import numpy as np
+import pypulseq as pp
+import sigpy.mri.rf as rf_ext
+
+import sigpy2pulseq as sp
+
+# ======
+# SETUP
+# ======
+# Create a new sequence object and library options
+seq = pp.Sequence()
+
+# Define FOV and resolution
+fov = 256e-3
+Nx = 256
+Ny = 256
+alpha = 90  # flip angle
+slice_thickness = 3e-3  # slice
+TE = np.array([4.3e-3])
+TR = 300e-3
+rf_spoiling_inc = 117  # RF spoiling increment
+
+ext_pulse_library = 1  # 0 - pypulseq, 1 - sigpy
+ext_pulse_type = 'slr'
+disp_pulse = 0
+
+system = pp.Opts(max_grad=28, grad_unit='mT/m', max_slew=150, slew_unit='T/m/s', rf_ringdown_time=20e-6,
+                 rf_dead_time=100e-6, adc_dead_time=10e-6, rf_raster_time=1e-6)
+
+# Common pulse design parameters
+duration = 3e-3  # seconds
+tb = 4
+
+# ======
+# CREATE EVENTS
+# ======
+if (ext_pulse_library):
+    # Declare pulse configuration inputs here
+    pulse_freq_offset = 0.0  # if designing rf pulse without a FM
+    if ext_pulse_type == 'slr':
+        pulse = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls', d1=0.01,
+                                d2=0.01,
+                                cancel_alpha_phs=True)
+
+    if ext_pulse_type == 'sms':
+        pulse_in = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls',
+                                   d1=0.01, d2=0.01,
+                                   cancel_alpha_phs=False)
+        pulse = rf_ext.multiband.mb_rf(pulse_in, n_bands=3, band_sep=25, phs_0_pt='None')
+
+    if ext_pulse_type == 'hypsec':
+        # Let us design two pulses here
+        # - STA for excitation
+        pulse = rf_ext.slr.dzrf(n=int(round(duration / system.rf_raster_time)), tb=tb, ptype='st', ftype='ls', d1=0.01,
+                                d2=0.01,
+                                cancel_alpha_phs=True)
+
+        # - Hyperbolic secant pulse for inversion as part of magnetization prep
+        duration_hypsec = 4 * duration  # adiabatic pulses need more time
+        pulse_hypsec, pulse_freq_offset = rf_ext.adiabatic.hypsec(n=2002,
+                                                                  beta=500,
+                                                                  mu=60, dur=duration_hypsec)
+
+        rf_inv = sp.sig_2_seq(pulse=pulse_hypsec, flip_angle=alpha * math.pi / 180, system=system,
+                              duration=duration_hypsec,
+                              slice_thickness=slice_thickness,
+                              return_gz=False, time_bw_product=tb, rf_freq=pulse_freq_offset)
+
+    # Display if needed
+    if (disp_pulse):
+        if (ext_pulse_type == 'hypsec'):
+            sp.disp_pulse(pulse_hypsec, pulse_freq_offset, tb, duration_hypsec, system)
+        else:
+            sp.disp_pulse(pulse, 0, tb, duration, system)
+
+    # Convert and integrate with pypulseq
+    rf, gz, gzr, _ = sp.sig_2_seq(pulse=pulse, flip_angle=alpha * math.pi / 180, system=system, duration=duration,
+                                  slice_thickness=slice_thickness,
+                                  return_gz=True, time_bw_product=tb, rf_freq=0)
+else:
+    rf, gz, gzr = pp.make_sinc_pulse(flip_angle=alpha * math.pi / 180, duration=duration,
+                                     slice_thickness=slice_thickness,
+                                     apodization=0.5, time_bw_product=tb, system=system, return_gz=True)
+
+# Define other gradients and ADC events
+delta_k = 1 / fov
+gx = pp.make_trapezoid(channel='x', flat_area=Nx * delta_k, flat_time=3.2e-3, system=system)
+adc = pp.make_adc(num_samples=Nx, duration=gx.flat_time, delay=gx.rise_time, system=system)
+gx_pre = pp.make_trapezoid(channel='x', area=-gx.area / 2, duration=1e-3, system=system)
+gz_reph = pp.make_trapezoid(channel='z', area=-gz.area / 2, duration=1e-3, system=system)
+phase_areas = (np.arange(Ny) - Ny / 2) * delta_k
+
+# gradient spoiling
+gx_spoil = pp.make_trapezoid(channel='x', area=2 * Nx * delta_k, system=system)
+gz_spoil = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system)
+
+# Calculate timing
+delay_TE = np.ceil((TE - pp.calc_duration(gx_pre) - gz.fall_time - gz.flat_time / 2 - pp.calc_duration(
+    gx) / 2) / seq.grad_raster_time) * seq.grad_raster_time
+delay_TR = np.ceil((TR - pp.calc_duration(gz) - pp.calc_duration(gx_pre) - pp.calc_duration(
+    gx) - delay_TE) / seq.grad_raster_time) * seq.grad_raster_time
+
+assert np.all(delay_TE >= 0)
+assert np.all(delay_TR >= pp.calc_duration(gx_spoil, gz_spoil))
+
+rf_phase = 0
+rf_inc = 0
+# demo slice profile
+gz.channel = 'x'
+# ======
+# CONSTRUCT SEQUENCE
+# ======
+# Loop over phase encodes and define sequence blocks
+# Insert adiabatic inversion here if chosen - 12 ms - for demo purposes only
+
+for i in range(Ny):
+
+    for j in range(len(TE)):
+        rf.phase_offset = rf_phase / 180 * np.pi
+        adc.phase_offset = rf_phase / 180 * np.pi
+        rf_inc = divmod(rf_inc + rf_spoiling_inc, 360.0)[1]
+        rf_phase = divmod(rf_phase + rf_inc, 360.0)[1]
+
+        if ext_pulse_type == 'hypsec':
+            if np.mod(i, 64) == 0:
+                print('Mag prep now')
+                for n in range(len(rf_inv)):
+                    seq.add_block(rf_inv[n])
+                seq.add_block(pp.make_delay(300e-3))
+
+        # seq.plot()
+        seq.add_block(rf, gz)
+
+        gy_pre = pp.make_trapezoid(channel='y', area=phase_areas[i], duration=pp.calc_duration(gx_pre), system=system)
+        seq.add_block(gx_pre, gy_pre, gz_reph)
+        seq.add_block(pp.make_delay(delay_TE[j]))
+        seq.add_block(gx, adc)
+        gy_pre.amplitude = -gy_pre.amplitude
+        seq.add_block(pp.make_delay(delay_TR[j]), gx_spoil, gy_pre, gz_spoil)
+
+ok, error_report = seq.check_timing()  # Check whether the timing of the sequence is correct
+if ok:
+    print('Timing check passed successfully')
+else:
+    print('Timing check failed. Error listing follows:')
+    [print(e) for e in error_report]
+
+# ======
+# VISUALIZATION
+# ======
+# seq.plot()
+
+# Trajectory calculation and plotting
+# ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc = seq.calculate_kspace()
+# time_axis = np.arange(1, ktraj.shape[1] + 1) * system.grad_raster_time
+# plt.figure()
+# plt.plot(time_axis, ktraj.T)  # Plot the entire k-space trajectory
+# plt.plot(t_adc, ktraj_adc[0], '.')  # Plot sampling points on the kx-axis
+# plt.figure()
+# plt.plot(ktraj[0], ktraj[1], 'b')  # 2D plot
+# plt.axis('equal')  # Enforce aspect ratio for the correct trajectory display
+# plt.plot(ktraj_adc[0], ktraj_adc[1], 'r.')  # Plot  sampling points
+# plt.show()
+
+# Prepare the sequence output for the scanner
+seq.set_definition('FOV', [fov, fov, slice_thickness])
+seq.set_definition('Name', 'gre')
+
+seq.write('gre_pypulseq_ext_profile_' + str(alpha) + '_' + str(ext_pulse_type) + '.seq')
+
+# Very optional slow step, but useful for testing during development e.g. for the real TE, TR or for staying within
+# slew-rate limits
+# rep = seq.test_report()
+# print(rep)
diff --git a/recon_cart.py b/recon_cart.py
new file mode 100644
index 0000000..b78a9d0
--- /dev/null
+++ b/recon_cart.py
@@ -0,0 +1,18 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from rawdatarinator import twixread
+
+data = twixread('./Data/110121/slr_90.dat', A=True)
+im_data = np.zeros(data.shape)
+
+for ch in range(data.shape[3]):  #
+    temp = np.squeeze(data[:, :, 0, ch])
+    temp2 = temp
+    im_temp = np.fft.fftshift(np.fft.fft2(temp2))
+    im_data[:, :, 0, ch] = np.abs(im_temp)
+
+im_data_sos = np.sum(((np.squeeze(im_data))), axis=-1)
+print(im_data_sos.shape)
+im_data_sos = np.rot90(im_data_sos, k=1)
+plt.imshow(im_data_sos, cmap='gray')
+plt.show()
diff --git a/sigpy2pulseq.py b/sigpy2pulseq.py
new file mode 100644
index 0000000..72697bf
--- /dev/null
+++ b/sigpy2pulseq.py
@@ -0,0 +1,241 @@
+import math
+from types import SimpleNamespace
+from typing import Tuple, Union
+
+import numpy as np
+import pypulseq as pp
+import sigpy.mri.rf as rf_ext
+import sigpy.plot as pl
+from matplotlib import pyplot as plt
+from pypulseq.make_trap_pulse import make_trapezoid
+from pypulseq.opts import Opts
+from scipy import integrate
+
+
+def sig_2_seq(pulse: float, flip_angle: float, delay: float = 0, duration: float = 0,
+              freq_offset: float = 0, center_pos: float = 0.5, max_grad: float = 0, max_slew: float = 0,
+              phase_offset: float = 0, return_gz: bool = True, slice_thickness: float = 0, system: Opts = Opts(),
+              time_bw_product: float = 4, rf_freq: float = 0,
+              use: str = str()) -> Union[SimpleNamespace,
+                                         Tuple[SimpleNamespace, SimpleNamespace,
+                                               SimpleNamespace]]:
+    """
+    Creates a radio-frequency sinc pulse event using the sigpy rf pulse library and optionally accompanying slice select, slice select rephasing
+    trapezoidal gradient events.
+
+    Parameters
+    ----------
+    pulse : float
+        Sigpy designed RF pulse
+    flip_angle : float
+        Flip angle in radians.
+    apodization : float, optional, default=0
+        Apodization.
+    center_pos : float, optional, default=0.5
+        Position of peak.5 (midway).
+    delay : float, optional, default=0
+        Delay in milliseconds (ms).
+    duration : float, optional, default=0
+        Duration in milliseconds (ms).
+    freq_offset : float, optional, default=0
+        Frequency offset in Hertz (Hz).
+    max_grad : float, optional, default=0
+        Maximum gradient strength of accompanying slice select trapezoidal event.
+    max_slew : float, optional, default=0
+        Maximum slew rate of accompanying slice select trapezoidal event.
+    phase_offset : float, optional, default=0
+        Phase offset in Hertz (Hz).
+    return_gz:bool, default=False
+        Boolean flag to indicate if slice-selective gradient has to be returned.
+    slice_thickness : float, optional, default=0
+        Slice thickness of accompanying slice select trapezoidal event. The slice thickness determines the area of the
+        slice select event.
+    system : Opts, optional
+        System limits. Default is a system limits object initialised to default values.
+    time_bw_product : float, optional, default=4
+        Time-bandwidth product.
+    rf_freq : float, optional, default = 0
+        frequency offset
+    use : str, optional, default=str()
+        Use of radio-frequency sinc pulse. Must be one of 'excitation', 'refocusing' or 'inversion'.
+
+
+    Returns
+    -------
+    rf : SimpleNamespace
+        Radio-frequency sinc pulse event.
+    gz : SimpleNamespace, optional
+        Accompanying slice select trapezoidal gradient event. Returned only if `slice_thickness` is provided.
+    gzr : SimpleNamespace, optional
+        Accompanying slice select rephasing trapezoidal gradient event. Returned only if `slice_thickness` is provided.
+
+    Raises
+    ------
+    ValueError
+        If invalid `use` parameter was passed. Must be one of 'excitation', 'refocusing' or 'inversion'.
+        If `return_gz=True` and `slice_thickness` was not provided.
+    """
+
+    valid_use_pulses = ['excitation', 'refocusing', 'inversion']
+    if use != '' and use not in valid_use_pulses:
+        raise ValueError(
+            f"Invalid use parameter. Must be one of 'excitation', 'refocusing' or 'inversion'. Passed: {use}")
+    if np.sum(rf_freq) == 0:
+        # Step 1: Get signal, t
+        [signal, t] = get_signal(pulse=pulse, duration=duration, flip_angle=flip_angle, system=system)
+
+        rfp = SimpleNamespace()
+        rfp.type = 'rf'
+        rfp.signal = signal
+        rfp.t = t
+        rfp.freq_offset = rf_freq
+        rfp.phase_offset = phase_offset
+        rfp.dead_time = system.rf_dead_time
+        rfp.ringdown_time = system.rf_ringdown_time
+        rfp.delay = delay
+
+        if use != '':
+            rfp.use = use
+
+        if rfp.dead_time > rfp.delay:
+            rfp.delay = rfp.dead_time
+
+        if return_gz:
+            if slice_thickness == 0:
+                raise ValueError('Slice thickness must be provided')
+
+            if max_grad > 0:
+                system.max_grad = max_grad
+
+            if max_slew > 0:
+                system.max_slew = max_slew
+            BW = time_bw_product / duration
+            amplitude = BW / slice_thickness
+            area = amplitude * duration
+            gz = make_trapezoid(channel='z', system=system, flat_time=duration, flat_area=area)
+            gzr = make_trapezoid(channel='z', system=system, area=-area * (1 - center_pos) - 0.5 * (gz.area - area))
+
+            if rfp.delay > gz.rise_time:
+                gz.delay = math.ceil((rfp.delay - gz.rise_time) / system.grad_raster_time) * system.grad_raster_time
+
+            if rfp.delay < (gz.rise_time + gz.delay):
+                rfp.delay = gz.rise_time + gz.delay
+
+        if rfp.ringdown_time > 0:
+            t_fill = np.arange(1, round(rfp.ringdown_time / 1e-6) + 1) * 1e-6
+            rfp.t = np.concatenate((rfp.t, rfp.t[-1] + t_fill))
+            rfp.signal = np.concatenate((rfp.signal, np.zeros(len(t_fill))))
+
+        # Following 2 lines of code are workarounds for numpy returning 3.14... for np.angle(-0.00...)
+        negative_zero_indices = np.where(rfp.signal == -0.0)
+        rfp.signal[negative_zero_indices] = 0
+
+        if return_gz:
+            return rfp, gz, gzr, pulse
+        else:
+            return rfp
+    else:
+        print('Working on a FM pulse')
+        rf_train = get_RF_train(am_pulse=pulse, fm_pulse=rf_freq, duration=duration, system=system,
+                                flip_angle=flip_angle)
+
+        return rf_train
+
+
+def get_signal(pulse: float, duration: float, flip_angle: float, system: Opts = Opts()):
+    # N = int(round(duration / 1e-6))
+    N = int(round(duration / system.rf_raster_time))
+    t = np.arange(1, N + 1) * system.rf_raster_time
+    flip = np.sum(pulse) * system.rf_raster_time * 2 * np.pi
+    signal = pulse * flip_angle / flip
+    return signal, t
+
+
+def disp_pulse(am_pulse: float, fm_pulse: float = 0, tb=4, duration=3e-3, system: Opts = Opts()):
+    if np.sum(fm_pulse) == 0:
+        pl.LinePlot(am_pulse)
+
+        # Simulate it
+        # x =np.arange(-20 * tb, 20 * tb, 40 * tb / 2000)
+        a = 20
+        x = np.arange(-a * tb, a * tb, 40 * tb / 2000)
+        [a, b] = rf_ext.sim.abrm(am_pulse,
+                                 x,
+                                 True)
+        Mxy = 2 * np.multiply(np.conj(a), b)
+        pl.LinePlot(Mxy)
+    else:
+        T = duration
+        Ns = am_pulse.shape[0]
+        t = np.linspace(-T / 2, T / 2, Ns, endpoint=True) * 1000
+        plt.figure()
+        plt.plot(t, np.abs(am_pulse))
+        plt.xlabel('ms')
+        plt.ylabel('a.u.')
+        plt.title('|AM|')
+        plt.figure()
+        plt.plot(t, fm_pulse / 1000)
+        plt.xlabel('ms')
+        plt.ylabel('kHz')
+        plt.title('FM')
+
+        # Simulate it
+        b1 = np.arange(0, 1, 0.01)  # b1 grid we simulate the pulse over, Gauss
+        b1 = np.reshape(b1, (np.size(b1), 1))
+        a = np.zeros(np.shape(b1), dtype='complex')
+        b = np.zeros(np.shape(b1), dtype='complex')
+        for ii in range(0, np.size(b1)):
+            [a[ii], b[ii]] = rf_ext.sim.abrm_nd(2 * np.pi * (T / t.shape[0]) * 4258 * b1[ii] * am_pulse, np.ones(1),
+                                                T / t.shape[0] * np.reshape(fm_pulse, (np.size(fm_pulse), 1)))
+        Mz = 1 - 2 * np.abs(b) ** 2
+        plt.figure()
+        plt.plot(b1, Mz)
+        plt.xlabel('Gauss')
+        plt.ylabel('Mz')
+
+
+def get_RF_train(am_pulse, fm_pulse, duration, system, flip_angle):
+    #scale the amplitude according to pulseq requirements
+
+
+    RF_train = []
+    Ns = am_pulse.shape[0]
+    Tseg = np.round(duration / Ns, 6)
+    t_rf_raster = system.rf_raster_time
+    system.rf_raster_time = Tseg
+    [signal, _] = get_signal(pulse=am_pulse, duration=duration, flip_angle=flip_angle, system=system)
+
+
+    t_block = np.round(0.5 * Tseg, 6)  # Fixed at 50% duty cycle
+    system.rf_raster_time = t_rf_raster
+    system.rf_ringdown_time = np.round(0.5 * Tseg, 6)  # To ensure 50% duty cycle
+
+    ph = 0
+    alpha = 0
+
+    for n in range(1, Ns):
+        flip_angle_block = 2 * math.pi * signal[n - 1] * t_block
+        pulse = pp.make_block_pulse(flip_angle=flip_angle_block, duration=t_block,
+                                    phase_offset=ph, system=system)
+
+
+        rfp = SimpleNamespace()
+        rfp.type = 'rf'
+        rfp.signal = pulse.signal
+        rfp.t = pulse.t
+        rfp.freq_offset = pulse.freq_offset
+        rfp.phase_offset = pulse.phase_offset
+        rfp.dead_time = 0
+        rfp.ringdown_time = system.rf_ringdown_time
+        rfp.delay = 0
+
+        RF_train.append(rfp)
+
+        alpha = alpha + 2 * math.pi * signal[n] * (Tseg)
+        del_ph = integrate.trapz([fm_pulse[n], fm_pulse[n - 1]], x=None, dx=t_block)
+
+        ph = ph + del_ph * 2 * math.pi
+
+    print(flip_angle)
+    print(alpha)
+    return RF_train