diff --git a/pypulseq/Sequence/__init__.py b/pypulseq/Sequence/__init__.py
index d2356ad..e69de29 100755
--- a/pypulseq/Sequence/__init__.py
+++ b/pypulseq/Sequence/__init__.py
@@ -1,5 +0,0 @@
-__all__ = ['block',
-           'read_seq',
-           'sequence',
-           'test_report',
-           'write_seq']
diff --git a/pypulseq/Sequence/sequence.py b/pypulseq/Sequence/sequence.py
index be7c33a..c4a073e 100755
--- a/pypulseq/Sequence/sequence.py
+++ b/pypulseq/Sequence/sequence.py
@@ -1,473 +1,472 @@
 import math
 from types import SimpleNamespace
 
-import numpy as np
 import matplotlib as mpl
+import numpy as np
+
 mpl.use('TkAgg')
 from matplotlib import pyplot as plt
 
-import pypulseq.Sequence.test_report
+import pypulseq.sequence.test_report
 import pypulseq.check_timing
-from pypulseq.Sequence import block
-from pypulseq.Sequence.read_seq import read
-from pypulseq.Sequence.write_seq import write
+from pypulseq.sequence import block
+from pypulseq.sequence.read_seq import read
+from pypulseq.sequence.write_seq import write
 from pypulseq.calc_duration import calc_duration
 from pypulseq.calc_rf_center import calc_rf_center
 from pypulseq.decompress_shape import decompress_shape
 from pypulseq.event_lib import EventLibrary
 from pypulseq.opts import Opts
 from pypulseq.points_to_waveform import points_to_waveform
 
 
 class Sequence:
     def __init__(self, system=Opts()):
         self.version_major = 1
         self.version_minor = 2
         self.version_revision = 0
         self.system = system
         self.definitions = dict()
         self.grad_library = EventLibrary()
         self.shape_library = EventLibrary()
         self.rf_library = EventLibrary()
         self.adc_library = EventLibrary()
         self.delay_library = EventLibrary()
         self.block_events = {}
         self.rf_raster_time = self.system.rf_raster_time
         self.grad_raster_time = self.system.grad_raster_time
 
     def __str__(self):
         s = "Sequence:"
         s += "\nshape_library: " + str(self.shape_library)
         s += "\nrf_library: " + str(self.rf_library)
         s += "\ngrad_library: " + str(self.grad_library)
         s += "\nadc_library: " + str(self.adc_library)
         s += "\ndelay_library: " + str(self.delay_library)
         s += "\nrf_raster_time: " + str(self.rf_raster_time)
         s += "\ngrad_raster_time: " + str(self.grad_raster_time)
         s += "\nblock_events: " + str(len(self.block_events))
         return s
 
     def duration(self):
         """
         Get duration of this sequence.
 
         Returns
         -------
         duration : float
             Duration of this sequence in millis.
         num_blocks : int
             Number of blocks in this sequence.
         event_count : int
             Number of events in this sequence.
         """
 
         num_blocks = len(self.block_events)
         event_count = np.zeros(len(self.block_events[1]))
         duration = 0
         for i in range(num_blocks):
             block = self.get_block(i + 1)
             event_count += self.block_events[i + 1] > 0
             duration += calc_duration(block)
 
         return duration, num_blocks, event_count
 
     def check_timing(self):
         """
         Check timing of the events in each block based on grad raster time system limit.
 
         Returns
         -------
         is_ok : bool
             Boolean flag indicating timing errors.
         error_report : str
             Error report in case of timing errors.
         """
         num_blocks = len(self.block_events)
         is_ok = True
         error_report = []
         for i in range(num_blocks):
             block = self.get_block(i + 1)
             event_names = ['rf', 'gx', 'gy', 'gz', 'adc', 'delay']
             ind = [hasattr(block, 'rf'), hasattr(block, 'gx'), hasattr(block, 'gy'), hasattr(block, 'gz'),
                    hasattr(block, 'adc'), hasattr(block, 'delay')]
             ev = [getattr(block, event_names[i]) for i in range(len(event_names)) if ind[i] == 1]
-            res, rep = pypulseq.check_timing.check_timing(self, *ev)
+            res, rep = pypulseq.check_timing(self, *ev)
             is_ok = is_ok and res
             if len(rep) != 0:
                 error_report.append(f'Block: {i} - {rep}\n')
 
         return is_ok, error_report
 
     def test_report(self):
         """
         Analyze the sequence and return a text report.
         """
-        return pypulseq.Sequence.test_report.test_report(self)
+        return pypulseq.sequence.test_report.main(self)
 
     def set_definition(self, key: str, val: str):
         """
-        Sets custom definition to the `Sequence`.
+        Sets custom definition to the `sequence`.
 
         Parameters
         ----------
         key : str
             Definition key.
         val : str
             Definition value.
         """
         self.definitions[key] = val
 
     def get_definition(self, key: str) -> str:
         """
         Retrieves definition identified by `key` from `Sequence`. Returns `None` if no matching definition is found.
 
         Parameters
         ----------
         key : str
             Key of definition to retrieve.
 
         Returns
         -------
         str
             Definition identified by `key` if found, else returns `None`.
         """
         if key in self.definitions:
             return self.definitions[key]
         else:
             return None
 
     def add_block(self, *args):
         """
         Adds event(s) as a block to `Sequence`.
 
         Parameters
         ----------
         args
             Event or list of events to be added as a block to `Sequence`.
         """
         block.add_block(self, len(self.block_events) + 1, *args)
 
     def get_block(self, block_index: int) -> SimpleNamespace:
         """
         Retrieves block of events identified by `block_index` from `Sequence`.
 
         Parameters
         ----------
         block_index : int
             Index of block to be retrieved from `Sequence`.
 
         Returns
         -------
         SimpleNamespace
             Block identified by `block_index`.
         """
         return block.get_block(self, block_index)
 
     def read(self, file_path: str):
         """
         Read `.seq` file from `file_path`.
 
         Parameters
         ----------
         file_path : str
             Path to `.seq` file to be read.
         """
         read(self, file_path)
 
     def write(self, name: str):
         """
         Writes the calling `Sequence` object as a `.seq` file with filename `name`.
 
         Parameters
         ----------
         name :str
             Filename of `.seq` file to be written to disk.
         """
         write(self, name)
 
     def rf_from_lib_data(self, lib_data):
         """
         Construct RF object from `lib_data`.
 
         Parameters
         ----------
         lib_data : list
             RF envelope.
 
         Returns
         -------
         rf : SimpleNamespace
             RF object constructed from lib_data.
         """
         rf = SimpleNamespace()
         rf.type = 'rf'
 
         amplitude, mag_shape, phase_shape = lib_data[0], lib_data[1], lib_data[2]
         shape_data = self.shape_library.data[mag_shape]
         compressed = SimpleNamespace()
         compressed.num_samples = shape_data[0]
         compressed.data = shape_data[1:]
         mag = decompress_shape(compressed)
         shape_data = self.shape_library.data[phase_shape]
         compressed.num_samples = shape_data[0]
         compressed.data = shape_data[1:]
         phase = decompress_shape(compressed)
         rf.signal = 1j * 2 * np.pi * phase
         rf.signal = amplitude * mag * np.exp(rf.signal)
         rf.t = np.arange(1, max(mag.shape) + 1) * self.rf_raster_time
 
-
         if max(lib_data.shape) < 6:
             rf.delay = 0
             rf.freq_offset = lib_data[3]
             rf.phase_offset = lib_data[4]
             lib_data = np.append(lib_data, 0)
         else:
             rf.delay = lib_data[3]
             rf.freq_offset = lib_data[4]
             rf.phase_offset = lib_data[5]
 
-
         if max(lib_data.shape) < 7:
             lib_data = np.append(lib_data, 0)
         rf.dead_time = lib_data[6]
 
         if max(lib_data.shape) < 8:
             lib_data = np.append(lib_data, 0)
         rf.ringdown_time = lib_data[7]
 
         if max(lib_data.shape) < 9:
             lib_data = np.append(lib_data, 0)
 
         use_cases = {1: 'excitation', 2: 'refocusing', 3: 'inversion'}
         if lib_data[8] in use_cases:
             rf.use = use_cases[lib_data[8]]
 
         return rf
 
     def calculate_kspace(self, trajectory_delay: int = 0):
         """
         Calculates the k-space trajectory of the entire pulse sequence.
 
         Parameters
         ----------
         trajectory_delay : int
             Compensation factor in millis to align ADC and gradients in the reconstruction.
 
         Returns
         -------
         k_traj_adc : numpy.ndarray
             K-space trajectory sampled at `t_adc` timepoints.
         k_traj : numpy.ndarray
             K-space trajectory of the entire pulse sequence.
         t_excitation : numpy.ndarray
             Excitation timepoints.
         t_refocusing : numpy.ndarray
             Refocusing timepoints.
         t_adc : numpy.ndarray
             Sampling timepoints.
         """
         c_excitation = 0
         c_refocusing = 0
         c_adc_samples = 0
 
         for i in range(len(self.block_events)):
             block = self.get_block(i + 1)
             if hasattr(block, 'rf'):
                 if not hasattr(block.rf, 'use') or block.rf.use != 'refocusing':
                     c_excitation += 1
                 else:
                     c_refocusing += 1
 
             if hasattr(block, 'adc'):
                 c_adc_samples += int(block.adc.num_samples)
 
         t_excitation = np.zeros(c_excitation)
         t_refocusing = np.zeros(c_refocusing)
         k_time = np.zeros(c_adc_samples)
         current_dur = 0
         c_excitation = 0
         c_refocusing = 0
         k_counter = 0
         traj_recon_delay = trajectory_delay
 
         for i in range(len(self.block_events)):
             block = self.get_block(i + 1)
             if hasattr(block, 'rf'):
                 rf = block.rf
                 rf_center, _ = calc_rf_center(rf)
                 t = rf.delay + rf_center
                 if not hasattr(block.rf, 'use') or block.rf.use != 'refocusing':
                     t_excitation[c_excitation] = current_dur + t
                     c_excitation += 1
                 else:
                     t_refocusing[c_refocusing] = current_dur + t
                     c_refocusing += 1
             if hasattr(block, 'adc'):
                 k_time[k_counter:k_counter + block.adc.num_samples] = np.arange(
                     block.adc.num_samples) * block.adc.dwell + block.adc.delay + current_dur + traj_recon_delay
                 k_counter += block.adc.num_samples
             current_dur += calc_duration(block)
 
         gw = self.gradient_waveforms()
         i_excitation = np.round(t_excitation / self.grad_raster_time)
         i_refocusing = np.round(t_refocusing / self.grad_raster_time)
         k_traj = np.zeros(gw.shape)
         k = [0, 0, 0]
         for i in range(gw.shape[1]):
             k += gw[:, i] * self.grad_raster_time
             k_traj[:, i] = k
             if len(np.where(i_excitation == i + 1)[0]) >= 1:
                 k = 0
                 k_traj[:, i] = np.nan
             if len(np.where(i_refocusing == i + 1)[0]) >= 1:
                 k = -k
 
         k_traj_adc = []
         for i in range(k_traj.shape[0]):
             k_traj_adc.append(np.interp(k_time, np.arange(1, k_traj.shape[1] + 1) * self.grad_raster_time, k_traj[i]))
         k_traj_adc = np.asarray(k_traj_adc)
         t_adc = k_time
 
         return k_traj_adc, k_traj, t_excitation, t_refocusing, t_adc
 
     def gradient_waveforms(self) -> np.ndarray:
         """
         Decompress the entire gradient waveform. Returns an array of shape `gradient_axesxtimepoints`.
         `gradient_axes` is typically 3.
 
         Returns
         -------
         grad_waveforms : numpy.ndarray
             Decompressed gradient waveform.
         """
         duration, num_blocks, _ = self.duration()
 
         wave_length = math.ceil(duration / self.grad_raster_time)
         grad_channels = 3
         grad_waveforms = np.zeros((grad_channels, wave_length))
         grad_channels = ['gx', 'gy', 'gz']
 
         t0 = 0
         t0_n = 0
         for i in range(num_blocks):
             block = self.get_block(i + 1)
             for j in range(len(grad_channels)):
                 if hasattr(block, grad_channels[j]):
                     grad = getattr(block, grad_channels[j])
                     if grad.type == 'grad':
                         nt_start = round((grad.delay + grad.t[0]) / self.grad_raster_time)
                         waveform = grad.waveform
                     else:
                         nt_start = round(grad.delay / self.grad_raster_time)
                         if abs(grad.flat_time) > np.finfo(float).eps:
                             t = np.cumsum([0, grad.rise_time, grad.flat_time, grad.fall_time])
                             trap_form = np.multiply([0, 1, 1, 0], grad.amplitude)
                         else:
                             t = np.cumsum([0, grad.rise_time, grad.fall_time])
                             trap_form = np.multiply([0, 1, 0], grad.amplitude)
 
                         tn = math.floor(t[-1] / self.grad_raster_time)
                         t = np.append(t, t[-1] + self.grad_raster_time)
                         trap_form = np.append(trap_form, 0)
 
                         if abs(grad.amplitude) > np.finfo(float).eps:
                             waveform = points_to_waveform(t, trap_form, self.grad_raster_time)
                         else:
                             waveform = np.zeros(tn + 1)
 
                     if waveform.size != np.sum(np.isfinite(waveform)):
                         raise Warning('Not all elements of the generated waveform are finite')
 
                     """
                     Matlab dynamically resizes arrays during slice assignment operation if assignment is out of bounds
                     Numpy does not
                     Following lines are a workaround
                     """
                     l1, l2 = int(t0_n + nt_start), int(t0_n + nt_start + max(waveform.shape))
                     if l2 > grad_waveforms.shape[1]:
                         grad_waveforms.resize((len(grad_channels), l2))
                     grad_waveforms[j, l1:l2] = waveform
 
             t0 += calc_duration(block)
             t0_n = round(t0 / self.grad_raster_time)
 
         return grad_waveforms
 
     def plot(self, type: str = 'Gradient', time_range=(0, np.inf), time_disp: str = 's'):
         """
         Plot `Sequence`.
 
         Parameters
         ----------
         type : str
             Gradients display type, must be one of either 'Gradient' or 'Kspace'.
         time_range : List
             Time range (x-axis limits) for plotting the sequence. Default is 0 to infinity (entire sequence).
         time_disp : str
             Time display type, must be one of `s`, `ms` or `us`.
         """
 
         valid_plot_types = ['Gradient', 'Kspace']
         valid_time_units = ['s', 'ms', 'us']
         if type not in valid_plot_types:
             raise Exception()
         if time_disp not in valid_time_units:
             raise Exception()
 
         fig1, fig2 = plt.figure(1), plt.figure(2)
         sp11 = fig1.add_subplot(311)
         sp12, sp13 = fig1.add_subplot(312, sharex=sp11), fig1.add_subplot(313, sharex=sp11)
         fig2_sp_list = [fig2.add_subplot(311, sharex=sp11), fig2.add_subplot(312, sharex=sp11),
                         fig2.add_subplot(313, sharex=sp11)]
 
         t_factor_list = [1, 1e3, 1e6]
         t_factor = t_factor_list[valid_time_units.index(time_disp)]
         t0 = 0
         for iB in range(1, len(self.block_events) + 1):
             block = self.get_block(iB)
             is_valid = time_range[0] <= t0 <= time_range[1]
             if is_valid:
                 if hasattr(block, 'adc'):
                     adc = block.adc
                     t = adc.delay + [(x * adc.dwell) for x in range(0, int(adc.num_samples))]
                     sp11.plot((t0 + t), np.zeros(len(t)), 'rx')
                 if hasattr(block, 'rf'):
                     rf = block.rf
                     tc, ic = calc_rf_center(rf)
                     t = rf.t + rf.delay
                     tc = tc + rf.delay
                     sp12.plot(t_factor * (t0 + t), abs(rf.signal))
                     sp13.plot(t_factor * (t0 + t), np.angle(
                         rf.signal * np.exp(1j * rf.phase_offset) * np.exp(1j * 2 * math.pi * rf.t * rf.freq_offset)),
                               t_factor * (t0 + tc), np.angle(rf.signal[ic] * np.exp(1j * rf.phase_offset) * np.exp(
                             1j * 2 * math.pi * rf.t[ic] * rf.freq_offset)), 'xb')
                 grad_channels = ['gx', 'gy', 'gz']
                 for x in range(0, len(grad_channels)):
                     if hasattr(block, grad_channels[x]):
                         grad = getattr(block, grad_channels[x])
                         if grad.type == 'grad':
                             # In place unpacking of grad.t with the starred expression
                             t = grad.delay + [0, *(grad.t + (grad.t[1] - grad.t[0]) / 2),
                                               grad.t[-1] + grad.t[1] - grad.t[0]]
                             waveform = np.array([grad.first, grad.last])
                             waveform = 1e-3 * np.insert(waveform, 1, grad.waveform)
                         else:
                             t = np.cumsum([0, grad.delay, grad.rise_time, grad.flat_time, grad.fall_time])
                             waveform = 1e-3 * grad.amplitude * np.array([0, 0, 1, 1, 0])
                         fig2_sp_list[x].plot(t_factor * (t0 + t), waveform)
             t0 += calc_duration(block)
 
         grad_plot_labels = ['x', 'y', 'z']
         sp11.set_ylabel('ADC')
         sp12.set_ylabel('RF mag (Hz)')
         sp13.set_ylabel('RF phase (rad)')
         [fig2_sp_list[x].set_ylabel(f'G{grad_plot_labels[x]} (kHz/m)') for x in range(3)]
         # Setting display limits
         disp_range = t_factor * np.array([time_range[0], min(t0, time_range[1])])
         sp11.set_xlim(disp_range)
         sp12.set_xlim(disp_range)
         sp13.set_xlim(disp_range)
         [x.set_xlim(disp_range) for x in fig2_sp_list]
 
         plt.show()
diff --git a/pypulseq/Sequence/test_report.py b/pypulseq/Sequence/test_report.py
index fdc353a..8b24880 100644
--- a/pypulseq/Sequence/test_report.py
+++ b/pypulseq/Sequence/test_report.py
@@ -1,141 +1,141 @@
 import numpy as np
 from pypulseq.convert import convert
 
 
-def test_report(self):
+def main(self):
     """
     Analyze the sequence and return a text report.
     """
     # Find RF pulses and list flip angles
     flip_angles_deg = []
     for k in self.rf_library.keys:
         lib_data = self.rf_library.data[k]
         rf = self.rf_from_lib_data(lib_data)
         flip_angles_deg.append(np.abs(np.sum(rf.signal) * rf.t[0] * 360))
 
     flip_angles_deg = np.unique(flip_angles_deg)
 
     # Calculate TE, TR
     duration, num_blocks, event_count = self.duration()
 
     k_traj_adc, k_traj, t_excitation, t_refocusing, t_adc = self.calculate_kspace()
 
     k_abs_adc = np.sqrt(np.sum(np.square(k_traj_adc), axis=0))
     k_abs_echo, index_echo = np.min(k_abs_adc), np.argmin(k_abs_adc)
     t_echo = t_adc[index_echo]
     t_ex_tmp = t_excitation[t_excitation < t_echo]
     TE = t_echo - t_ex_tmp[-1]
 
     if len(t_excitation) < 2:
         TR = duration
     else:
         t_ex_tmp1 = t_excitation[t_excitation > t_echo]
         if len(t_ex_tmp1) == 0:
             TR = t_ex_tmp[-1] - t_ex_tmp[-2]
         else:
             TR = t_ex_tmp1[0] - t_ex_tmp[-1]
 
     # Check sequence dimensionality and spatial resolution
     k_extent = np.max(np.abs(k_traj_adc), axis=1)
     k_scale = np.max(k_extent)
     if k_scale != 0:
         k_bins = 4e6
         k_threshold = k_scale / k_bins
 
         # Detect unused dimensions and delete them
         if np.any(k_extent < k_threshold):
             k_traj_adc = np.delete(k_traj_adc, np.where(k_extent < k_threshold), axis=0)
             k_extent = np.delete(k_extent, np.where(k_extent < k_threshold), axis=0)
 
         # Bin the k-space trajectory to detect repetitions / slices
         k_len = k_traj_adc.shape[1]
         k_repeat = np.zeros(k_len)
 
         k_map = dict()
         for i in range(k_len):
             l = k_bins + np.round(k_traj_adc[:, i] / k_threshold)
             key_string = ('{:.0f} ' * len(l)).format(*l)
             if key_string in k_map:
                 k_repeat[i] = k_map[key_string] + 1
             else:
                 k_repeat[i] = 1
             k_map[key_string] = k_repeat[i]
 
         repeats = np.max(k_repeat)
 
         k_traj_rep1 = k_traj_adc[:, k_repeat == 1]
 
         k_counters = np.zeros(k_traj_rep1.shape)
         dims = k_traj_rep1.shape[0]
         ordering = dict()
         for j in range(dims):
             c = 1
             k_map = dict()
             for i in range(k_traj_rep1.shape[1]):
                 key = round(k_traj_rep1[j, i] / k_threshold)
                 if key in k_map:
                     k_counters[j, i] = k_map[key]
                 else:
                     k_counters[j, i] = c
                     k_map[key] = c
                     c += 1
             ordering[j] = k_map.values()
 
         unique_k_positions = np.max(k_counters, axis=1)
         is_cartesian = np.prod(unique_k_positions) == k_traj_rep1.shape[1]
     else:
         unique_k_positions = 1
 
     gw = self.gradient_waveforms()
     gws = (gw[:, 1:] - gw[:, :-1]) / self.system.grad_raster_time
     ga = np.max(np.abs(gw), axis=1)
     gs = np.max(np.abs(gws), axis=1)
 
     ga_abs = np.max(np.sqrt(np.sum(np.square(gw), axis=0)))
     gs_abs = np.max(np.sqrt(np.sum(np.square(gws), axis=0)))
 
     timing_ok, timing_error_report = self.check_timing()
 
     report = f'Number of blocks: {num_blocks}\n' \
         f'Number of events:\n' \
         f'RF: {event_count[1]:6.0f}\n' \
         f'Gx: {event_count[2]:6.0f}\n' \
         f'Gy: {event_count[3]:6.0f}\n' \
         f'Gz: {event_count[4]:6.0f}\n' \
         f'ADC: {event_count[5]:6.0f}\n' \
         f'Delay: {event_count[0]:6.0f}\n' \
         f'Sequence duration: {duration:.6f} s\n' \
         f'TE: {TE:.6f} s\n' \
         f'TR: {TR:.6f} s\n'
     report += 'Flip angle: ' + ('{:.02f} ' * len(flip_angles_deg)).format(*flip_angles_deg) + 'deg\n'
     report += 'Unique k-space positions (aka cols, rows, etc.): ' + ('{:.0f} ' * len(unique_k_positions)).format(
         *unique_k_positions) + '\n'
 
     if len(np.where(unique_k_positions > 1)):
         report += f'Dimensions: {len(k_extent)}\n'
         report += ('Spatial resolution: {:.02f} mm\n' * len(k_extent)).format(*(0.5 / k_extent * 1e3))
         report += f'Repetitions/slices/contrasts: {repeats}\n'
 
         if is_cartesian:
             report += 'Cartesian encoding trajectory detected\n'
         else:
             report += 'Non-cartesian/irregular encoding trajectory detected (eg: EPI, spiral, radial, etc.)\n'
 
     if timing_ok:
         report += 'Block timing check passed successfully\n'
     else:
         report += f'Block timing check failed. Error listing follows:\n {timing_error_report}'
 
     ga_converted = convert(from_value=ga, from_unit='Hz/m', to_unit='mT/m')
     gs_converted = convert(from_value=gs, from_unit='Hz/m/s', to_unit='T/m/s')
     report += 'Max gradient: ' + ('{:.0f} ' * len(ga)).format(*ga) + 'Hz/m == ' + (
             '{:.02f} ' * len(ga_converted)).format(*ga_converted) + 'mT/m\n'
     report += 'Max slew rate: ' + ('{:.0f} ' * len(gs)).format(*gs) + 'Hz/m/s == ' + (
             '{:.02f} ' * len(gs_converted)).format(*gs_converted) + 'mT/m/s\n'
 
     ga_abs_converted = convert(from_value=ga_abs, from_unit='Hz/m', to_unit='mT/m')
     gs_abs_converted = convert(from_value=gs_abs, from_unit='Hz/m/s', to_unit='T/m/s')
     report += f'Max absolute gradient: {ga_abs:.0f} Hz/m == {ga_abs_converted:.2f} mT/m\n'
     report += f'Max absolute slew rate: {gs_abs:g} Hz/m/s == {gs_abs_converted:.2f} T/m/s'
 
     return report
diff --git a/pypulseq/__init__.py b/pypulseq/__init__.py
index 9b6c390..cf5e3c6 100755
--- a/pypulseq/__init__.py
+++ b/pypulseq/__init__.py
@@ -1,26 +1,26 @@
-from pypulseq.Sequence.sequence import Sequence
+from pypulseq.sequence.sequence import Sequence
 from pypulseq.add_gradients import add_gradients
 from pypulseq.add_ramps import add_ramps
 from pypulseq.align import align
 from pypulseq.calc_duration import calc_duration
 from pypulseq.calc_ramp import calc_ramp
 from pypulseq.calc_rf_center import calc_rf_center
 from pypulseq.check_timing import check_timing
 from pypulseq.compress_shape import compress_shape
 from pypulseq.decompress_shape import decompress_shape
 from pypulseq.event_lib import EventLibrary
 from pypulseq.make_adc import make_adc
 from pypulseq.make_arbitrary_grad import make_arbitrary_grad
 from pypulseq.make_arbitrary_rf import make_trapezoid
 from pypulseq.make_block_pulse import make_block_pulse
 from pypulseq.make_delay import make_delay
 from pypulseq.make_extended_trapezoid import make_arbitrary_grad
 from pypulseq.make_gauss_pulse import make_gauss_pulse
 from pypulseq.make_sinc_pulse import make_sinc_pulse
 from pypulseq.make_trap_pulse import make_trapezoid
 from pypulseq.opts import Opts
 from pypulseq.points_to_waveform import points_to_waveform
 from pypulseq.split_gradient_at import split_gradient_at
 from pypulseq.utils.SAR.SAR_calc import calc_SAR
 from pypulseq.utils.k2g import k2g
 from pypulseq.utils.vds_2d import vds_2d
diff --git a/pypulseq/calc_duration.py b/pypulseq/calc_duration.py
index 1edeb9d..78765f3 100755
--- a/pypulseq/calc_duration.py
+++ b/pypulseq/calc_duration.py
@@ -1,41 +1,41 @@
 from types import SimpleNamespace
 
 
 def calc_duration(*events: list) -> float:
     """
     Calculate the cumulative duration of Events.
 
     Parameters
     ----------
     events : list
         List of `SimpleNamespace` events. Can also be a list containing a single block (see
-        `pypulseq.Sequence.sequence.plot()`).
+        `pypulseq.sequence.sequence.Sequence.plot()`).
 
     Returns
     -------
     duration : float
         The cumulative duration of the pulse events in `events`.
     """
 
     for e in events:
         if isinstance(e, SimpleNamespace) and isinstance(getattr(e, list(e.__dict__.keys())[0]), SimpleNamespace):
             events = list(e.__dict__.values())
             break
 
     duration = 0
     for event in events:
         if not isinstance(event, SimpleNamespace):
             raise TypeError("input(s) should be of type SimpleNamespace")
 
         if event.type == 'delay':
             duration = max(duration, event.delay)
         elif event.type == 'rf':
             duration = max(duration, event.delay + event.t[-1])
         elif event.type == 'grad':
             duration = max(duration, event.t[-1] + event.t[1] - event.t[0] + event.delay)
         elif event.type == 'adc':
             duration = max(duration, event.delay + event.num_samples * event.dwell + event.dead_time)
         elif event.type == 'trap':
             duration = max(duration, event.delay + event.rise_time + event.flat_time + event.fall_time)
 
     return round(duration, ndigits=6)
diff --git a/pypulseq/seq_examples/scripts/write_SE.py b/pypulseq/seq_examples/scripts/write_SE.py
deleted file mode 100644
index e69de29..0000000
diff --git a/pypulseq/seq_examples/scripts/write_epi.py b/pypulseq/seq_examples/scripts/write_epi.py
index b7134fe..8742ee2 100644
--- a/pypulseq/seq_examples/scripts/write_epi.py
+++ b/pypulseq/seq_examples/scripts/write_epi.py
@@ -1,64 +1,64 @@
 import math
 
 import numpy as np
 
-from pypulseq.Sequence.sequence import Sequence
+from pypulseq.sequence.sequence import Sequence
 from pypulseq.make_adc import make_adc
 from pypulseq.make_sinc_pulse import make_sinc_pulse
 from pypulseq.make_trap_pulse import make_trapezoid
 from pypulseq.opts import Opts
 import matplotlib.pyplot as plt
 
 seq = Sequence()
 fov = 220e-3
 Nx = 64
 Ny = 64
 slice_thickness = 3e-3
 n_slices = 1
 
 system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s', rf_ringdown_time=30e-6,
               rf_dead_time=100e-6)
 
 rf, gz, _ = make_sinc_pulse(flip_angle=math.pi / 2, system=system, duration=3e-3, slice_thickness=slice_thickness,
                             apodization=0.5, time_bw_product=4)
 
 delta_k = 1 / fov
 k_width = Nx * delta_k
 dwell_time = 4e-6
 readout_time = Nx * dwell_time
 flat_time = math.ceil(readout_time * 1e5) * 1e-5
 gx = make_trapezoid(channel='x', system=system, amplitude=k_width / readout_time, flat_time=flat_time)
 adc = make_adc(num_samples=Nx, duration=readout_time,
                delay=gx.rise_time + flat_time / 2 - (readout_time - dwell_time) / 2)
 
 pre_time = 8e-4
 gx_pre = make_trapezoid(channel='x', system=system, area=-gx.area / 2, duration=pre_time)
 gz_reph = make_trapezoid(channel='z', system=system, area=-gz.area / 2, duration=pre_time)
 gy_pre = make_trapezoid(channel='y', system=system, area=-Ny / 2 * delta_k, duration=pre_time)
 
 dur = math.ceil(2 * math.sqrt(delta_k / system.max_slew) / 10e-6) * 10e-6
 gy = make_trapezoid(channel='y', system=system, area=delta_k, duration=dur)
 
 for s in range(n_slices):
     rf.freq_offset = gz.amplitude * slice_thickness * (s - (n_slices - 1) / 2)
     seq.add_block(rf, gz)
     seq.add_block(gx_pre, gy_pre, gz_reph)
     for i in range(Ny):
         seq.add_block(gx, adc)
         seq.add_block(gy)
         gx.amplitude = -gx.amplitude
 
 seq.plot()
 
 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.plot(time_axis, ktraj.T)
 plt.plot(t_adc, ktraj_adc[0, :], '.')
 plt.figure()
 plt.plot(ktraj[0, :], ktraj[1, :], 'b')
 plt.axis('equal')
 plt.plot(ktraj_adc[0, :], ktraj_adc[1, :], 'r.')
 plt.show()
 
 seq.write('epi_pypulseq.seq')
diff --git a/pypulseq/seq_examples/scripts/write_epi_rs.py b/pypulseq/seq_examples/scripts/write_epi_rs.py
index ecadde2..e3e7d53 100644
--- a/pypulseq/seq_examples/scripts/write_epi_rs.py
+++ b/pypulseq/seq_examples/scripts/write_epi_rs.py
@@ -1,116 +1,116 @@
 import math
 from warnings import warn
 
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pypulseq.Sequence.sequence import Sequence
+from pypulseq.sequence.sequence import Sequence
 from pypulseq.add_gradients import add_gradients
 from pypulseq.align import align
 from pypulseq.calc_duration import calc_duration
 from pypulseq.make_adc import make_adc
 from pypulseq.make_gauss_pulse import make_gauss_pulse
 from pypulseq.make_sinc_pulse import make_sinc_pulse
 from pypulseq.make_trap_pulse import make_trapezoid
 from pypulseq.opts import Opts
 from pypulseq.split_gradient_at import split_gradient_at
 
 seq = Sequence()
 fov = 250e-3
 Nx = 64
 Ny = 64
 slice_thickness = 3e-3
 n_slices = 1
 
 pe_enable = 1
 
 system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s', rf_ringdown_time=30e-6,
               rf_dead_time=100e-6)
 
 b0 = 2.89
 sat_ppm = -3.45
 sat_freq = sat_ppm * 1e-6 * b0 * system.gamma
 rf_fs, _, _ = make_gauss_pulse(flip_angle=110 * math.pi / 180, system=system, duration=8e-3, bandwidth=abs(sat_freq),
                                freq_offset=sat_freq)
 gz_fs = make_trapezoid(channel='z', system=system, delay=calc_duration(rf_fs), area=1 / 1e-4)
 
 rf, gz, gz_reph = make_sinc_pulse(flip_angle=math.pi / 2, system=system, duration=3e-3, slice_thickness=slice_thickness,
                                   apodization=0.5, time_bw_product=4)
 
 delta_k = 1 / fov
 k_width = Nx * delta_k
 readout_time = 4.2e-4
 
 blip_dur = math.ceil(2 * math.sqrt(delta_k / system.max_slew) / 10e-6 / 2) * 10e-6 * 2
 gy = make_trapezoid(channel='y', system=system, area=delta_k, duration=blip_dur)
 
 extra_area = blip_dur / 2 * blip_dur / 2 * system.max_slew
 gx = make_trapezoid(channel='x', system=system, area=k_width + extra_area, duration=readout_time + blip_dur)
 actual_area = gx.area - gx.amplitude / gx.rise_time * blip_dur / 2 * blip_dur / 2 / 2 - gx.amplitude / gx.fall_time * blip_dur / 2 * blip_dur / 2 / 2
 gx.amplitude = gx.amplitude / actual_area * k_width
 gx.area = gx.amplitude * (gx.flat_time + gx.rise_time / 2 + gx.fall_time / 2)
 gx.flat_area = gx.amplitude * gx.flat_time
 
 adc_dwell_nyquist = delta_k / gx.amplitude
 adc_dwell = math.floor(adc_dwell_nyquist * 1e7) * 1e-7
 adc_samples = math.floor(readout_time / adc_dwell / 4) * 4
 
 adc = make_adc(num_samples=adc_samples, dwell=adc_dwell, delay=blip_dur / 2)
 
 time_to_center = adc.dwell * (adc_samples - 1) / 2
 adc.delay = round((gx.rise_time + gx.flat_time / 2 - time_to_center) * 1e6) * 1e-6
 
 gy_parts = split_gradient_at(grad=gy, time_point=blip_dur / 2, system=system)
 gy_blipup, gy_blipdown, _ = align('right', gy_parts[0], 'left', gy_parts[1], gx)
 gy_blipdownup = add_gradients([gy_blipdown, gy_blipup], system=system)
 
 gy_blipup.waveform = gy_blipup.waveform * pe_enable
 gy_blipdown.waveform = gy_blipdown.waveform * pe_enable
 gy_blipdownup.waveform = gy_blipdownup.waveform * pe_enable
 
 gx_pre = make_trapezoid(channel='x', system=system, area=-gx.area / 2)
 gy_pre = make_trapezoid(channel='y', system=system, area=-Ny / 2 * delta_k)
 gx_pre, gy_pre, gz_reph = align('right', gx_pre, 'left', gy_pre, gz_reph)
 gy_pre = make_trapezoid(channel='y', system=system, area=gy_pre.area, duration=calc_duration(gx_pre, gy_pre, gz_reph))
 gy_pre.amplitude = gy_pre.amplitude * pe_enable
 
 for s in range(n_slices):
     seq.add_block(rf_fs, gz_fs)
     rf.freq_offset = gz.amplitude * slice_thickness * (s - (n_slices - 1) / 2)
     seq.add_block(rf, gz)
     seq.add_block(gx_pre, gy_pre, gz_reph)
     for i in range(1, Ny + 1):
         if i == 1:
             seq.add_block(gx, gy_blipup, adc)
         elif i == Ny:
             seq.add_block(gx, gy_blipdown, adc)
         else:
             seq.add_block(gx, gy_blipdownup, adc)
         gx.amplitude = -gx.amplitude
 
 ok, error_report = seq.check_timing()
 
 if ok:
     print('Timing check passed succesfully.')
 else:
     warn('Timing check failed! Error listing follows:\n')
     print(error_report)
     print('\n')
 
 seq.plot()
 
 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.plot(time_axis, ktraj.T)
 plt.plot(t_adc, ktraj_adc[0], '.')
 plt.figure()
 plt.plot(ktraj_adc[0], ktraj_adc[1], 'b')
 plt.axis('equal')
 plt.plot(ktraj_adc[0], ktraj_adc[1], 'r.')
 plt.show()
 
 seq.set_definition('FOV', np.array([fov, fov, slice_thickness]) * 1e3)
 seq.set_definition('Name', 'epi')
 
 seq.write('epi_rs_pypulseq.seq')
diff --git a/pypulseq/seq_examples/scripts/write_gre.py b/pypulseq/seq_examples/scripts/write_gre.py
index b6089bd..a7a657b 100644
--- a/pypulseq/seq_examples/scripts/write_gre.py
+++ b/pypulseq/seq_examples/scripts/write_gre.py
@@ -1,69 +1,69 @@
 import math
 
 import numpy as np
 
-from pypulseq.Sequence.sequence import Sequence
+from pypulseq.sequence.sequence import Sequence
 from pypulseq.calc_duration import calc_duration
 from pypulseq.make_adc import make_adc
 from pypulseq.make_delay import make_delay
 from pypulseq.make_sinc_pulse import make_sinc_pulse
 from pypulseq.make_trap_pulse import make_trapezoid
 from pypulseq.opts import Opts
 
 seq = Sequence()
 fov = 250e-3
 Nx = 256
 Ny = 256
 alpha = 10
 slice_thickness = 3e-3
 TE = np.array([7.38, 9.84]) * 1e-3
 TR = 100e-3
 
 rf_spoiling_inc = 117
 
 sys = 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, gz, gzr = make_sinc_pulse(flip_angle=alpha * math.pi / 180, duration=4e-3, slice_thickness=slice_thickness,
                               apodization=0.5, time_bw_product=4, system=sys)
 
 delta_k = 1 / fov
 gx = make_trapezoid(channel='x', flat_area=Nx * delta_k, flat_time=6.4e-3, system=sys)
 adc = make_adc(num_samples=Nx, duration=gx.flat_time, delay=gx.rise_time, system=sys)
 gx_pre = make_trapezoid(channel='x', area=-gx.area / 2, duration=2e-3, system=sys)
 gz_reph = make_trapezoid(channel='z', area=-gz.area / 2, duration=2e-3, system=sys)
 phase_areas = (np.arange(Ny) - Ny / 2) * delta_k
 
 gx_spoil = make_trapezoid(channel='x', area=2 * Nx * delta_k, system=sys)
 gz_spoil = make_trapezoid(channel='z', area=4 / slice_thickness, system=sys)
 
 delay_TE = np.ceil((TE - calc_duration(gx_pre) - gz.fall_time - gz.flat_time / 2 - calc_duration(
     gx) / 2) / seq.grad_raster_time) * seq.grad_raster_time
 delay_TR = np.ceil((TR - calc_duration(gx_pre) - calc_duration(gz) - calc_duration(
     gx) - delay_TE) / seq.grad_raster_time) * seq.grad_raster_time
 
 assert np.all(delay_TR >= calc_duration(gx_spoil, gz_spoil))
 
 rf_phase = 0
 rf_inc = 0
 
 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]
 
         seq.add_block(rf, gz)
         gy_pre = make_trapezoid(channel='y', area=phase_areas[i], duration=2e-3, system=sys)
         seq.add_block(gx_pre, gy_pre, gz_reph)
         seq.add_block(make_delay(delay_TE[j]))
         seq.add_block(gx, adc)
         gy_pre.amplitude = -gy_pre.amplitude
         seq.add_block(make_delay(delay_TR[j]), gx_spoil, gy_pre, gz_spoil)
 
 report = seq.test_report()
 print(report)
 seq.calculate_kspace()
 seq.plot()
 seq.write('gre_pypulseq.seq')
diff --git a/pypulseq/seq_examples/scripts/write_haste.py b/pypulseq/seq_examples/scripts/write_haste.py
index 152b119..6fb7ed3 100644
--- a/pypulseq/seq_examples/scripts/write_haste.py
+++ b/pypulseq/seq_examples/scripts/write_haste.py
@@ -1,180 +1,180 @@
 import math
 import warnings
 
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pypulseq.Sequence.sequence import Sequence
+from pypulseq.sequence.sequence import Sequence
 from pypulseq.calc_rf_center import calc_rf_center
 from pypulseq.make_adc import make_adc
 from pypulseq.make_delay import make_delay
 from pypulseq.make_extended_trapezoid import make_extended_trapezoid
 from pypulseq.make_sinc_pulse import make_sinc_pulse
 from pypulseq.make_trap_pulse import make_trapezoid
 from pypulseq.opts import Opts
 
 dG = 250e-6
 system = Opts(max_grad=30, grad_unit='mT/m', max_slew=170, slew_unit='T/m/s', rf_ringdown_time=100e-6,
               rf_dead_time=100e-6, adc_dead_time=10e-6)
 
 seq = Sequence(system=system)
 
 fov = 256e-3
 Ny_pre = 8
 Nx, Ny = 128, 128
 n_echo = int(Ny / 2 + Ny_pre)
 n_slices = 1
 rf_flip = 180
 if isinstance(rf_flip, int):
     rf_flip = np.zeros(n_echo) + rf_flip
 slice_thickness = 5e-3
 TE = 12e-3
 TR = 2000e-3
 TE_eff = 60e-3
 k0 = round(TE_eff / TE)
 PE_type = 'linear'
 
 readout_time = 6.4e-3 + 2 * system.adc_dead_time
 t_ex = 2.5e-3
 t_ex_wd = t_ex + system.rf_ringdown_time + system.rf_dead_time
 t_ref = 2e-3
 tf_ref_wd = t_ref + system.rf_ringdown_time + system.rf_dead_time
 t_sp = 0.5 * (TE - readout_time - tf_ref_wd)
 t_sp_ex = 0.5 * (TE - t_ex_wd - tf_ref_wd)
 fspR = 1.0
 fspS = 0.5
 
 rfex_phase = math.pi / 2
 rfref_phase = 0
 
 flipex = 90 * math.pi / 180
 rfex, gz, _ = make_sinc_pulse(flip_angle=flipex, system=system, duration=t_ex, slice_thickness=slice_thickness,
                               apodization=0.5, time_bw_product=4, phase_offset=rfex_phase)
 GS_ex = make_trapezoid(channel='z', system=system, amplitude=gz.amplitude, flat_time=t_ex_wd, rise_time=dG)
 
 flipref = rf_flip[0] * math.pi / 180
 rfref, gz, _ = make_sinc_pulse(flip_angle=flipref, system=system, duration=t_ref, slice_thickness=slice_thickness,
                                apodization=0.5, time_bw_product=4, phase_offset=rfref_phase, use='refocusing')
 GS_ref = make_trapezoid(channel='z', system=system, amplitude=GS_ex.amplitude, flat_time=tf_ref_wd, rise_time=dG)
 
 AGS_ex = GS_ex.area / 2
 GS_spr = make_trapezoid(channel='z', system=system, area=AGS_ex * (1 + fspS), duration=t_sp, rise_time=dG)
 GS_spex = make_trapezoid(channel='z', system=system, area=AGS_ex * fspS, duration=t_sp_ex, rise_time=dG)
 
 delta_k = 1 / fov
 k_width = Nx * delta_k
 
 GR_acq = make_trapezoid(channel='x', system=system, flat_area=k_width, flat_time=readout_time, rise_time=dG)
 adc = make_adc(num_samples=Nx, duration=GR_acq.flat_time - 40e-6, delay=20e-6)
 GR_spr = make_trapezoid(channel='x', system=system, area=GR_acq.area * fspR, duration=t_sp, rise_time=dG)
 GR_spex = make_trapezoid(channel='x', system=system, area=GR_acq.area * (1 + fspR), duration=t_sp_ex, rise_time=dG)
 
 AGR_spr = GR_spr.area
 AGR_preph = GR_acq.area / 2 + AGR_spr
 GR_preph = make_trapezoid(channel='x', system=system, area=AGR_preph, duration=t_sp_ex, rise_time=dG)
 
 n_ex = 1
 PE_order = np.arange(-Ny_pre, Ny + 1).T
 phase_areas = PE_order * delta_k
 
 GS1_times = [0, GS_ex.rise_time]
 GS1_amp = [0, GS_ex.amplitude]
 GS1 = make_extended_trapezoid(channel='z', times=GS1_times, amplitudes=GS1_amp)
 
 GS2_times = [0, GS_ex.flat_time]
 GS2_amp = [GS_ex.amplitude, GS_ex.amplitude]
 GS2 = make_extended_trapezoid(channel='z', times=GS2_times, amplitudes=GS2_amp)
 
 GS3_times = [0, GS_spex.rise_time, GS_spex.rise_time + GS_spex.flat_time,
              GS_spex.rise_time + GS_spex.flat_time + GS_spex.fall_time]
 GS3_amp = [GS_ex.amplitude, GS_spex.amplitude, GS_spex.amplitude, GS_ref.amplitude]
 GS3 = make_extended_trapezoid(channel='z', times=GS3_times, amplitudes=GS3_amp)
 
 GS4_times = [0, GS_ref.flat_time]
 GS4_amp = [GS_ref.amplitude, GS_ref.amplitude]
 GS4 = make_extended_trapezoid(channel='z', times=GS4_times, amplitudes=GS4_amp)
 
 GS5_times = [0, GS_spr.rise_time, GS_spr.rise_time + GS_spr.flat_time,
              GS_spr.rise_time + GS_spr.flat_time + GS_spr.fall_time]
 GS5_amp = [GS_ref.amplitude, GS_spr.amplitude, GS_spr.amplitude, 0]
 GS5 = make_extended_trapezoid(channel='z', times=GS5_times, amplitudes=GS5_amp)
 
 GS7_times = [0, GS_spr.rise_time, GS_spr.rise_time + GS_spr.flat_time,
              GS_spr.rise_time + GS_spr.flat_time + GS_spr.fall_time]
 GS7_amp = [0, GS_spr.amplitude, GS_spr.amplitude, GS_ref.amplitude]
 GS7 = make_extended_trapezoid(channel='z', times=GS7_times, amplitudes=GS7_amp)
 
 GR3 = GR_preph
 
 GR5_times = [0, GR_spr.rise_time, GR_spr.rise_time + GR_spr.flat_time,
              GR_spr.rise_time + GR_spr.flat_time + GR_spr.fall_time]
 GR5_amp = [0, GR_spr.amplitude, GR_spr.amplitude, GR_acq.amplitude]
 GR5 = make_extended_trapezoid(channel='x', times=GR5_times, amplitudes=GR5_amp)
 
 GR6_times = [0, readout_time]
 GR6_amp = [GR_acq.amplitude, GR_acq.amplitude]
 GR6 = make_extended_trapezoid(channel='x', times=GR6_times, amplitudes=GR6_amp)
 
 GR7_times = [0, GR_spr.rise_time, GR_spr.rise_time + GR_spr.flat_time,
              GR_spr.rise_time + GR_spr.flat_time + GR_spr.fall_time]
 GR7_amp = [GR_acq.amplitude, GR_spr.amplitude, GR_spr.amplitude, 0]
 GR7 = make_extended_trapezoid(channel='x', times=GR7_times, amplitudes=GR7_amp)
 
 tex = GS1.t[-1] + GS2.t[-1] + GS3.t[-1]
 tref = GS4.t[-1] + GS5.t[-1] + GS7.t[-1] + readout_time
 tend = GS4.t[-1] + GS5.t[-1]
 TE_train = tex + n_echo * tref + tend
 TR_fill = (TR - n_slices * TE_train) / n_slices
 
 TR_fill = system.grad_raster_time * round(TR_fill / system.grad_raster_time)
 if TR_fill < 0:
     TR_fill = 1e-3
     warnings.warn(f'TR too short, adapted to include all slices to: {1000 * n_slices * (TE_train + TR_fill)} ms')
 else:
     print(f'TR fill: {1000 * TR_fill} ms')
 delay_TR = make_delay(TR_fill)
 delay_end = make_delay(5)
 
 for k_ex in range(n_ex):
     for s in range(n_slices):
         rfex.freq_offset = GS_ex.amplitude * slice_thickness * (s - (n_slices - 1) / 2)
         rfref.freq_offset = GS_ref.amplitude * slice_thickness * (s - (n_slices - 1) / 2)
         rfex.phase_offset = rfex_phase - 2 * math.pi * rfex.freq_offset * calc_rf_center(rfex)[0]
         rfref.phase_offset = rfref_phase - 2 * math.pi * rfref.freq_offset * calc_rf_center(rfref)[0]
 
         seq.add_block(GS1)
         seq.add_block(GS2, rfex)
         seq.add_block(GS3, GR3)
 
         for k_ech in range(n_echo):
             if k_ex >= 0:
                 phase_area = phase_areas[k_ech]
             else:
                 phase_area = 0
 
             GP_pre = make_trapezoid(channel='y', system=system, area=phase_area, duration=t_sp, rise_time=dG)
             GP_rew = make_trapezoid(channel='y', system=system, area=-phase_area, duration=t_sp, rise_time=dG)
 
             seq.add_block(GS4, rfref)
             seq.add_block(GS5, GR5, GP_pre)
 
             if k_ex >= 0:
                 seq.add_block(GR6, adc)
             else:
                 seq.add_block(GR6)
 
             seq.add_block(GS7, GR7, GP_rew)
 
         seq.add_block(GS4)
         seq.add_block(GS5)
         seq.add_block(delay_TR)
 
 seq.add_block(delay_end)
 
 ktraj_adc, ktraj, t_excitation, t_refocusing, _ = seq.calculate_kspace()
 plt.plot(ktraj.T)
 plt.figure()
 plt.plot(ktraj[0], ktraj[1], 'b', ktraj_adc[0], ktraj_adc[1], 'r.')
 plt.show()
 
 seq.write('haste_pypulseq.seq')
diff --git a/pypulseq/seq_examples/scripts/write_tse.py b/pypulseq/seq_examples/scripts/write_tse.py
index fa27eea..d17bcd4 100644
--- a/pypulseq/seq_examples/scripts/write_tse.py
+++ b/pypulseq/seq_examples/scripts/write_tse.py
@@ -1,175 +1,175 @@
 import math
 import warnings
 
 import numpy as np
 
-from pypulseq.Sequence.sequence import Sequence
+from pypulseq.sequence.sequence import Sequence
 from pypulseq.calc_rf_center import calc_rf_center
 from pypulseq.make_adc import make_adc
 from pypulseq.make_delay import make_delay
 from pypulseq.make_extended_trapezoid import make_extended_trapezoid
 from pypulseq.make_sinc_pulse import make_sinc_pulse
 from pypulseq.make_trap_pulse import make_trapezoid
 from pypulseq.opts import Opts
 
 dG = 250e-6
 
 system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s', rf_ringdown_time=100e-6,
               rf_dead_time=100e-6, adc_dead_time=10e-6)
 seq = Sequence(system)
 
 fov = 256e-3
 Nx, Ny = 128, 128
 n_echo = 16
 n_slices = 1
 rf_flip = 180
 
 if isinstance(rf_flip, int):
     rf_flip = np.zeros(n_echo) + rf_flip
 slice_thickness = 5e-3
 TE = 12e-3
 TR = 2000e-3
 TE_eff = 60e-3
 k0 = round(TE_eff / TE)
 pe_type = 'linear'
 
 readout_time = 6.4e-3 + 2 * system.adc_dead_time
 t_ex = 2.5e-3
 t_exwd = t_ex + system.rf_ringdown_time + system.rf_dead_time
 t_ref = 2e-3
 t_refwd = t_ref + system.rf_ringdown_time + system.rf_dead_time
 t_sp = 0.5 * (TE - readout_time - t_refwd)
 t_spex = 0.5 * (TE - t_exwd - t_refwd)
 fsp_r = 1
 fsp_s = 0.5
 
 rf_ex_phase = np.pi / 2
 rf_ref_phase = 0
 
 flip_ex = 90 * np.pi / 180
 rf_ex, gz, _ = make_sinc_pulse(flip_angle=flip_ex, system=system, duration=t_ex, slice_thickness=slice_thickness,
                                apodization=0.5, time_bw_product=4, phase_offset=rf_ex_phase)
 gs_ex = make_trapezoid(channel='z', system=system, amplitude=gz.amplitude, flat_time=t_exwd, rise_time=dG)
 
 flip_ref = rf_flip[0] * np.pi / 180
 rf_ref, gz, _ = make_sinc_pulse(flip_angle=flip_ref, system=system, duration=t_ref, slice_thickness=slice_thickness,
                                 apodization=0.5, time_bw_product=4, phase_offset=rf_ref_phase, use='refocusing')
 gs_ref = make_trapezoid(channel='z', system=system, amplitude=gs_ex.amplitude, flat_time=t_refwd, rise_time=dG)
 
 ags_ex = gs_ex.area / 2
 gs_spr = make_trapezoid(channel='z', system=system, area=ags_ex * (1 + fsp_s), duration=t_sp, rise_time=dG)
 gs_spex = make_trapezoid(channel='z', system=system, area=ags_ex * fsp_s, duration=t_spex, rise_time=dG)
 
 delta_k = 1 / fov
 k_width = Nx * delta_k
 
 gr_acq = make_trapezoid(channel='x', system=system, flat_area=k_width, flat_time=readout_time, rise_time=dG)
 adc = make_adc(num_samples=Nx, duration=gr_acq.flat_time - 40e-6, delay=20e-6)
 gr_spr = make_trapezoid(channel='x', system=system, area=gr_acq.area * fsp_r, duration=t_sp, rise_time=dG)
 gr_spex = make_trapezoid(channel='x', system=system, area=gr_acq.area * (1 + fsp_r), duration=t_spex, rise_time=dG)
 
 agr_spr = gr_spr.area
 agr_preph = gr_acq.area / 2 + agr_spr
 gr_preph = make_trapezoid(channel='x', system=system, area=agr_preph, duration=t_spex, rise_time=dG)
 
 n_ex = math.floor(Ny / n_echo)
 pe_steps = np.arange(1, n_echo * n_ex + 1) - 0.5 * n_echo * n_ex - 1
 if divmod(n_echo, 2)[1] == 0:
     pe_steps = np.roll(pe_steps, -round(n_ex / 2))
 pe_order = pe_steps.reshape((n_ex, n_echo), order='F').T
 phase_areas = pe_order * delta_k
 
 # Split gradients and recombine into blocks
 gs1_times = [0, gs_ex.rise_time]
 gs1_amp = [0, gs_ex.amplitude]
 gs1 = make_extended_trapezoid(channel='z', times=gs1_times, amplitudes=gs1_amp)
 
 gs2_times = [0, gs_ex.flat_time]
 gs2_amp = [gs_ex.amplitude, gs_ex.amplitude]
 gs2 = make_extended_trapezoid(channel='z', times=gs2_times, amplitudes=gs2_amp)
 
 gs3_times = [0, gs_spex.rise_time, gs_spex.rise_time + gs_spex.flat_time,
              gs_spex.rise_time + gs_spex.flat_time + gs_spex.fall_time]
 gs3_amp = [gs_ex.amplitude, gs_spex.amplitude, gs_spex.amplitude, gs_ref.amplitude]
 gs3 = make_extended_trapezoid(channel='z', times=gs3_times, amplitudes=gs3_amp)
 
 gs4_times = [0, gs_ref.flat_time]
 gs4_amp = [gs_ref.amplitude, gs_ref.amplitude]
 gs4 = make_extended_trapezoid(channel='z', times=gs4_times, amplitudes=gs4_amp)
 
 gs5_times = [0, gs_spr.rise_time, gs_spr.rise_time + gs_spr.flat_time,
              gs_spr.rise_time + gs_spr.flat_time + gs_spr.fall_time]
 gs5_amp = [gs_ref.amplitude, gs_spr.amplitude, gs_spr.amplitude, 0]
 gs5 = make_extended_trapezoid(channel='z', times=gs5_times, amplitudes=gs5_amp)
 
 gs7_times = [0, gs_spr.rise_time, gs_spr.rise_time + gs_spr.flat_time,
              gs_spr.rise_time + gs_spr.flat_time + gs_spr.fall_time]
 gs7_amp = [0, gs_spr.amplitude, gs_spr.amplitude, gs_ref.amplitude]
 gs7 = make_extended_trapezoid(channel='z', times=gs7_times, amplitudes=gs7_amp)
 
 gr3 = gr_preph
 
 gr5_times = [0, gr_spr.rise_time, gr_spr.rise_time + gr_spr.flat_time,
              gr_spr.rise_time + gr_spr.flat_time + gr_spr.fall_time]
 gr5_amp = [0, gr_spr.amplitude, gr_spr.amplitude, gr_acq.amplitude]
 gr5 = make_extended_trapezoid(channel='x', times=gr5_times, amplitudes=gr5_amp)
 
 gr6_times = [0, readout_time]
 gr6_amp = [gr_acq.amplitude, gr_acq.amplitude]
 gr6 = make_extended_trapezoid(channel='x', times=gr6_times, amplitudes=gr6_amp)
 
 gr7_times = [0, gr_spr.rise_time, gr_spr.rise_time + gr_spr.flat_time,
              gr_spr.rise_time + gr_spr.flat_time + gr_spr.fall_time]
 gr7_amp = [gr_acq.amplitude, gr_spr.amplitude, gr_spr.amplitude, 0]
 gr7 = make_extended_trapezoid(channel='x', times=gr7_times, amplitudes=gr7_amp)
 
 t_ex = gs1.t[-1] + gs2.t[-1] + gs3.t[-1]
 t_ref = gs4.t[-1] + gs5.t[-1] + gs7.t[-1] + readout_time
 t_end = gs4.t[-1] + gs5.t[-1]
 TE_train = t_ex + n_echo * t_ref + t_end
 TR_fill = (TR - n_slices * TE_train) / n_slices
 TR_fill = system.grad_raster_time * round(TR_fill / system.grad_raster_time)
 if TR_fill < 0:
     TR_fill = 1e-3
     warnings.warn(f'TR too short, adapted to include all slices to: {1000 * n_slices * (TE_train + TR_fill)} ms')
 else:
     print(f'TR fill: {1000 * TR_fill} ms')
 delay_TR = make_delay(TR_fill)
 
 for k_ex in range(n_ex + 1):
     for s in range(n_slices):
         rf_ex.freq_offset = gs_ex.amplitude * slice_thickness * (s - (n_slices - 1) / 2)
         rf_ref.freq_offset = gs_ref.amplitude * slice_thickness * (s - (n_slices - 1) / 2)
         rf_ex.phase_offset = rf_ex_phase - 2 * np.pi * rf_ex.freq_offset * calc_rf_center(rf_ex)[0]
         rf_ref.phase_offset = rf_ref_phase - 2 * np.pi * rf_ref.freq_offset * calc_rf_center(rf_ref)[0]
 
         seq.add_block(gs1)
         seq.add_block(gs2, rf_ex)
         seq.add_block(gs3, gr3)
 
         for k_echo in range(n_echo):
             if k_ex > 0:
                 phase_area = phase_areas[k_echo, k_ex - 1]
             else:
                 phase_area = 0.0  # 0.0 and not 0 because -phase_area should successfully result in negative zero
 
             gp_pre = make_trapezoid(channel='y', system=system, area=phase_area, duration=t_sp, rise_time=dG)
             gp_rew = make_trapezoid(channel='y', system=system, area=-phase_area, duration=t_sp, rise_time=dG)
             seq.add_block(gs4, rf_ref)
             seq.add_block(gs5, gr5, gp_pre)
             if k_ex > 0:
                 seq.add_block(gr6, adc)
             else:
                 seq.add_block(gr6)
 
             seq.add_block(gs7, gr7, gp_rew)
 
         seq.add_block(gs4)
         seq.add_block(gs5)
         seq.add_block(delay_TR)
 
 report = seq.test_report()
 # print(report)
 # seq.calculate_kspace()
 # seq.plot()
 seq.write('tse_pypulseq.seq')
diff --git a/pypulseq/utils/SAR/SAR_calc.py b/pypulseq/utils/SAR/SAR_calc.py
index 1292171..2d88355 100644
--- a/pypulseq/utils/SAR/SAR_calc.py
+++ b/pypulseq/utils/SAR/SAR_calc.py
@@ -1,252 +1,252 @@
 # Copyright of the Board of Trustees of Columbia University in the City of New York
 
 from pathlib import Path
 
 import matplotlib.pyplot as plt
 import numpy as np
 import numpy.matlib
 import scipy.io as sio
 from scipy import interpolate
 
-from pypulseq.Sequence.sequence import Sequence
+from pypulseq.sequence.sequence import Sequence
 from pypulseq.calc_duration import calc_duration
 
 
 def __calc_SAR(Q, I):
     """
         Compute the SAR output for a given Q matrix and I current values.
 
         Parameters
         ----------
         Q : numpy.ndarray
         I : numpy.ndarray
 
         Returns
         -------
         SAR : numpy.ndarray
            contains the SAR value for a particular Q matrix
         """
 
     if len(I.shape) == 1:  # Just to fit the multi-transmit case for now, ToDo
         I = np.matlib.repmat(I, Q.shape[0], 1)  # Nc x Nt
 
     Ifact = np.divide(np.matmul(I, np.matrix(I).getH()), I.shape[1])
     SAR_temp = np.multiply(Q, Ifact)
     SAR = np.abs(np.sum(SAR_temp[:]))
     return SAR
 
 
 def __loadQ():
     """
     Load Q matrix that is precomputed based on the VHM model for 8 channels.
 
     Returns
     -------
     Qtmf, Qhmf : numpy.ndarray
         Contains the Q-matrix, GSAR head and body for now.
     """
 
     # Load relevant Q matrices computed from the model - this code will be integrated later - starting from E fields
     Qpath = str(Path(__file__).parent / 'QGlobal.mat')
     Qmat = sio.loadmat(Qpath)
     Q = Qmat['Q']
     val = Q[0, 0]
 
     Qtmf = val['Qtmf']
     Qhmf = val['Qhmf']
     return Qtmf, Qhmf
 
 
 def __SAR_from_seq(seq, Qtmf, Qhmf):
     """
     Compute global whole body and head only SAR values.
 
     Parameters
     ----------
     seq : Sequence
     Qtmf : numpy.ndarray
     Qhmf : numpy.ndarray
 
     Returns
     -------
     SAR_wbg_vec : numpy.ndarray
     SAR_hg_vec : numpy.ndarray
     t_vec : numpy.ndarray
         Contains the Q-matrix, GSAR head and body for now.
     """
 
     # Identify RF blocks and compute SAR - 10 seconds must be less than twice and 6 minutes must be less than
     # 4 (WB) and 3.2 (head-20)
     block_events = seq.block_events
     num_events = len(block_events)
     t_vec = np.zeros(num_events)
     SAR_wbg_vec = np.zeros(t_vec.shape)
     SAR_hg_vec = np.zeros(t_vec.shape)
     t_prev = 0
 
     for iB in block_events:
         block = seq.get_block(iB)
         block_dur = calc_duration(block)
         t_vec[iB - 1] = t_prev + block_dur
         t_prev = t_vec[iB - 1]
         if hasattr(block, 'rf'):  # has rf
             rf = block.rf
             t = rf.t
             signal = rf.signal
             # This rf could be parallel transmit as well
             SAR_wbg_vec[iB] = __calc_SAR(Qtmf, signal)
             SAR_hg_vec[iB] = __calc_SAR(Qhmf, signal)
 
     return SAR_wbg_vec, SAR_hg_vec, t_vec
 
 
 def __SAR_interp(SAR, t):
     """
     Interpolate SAR values for one second resolution.
 
     Parameters
     ----------
     SAR : numpy.ndarray
     t : numpy.ndarray
 
     Returns
     -------
     __SAR_interp : numpy.ndarray
     tsec : numpy.ndarray
         Interpolated values of SAR for a temporal resolution of 1 second
     """
     tsec = np.arange(1, np.floor(t[-1]) + 1, 1)
     f = interpolate.interp1d(t, SAR)
     SARinterp = f(tsec)
     return SARinterp, tsec
 
 
 def __SAR_lims_check(SARwbg_lim_s, SARhg_lim_s, tsec):
     """
     Check for SAR violations as compared to IEC 10 second and 6 minute averages;
     returns SAR values that are interpolated for the fixed IEC time intervals.
 
     Parameters
     ----------
     SARwbg_lim_s : numpy.ndarray
     SARhg_lim_s : numpy.ndarray
     tsec : numpy.ndarray
 
     Returns
     -------
     SAR_wbg_tensec : numpy.ndarray
     SAR_wbg_sixmin : numpy.ndarray
     SAR_hg_tensec : numpy.ndarray
     SAR_hg_sixmin : numpy.ndarray
     SAR_wbg_sixmin_peak : numpy.ndarray
     SAR_hg_sixmin_peak : numpy.ndarray
     SAR_wbg_tensec_peak : numpy.ndarray
     SAR_hg_tensec_peak : numpy.ndarray
     """
     if tsec[-1] > 10:
         SixMinThresh_wbg = 4
         TenSecThresh_wbg = 8
 
         SixMinThresh_hg = 3.2
         TenSecThresh_hg = 6.4
 
         SARwbg_lim_app = np.concatenate((np.zeros(5), SARwbg_lim_s, np.zeros(5)), axis=0)
         SARhg_lim_app = np.concatenate((np.zeros(5), SARhg_lim_s, np.zeros(5)), axis=0)
 
         SAR_wbg_tensec = __do_sw_sar(SARwbg_lim_app, tsec, 10)  # < 2  SARmax
         SAR_hg_tensec = __do_sw_sar(SARhg_lim_app, tsec, 10)  # < 2 SARmax
         SAR_wbg_tensec_peak = np.round(np.max(SAR_wbg_tensec), 2)
         SAR_hg_tensec_peak = np.round(np.max(SAR_hg_tensec), 2)
 
         if ((np.max(SAR_wbg_tensec) > TenSecThresh_wbg) or (np.max(SAR_hg_tensec) > TenSecThresh_hg)):
             print('Pulse exceeding 10 second Global SAR limits, increase TR')
         SAR_wbg_sixmin = 'NA'
         SAR_hg_sixmin = 'NA'
         SAR_wbg_sixmin_peak = 'NA'
         SAR_hg_sixmin_peak = 'NA'
 
         if tsec[-1] > 600:
             SARwbg_lim_app = np.concatenate((np.zeros(300), SARwbg_lim_s, np.zeros(300)), axis=0)
             SARhg_lim_app = np.concatenate((np.zeros(300), SARhg_lim_s, np.zeros(300)), axis=0)
 
             SAR_hg_sixmin = __do_sw_sar(SARhg_lim_app, tsec, 600)
             SAR_wbg_sixmin = __do_sw_sar(SARwbg_lim_app, tsec, 600)
             SAR_wbg_sixmin_peak = np.round(np.max(SAR_wbg_sixmin), 2)
             SAR_hg_sixmin_peak = np.round(np.max(SAR_hg_sixmin), 2)
 
             if ((np.max(SAR_hg_sixmin) > SixMinThresh_wbg) or (np.max(SAR_hg_sixmin) > SixMinThresh_hg)):
                 print('Pulse exceeding 10 second Global SAR limits, increase TR')
     else:
         print('Need at least 10 seconds worth of sequence to calculate SAR')
         SAR_wbg_tensec = 'NA'
         SAR_wbg_sixmin = 'NA'
         SAR_hg_tensec = "NA"
         SAR_hg_sixmin = "NA"
         SAR_wbg_sixmin_peak = 'NA'
         SAR_hg_sixmin_peak = 'NA'
         SAR_wbg_tensec_peak = 'NA'
         SAR_hg_tensec_peak = 'NA'
 
     return SAR_wbg_tensec, SAR_wbg_sixmin, SAR_hg_tensec, SAR_hg_sixmin, SAR_wbg_sixmin_peak, SAR_hg_sixmin_peak, SAR_wbg_tensec_peak, SAR_hg_tensec_peak
 
 
 def __do_sw_sar(SAR, tsec, t):
     """
     Compute a sliding window average of SAR values.
 
     Parameters
     ----------
     SAR : numpy.ndarray
     tsec : numpy.ndarray
     t : numpy.ndarray
 
     Returns
     -------
     SAR_timeavag : numpy.ndarray
         Sliding window time average of SAR values
     """
     SAR_timeavg = np.zeros(len(tsec) + int(t))
     for instant in range(int(t / 2), int(t / 2) + (int(tsec[-1]))):  # better to go from  -sw / 2: sw / 2
         SAR_timeavg[instant] = sum(SAR[range(instant - int(t / 2), instant + int(t / 2) - 1)]) / t
     SAR_timeavg = SAR_timeavg[int(t / 2):int(t / 2) + (int(tsec[-1]))]
     return SAR_timeavg
 
 
 def calc_SAR(seq):
     """
     Compute Global SAR values on the `seq` object for head and whole body over the specified time averages.
 
     Parameters
     ----------
     seq : Sequence
         pypulseq `Sequence` object for which global SAR values will be computed.
     """
 
     if isinstance(seq, str):
         seq_obj = Sequence()
         seq_obj.read(seq)
         seq = seq_obj
 
     Qtmf, Qhmf = __loadQ()
     SARwbg, SARhg, t_vec = __SAR_from_seq(seq, Qtmf, Qhmf)
     SARwbg_lim, tsec = __SAR_interp(SARwbg, t_vec)
     SARhg_lim, tsec = __SAR_interp(SARhg, t_vec)
     SAR_wbg_tensec, SAR_wbg_sixmin, SAR_hg_tensec, SAR_hg_sixmin, SAR_wbg_sixmin_peak, SAR_hg_sixmin_peak, SAR_wbg_tensec_peak, SAR_hg_tensec_peak = __SAR_lims_check(
         SARwbg_lim, SARhg_lim, tsec)
 
     # Plot 10 sec average SAR
     if (tsec[-1] > 10):
         plt.plot(tsec, SAR_wbg_tensec, label='Whole Body: 10sec')
         plt.plot(tsec, SAR_hg_tensec, label='Head only: 10sec')
 
         # plt.plot(t_vec, SARwbg, label='Whole Body - instant')
         # plt.plot(t_vec, SARhg, label='Whole Body - instant')
 
         plt.xlabel('Time (s)')
         plt.ylabel('SAR (W/kg)')
         plt.title('Global SAR  - Mass Normalized -  Whole body and head only')
         ax = plt.subplot(111)
 
         # ax.legend()
         plt.grid(True)
         plt.show()