diff --git a/.gitignore b/.gitignore
index c1ac5a914..ed5164d57 100644
--- a/.gitignore
+++ b/.gitignore
@@ -34,4 +34,5 @@ examples/5.koma_paper/comparison_accuracy/**/*.mrd
 examples/5.koma_paper/mrf/**/*.mrd
 examples/5.koma_paper/comparison_accuracy/**/*.mat
 examples/5.koma_paper/comparison_accuracy/**/*.h5
-!examples/5.koma_paper/Manifest.toml
\ No newline at end of file
+!examples/5.koma_paper/Manifest.toml
+*_w.seq
\ No newline at end of file
diff --git a/KomaMRIBase/src/KomaMRIBase.jl b/KomaMRIBase/src/KomaMRIBase.jl
index 3fd03d493..393ad1b71 100644
--- a/KomaMRIBase/src/KomaMRIBase.jl
+++ b/KomaMRIBase/src/KomaMRIBase.jl
@@ -39,7 +39,7 @@ export Grad, RF, ADC, Delay
 export dur, get_block_start_times, get_samples
 export DiscreteSequence
 export discretize, get_adc_phase_compensation, get_adc_sampling_times
-export is_Gx_on, is_Gy_on, is_Gz_on, is_RF_on, is_ADC_on
+export is_Gx_on, is_Gy_on, is_Gz_on, is_RF_on, is_ADC_on, is_on
 export times, ampls, freqs
 # This are also used for simulation
 export kfoldperm, trapz, cumtrapz
@@ -64,7 +64,7 @@ export get_Mk, get_kspace, get_M0, get_M1, get_M2
 include("sequences/PulseDesigner.jl")
 export PulseDesigner
 
-#Package version, KomaMRIBase.__VERSION__
+# Package version, KomaMRIBase.__VERSION__
 using Pkg
 __VERSION__ = VersionNumber(
     Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
diff --git a/KomaMRIBase/src/datatypes/Sequence.jl b/KomaMRIBase/src/datatypes/Sequence.jl
index 5501040a2..653a7cd36 100644
--- a/KomaMRIBase/src/datatypes/Sequence.jl
+++ b/KomaMRIBase/src/datatypes/Sequence.jl
@@ -121,22 +121,35 @@ recursive_merge(x...) = x[end]
 -(x::Sequence) = Sequence(-x.GR, x.RF, x.ADC, x.DUR, x.DEF)
 *(x::Sequence, α::Real) = Sequence(α .* x.GR, x.RF, x.ADC, x.DUR, x.DEF)
 *(α::Real, x::Sequence) = Sequence(α .* x.GR, x.RF, x.ADC, x.DUR, x.DEF)
-*(x::Sequence, α::ComplexF64) = Sequence(x.GR, α.*x.RF, x.ADC, x.DUR, x.DEF)
-*(α::ComplexF64, x::Sequence) = Sequence(x.GR, α.*x.RF, x.ADC, x.DUR, x.DEF)
-*(x::Sequence, A::Matrix{Float64}) = Sequence(A*x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
-*(A::Matrix{Float64}, x::Sequence) = Sequence(A*x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
-/(x::Sequence, α::Real) = Sequence(x.GR/α, x.RF, x.ADC, x.DUR, x.DEF)
+*(x::Sequence, α::ComplexF64) = Sequence(x.GR, α .* x.RF, x.ADC, x.DUR, x.DEF)
+*(α::ComplexF64, x::Sequence) = Sequence(x.GR, α .* x.RF, x.ADC, x.DUR, x.DEF)
+*(x::Sequence, A::Matrix{Float64}) = Sequence(A * x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
+*(A::Matrix{Float64}, x::Sequence) = Sequence(A * x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
+/(x::Sequence, α::Real) = Sequence(x.GR / α, x.RF, x.ADC, x.DUR, x.DEF)
 #Grad operations
-+(s::Sequence, g::Grad) = s + Sequence(reshape([g],1,1)) #Changed [a;;] for reshape(a,1,1) for Julia 1.6
-+(g::Grad, s::Sequence) = Sequence(reshape([g],1,1)) + s #Changed [a;;] for reshape(a,1,1) for Julia 1.6
++(s::Sequence, g::Grad) = s + Sequence([g;;])
++(g::Grad, s::Sequence) = Sequence([g;;]) + s
 #RF operations
-+(s::Sequence, r::RF) = s + Sequence(reshape([Grad(0.0,0.0)],1,1),reshape([r],1,1)) #Changed [a;;] for reshape(a,1,1) for Julia 1.6
-+(r::RF, s::Sequence) = Sequence(reshape([Grad(0.0,0.0)],1,1),reshape([r],1,1)) + s #Changed [a;;] for reshape(a,1,1) for Julia 1.6
++(s::Sequence, r::RF) = s + Sequence([Grad(0.0, 0.0);;], [r;;])
++(r::RF, s::Sequence) = Sequence([Grad(0.0, 0.0);;], [r;;]) + s
 #ADC operations
-+(s::Sequence, a::ADC) = s + Sequence(reshape([Grad(0.0,0.0)],1,1),reshape([RF(0.0,0.0)],1,1),[a]) #Changed [a;;] for reshape(a,1,1) for Julia 1.6
-+(a::ADC, s::Sequence) = Sequence(reshape([Grad(0.0,0.0)],1,1),reshape([RF(0.0,0.0)],1,1),[a]) + s #Changed [a;;] for reshape(a,1,1) for Julia 1.6
++(s::Sequence, a::ADC) = s + Sequence([Grad(0.0, 0.0);;], [RF(0.0, 0.0);;], [a])
++(a::ADC, s::Sequence) = Sequence([Grad(0.0, 0.0);;], [RF(0.0, 0.0);;], [a]) + s
 #Sequence object functions
-size(x::Sequence) = size(x.GR[1,:])
+size(x::Sequence) = size(x.GR[1, :])
+
+"""
+For comparing two `Sequence`s custom types
+"""
+Base.isapprox(s1::Sequence, s2::Sequence) = begin
+    (length(s1) != length(s2)) && return false
+    return all([
+        s1.ADC[i] ≈ s2.ADC[i] &&
+        s1.RF[i] ≈ s2.RF[i] &&
+        s1.GR[i] ≈ s2.GR[i] &&
+        s1.DUR[i] ≈ s2.DUR[i] for i in 1:length(s1)
+    ])
+end
 
 """
     y = is_ADC_on(x::Sequence)
diff --git a/KomaMRIBase/src/datatypes/sequence/ADC.jl b/KomaMRIBase/src/datatypes/sequence/ADC.jl
index 8b46be643..6958a7697 100644
--- a/KomaMRIBase/src/datatypes/sequence/ADC.jl
+++ b/KomaMRIBase/src/datatypes/sequence/ADC.jl
@@ -39,10 +39,16 @@ mutable struct ADC
     end
 end
 
-# ADC comparison
-Base.isapprox(adc1::ADC, adc2::ADC) = begin
-    return all(length(getfield(adc1, k)) ≈ length(getfield(adc2, k)) for k ∈ fieldnames(ADC))
-        all(getfield(adc1, k) ≈ getfield(adc2, k) for k ∈ fieldnames(ADC))
+"""
+For comparing two `ADC`s custom types
+"""
+function Base.isapprox(adc1::ADC, adc2::ADC)
+    for k in fieldnames(ADC)
+        if length(getfield(adc1, k)) != length(getfield(adc2, k))
+            return false
+        end
+    end
+    return all(getfield(adc1, k) ≈ getfield(adc2, k) for k in fieldnames(ADC))
 end
 
 """
diff --git a/KomaMRIBase/src/datatypes/sequence/Grad.jl b/KomaMRIBase/src/datatypes/sequence/Grad.jl
index 456d6694f..19e055c7e 100644
--- a/KomaMRIBase/src/datatypes/sequence/Grad.jl
+++ b/KomaMRIBase/src/datatypes/sequence/Grad.jl
@@ -213,11 +213,16 @@ getproperty(x::Matrix{Grad}, f::Symbol) = begin
     end
 end
 
-# Gradient comparison
+"""
+For comparing two `Grad`s custom types
+"""
 function Base.isapprox(gr1::Grad, gr2::Grad)
-    return all(
-        length(getfield(gr1, k)) ≈ length(getfield(gr2, k)) for k in fieldnames(Grad)
-    ) && all(getfield(gr1, k) ≈ getfield(gr2, k) for k in fieldnames(Grad))
+    for k in fieldnames(Grad)
+        if length(getfield(gr1, k)) != length(getfield(gr2, k))
+            return false
+        end
+    end
+    return all(getfield(gr1, k) ≈ getfield(gr2, k) for k in fieldnames(Grad))
 end
 
 # Gradient operations
diff --git a/KomaMRIBase/src/datatypes/sequence/RF.jl b/KomaMRIBase/src/datatypes/sequence/RF.jl
index 9d303e5bb..494b856ee 100644
--- a/KomaMRIBase/src/datatypes/sequence/RF.jl
+++ b/KomaMRIBase/src/datatypes/sequence/RF.jl
@@ -105,10 +105,16 @@ getproperty(x::Matrix{RF}, f::Symbol) = begin
     end
 end
 
-# RF comparison
+"""
+For comparing two `RF`s custom types
+"""
 function Base.isapprox(rf1::RF, rf2::RF)
-    return all(length(getfield(rf1, k)) == length(getfield(rf2, k)) for k in fieldnames(RF))
-    return all(≈(getfield(rf1, k), getfield(rf2, k); atol=1e-9) for k in fieldnames(RF))
+    for k in fieldnames(RF)
+        if length(getfield(rf1, k)) != length(getfield(rf2, k))
+            return false
+        end
+    end
+    return all(≈(getfield(rf1, k), getfield(rf2, k), atol=1e-9) for k in fieldnames(RF))
 end
 
 # Properties
diff --git a/KomaMRICore/src/KomaMRICore.jl b/KomaMRICore/src/KomaMRICore.jl
index b61f4b125..be239632c 100644
--- a/KomaMRICore/src/KomaMRICore.jl
+++ b/KomaMRICore/src/KomaMRICore.jl
@@ -29,8 +29,10 @@ export simulate, simulate_slice_profile
 # Spinors
 export Spinor, Rx, Ry, Rz, Q, Un
 
-#Package version, KomaMRICore.__VERSION__
+# Package version, KomaMRICore.__VERSION__
 using Pkg
-__VERSION__ = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
+__VERSION__ = VersionNumber(
+    Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
+)
 
 end
diff --git a/KomaMRIFiles/Project.toml b/KomaMRIFiles/Project.toml
index 240c9da81..22720d30f 100644
--- a/KomaMRIFiles/Project.toml
+++ b/KomaMRIFiles/Project.toml
@@ -9,8 +9,10 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 KomaMRIBase = "d0bc0b20-b151-4d03-b2a4-6ca51751cb9c"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+MD5 = "6ac74813-4b46-53a4-afec-0b5dc9d7885c"
 MRIFiles = "5a6f062f-bf45-497d-b654-ad17aae2a530"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Scanf = "6ef1bc8b-493b-44e1-8d40-549aa65c4b41"
 
@@ -20,6 +22,7 @@ HDF5 = "0.16, 0.17"
 InteractiveUtils = "1"
 KomaMRIBase = "0.8"
 MAT = "0.10"
+MD5 = "0.2"
 MRIFiles = "0.1, 0.2, 0.3"
 Pkg = "1.4"
 Reexport = "1"
diff --git a/KomaMRIFiles/src/KomaMRIFiles.jl b/KomaMRIFiles/src/KomaMRIFiles.jl
index 259639a62..a48cfc9f2 100644
--- a/KomaMRIFiles/src/KomaMRIFiles.jl
+++ b/KomaMRIFiles/src/KomaMRIFiles.jl
@@ -1,24 +1,37 @@
 module KomaMRIFiles
 
 using KomaMRIBase
-using Scanf, FileIO, HDF5, MAT, InteractiveUtils # IO related
+# .mrd .jld2
+using FileIO
+# For reading .seq
+using Scanf
+# For reading .h5/.mat phantoms
+using HDF5, MAT
+# For writing .seq
+using Printf, MD5
+# For writing .phantom
+using InteractiveUtils
 
 using Reexport
 using MRIFiles
-import MRIFiles: insertNode
 @reexport using MRIFiles: ISMRMRDFile
 @reexport using FileIO: save
 
-include("Sequence/Pulseq.jl")
+include("Sequence/ReadPulseq.jl")
+include("Sequence/WritePulseq.jl")
 include("Phantom/JEMRIS.jl")
 include("Phantom/MRiLab.jl")
 include("Phantom/Phantom.jl")
 
-export read_seq                                                                     # Pulseq
-export read_phantom_jemris, read_phantom_MRiLab, read_phantom, write_phantom        # Phantom
+# Pulseq
+export read_seq, write_seq
+# Phantom
+export read_phantom_jemris, read_phantom_MRiLab, read_phantom, write_phantom
 
 # Package version: KomaMRIFiles.__VERSION__
 using Pkg
-__VERSION__ = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
+__VERSION__ = VersionNumber(
+    Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
+)
 
 end # module KomaMRIFiles
diff --git a/KomaMRIFiles/src/Sequence/Pulseq.jl b/KomaMRIFiles/src/Sequence/ReadPulseq.jl
similarity index 99%
rename from KomaMRIFiles/src/Sequence/Pulseq.jl
rename to KomaMRIFiles/src/Sequence/ReadPulseq.jl
index 2c41575b6..b287b0f96 100644
--- a/KomaMRIFiles/src/Sequence/Pulseq.jl
+++ b/KomaMRIFiles/src/Sequence/ReadPulseq.jl
@@ -520,7 +520,7 @@ function read_Grad(gradLibrary, shapeLibrary, Δt_gr, i)
             G = Grad(gA,gT,0.0,0.0,delay)
         end
     end
-    G
+    return G
 end
 
 """
diff --git a/KomaMRIFiles/src/Sequence/WritePulseq.jl b/KomaMRIFiles/src/Sequence/WritePulseq.jl
new file mode 100644
index 000000000..6988f45c3
--- /dev/null
+++ b/KomaMRIFiles/src/Sequence/WritePulseq.jl
@@ -0,0 +1,481 @@
+"""
+Returns a boolean value indicating whether two vector of angles are equal.
+"""
+function safe_equal_angles(a1, a2)
+    if length(a1) != length(a2)
+        return false
+    end
+    return abs(sum(exp.(1im * 2π * a1) .* exp.(-1im * 2π * a2)) / length(a1)) == 1
+end
+
+"""
+Returns a boolean value indicating whether a vector is present in a list of vectors.
+"""
+function not_in_list_angles(angle, angle_list)
+    for phase in angle_list
+        if safe_equal_angles(angle, phase)
+            return false
+        end
+    end
+    return true
+end
+
+"""
+Returns the events and IDs which are unique given an input vector of events.
+"""
+function get_unique_events(events::Vector)
+    events_obj, events_id, id_cnt = [], [], 1
+    for obj in events
+        if !(obj in events_obj)
+            push!(events_obj, obj)
+            push!(events_id, id_cnt)
+            id_cnt += 1
+        end
+    end
+    return (obj = events_obj , id = events_id)
+end
+
+"""
+Returns the vector of IDs for all the blocks of an array of events.
+It is neccessary to add the input vector `events` which contains the uniques objects
+(RF, Grad, or ADC) with its repective ID.
+"""
+function get_events_id(event_array, events)
+    events_id = fill(0, length(event_array))
+    for (blk, eve) in enumerate(event_array)
+        if is_on(eve)
+            for (obj, id) in zip(events.obj, events.id)
+                if eve ≈ obj
+                    events_id[blk] = id
+                    break
+                end
+            end
+        end
+    end
+    return events_id
+end
+
+"""
+Separates the `grs` into the gradients with arbitrary waveform and trapezoidal
+"""
+function get_unique_gradients(grs)
+    grads = (obj = [], id = [])
+    traps = (obj = [], id = [])
+    for (obj, id) in zip(grs.obj, grs.id)
+        if length(obj.A) == 1
+            push!(traps.obj, obj)
+            push!(traps.id, id)
+        else
+            push!(grads.obj, obj)
+            push!(grads.id, id)
+        end
+    end
+    return grads, traps
+end
+
+"""
+Returns the unique shapes for the magnitude, angle and time of the `rfs` vector.
+Requires an initial integer counter `id_shape_cnt` to asign IDs incrementally.
+"""
+function get_rf_shapes(rfs, id_shape_cnt, Δt_rf)
+    # Find the unique shapes (magnitude, phase and time shapes) and assign IDs
+    rfs_abs, rfs_ang, rfs_tim = [], [], []
+    rfs_abs_id, rfs_ang_id, rfs_tim_id = [], [], []
+    for obj in rfs.obj
+        shape_abs = abs.(obj.A) / maximum(abs.(obj.A))
+        if !(shape_abs in rfs_abs)
+            push!(rfs_abs, shape_abs)
+            push!(rfs_abs_id, id_shape_cnt)
+            id_shape_cnt += 1
+        end
+        shape_ang = mod.(angle.(obj.A), 2π) / 2π
+        ang = shape_ang .- shape_ang[1]
+        list_ang_unique = [shape .- shape[1] for shape in rfs_ang]
+        if not_in_list_angles(ang, list_ang_unique)
+            push!(rfs_ang, shape_ang)
+            push!(rfs_ang_id, id_shape_cnt)
+            id_shape_cnt += 1
+        end
+        if isa(obj.T, Vector{<:Number})
+            shape_tim = cumsum([0; obj.T]) / Δt_rf
+            if !(shape_tim in rfs_tim)
+                push!(rfs_tim, shape_tim)
+                push!(rfs_tim_id, id_shape_cnt)
+                id_shape_cnt += 1
+            end
+        end
+    end
+    rfmags = (data = rfs_abs, id = rfs_abs_id)
+    rfangs = (data = rfs_ang, id = rfs_ang_id)
+    rftimes = (data = rfs_tim, id = rfs_tim_id)
+    return rfmags, rfangs, rftimes, id_shape_cnt
+end
+
+"""
+Returns the unique shapes for the amplitude and time of the `grads` vector.
+Requires an initial integer counter `id_shape_cnt` to asign IDs incrementally.
+"""
+function get_grad_shapes(grads, id_shape_cnt, Δt_gr)
+    # Find shapes for magnitude and time gradients
+    grads_amp, grads_tim = [], []
+    grads_amp_id, grads_tim_id = [], []
+    for obj in grads.obj
+        shape_amp = obj.A / maximum(abs.(obj.A))
+        if !(shape_amp in grads_amp)
+            push!(grads_amp, shape_amp)
+            push!(grads_amp_id, id_shape_cnt)
+            id_shape_cnt = id_shape_cnt + 1
+        end
+        if isa(obj.T, Vector{<:Number})
+            shape_tim = cumsum([0; obj.T]) / Δt_gr
+            if !(shape_tim in grads_tim)
+                push!(grads_tim, shape_tim)
+                push!(grads_tim_id, id_shape_cnt)
+                id_shape_cnt = id_shape_cnt + 1
+            end
+        end
+    end
+    gramps = (data = grads_amp, id = grads_amp_id)
+    grtimes = (data = grads_tim, id = grads_tim_id)
+    return gramps, grtimes, id_shape_cnt
+end
+
+"""
+Defines the library to be written in the [BLOCKS] section
+Columns of block_events: [blk, id_rf, id_gx, id_gy, id_gz, id_adc, id_ext]
+"""
+function get_block_events(seq, rfs, grs, adcs)
+    r = get_events_id(seq.RF, rfs)
+    x = get_events_id(seq.GR.x, grs)
+    y = get_events_id(seq.GR.y, grs)
+    z = get_events_id(seq.GR.z, grs)
+    a = get_events_id(seq.ADC, adcs)
+    block_events = [[b, 0, r[b], x[b], y[b], z[b], a[b], 0] for (b, _) in enumerate(seq)]
+    for row in block_events
+        blk = row[1]
+        bd = seq.DUR[blk] / seq.DEF["BlockDurationRaster"]
+        bdr = round(bd)
+        if abs(bdr - bd) > 1e-6
+            @warn "Block $blk duration rounded"
+        end
+        row[2] = bdr
+    end
+    return block_events
+end
+
+"""
+Defines the library to be written in the [RF] section
+Columns of rf_library: [id, amp, id_mag, id_phase, id_time, delay, freq, phase]
+"""
+function get_rf_library(rfs, rfmags, rfangs, rftimes, Δt_rf)
+    rf_library = [[id, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] for id in rfs.id]
+    for (i, row) in enumerate(rf_library)
+        obj = rfs.obj[i]
+        row[2] = γ * maximum(abs.(obj.A))
+        row[7] = obj.Δf
+        shape_abs = abs.(obj.A) / maximum(abs.(obj.A))
+        for (shape_abs_unique, id_abs) in zip(rfmags.data, rfmags.id)
+            if shape_abs == shape_abs_unique
+                row[3] = id_abs
+                break
+            end
+        end
+        shape_ang = mod.(angle.(obj.A), 2π) / 2π
+        ang = shape_ang .- shape_ang[1]
+        for (shape_ang_unique, id_ang) in zip(rfangs.data, rfangs.id)
+            ang_unique = shape_ang_unique .- shape_ang_unique[1]
+            if safe_equal_angles(ang, ang_unique)
+                row[4] = id_ang
+                row[8] = angle(
+                    sum(exp.(1im * 2π * shape_ang) .* exp.(-1im * 2π * shape_ang_unique)) /
+                    length(shape_ang),
+                )
+                break
+            end
+        end
+        if isa(obj.T, Vector{<:Number})
+            shape_tim = cumsum([0; obj.T]) / Δt_rf
+            for (shape_tim_unique, id_tim) in zip(rftimes.data, rftimes.id)
+                if shape_tim == shape_tim_unique
+                    row[5] = id_tim
+                    break
+                end
+            end
+        end
+        delay_compensation_rf_koma = (row[5] == 0) * Δt_rf / 2
+        row[6] = round((obj.delay - delay_compensation_rf_koma) / Δt_rf) * Δt_rf * 1e6
+    end
+    return rf_library
+end
+
+"""
+Defines the library to be written in the [GRADIENTS] section
+Columns of grad_library_arb: [id, amp, id_amp, id_time, delay]
+"""
+function get_grad_library_arb(grads, gramps, grtimes, Δt_gr)
+    grad_library_arb = [[id, 0.0, 0.0, 0.0, 0.0] for id in grads.id]
+    for (i, row) in enumerate(grad_library_arb)
+        obj = grads.obj[i]
+        row[2] = γ * maximum(abs.(obj.A))    # this always stores positive values, the waveform vector have the respective positive or negative values
+        row[5] = round(1e6 * obj.delay)
+        shape_amp = obj.A / maximum(abs.(obj.A))
+        for (shape_amp_unique, id_amp) in zip(gramps.data, gramps.id)
+            if shape_amp == shape_amp_unique
+                row[3] = id_amp
+                break
+            end
+        end
+        if isa(obj.T, Vector{<:Number})
+            shape_tim = cumsum([0; obj.T]) / Δt_gr
+            for (shape_tim_unique, id_tim) in zip(grtimes.data, grtimes.id)
+                if shape_tim == shape_tim_unique
+                    row[4] = id_tim
+                    break
+                end
+            end
+        end
+    end
+    return grad_library_arb
+end
+
+"""
+Defines the library to be written in the [TRAP] section
+Columns of grad_library_trap: [id, amp, rise, flat, fall, delay]
+"""
+function get_grad_library_trap(traps)
+    grad_library_trap = [[id, 0.0, 0.0, 0.0, 0.0, 0.0] for id in traps.id]
+    for (i, row) in enumerate(grad_library_trap)
+        obj = traps.obj[i]
+        row[2] = γ * obj.A
+        row[3] = 1e6 * obj.rise
+        row[4] = 1e6 * obj.T
+        row[5] = 1e6 * obj.fall
+        row[6] = 1e6 * obj.delay
+    end
+    return grad_library_trap
+end
+
+"""
+Defines the library to be written in the [ADC] section
+Columns of adc_library: [id, num, dwell, delay, freq, phase]
+"""
+function get_adc_library(adcs)
+    adc_library = [[id, 0.0, 0.0, 0.0, 0.0, 0.0] for id in adcs.id]
+    for (i, row) in enumerate(adc_library)
+        obj = adcs.obj[i]
+        row[2] = obj.N
+        row[3] = obj.T * 1e9 / (obj.N - 1)
+        row[4] = (obj.delay - 0.5 * obj.T / (obj.N - 1)) * 1e6
+        row[5] = obj.Δf
+        row[6] = obj.ϕ
+    end
+    return adc_library
+end
+
+"""
+Defines the library to be written for the [SHAPES] section
+Elements of shape_library: [id, num, data]
+"""
+function get_shape_library(rfmags, rfangs, rftimes, gramps, grtimes)
+    shape_library = []
+    for shapes in [rfmags, rfangs, rftimes, gramps, grtimes]
+        for (data, id) in zip(shapes.data, shapes.id)
+            num, data = compress_shape(data)
+            push!(shape_library, (id = id, num = num, data = data))
+        end
+    end
+    return shape_library
+end
+
+"""
+Returns the with its libraries to be written in the obj file
+"""
+function get_pulseq_object(seq::Sequence)
+    # Get the raster times
+    Δt_rf = seq.DEF["RadiofrequencyRasterTime"]
+    Δt_gr = seq.DEF["GradientRasterTime"]
+    # Get the unique events and its IDs
+    rfs = get_unique_events(seq.RF[is_on.(seq.RF)])
+    grs = get_unique_events(seq.GR[is_on.(seq.GR)])
+    adcs = get_unique_events(seq.ADC[is_on.(seq.ADC)])
+    grads, traps = get_unique_gradients(grs)
+    rfmags, rfangs, rftimes, cnt = get_rf_shapes(rfs, 1, Δt_rf)
+    gramps, grtimes, _ = get_grad_shapes(grads, cnt, Δt_gr)
+    # Return the pulseq object
+    return (
+        blockEvents = get_block_events(seq, rfs, grs, adcs),
+        rfLibrary = get_rf_library(rfs, rfmags, rfangs, rftimes, Δt_rf),
+        gradLibrary = (
+            arb = get_grad_library_arb(grads, gramps, grtimes, Δt_gr),
+            trap = get_grad_library_trap(traps),
+        ),
+        adcLibrary = get_adc_library(adcs),
+        shapeLibrary = get_shape_library(rfmags, rfangs, rftimes, gramps, grtimes),
+    )
+end
+
+"""
+    write_seq(seq::Sequence, filename::String)
+
+Writes a .seq file for a given sequence `seq` y the location `filename`
+"""
+function write_seq(seq::Sequence, filename)
+    # Just a warning message for extensions
+    #@warn "EXTENSIONS will not be handled"
+    # Get the pulseq object
+    obj = get_pulseq_object(seq)
+    # Write the .seq file
+    open(filename, "w") do fid
+        @printf(
+            fid,
+            """
+            # Pulseq sequence file
+            # Created with KomaMRIFiles.jl %s
+
+            [VERSION]
+            major 1
+            minor 4
+            revision 2
+            """,
+            string(KomaMRIFiles.__VERSION__)
+        )
+        if !isempty(seq.DEF)
+            @printf(
+                fid,
+                """
+
+                [DEFINITIONS]
+                """
+            )
+            sorted_keys = sort(collect(keys(seq.DEF)))
+            filter!(x -> x != "extension" && x != "additional_text", sorted_keys) # For now remove extensions
+            for key in sorted_keys
+                val = seq.DEF[key]
+                @printf(fid, "%s ", key)
+                if isa(val, String)
+                    @printf(fid, "%s ", val)
+                else
+                    if isa(val, Vector{<:Number})
+                        for v in val
+                            @printf(fid, "%.9g ", v)
+                        end
+                    else
+                        @printf(fid, "%.9g ", val)
+                    end
+                end
+                @printf(fid, "\n")
+            end
+        end
+        if !isempty(obj.blockEvents)
+            @printf(
+                fid,
+                """
+
+                # Format of blocks:
+                # NUM DUR RF  GX  GY  GZ  ADC  EXT
+                [BLOCKS]
+                """
+            )
+            id_format_str = "%" * string(length(string(length(seq)))) * "d "
+            fmt = Printf.Format(id_format_str * "%3d %3d %3d %3d %3d %2d %2d\n")
+            for row in obj.blockEvents
+                Printf.format(fid, fmt, row...)
+            end
+        end
+        if !isempty(obj.rfLibrary)
+            @printf(
+                fid,
+                """
+
+                # Format of RF events:
+                # id amplitude mag_id phase_id time_shape_id delay freq phase
+                # ..        Hz   ....     ....          ....    us   Hz   rad
+                [RF]
+                """
+            )
+            fmt = Printf.Format("%d %12g %d %d %d %g %g %g\n")
+            for row in obj.rfLibrary
+                Printf.format(fid, fmt, row...)
+            end
+        end
+        if !isempty(obj.gradLibrary.arb)
+            @printf(
+                fid,
+                """
+
+                # Format of arbitrary gradients:
+                #   time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)
+                # id amplitude amp_shape_id time_shape_id delay
+                # ..      Hz/m       ..         ..          us
+                [GRADIENTS]
+                """
+            )
+            for row in obj.gradLibrary.arb
+                @printf(fid, "%d %12g %d %d %d\n", row...)
+            end
+        end
+        if !isempty(obj.gradLibrary.trap)
+            @printf(
+                fid,
+                """
+
+                # Format of trapezoid gradients:
+                # id amplitude rise flat fall delay
+                # ..      Hz/m   us   us   us    us
+                [TRAP]
+                """
+            )
+            for row in obj.gradLibrary.trap
+                @printf(fid, "%2d %12g %3d %4d %3d %3d\n", row...)
+            end
+        end
+        if !isempty(obj.adcLibrary)
+            @printf(
+                fid,
+                """
+
+                # Format of ADC events:
+                # id num dwell delay freq phase
+                # ..  ..    ns    us   Hz   rad
+                [ADC]
+                """
+            )
+            for row in obj.adcLibrary
+                @printf(fid, "%d %d %.0f %.0f %g %g\n", row...)
+            end
+        end
+        if !isempty(obj.shapeLibrary)
+            @printf(
+                fid,
+                """
+
+                # Sequence Shapes
+                [SHAPES]
+
+                """
+            )
+            for row in obj.shapeLibrary
+                @printf(fid, "shape_id %d\n", row.id)
+                @printf(fid, "num_samples %d\n", row.num)
+                [@printf(fid, "%.9g\n", i) for i in row.data]
+                @printf(fid, "\n")
+            end
+        end
+    end
+    md5hash = bytes2hex(open(md5, filename))
+    open(filename, "a") do fid
+        @printf(
+            fid,
+            """
+
+            [SIGNATURE]
+            # This is the hash of the Pulseq file, calculated right before the [SIGNATURE] section was added
+            # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]
+            # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be sripped away for recalculating/verification)
+            Type md5
+            Hash %s
+            """,
+            md5hash
+        )
+    end
+end
diff --git a/KomaMRIFiles/test/runtests.jl b/KomaMRIFiles/test/runtests.jl
index deccbf82d..7cee3da36 100644
--- a/KomaMRIFiles/test/runtests.jl
+++ b/KomaMRIFiles/test/runtests.jl
@@ -1,14 +1,14 @@
 using TestItems, TestItemRunner
 
-@run_package_tests filter=t_start->!(:skipci in t_start.tags)&&(:files in t_start.tags) #verbose=true
+@run_package_tests filter = ti -> !(:skipci in ti.tags) && (:files in ti.tags) #verbose=true
 
-@testitem "Files" tags=[:files] begin
+@testitem "Files" tags = [:files] begin
     using Suppressor
 
-    # Test Pulseq
-    @testset "Pulseq" begin
+    # Test ReadPulseq
+    @testset "ReadPulseq" begin
         path = @__DIR__
-        seq = @suppress read_seq(path*"/test_files/epi.seq") #Pulseq v1.4.0, RF arbitrary
+        seq = @suppress read_seq(path * "/test_files/epi.seq") #Pulseq v1.4.0, RF arbitrary
         @test seq.DEF["FileName"] == "epi.seq"
         @test seq.DEF["PulseqVersion"] ≈ 1004000
         @test seq.DEF["signature"] == "67ebeffe6afdf0c393834101c14f3990"
@@ -34,10 +34,11 @@ using TestItems, TestItemRunner
         shape2 = KomaMRIFiles.decompress_shape(num_samples, compressed_data)
         @test shape == shape2
     end
+
     # Test JEMRIS
     @testset "JEMRIS" begin
         path = @__DIR__
-        obj = read_phantom_jemris(path*"/test_files/column1d.h5")
+        obj = read_phantom_jemris(path * "/test_files/column1d.h5")
         @test obj.name == "column1d.h5"
     end
     # Test MRiLab
@@ -93,7 +94,7 @@ using TestItems, TestItemRunner
     end
 end
 
-@testitem "Pulseq compat" tags=[:files, :pulseq] begin
+@testitem "Pulseq compat" tags = [:files, :pulseq] begin
     using MAT, KomaMRIBase, Suppressor
 
     # Aux functions
@@ -103,7 +104,7 @@ end
     not_empty = ((ek, ep),) -> !isempty(ep.t)
 
     # Reading files
-    path         = joinpath(@__DIR__, "test_files/pulseq_read_comparison")
+    path = joinpath(@__DIR__, "test_files/pulseq_examples")
     pulseq_files = filter(endswith(".seq"), readdir(path)) .|> x -> splitext(x)[1]
     for pulseq_file in pulseq_files
         #@show pulseq_file
@@ -123,3 +124,17 @@ end
         end
     end
 end
+
+@testitem "WritePulseq" tags = [:files, :pulseq] begin
+    path = joinpath(@__DIR__, "test_files/pulseq_examples")
+    pulseq_files =
+        filter(endswith(r"^(?!.*_w\.seq$).*\.seq$"), readdir(path)) .|> x -> splitext(x)[1]
+    for pulseq_file in pulseq_files
+        seq_original = read_seq("$path/$(pulseq_file).seq")
+        write_seq(seq_original, "$path/$(pulseq_file)_w.seq")
+        seq_written = read_seq("$path/$(pulseq_file)_w.seq")
+        @testset "$pulseq_file" begin
+            @test seq_original ≈ seq_written
+        end
+    end
+end
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/epi.mat b/KomaMRIFiles/test/test_files/pulseq_examples/epi.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/epi.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/epi.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/epi.seq b/KomaMRIFiles/test/test_files/pulseq_examples/epi.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/epi.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/epi.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/fid.mat b/KomaMRIFiles/test/test_files/pulseq_examples/fid.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/fid.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/fid.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/fid.seq b/KomaMRIFiles/test/test_files/pulseq_examples/fid.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/fid.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/fid.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-time-shaped.mat b/KomaMRIFiles/test/test_files/pulseq_examples/gr-time-shaped.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-time-shaped.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/gr-time-shaped.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-time-shaped.seq b/KomaMRIFiles/test/test_files/pulseq_examples/gr-time-shaped.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-time-shaped.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/gr-time-shaped.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-trapezoidal.mat b/KomaMRIFiles/test/test_files/pulseq_examples/gr-trapezoidal.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-trapezoidal.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/gr-trapezoidal.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-trapezoidal.seq b/KomaMRIFiles/test/test_files/pulseq_examples/gr-trapezoidal.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-trapezoidal.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/gr-trapezoidal.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-uniformly-shaped.mat b/KomaMRIFiles/test/test_files/pulseq_examples/gr-uniformly-shaped.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-uniformly-shaped.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/gr-uniformly-shaped.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-uniformly-shaped.seq b/KomaMRIFiles/test/test_files/pulseq_examples/gr-uniformly-shaped.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/gr-uniformly-shaped.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/gr-uniformly-shaped.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/gre.mat b/KomaMRIFiles/test/test_files/pulseq_examples/gre.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/gre.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/gre.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/gre.seq b/KomaMRIFiles/test/test_files/pulseq_examples/gre.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/gre.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/gre.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-pulse.mat b/KomaMRIFiles/test/test_files/pulseq_examples/rf-pulse.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-pulse.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/rf-pulse.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-pulse.seq b/KomaMRIFiles/test/test_files/pulseq_examples/rf-pulse.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-pulse.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/rf-pulse.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-time-shaped.mat b/KomaMRIFiles/test/test_files/pulseq_examples/rf-time-shaped.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-time-shaped.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/rf-time-shaped.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-time-shaped.seq b/KomaMRIFiles/test/test_files/pulseq_examples/rf-time-shaped.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-time-shaped.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/rf-time-shaped.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-uniformly-shaped.mat b/KomaMRIFiles/test/test_files/pulseq_examples/rf-uniformly-shaped.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-uniformly-shaped.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/rf-uniformly-shaped.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-uniformly-shaped.seq b/KomaMRIFiles/test/test_files/pulseq_examples/rf-uniformly-shaped.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/rf-uniformly-shaped.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/rf-uniformly-shaped.seq
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/spiral.mat b/KomaMRIFiles/test/test_files/pulseq_examples/spiral.mat
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/spiral.mat
rename to KomaMRIFiles/test/test_files/pulseq_examples/spiral.mat
diff --git a/KomaMRIFiles/test/test_files/pulseq_read_comparison/spiral.seq b/KomaMRIFiles/test/test_files/pulseq_examples/spiral.seq
similarity index 100%
rename from KomaMRIFiles/test/test_files/pulseq_read_comparison/spiral.seq
rename to KomaMRIFiles/test/test_files/pulseq_examples/spiral.seq
diff --git a/KomaMRIPlots/src/KomaMRIPlots.jl b/KomaMRIPlots/src/KomaMRIPlots.jl
index 8594482e9..0542f1d9a 100644
--- a/KomaMRIPlots/src/KomaMRIPlots.jl
+++ b/KomaMRIPlots/src/KomaMRIPlots.jl
@@ -22,7 +22,7 @@ export plot_seq,
     plot_image,
     plot_dict
 
-#Package version, KomaMRIPlots.__VERSION__
+# Package version, KomaMRIPlots.__VERSION__
 using Pkg
 __VERSION__ = VersionNumber(
     Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
diff --git a/docs/src/reference/4-koma-files.md b/docs/src/reference/4-koma-files.md
index 2640f7205..c9b927a3f 100644
--- a/docs/src/reference/4-koma-files.md
+++ b/docs/src/reference/4-koma-files.md
@@ -8,6 +8,7 @@ CurrentModule = KomaMRIFiles
 
 ```@docs
 read_seq
+write_seq
 ```
 
 ## Phantom
diff --git a/src/KomaMRI.jl b/src/KomaMRI.jl
index b4412892e..32929e94a 100644
--- a/src/KomaMRI.jl
+++ b/src/KomaMRI.jl
@@ -28,8 +28,10 @@ include("KomaUI.jl")
 export KomaUI
 export sys_ui, seq_ui, obj_ui, raw_ui, img_ui
 
-#Package version, KomaMRI.__VERSION__
+# Package version, KomaMRI.__VERSION__
 using Pkg
-__VERSION__ = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
+__VERSION__ = VersionNumber(
+    Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
+)
 
 end