diff --git a/src/FloquetUtils.jl b/src/FloquetUtils.jl
new file mode 100644
index 0000000..ad2ca3f
--- /dev/null
+++ b/src/FloquetUtils.jl
@@ -0,0 +1,123 @@
+module FloquetUtils
+
+import LinearAlgebra: I, eigvecs, eigen, Diagonal
+import DifferentialEquations as de
+
+export time_evolution_op
+export floquet_hamiltonian
+export KickOperator
+
+using ..Utilities
+
+"""
+    time_evolution_op(H, T, args...; kwargs...)
+
+Calculate the time evolution operator for a Hamiltonian `H` up to a
+total time `T`. The rest arguments are passed on to
+`DifferentialEquations.solve`.
+"""
+function time_evolution_op(H::Function, T::Real, args...; kwargs...)
+    function u_step!(dU, U, _, t)
+        dU .= -1im * H(t) * U
+    end
+
+    N = size(H(0), 1)
+    u0 = Matrix{ComplexF64}(I, N, N)
+
+    tspan = (0.0, T)
+    prob = de.ODEProblem(u_step!, u0, tspan)
+    de.solve(prob, args...; kwargs...)
+end
+
+
+"""
+    floquet_hamiltonian(U, T)
+
+Returns the Floquet Hamiltonian given a unitary matrix `U` and a time
+`T`.
+"""
+function floquet_hamiltonian(U::Matrix, T::Real)
+    1im / T * log(U)
+end
+
+
+
+"""
+    floquet_hamiltonian(H(t), T)
+
+Returns the Floquet Hamiltonian given a Hamiltionian `H` and a time
+`T`. The rest arguments are passed on to
+`DifferentialEquations.solve`.
+"""
+function floquet_hamiltonian(H::Function, T::Real, args...; kwargs...)
+    U = time_evolution_op(H, T, args...; kwargs...)
+    floquet_hamiltonian(U(T), T)
+end
+
+"""
+    trivial_floquet_hamiltonian(H, T, restrict_energies=false)
+
+Returns the Floquet Hamiltonian for a Hamiltonain `H` that commutes
+with itself at different times.
+"""
+function trivial_floquet_hamiltonian(H::Function, T::Real, restrict_energies::Bool=false)
+    function h_step!(dU, _, _, t)
+        dU .= H(t)
+    end
+
+    dimension = size(H(0), 1)
+    u0 = zeros(ComplexF64, dimension, dimension)
+    prob = de.ODEProblem(h_step!, u0, (0.0, T))
+
+    integrated_hamiltonian = de.solve(prob)
+
+    H_F = integrated_hamiltonian(T) / T
+
+    if restrict_energies
+        energies, vecs = eigen(H_F)
+        energies .= restrict_to_range.(energies, π / T)
+        H_F .= vecs * Diagonal(energies) * vecs'
+    end
+
+    H_F
+end
+
+
+struct KickOperator
+    U
+    H_F::Matrix
+
+
+    """
+        KickOperator(U(t), H_F)
+
+    Return the Kick operator Given the time evolution operator `U(t)` and the Floquet Hamiltonian
+    `H_F`.
+    """
+    function KickOperator(U, H_F::Matrix)
+        if size(U(0)) != size(H_F)
+            throw(DimensionMismatch("`U` and `H_F` should have the same dimension."))
+        end
+        new(U, H_F)
+    end
+end
+
+
+function KickOperator(H::Function, T::Real)
+    U = time_evolution_op(H, T)
+    H_F = floquet_hamiltonian(U(T), T)
+    KickOperator(U, H_F)
+end
+
+"""
+    K(t)
+
+Return the Kick operator at time t.
+"""
+function (K::KickOperator)(t::Real)
+    K.U(t) * exp(1im * K.H_F * t)
+end
+
+
+
+end
diff --git a/src/Looping.jl b/src/Looping.jl
index eef1937..0a7677c 100644
--- a/src/Looping.jl
+++ b/src/Looping.jl
@@ -1,5 +1,8 @@
 module Looping
+export hi
 
-# Write your package code here.
+hi() = print("hi")
 
+include("Utilities.jl")
+include("FloquetUtils.jl")
 end
diff --git a/src/Utilities.jl b/src/Utilities.jl
new file mode 100644
index 0000000..5bd3a22
--- /dev/null
+++ b/src/Utilities.jl
@@ -0,0 +1,47 @@
+module Utilities
+
+export periodic_distance
+export restrict_to_range
+
+"""
+    periodic_distance(m, n, N)
+
+Calculate the distance `(m-n)`, but assuming periodic boundary
+conditions for a chain of length `N`.
+
+# Examples
+```jldoctest
+julia> periodic_difference(1,4,5)
+2
+julia> periodic_difference(1,3,5)
+-2
+```
+"""
+function periodic_distance(m::Integer, n::Integer, N::Integer)
+    diff = m - n
+    adiff = abs(diff)
+    if abs(diff + N) < adiff
+        diff + N
+    elseif abs(diff - N) < adiff
+        diff - N
+    else
+        diff
+    end
+end
+
+"""
+    restrict_to_range(x, border)
+
+Returns the value `x` restricted to the range `[-border, border]`.
+"""
+function restrict_to_range(x::Real, border::Real)
+    sign = x ≥ border ? -1 : 1
+
+    while abs(x) > border
+        x += sign * border * 2
+    end
+
+    x
+end
+
+end