diff --git a/docs/src/guide.md b/docs/src/guide.md
index 6d05d9d..3f8bd46 100644
--- a/docs/src/guide.md
+++ b/docs/src/guide.md
@@ -81,13 +81,16 @@ N=301; # discretization size
 model=OneSpinModel(T,L,N,B)
 println(model)
 ```
+This provides us an overview of the model. It's a single spin shuttling problem with initial state `Ψ₀` and an Ornstein-Uhlenbeck noise. The total time of simulation is `T`, which is discretized into `N` steps.
 
-
-The fidelity of the spin state after shuttling can be calculated using numerical integration of the covariance matrix.  
-
-This provides us an overview of the model. It's a single spin shuttling problem with initial state `Ψ₀` and an Ornstein-Uhlenbeck noise. The total time of simulation is `T`, which is discretized into `N` steps.  
-
-The state fidelity after such a quantum process can be obtained by different numerical methods. 
+The effective noise of this spin qubit is completely characterized by its covariance matrix.  
+```@example quickstart
+heatmap(collect(sqrt.(model.R.Σ)), title="sqrt cov, 1-spin one-way shuttling", 
+size=(400,300), 
+xlabel="t1", ylabel="t2", dpi=300,
+right_margin=5Plots.mm)
+```
+The state fidelity after such a quantum process can be obtained using numerical integration of the covariance matrix.  
 ```@example quickstart
 f1=averagefidelity(model); # direct integration
 
@@ -119,7 +122,6 @@ We can check that the fidelity between the initial and final state is consistent
 f=(Ψ'*ρt*Ψ)
 ```
 
-
 ## Dephasing of entangled spin pairs during shuttling. 
 Following the approach above, we can further explore the multi-spin system. 
 The general abstraction on such a problem is given by the data type `ShuttlingModel`.  
@@ -148,7 +150,7 @@ plot!(model.R.P[N+1:2N,1], label="x2(t)")
 
 
 ```@example quickstart
-heatmap(collect(model.R.Σ)*1e3, title="covariance matrix, two spin EPR", 
+heatmap(collect(model.R.Σ)*1e3, title="covariance, 2-spin sequential shuttling", 
 size=(400,300), 
 xlabel="t1", ylabel="t2", dpi=300,
 right_margin=5Plots.mm)
diff --git a/docs/src/manual.md b/docs/src/manual.md
index 7283a96..8e0a4be 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -31,7 +31,6 @@ TwoSpinParallelModel
 dephasingmatrix
 ```
 
-
 ```@docs
 fidelity
 ```
diff --git a/src/SpinShuttling.jl b/src/SpinShuttling.jl
index 9eed241..2b43873 100644
--- a/src/SpinShuttling.jl
+++ b/src/SpinShuttling.jl
@@ -11,9 +11,9 @@ include("integration.jl")
 include("analytics.jl")
 include("stochastics.jl")
 
-export ShuttlingModel, OneSpinModel, TwoSpinModel, 
-OneSpinForthBackModel, TwoSpinParallelModel, RandomFunction, CompositeRandomFunction,
-OrnsteinUhlenbeckField, PinkBrownianField
+export ShuttlingModel, OneSpinModel, TwoSpinModel,
+    OneSpinForthBackModel, TwoSpinParallelModel, RandomFunction, CompositeRandomFunction,
+    OrnsteinUhlenbeckField, PinkBrownianField
 export averagefidelity, fidelity, sampling, characteristicfunction, characteristicvalue
 export dephasingmatrix, covariance, covariancematrix
 export W
@@ -49,9 +49,9 @@ function Base.show(io::IO, model::ShuttlingModel)
     println(io, "Time Discretization: N=$(model.N)")
     println(io, "Process Time: T=$(model.T)")
     println(io, "Shuttling Paths:")
-    t=range(0, model.T, model.N)
-    fig=lineplot(t, model.X[1].(t); width=30, height=9,
-    name="x1(t)")
+    t = range(0, model.T, model.N)
+    fig = lineplot(t, model.X[1].(t); width=30, height=9,
+        name="x1(t)")
     for i in 2:model.n
         lineplot!(fig, t, model.X[i].(t), name="x$i(t)")
     end
@@ -101,12 +101,12 @@ with total time `T` in `μs` and length `L` in `μm`.
 - `B::GaussianRandomField`: Noise field
 - `v::Real`: Velocity of the shuttling
 """
-function OneSpinForthBackModel(t::Real, T::Real, L::Real, N::Int, B::GaussianRandomField) 
-    x(t::Real, v::Real, L::Real)::Real = (t=t%(2L/v); v*t < L ? v*t : 2L-v*t)
-    return OneSpinModel(1 / √2 * [1, 1], t, N, B, τ -> x(τ, 2L/T, L))
+function OneSpinForthBackModel(t::Real, T::Real, L::Real, N::Int, B::GaussianRandomField)
+    x(t::Real, v::Real, L::Real)::Real = (t = t % (2L / v); v * t < L ? v * t : 2L - v * t)
+    return OneSpinModel(1 / √2 * [1, 1], t, N, B, τ -> x(τ, 2L / T, L))
 end
 
-function OneSpinForthBackModel(T::Real, L::Real, N::Int, B::GaussianRandomField) 
+function OneSpinForthBackModel(T::Real, L::Real, N::Int, B::GaussianRandomField)
     return OneSpinForthBackModel(T, T, L, N, B)
 end
 
@@ -174,7 +174,7 @@ The total shuttling time is `T` and the length of the path is `L` in `μm`.
 """
 function TwoSpinParallelModel(T::Real, D::Real, L::Real, N::Int,
     B::GaussianRandomField)
-    @assert length(B.θ)>=3 
+    @assert length(B.θ) >= 3
     x₁(t::Real)::Tuple{Real,Real} = (L / T * t, 0)
     x₂(t::Real)::Tuple{Real,Real} = (L / T * t, D)
     Ψ = 1 / √2 .* [0, 1, -1, 0]
@@ -233,8 +233,8 @@ function sampling(samplingfunction::Function, M::Int)::Union{Tuple{Real,Real},Tu
     Q = copy(A)
     for k in 1:M
         x = samplingfunction(k)::Union{Real,Vector{<:Real}}
-        Q = Q +(k-1)/k*(x-A).^ 2
-        A = A + (x-A)/k
+        Q = Q + (k - 1) / k * (x - A) .^ 2
+        A = A + (x - A) / k
     end
     return A, Q / (M - 1)
 end
@@ -248,7 +248,7 @@ function parallelsampling(samplingfunction::Function, M::Int)::Union{Tuple{Real,
         end
         A = mean(cache, dims=2)
         Q = var(cache, dims=2)
-        return A, Q 
+        return A, Q
     else
         cache = zeros(M)
         @threads for i in 1:M
@@ -256,7 +256,7 @@ function parallelsampling(samplingfunction::Function, M::Int)::Union{Tuple{Real,
         end
         A = mean(cache)
         Q = var(cache)
-        return A, Q 
+        return A, Q
     end
 end
 
@@ -282,34 +282,34 @@ end
 - `p::Int`: index of spin state, range from (1,2^n)
 - `n::Int`: number of spins
 """
-function m(i::Int,p::Int,n::Int)
-    1/2-digits(p, base=2, pad=n)[i]
+function m(i::Int, p::Int, n::Int)
+    1 / 2 - digits(p, base=2, pad=n)[i]
 end
 
 """
 Calculate the dephasing matrix of a given spin shuttling model.
 """
 function dephasingmatrix(model::ShuttlingModel)::Symmetric{<:Real}
-    n=model.n
-    W=zeros(2^n,2^n)
+    n = model.n
+    W = zeros(2^n, 2^n)
     for j in 1:2^n
-        W[j,j]=1
+        W[j, j] = 1
         for k in 1:j-1
-            c=[trunc(Int,m(i,j-1,n)-m(i,k-1,n)) for i in 1:n]
+            c = [trunc(Int, m(i, j - 1, n) - m(i, k - 1, n)) for i in 1:n]
             R = CompositeRandomFunction(model.R, c)
-            W[j,k] = characteristicvalue(R)
-            W[k,j] = W[j,k]
+            W[j, k] = characteristicvalue(R)
+            W[k, j] = W[j, k]
         end
     end
     return Symmetric(W)
 end
 
 function dephasingcoeffs(n::Int)::Array{Real,3}
-    M=zeros(2^n,2^n, n)
+    M = zeros(2^n, 2^n, n)
     for j in 1:2^n
         for k in 1:2^n
-            c=[m(i,j-1,n)-m(i,k-1,n) for i in 1:n]
-            M[j,k, :] = c
+            c = [m(i, j - 1, n) - m(i, k - 1, n) for i in 1:n]
+            M[j, k, :] = c
         end
     end
     return M
@@ -330,11 +330,11 @@ function fidelity(model::ShuttlingModel, randseq::Vector{<:Real}; vector::Bool=f
         Z = A
     elseif model.n == 2
         # only valid for two-spin EPR pair, ψ=1/√2(|↑↓⟩-|↓↑⟩)
-        Z = A[1:N] - A[N+1:end] 
+        Z = A[1:N] - A[N+1:end]
     else
         Z = missing
     end
-    ϕ = vector ? cumsum(Z)* dt : sum(Z) * dt
+    ϕ = vector ? cumsum(Z) * dt : sum(Z) * dt
     return (1 .+ cos.(ϕ)) / 2
 end
 
@@ -349,41 +349,41 @@ Analytical dephasing factor of a one-spin shuttling model.
 - `B<:GaussianRandomField`: Noise field, Ornstein-Uhlenbeck or Pink-Brownian
 - `path::Symbol`: Path of the shuttling model, `:straight` or `:forthback`
 """
-function W(T::Real,L::Real,B::OrnsteinUhlenbeckField; path=:straight)::Real
-    κₜ=B.θ[1]
-    κₓ=B.θ[2]
-    σ =B.σ
-    β = κₜ*T
-    γ = κₓ*L
+function W(T::Real, L::Real, B::OrnsteinUhlenbeckField; path=:straight)::Real
+    κₜ = B.θ[1]
+    κₓ = B.θ[2]
+    σ = B.σ
+    β = κₜ * T
+    γ = κₓ * L
     if path == :straight
-        return exp(- σ^2/(4*κₜ*κₓ)/κₜ^2*P1(β, γ)/2)
+        return exp(-σ^2 / (4 * κₜ * κₓ) / κₜ^2 * P1(β, γ) / 2)
     elseif path == :forthback
-        β/=2
-        return exp(- σ^2/(4*κₜ*κₓ)/κₜ^2*(P1(β, γ)+P4(β,γ)))
+        β /= 2
+        return exp(-σ^2 / (4 * κₜ * κₓ) / κₜ^2 * (P1(β, γ) + P4(β, γ)))
     else
         error("Path not recognized. Use :straight or :forthback for one-spin shuttling model.")
     end
 end
 
 function W(T::Real, L::Real, B::PinkBrownianField)::Real
-    β= T.*B.γ
-    γ= L*B.θ[1]
-    return exp(-B.σ^2*T^2*F3(β,γ))
+    β = T .* B.γ
+    γ = L * B.θ[1]
+    return exp(-B.σ^2 * T^2 * F3(β, γ))
 end
 
 
 """
 Analytical dephasing factor of a sequenced two-spin EPR pair shuttling model.
 """
-function W(T0::Real,T1::Real,L::Real,B::OrnsteinUhlenbeckField; path=:sequenced)::Real
-    κₜ=B.θ[1]
-    κₓ=B.θ[2]
-    σ =B.σ
-    τ = κₜ*T0
-    β = κₜ*T1
-    γ = κₓ*L
+function W(T0::Real, T1::Real, L::Real, B::OrnsteinUhlenbeckField; path=:sequenced)::Real
+    κₜ = B.θ[1]
+    κₓ = B.θ[2]
+    σ = B.σ
+    τ = κₜ * T0
+    β = κₜ * T1
+    γ = κₓ * L
     if path == :sequenced
-        return exp(-σ^2/(4*κₜ*κₓ)/κₜ^2*(F1(β, γ, τ)-F2(β, γ, τ)))
+        return exp(-σ^2 / (4 * κₜ * κₓ) / κₜ^2 * (F1(β, γ, τ) - F2(β, γ, τ)))
     elseif path == :parallel
         missing("Parallel path not implemented yet.")
     else
diff --git a/src/analytics.jl b/src/analytics.jl
index 6122fc9..0958cae 100644
--- a/src/analytics.jl
+++ b/src/analytics.jl
@@ -1,73 +1,72 @@
-function C1(β::Real,γ::Real,τ::Real)::Real
-    if β<=τ
-        return ℯ^(-β-γ-τ)*(1+ℯ^(2β)-2ℯ^τ+2ℯ^(β+τ)*(-β+τ))
+function C1(β::Real, γ::Real, τ::Real)::Real
+    if β <= τ
+        return ℯ^(-β - γ - τ) * (1 + ℯ^(2β) - 2ℯ^τ + 2ℯ^(β + τ) * (-β + τ))
     else
-        return ℯ^(-β-γ-τ)*(-1+ℯ^τ)^2
+        return ℯ^(-β - γ - τ) * (-1 + ℯ^τ)^2
     end
 end
 
-function C2(β::Real,γ::Real,τ::Real)::Real
-    if β<=τ
-        return ℯ^(-γ)*β*((ℯ^(-τ)*(-ℯ^β+ℯ^γ))/(β-γ)+(ℯ^(-β)*(γ-2ℯ^β*(β+γ)+ℯ^(β+γ)*(2β+γ)))/(γ*(β+γ)))
+function C2(β::Real, γ::Real, τ::Real)::Real
+    if β <= τ
+        return ℯ^(-γ) * β * ((ℯ^(-τ) * (-ℯ^β + ℯ^γ)) / (β - γ) + (ℯ^(-β) * (γ - 2ℯ^β * (β + γ) + ℯ^(β + γ) * (2β + γ))) / (γ * (β + γ)))
     else
-        return (ℯ^(-(((β+γ)*(β+3τ))/β))*β*(2ℯ^(β+γ+3τ+(2*γ*τ)/β)*(-1+ℯ^((γ*τ)/β))*β^2-ℯ^((2+(3*γ)/β)*τ)*(-1+ℯ^τ)*γ*(ℯ^τ*(β-γ)+ℯ^(β+γ)*(β+γ))))/((β-γ)*γ*(β+γ))
+        return (ℯ^(-(((β + γ) * (β + 3τ)) / β)) * β * (2ℯ^(β + γ + 3τ + (2 * γ * τ) / β) * (-1 + ℯ^((γ * τ) / β)) * β^2 - ℯ^((2 + (3 * γ) / β) * τ) * (-1 + ℯ^τ) * γ * (ℯ^τ * (β - γ) + ℯ^(β + γ) * (β + γ)))) / ((β - γ) * γ * (β + γ))
     end
 end
 
-function C3(β::Real,γ::Real,τ::Real)::Real
-    if β<=τ
-        return (ℯ^(-β-γ-τ)*β^2*(2(-1+ℯ^(2β))*β*γ+(1+ℯ^(2β)+2ℯ^(β+γ)*(-1+γ))*γ^2+β^2*(1+ℯ^(2β)-2ℯ^(β+γ)*(1+γ))))/(β^2-γ^2)^2
+function C3(β::Real, γ::Real, τ::Real)::Real
+    if β <= τ
+        return (ℯ^(-β - γ - τ) * β^2 * (2(-1 + ℯ^(2β)) * β * γ + (1 + ℯ^(2β) + 2ℯ^(β + γ) * (-1 + γ)) * γ^2 + β^2 * (1 + ℯ^(2β) - 2ℯ^(β + γ) * (1 + γ)))) / (β^2 - γ^2)^2
     else
-        return (1/((β^2-γ^2)^2))ℯ^(-(((β+γ)*(β+τ))/β))*β^2*(ℯ^((γ*τ)/β)*(β-γ)^2+ℯ^((2+γ/β)*τ)*(β-γ)^2-2ℯ^(β+γ+(γ*τ)/β)*(-((-1+γ)*γ^2)+β^2*(1+γ))+2ℯ^(β+γ+τ)*(β^3-β*(-2+γ)*γ-β^2*τ+γ^2*τ))
+        return (1 / ((β^2 - γ^2)^2))ℯ^(-(((β + γ) * (β + τ)) / β)) * β^2 * (ℯ^((γ * τ) / β) * (β - γ)^2 + ℯ^((2 + γ / β) * τ) * (β - γ)^2 - 2ℯ^(β + γ + (γ * τ) / β) * (-((-1 + γ) * γ^2) + β^2 * (1 + γ)) + 2ℯ^(β + γ + τ) * (β^3 - β * (-2 + γ) * γ - β^2 * τ + γ^2 * τ))
     end
 end
 
-function C4(β::Real,γ::Real,τ::Real)::Real
-    if β<=τ
-        return ℯ^(-γ)*β*((ℯ^-τ*(-ℯ^β+ℯ^γ))/(β-γ)+(ℯ^(-β)*(γ-2ℯ^β*(β+γ)+ℯ^(β+γ)*(2β+γ)))/(γ*(β+γ)))
+function C4(β::Real, γ::Real, τ::Real)::Real
+    if β <= τ
+        return ℯ^(-γ) * β * ((ℯ^-τ * (-ℯ^β + ℯ^γ)) / (β - γ) + (ℯ^(-β) * (γ - 2ℯ^β * (β + γ) + ℯ^(β + γ) * (2β + γ))) / (γ * (β + γ)))
     else
-        return (ℯ^(-(((β+γ)*(β+τ))/β))*β*(2ℯ^(β+γ+τ)*(-1+ℯ^((γ*τ)/β))*β^2-ℯ^((γ*τ)/β)*(-1+ℯ^τ)*(ℯ^(β+γ)+ℯ^τ)*β*γ+ℯ^((γ*τ)/β)*(-1+ℯ^τ)*(-ℯ^(β+γ)+ℯ^τ)*γ^2))/(β^2*γ-γ^3)
+        return (ℯ^(-(((β + γ) * (β + τ)) / β)) * β * (2ℯ^(β + γ + τ) * (-1 + ℯ^((γ * τ) / β)) * β^2 - ℯ^((γ * τ) / β) * (-1 + ℯ^τ) * (ℯ^(β + γ) + ℯ^τ) * β * γ + ℯ^((γ * τ) / β) * (-1 + ℯ^τ) * (-ℯ^(β + γ) + ℯ^τ) * γ^2)) / (β^2 * γ - γ^3)
     end
 end
 
-function P1(β::Real,γ::Real)::Real
-    return -((2β^2 *(1-ℯ^(-β-γ)-β-γ))/(β+γ)^2)
+function P1(β::Real, γ::Real)::Real
+    return -((2β^2 * (1 - ℯ^(-β - γ) - β - γ)) / (β + γ)^2)
 end
 
-function P2(β::Real,γ::Real,τ::Real)::Real
-    return 2*(-1+ℯ^-τ+τ)
+function P2(β::Real, γ::Real, τ::Real)::Real
+    return 2 * (-1 + ℯ^-τ + τ)
 end
 
-function P3(β::Real,γ::Real,τ::Real)::Real
-    return (ℯ^(-β-γ-τ)*(-1+ℯ^(β+γ))*(-1+ℯ^τ)*β)/(β+γ)
+function P3(β::Real, γ::Real, τ::Real)::Real
+    return (ℯ^(-β - γ - τ) * (-1 + ℯ^(β + γ)) * (-1 + ℯ^τ) * β) / (β + γ)
 end
 
-function P4(β::Real,γ::Real)::Real
-    return ((ℯ^(-2β)-1)*(γ/β)-2*ℯ^(-β-γ)+ℯ^(-2β)+1)/(1-γ^2/β^2)
+function P4(β::Real, γ::Real)::Real
+    return ((ℯ^(-2β) - 1) * (γ / β) - 2 * ℯ^(-β - γ) + ℯ^(-2β) + 1) / (1 - γ^2 / β^2)
 end
 
 """
 Ancillary function for the dephasing of the sequential shuttling pf 
 Bell state under OU sheets of noise
 """
-function F1(β::Real,γ::Real,τ::Real)::Real
-    return P1(β,γ)+P2(β,γ,τ)+2*P3(β,γ,τ)
+function F1(β::Real, γ::Real, τ::Real)::Real
+    return P1(β, γ) + P2(β, γ, τ) + 2 * P3(β, γ, τ)
 end
 
-function F2(β::Real,γ::Real,τ::Real)::Real
-    return C1(β,γ,τ)+C2(β,γ,τ)+C3(β,γ,τ)+C4(β,γ,τ)
+function F2(β::Real, γ::Real, τ::Real)::Real
+    return C1(β, γ, τ) + C2(β, γ, τ) + C3(β, γ, τ) + C4(β, γ, τ)
 end
 
 """
 Ancillary function for the dephasing of the Pink-Brownian noise.
 """
-function F3(β::Tuple{Real,Real},γ::Real)::Real
-    F(β::Real)=1/2*(expinti(-β)+(1-exp(-β))/β^2+(exp(-β)-2)/β)
-    F(β::Real,γ::Real)::Real=1/γ^2*(exp(-γ)*expinti(-β)+(γ-1)*(expinti(-β-γ)-log((β+γ)/β))-γ*((1-exp(-β-γ))/(β+γ)))
-    if γ==0 # pure 1/f noise
-        return (F(β[2])-F(β[1]))/log(β[2]/β[1])
+function F3(β::Tuple{Real,Real}, γ::Real)::Real
+    F(β::Real) = 1 / 2 * (expinti(-β) + (1 - exp(-β)) / β^2 + (exp(-β) - 2) / β)
+    F(β::Real, γ::Real)::Real = 1 / γ^2 * (exp(-γ) * expinti(-β) + (γ - 1) * (expinti(-β - γ) - log((β + γ) / β)) - γ * ((1 - exp(-β - γ)) / (β + γ)))
+    if γ == 0 # pure 1/f noise
+        return (F(β[2]) - F(β[1])) / log(β[2] / β[1])
     else
-        return (F(β[2],γ)-F(β[1],γ))/log(β[2]/β[1])
+        return (F(β[2], γ) - F(β[1], γ)) / log(β[2] / β[1])
     end
 end
-   
\ No newline at end of file
diff --git a/src/integration.jl b/src/integration.jl
index 27a8219..74c8418 100644
--- a/src/integration.jl
+++ b/src/integration.jl
@@ -1,4 +1,4 @@
-const ArrayOrSubArray{T,N} = Union{Array{T,N}, SubArray{T,N}}
+const ArrayOrSubArray{T,N} = Union{Array{T,N},SubArray{T,N}}
 
 """
 1D Simpson integral of function `f(x)` on a given array of `y=f(x₁), f(x₂)...` with constant
@@ -11,13 +11,13 @@ increment `h`
 
 """
 function integrate(y::ArrayOrSubArray{<:Real,1}, h::Real; method::Symbol=:simpson)::Real
-    n = length(y)-1
-    if method==:simpson
+    n = length(y) - 1
+    if method == :simpson
         n % 2 == 0 || error("`y` length (number of intervals) must be odd")
-        s = sum(y[1:2:n] + 4*y[2:2:n] + y[3:2:n+1])
-        return h/3 * s
-    elseif method==:trapezoid
-        s = sum(y[1:n]+y[2:n+1])
+        s = sum(y[1:2:n] + 4 * y[2:2:n] + y[3:2:n+1])
+        return h / 3 * s
+    elseif method == :trapezoid
+        s = sum(y[1:n] + y[2:n+1])
         return h / 2 * s
     else
         error("invalid integration method specified!")
@@ -42,10 +42,10 @@ end
 result = integrate(sin, 0, π, 101, method=:simpson)
 ```
 """
-function integrate(f::Function, a::Real, b::Real, n::Int; 
+function integrate(f::Function, a::Real, b::Real, n::Int;
     method::Symbol=:simpson)::Real
-    h = (b - a) / (n-1)
-    y = f.(range(a,b,n))
+    h = (b - a) / (n - 1)
+    y = f.(range(a, b, n))
     return integrate(y, h, method=method)
 end
 
@@ -70,10 +70,10 @@ end
 result = integrate((x, y) -> x * y, (0, 1, 50), (0, 1, 50), method=:simpson)
 ```
 """
-function integrate(f::Function, x_range::Tuple{Real,Real,Int}, y_range::Tuple{Real,Real,Int}; 
+function integrate(f::Function, x_range::Tuple{Real,Real,Int}, y_range::Tuple{Real,Real,Int};
     method::Symbol=:simpson)::Real
-    g = y-> integrate(x->f(x,y), x_range, method = method)
-    return integrate(g, y_range, method = method)
+    g = y -> integrate(x -> f(x, y), x_range, method=method)
+    return integrate(g, y_range, method=method)
 end
 
 """
@@ -95,21 +95,21 @@ function integrate(z::ArrayOrSubArray{<:Real,2}, h_x::Real, h_y::Real; method::S
     return integrate(integral_x_direction, h_y, method=method)
 end
 
-integrate(z::ArrayOrSubArray{<:Real,2}, h; method::Symbol=:simpson)=integrate(z, h, h; method=method)
+integrate(z::ArrayOrSubArray{<:Real,2}, h; method::Symbol=:simpson) = integrate(z, h, h; method=method)
 
 """
 Special methods for the double integral on symmetric matrix with singularity on diagonal entries.
 """
 function integrate(z::Symmetric, h::Real)::Real
     n, _ = size(z)
-    @assert n%2 == 1 "The matrix must be of odd size"
-    m=(n+1)÷2
-    _integrate(x::ArrayOrSubArray) = integrate(x, h; method = :trapezoid)
+    @assert n % 2 == 1 "The matrix must be of odd size"
+    m = (n + 1) ÷ 2
+    _integrate(x::ArrayOrSubArray) = integrate(x, h; method=:trapezoid)
     # Integrate along upper half trapzoid
-    int_upper = [_integrate((@view z[j:n, j])) for j = 1:m]|> _integrate
+    int_upper = [_integrate((@view z[j:n, j])) for j = 1:m] |> _integrate
     # Integrate along lower half trapzoid
-    int_lower = [_integrate((@view z[1:j, j])) for j = m:n]|> _integrate
+    int_lower = [_integrate((@view z[1:j, j])) for j = m:n] |> _integrate
     # Integrate the duplicated box
     int_box = _integrate((@view z[m:n, 1:m]))
-    return 2*(int_upper + int_lower - int_box)
+    return 2 * (int_upper + int_lower - int_box)
 end
\ No newline at end of file
diff --git a/src/stochastics.jl b/src/stochastics.jl
index 755b8e3..55a2287 100644
--- a/src/stochastics.jl
+++ b/src/stochastics.jl
@@ -44,8 +44,8 @@ struct RandomFunction
     Σ::Symmetric{<:Real} # covariance matrices
     C::Cholesky # decomposition
     function RandomFunction(P::Vector{<:Point}, process::GaussianRandomField)
-        μ=process.μ isa Function ? process.μ(P) : repeat([process.μ], length(P))
-        Σ=covariancematrix(P, process)
+        μ = process.μ isa Function ? process.μ(P) : repeat([process.μ], length(P))
+        Σ = covariancematrix(P, process)
         return new(μ, P, Σ, cholesky(Σ))
     end
 
@@ -66,12 +66,12 @@ Divide the covariance matrix of a direct summed random function into partitions.
 - `Matrix{Matrix{<:Real}}`: a matrix of covariance matrices
 """
 function covariancepartition(R::RandomFunction, n::Int)::Matrix{Matrix{<:Real}}
-    Λ=Matrix{Matrix{<:Real}}(undef, n, n)
-    N=length(R.P)÷n
-    Σ(i::Int,j::Int) = R.Σ[(i-1)*N+1: i*N , (j-1)*N+1: j*N]
+    Λ = Matrix{Matrix{<:Real}}(undef, n, n)
+    N = length(R.P) ÷ n
+    Σ(i::Int, j::Int) = R.Σ[(i-1)*N+1:i*N, (j-1)*N+1:j*N]
     for i in 1:n
         for j in 1:n
-            Λ[i,j]=Σ(i,j)
+            Λ[i, j] = Σ(i, j)
         end
     end
     return Λ
@@ -87,8 +87,8 @@ end
 - `Vector{Vector{<:Real}}`: a vector of mean vectors
 """
 function meanpartition(R::RandomFunction, n::Int)::Vector{Vector{<:Real}}
-    N=length(R.P)÷n
-    return [R.μ[(i-1)*N+1: i*N] for i in 1:n]
+    N = length(R.P) ÷ n
+    return [R.μ[(i-1)*N+1:i*N] for i in 1:n]
 end
 
 """
@@ -104,11 +104,11 @@ The output random function is a tensor contraction from the input.
 - `RandomFunction`: a new random function composed by a linear combination of random processes
 """
 function CompositeRandomFunction(R::RandomFunction, c::Vector{Int})::RandomFunction
-    n=length(c)
-    N=size(R.Σ,1)
-    μ = sum(c .*meanpartition(R, n))
-    Σ = Symmetric(sum((c*c') .* covariancepartition(R, n)))
-    t=[(p[1],) for p in R.P[1:(N÷n)]]
+    n = length(c)
+    N = size(R.Σ, 1)
+    μ = sum(c .* meanpartition(R, n))
+    Σ = Symmetric(sum((c * c') .* covariancepartition(R, n)))
+    t = [(p[1],) for p in R.P[1:(N÷n)]]
     return RandomFunction(μ, t, Σ, cholesky(Σ))
 end
 
@@ -124,11 +124,11 @@ Generate a random time series from a Gaussian random field.
 `R(randseq)` generates a random time series from a Gaussian random field `R` with a given random sequence `randseq`.
 """
 function (R::RandomFunction)(randseq::Vector{<:Real})
-    return R.μ .+ R.C.L*randseq
+    return R.μ .+ R.C.L * randseq
 end
 
-(R::RandomFunction)()=R(randn(size(R.Σ, 1)))
- 
+(R::RandomFunction)() = R(randn(size(R.Σ, 1)))
+
 
 """
 Covariance function of Gaussian random field.
@@ -139,7 +139,7 @@ Covariance function of Gaussian random field.
 - `process<:GaussianRandomField`: a Gaussian random field, e.g. `OrnsteinUhlenbeckField` or `PinkBrownianField`
 """
 function covariance(p₁::Point, p₂::Point, process::OrnsteinUhlenbeckField)::Real
-    process.σ^2 / prod(2* process.θ) * exp(-dot(process.θ, abs.(p₁ .- p₂)))
+    process.σ^2 / prod(2 * process.θ) * exp(-dot(process.θ, abs.(p₁ .- p₂)))
 end
 
 function covariance(p₁::Point, p₂::Point, process::PinkBrownianField)::Real
@@ -148,7 +148,7 @@ function covariance(p₁::Point, p₂::Point, process::PinkBrownianField)::Real
     x₁ = p₁[2:end]
     x₂ = p₂[2:end]
     γ = process.γ
-    cov_pink = t₁ != t₂ ? (expinti(-γ[2]abs(t₁ - t₂)) - expinti(-γ[1]abs(t₁ - t₂)))/log(γ[2]/γ[1]) : 1
+    cov_pink = t₁ != t₂ ? (expinti(-γ[2]abs(t₁ - t₂)) - expinti(-γ[1]abs(t₁ - t₂))) / log(γ[2] / γ[1]) : 1
     cov_brown = exp(-dot(process.θ, abs.(x₁ .- x₂)))
     return process.σ^2 * cov_pink * cov_brown
 end
@@ -202,13 +202,13 @@ Using Simpson's rule by default.
 """
 function characteristicfunction(R::RandomFunction)::Tuple{Vector{<:Real},Vector{<:Number}}
     # need further optimization
-    dt=R.P[2][1]-R.P[1][1]
-    N=size(R.Σ,1)
-    @assert N%2==1
-    χ(j::Int)=exp.(1im*integrate(view(R.μ, 1:j), dt))*exp.(-integrate(view(R.Σ, 1:j,1:j), dt, dt)/2)
-    t=[p[1] for p in R.P[2:2:N-1]]
-    f=[χ(j) for j in 3:2:N] # only for simpson's rule
-    return (t,f)
+    dt = R.P[2][1] - R.P[1][1]
+    N = size(R.Σ, 1)
+    @assert N % 2 == 1
+    χ(j::Int) = exp.(1im * integrate(view(R.μ, 1:j), dt)) * exp.(-integrate(view(R.Σ, 1:j, 1:j), dt, dt) / 2)
+    t = [p[1] for p in R.P[2:2:N-1]]
+    f = [χ(j) for j in 3:2:N] # only for simpson's rule
+    return (t, f)
 end
 
 """
@@ -217,6 +217,6 @@ numerical quadrature of the covariance matrix.
 Using Simpson's rule by default.
 """
 function characteristicvalue(R::RandomFunction)::Number
-    dt=R.P[2][1]-R.P[1][1]
-    return exp.(1im*integrate(R.μ, dt))*exp.(-integrate((@view R.Σ[:,:]), dt, dt)/2)
+    dt = R.P[2][1] - R.P[1][1]
+    return exp.(1im * integrate(R.μ, dt)) * exp.(-integrate((@view R.Σ[:, :]), dt, dt) / 2)
 end