diff --git a/Project.toml b/Project.toml
index 37c8cac..2283222 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ julia = "1.9"
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 
 [targets]
-test = ["Test", "ForwardDiff"]
+test = ["Test", "ForwardDiff","BenchmarkTools"]
diff --git a/docs/lit/01-autoDiff.jl b/docs/lit/01-autoDiff.jl
index fcafeee..9eaf3c8 100644
--- a/docs/lit/01-autoDiff.jl
+++ b/docs/lit/01-autoDiff.jl
@@ -62,6 +62,27 @@ lines(TE_vec,df,axis =(;title = "dS/dT2", xlabel="TE [ms]"))
 j = ForwardDiff.jacobian(x -> MESE_EPG(x,T1,TE,ETL,deltaB1),[T2])
 # ## Reproducibility
 
+# ## Derivation along multiple parameters
+function MESE_EPG2(T2T1,TE,ETL,delta)
+  T2,T1 = T2T1
+  T = complex(eltype(T2))
+  E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
+  echo_vec = Vector{Complex{eltype(T2)}}()
+
+  E = epgRotation(E,pi/2*delta, pi/2)
+  ##loop over refocusing-pulses
+  for i = 1:ETL
+    E = epgDephasing(E,1)
+    E = epgRelaxation(E,TE,T1,T2)
+    E  = epgRotation(E,pi*delta,0.0)
+    E  = epgDephasing(E,1)
+    push!(echo_vec,E.Fp[1])
+  end
+
+  return abs.(echo_vec)
+end
+
+j2 = ForwardDiff.jacobian(x -> MESE_EPG2(x,TE,ETL,deltaB1),[T2,T1])
 # This page was generated with the following version of Julia:
 using InteractiveUtils
 io = IOBuffer();
diff --git a/src/RegularEPG.jl b/src/RegularEPG.jl
index e5d4bb5..f806865 100644
--- a/src/RegularEPG.jl
+++ b/src/RegularEPG.jl
@@ -1,4 +1,5 @@
-export epgDephasing,epgRelaxation,epgRotation,rfRotation
+export epgDephasing!,epgRelaxation!,epgRotation!, epgThreshold!
+export rfRotation
 export EPGStates, getStates
 
 """
@@ -50,17 +51,22 @@ function getStates(E::EPGStates)
     return stack([E.Fp,E.Fn,E.Z],dims=1)
 end
 
+function Base.show(io::IO, E::EPGStates{T} ) where {T}
+  println(io, "EPGStates struct with fields : Fp, Fn, Z")
+  display(getStates(E))
+end
+
 """
-    epgDephasing(E::EPGStates, n=1) where T
+    epgDephasing!(E::EPGStates, n=1) where T
   
 shifts the transverse dephasing states `F` corresponding to n dephasing-cycles.
 n can be any integer
 """
-function epgDephasing(E::EPGStates, n::Int=1)
+function epgDephasing!(E::EPGStates, n::Int=1,threshold::Real=10e-6)
   
   if(abs(n)>1)
     for i in 1:abs(n)
-      E = epgDephasing(E, (n > 0 ? +1 : -1))
+      E = epgDephasing!(E, (n > 0 ? +1 : -1))
     end
   elseif(n == 1 || n == -1)
     push!(E.Fp,0)
@@ -68,15 +74,15 @@ function epgDephasing(E::EPGStates, n::Int=1)
     push!(E.Z,0)
 
     if n == 1   
-      E.Fp = circshift(E.Fp,+1)# Shift positive F states right
-      E.Fn = circshift(E.Fn,-1) # Shift negative F* states left
+      E.Fp[:] = circshift(E.Fp,+1)# Shift positive F states right
+      E.Fn[:] = circshift(E.Fn,-1) # Shift negative F* states left
 
       # Update extremal states: F_{+0} using F*_{-0}, F*_{-max+1}=0
       E.Fp[1] = conj(E.Fn[1]);
       E.Fn[end] = 0;
     else # 
-      E.Fp = circshift(E.Fp,-1)# Shift positive F states right
-      E.Fn = circshift(E.Fn,+1) # Shift negative F* states left
+      E.Fp[:] = circshift(E.Fp,-1)# Shift positive F states right
+      E.Fn[:] = circshift(E.Fn,+1) # Shift negative F* states left
 
       # Update extremal states: F_{+0} using F*_{-0}, F*_{-max+1}=0
       E.Fn[1] = conj(E.Fp[1]);
@@ -84,12 +90,32 @@ function epgDephasing(E::EPGStates, n::Int=1)
     end
 
   end
-  
+  E = epgThreshold!(E,threshold)
   return E
 end 
 
+#=
+function epgDephasing(E::EPGStates, n::Int,threshold::Real)
+  E = epgDephasing(E, n)
+  E = epgThreshold(E,threshold)
+end
+=#
+function epgThreshold!(E::EPGStates,threshold::Real)
+  threshold²=threshold^2
+  for i in length(E.Fp):-1:2
+      if abs.(E.Fp[i]^2 + E.Fn[i]^2 + E.Z[i]^2) < threshold²
+        pop!(E.Fp)
+        pop!(E.Fn)
+        pop!(E.Z)
+      else
+        return E
+      end
+  end
+  return E
+end
+
 """
-    epgRelaxation(E::EPGStates,t,T1, T2)
+    epgRelaxation!(E::EPGStates,t,T1, T2)
 
 applies relaxation matrices to a set of EPG states.
 
@@ -99,13 +125,11 @@ applies relaxation matrices to a set of EPG states.
 * `T1::AbstractFloat`   - T1
 * `T2::AbstractFloat`   - T2
 """
-function epgRelaxation(E::EPGStates,t,T1, T2)
+function epgRelaxation!(E::EPGStates,t,T1, T2)
   @. E.Fp = exp(-t/T2) * E.Fp
   @. E.Fn = exp(-t/T2) * E.Fn
-  relaxedZ = fill!(copy(E.Z), 0)
-  @. relaxedZ[2:end] = exp(-t/T1) * E.Z[2:end]
-  relaxedZ[1] = exp.(-t./T1) * (E.Z[1]-1.0) + 1.0
-  E.Z = relaxedZ
+  @. E.Z[2:end] = exp(-t/T1) * E.Z[2:end]
+  E.Z[1] = exp.(-t./T1) * (E.Z[1]-1.0) + 1.0
   return E
 end
 
@@ -126,7 +150,7 @@ end
 
 
 """
-    epgRotation(E::EPGStates, alpha::Float64, phi::Float64=0.0)
+    epgRotation!(E::EPGStates, alpha::Float64, phi::Float64=0.0)
 
 applies Bloch-rotation (<=> RF pulse) to a set of EPG states.
 
@@ -135,14 +159,31 @@ applies Bloch-rotation (<=> RF pulse) to a set of EPG states.
 * `alpha::Float64`           - flip angle of the RF pulse (rad)
 * `phi::Float64=0.0`         - phase of the RF pulse (rad)
 """
-function epgRotation(E::EPGStates, alpha::Real, phi::Real=0.0)
+function epgRotation!(E::EPGStates, alpha::Real, phi::Real=0.0)
+  R = rfRotation(alpha, phi)
+
+  epgRotation!(E, R)
+
+  return E
+end
+
+"""
+    epgRotation!(E::EPGStates, R::Matrix)
+
+applies rotation matrix from `rfRotation` function to the EPGStates
+
+# Arguments
+* `E::EPGStates``
+* `R::Matrix`           - rotation Matrix (rad)
+"""
+function epgRotation!(E::EPGStates, R::Matrix)
   # apply rotation to all states per default
   n = length(E.Z) # numStates
 
-  R = rfRotation(alpha, phi)
   for i = 1:n
     E.Fp[i],E.Fn[i],E.Z[i] = R*[E.Fp[i]; E.Fn[i]; E.Z[i]]
   end
 
   return E
-end
\ No newline at end of file
+end
+
diff --git a/test/epg/test_regular.jl b/test/epg/test_regular.jl
index 1b932ef..170a613 100644
--- a/test/epg/test_regular.jl
+++ b/test/epg/test_regular.jl
@@ -1,3 +1,22 @@
+function MESE_EPG_thresh(T2,T1,TE,ETL,delta,thresh)
+  T = eltype(complex(T2))
+  E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
+  echo_vec = Vector{Complex{eltype(T2)}}()
+
+  epgRotation!(E,pi/2*delta, pi/2)
+  # loop over refocusing-pulses
+  R = rfRotation(pi*delta,0.0)
+  for i = 1:ETL
+    epgDephasing!(E,1,thresh)
+    epgRelaxation!(E,TE,T1,T2)
+    epgRotation!(E,R)
+    epgDephasing!(E,1,thresh)
+    push!(echo_vec,E.Fp[1])
+  end
+
+  return abs.(echo_vec)
+end
+
 @testset "EPG" begin
   # test empty
   E=EPGStates()
@@ -10,30 +29,30 @@
  
   # test pulse
   E=EPGStates()
-  E = epgRotation(E,deg2rad(47),deg2rad(23))
+  epgRotation!(E,deg2rad(47),deg2rad(23))
   @test getStates(E) ≈ [
       0.2857626571584661 - im*0.6732146319308543, 
       0.2857626571584661 + im*0.6732146319308543, 
       0.6819983600624985]
 
   #test positive gradient
-  E = epgDephasing(E,1)
+  epgDephasing!(E,1)
   @test getStates(E) ≈ [[0, 0, 0.6819983600624985];;
   [0.2857626571584661 - im * 0.6732146319308543, 0, 0]]
 
   # test negative gradient
   E = EPGStates()
-  E = epgRotation(E,deg2rad(47),deg2rad(23))
-  E = epgDephasing(E,-1)
+  epgRotation!(E,deg2rad(47),deg2rad(23))
+  epgDephasing!(E,-1)
   @test getStates(E) ≈ [[0, 0, 0.6819983600624985];;
   [0, 0.2857626571584661 + im * 0.6732146319308543, 0]]
 
   # test multiple gradient
   E = EPGStates()
-  E = epgRotation(E,deg2rad(47),deg2rad(23))
-  E = epgDephasing(E,-2)
-  E = epgRotation(E,deg2rad(47),deg2rad(23))
-  E = epgDephasing(E,1)
+  epgRotation!(E,deg2rad(47),deg2rad(23))
+  epgDephasing!(E,-2)
+  epgRotation!(E,deg2rad(47),deg2rad(23))
+  epgDephasing!(E,1)
   @test getStates(E) ≈ [[0, 0, 0.4651217631279373];; 
   [0.19488966354917586-im*0.45913127494692113, 0.240326160353821+im*0.5661729534388877,0];;
   [0, 0, -0.26743911843603135];;
@@ -41,9 +60,35 @@
   
   # test relaxation
   E = EPGStates()
-  E = epgRotation(E,deg2rad(47),deg2rad(23))
-  E = epgDephasing(E,1)
-  E = epgRelaxation(E,10,1000,100)
+  epgRotation!(E,deg2rad(47),deg2rad(23))
+  epgDephasing!(E,1)
+  epgRelaxation!(E,10,1000,100)
   @test getStates(E) ≈ [[0, 0, 0.6851625292479138];;
   [0.2585687448743616 - im*0.6091497893403431, 0, 0]]
+
+  # test threshold
+  E = EPGStates([0+0*im,0+0.5im,0+0.01im],[0+0*im,0+0.5im,0+0.01im],[1+0*im,0+0im,0+0.0im])
+  epgDephasing!(E,1,10e-2)
+
+  @test getStates(E) ≈ [[0-0.5im, 0+0.5im, 1];; 
+  [0, 0+0.01im, 0];;
+  [0.5im, 0,0]]
+
+  # benchmark
+
+  b = @benchmark MESE_EPG_thresh(60.0,1000.0,7.0,50,0.9,0.0)
+  @info "without threshold (0.0) :\n 
+  time = $(median(b).time/1000) us\n
+  memory = $(median(b).memory)\n
+  allocs = $(median(b).allocs)\n
+  gctimes = $(median(b).gctime) ns\n"
+
+  b = @benchmark MESE_EPG_thresh(60.0,1000.0,7.0,50,0.9,10e-6)
+  @info "With threshold 10e-6 :\n 
+  time = $(median(b).time/1000) us\n
+  memory = $(median(b).memory)\n
+  allocs = $(median(b).allocs)\n
+  gctimes = $(median(b).gctime) ns\n"
+  
+  
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index f4d7453..de681fc 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,7 @@
 using EPGsim
 using Test
 using ForwardDiff
+using BenchmarkTools
 
 @testset "EPGsim.jl" begin
     include("epg/test_regular.jl")
diff --git a/test/test_AD.jl b/test/test_AD.jl
index 20e405a..5b9a9da 100644
--- a/test/test_AD.jl
+++ b/test/test_AD.jl
@@ -5,13 +5,13 @@ function MESE_EPG(T2,T1,TE,ETL,delta)
   E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
   echo_vec = Vector{Complex{eltype(T2)}}()
 
-  E = epgRotation(E,pi/2*delta, pi/2)
+  epgRotation!(E,pi/2*delta, pi/2)
   # loop over refocusing-pulses
   for i = 1:ETL
-    E = epgDephasing(E,1)
-    E = epgRelaxation(E,TE,T1,T2)
-    E  = epgRotation(E,pi*delta,0.0)
-    E  = epgDephasing(E,1)
+    epgDephasing!(E,1)
+    epgRelaxation!(E,TE,T1,T2)
+    epgRotation!(E,pi*delta,0.0)
+    epgDephasing!(E,1)
     push!(echo_vec,E.Fp[1])
   end
 
@@ -22,9 +22,9 @@ end
   #amp = MESE_EPG(60.0,1000.0,7,50,1) # Not used
   T2 = 60.0
   T1 = 1000.0
-  TE = 7
+  TE = 7.0
   ETL = 50
-  deltaB1 = 1
+  deltaB1 = 1.0
 
   # analytic gradient
   TE_vec = TE:TE:TE*50