diff --git a/Project.toml b/Project.toml
index f3f68857..3bc62f22 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,7 +25,7 @@ FFTW = "1"
 IterativeSolvers = "0.9"
 LinearMaps = "3"
 OrdinaryDiffEq = "5"
-QuantumOpticsBase = "0.2.7"
+QuantumOpticsBase = "0.3"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1.0"
 StochasticDiffEq = "6"
diff --git a/src/bloch_redfield_master.jl b/src/bloch_redfield_master.jl
index d239dd49..cb33f9cb 100644
--- a/src/bloch_redfield_master.jl
+++ b/src/bloch_redfield_master.jl
@@ -15,7 +15,7 @@ See QuTiP's documentation (http://qutip.org/docs/latest/guide/dynamics/dynamics-
 * `secular_cutoff=0.1`: Cutoff to allow a degree of partial secularization. Terms are discarded if they are greater than (dw\\_min * secular cutoff) where dw\\_min is the smallest (non-zero) difference between any two eigenenergies of H.
                         This argument is only taken into account if use_secular=true.
 """
-function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secular=true, secular_cutoff=0.1)
+function bloch_redfield_tensor(H::AbstractOperator, a_ops; J=SparseOpType[], use_secular=true, secular_cutoff=0.1)
 
     # Use the energy eigenbasis
     H_evals, transf_mat = eigen(Array(H.data)) #Array call makes sure H is a dense array
@@ -45,7 +45,7 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
     end
 
     #Transform interaction operators to Hamiltonian eigenbasis
-    A = Array{ComplexF64}(undef, N, N, K)
+    A = Array{eltype(transf_mat)}(undef, N, N, K)
     for k in 1:K
         A[:, :, k] = to_Heb(a_ops[k][1], transf_mat).data
     end
@@ -54,8 +54,8 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
     W = H_evals .- transpose(H_evals)
 
     #Array for spectral functions evaluated at transition frequencies
-    Jw = Array{Float64}(undef, N, N, K)
-    #Loop over all a_ops and calculate each spectral density at all transition frequencies
+    Jw = Array{eltype(transf_mat)}(undef, N, N, K)
+    # Jw = zeros(Complex{Float64}, N, N, K)
     for k in 1:K
         Jw[:, :, k] .= a_ops[k][2].(W)
     end
@@ -106,7 +106,6 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
 end #Function
 
 
-
 """
     timeevolution.master_bloch_redfield(tspan, rho0, R, H; <keyword arguments>)
 
@@ -127,9 +126,9 @@ Time-evolution according to a Bloch-Redfield master equation.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
 function master_bloch_redfield(tspan,
-        rho0::T, L::SuperOperator{Tuple{B,B},Tuple{B,B}},
+        rho0::Operator{B,B}, L::SuperOperator{Tuple{B,B},Tuple{B,B}},
         H::AbstractOperator{B,B}; fout::Union{Function,Nothing}=nothing,
-        kwargs...) where {B<:Basis,T<:Operator{B,B}}
+        kwargs...) where {B}
 
     #Prep basis transf
     evals, transf_mat = eigen(dense(H).data)
@@ -142,28 +141,28 @@ function master_bloch_redfield(tspan,
     L_ = isa(L, SparseSuperOpType) ? SparseOperator(basis_comp, L.data) : DenseOperator(basis_comp, L.data)
 
     # Derivative function
-    dmaster_br_(t, rho::T2, drho::T2) where T2<:Ket = dmaster_br(drho, rho, L_)
+    dmaster_br_(t, rho, drho) = dmaster_br(drho, rho, L_)
 
     return integrate_br(tspan, dmaster_br_, rho0_eb, transf_op, inv_transf_op, fout; kwargs...)
 end
 master_bloch_redfield(tspan, psi::Ket, args...; kwargs...) = master_bloch_redfield(tspan, dm(psi), args...; kwargs...)
 
 # Derivative ∂ₜρ = Lρ
-function dmaster_br(drho::T, rho::T, L::DataOperator{B,B}) where {B<:Basis,T<:Ket{B}}
+function dmaster_br(drho, rho, L)
     QuantumOpticsBase.mul!(drho,L,rho)
 end
 
 # Integrate if there is no fout specified
-function integrate_br(tspan, dmaster_br::Function, rho::T,
-                transf_op::T2, inv_transf_op::T2, ::Nothing;
-                kwargs...) where {T<:Ket,T2<:Operator}
+function integrate_br(tspan, dmaster_br, rho,
+                transf_op, inv_transf_op, ::Nothing;
+                kwargs...)
     # Pre-allocate for in-place back-transformation from eigenbasis
     rho_out = copy(transf_op)
     tmp = copy(transf_op)
     tmp2 = copy(transf_op)
 
     # Define fout
-    function fout(t, rho::T)
+    function fout(t, rho)
         tmp.data[:] = rho.data
         QuantumOpticsBase.mul!(tmp2,transf_op,tmp)
         QuantumOpticsBase.mul!(rho_out,tmp2,inv_transf_op)
@@ -174,9 +173,9 @@ function integrate_br(tspan, dmaster_br::Function, rho::T,
 end
 
 # Integrate with given fout
-function integrate_br(tspan, dmaster_br::Function, rho::T,
-                transf_op::T2, inv_transf_op::T2, fout::Function;
-                kwargs...) where {T<:Ket,T2<:Operator}
+function integrate_br(tspan, dmaster_br, rho,
+                transf_op, inv_transf_op, fout::Function;
+                kwargs...)
     # Pre-allocate for in-place back-transformation from eigenbasis
     rho_out = copy(transf_op)
     tmp = copy(transf_op)
@@ -185,7 +184,7 @@ function integrate_br(tspan, dmaster_br::Function, rho::T,
     tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
 
     # Perform back-transfomration before calling fout
-    function fout_(t, rho::T)
+    function fout_(t, rho)
         tmp.data[:] = rho.data
         QuantumOpticsBase.mul!(tmp2,transf_op,tmp)
         QuantumOpticsBase.mul!(rho_out,tmp2,inv_transf_op)
diff --git a/src/master.jl b/src/master.jl
index f65613a0..dd7673a9 100644
--- a/src/master.jl
+++ b/src/master.jl
@@ -1,5 +1,3 @@
-const DecayRates = Union{Vector, Matrix, Nothing}
-
 """
     timeevolution.master_h(tspan, rho0, H, J; <keyword arguments>)
 
@@ -7,14 +5,14 @@ Integrate the master equation with dmaster_h as derivative function.
 
 Further information can be found at [`master`](@ref).
 """
-function master_h(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector;
-                rates::DecayRates=nothing,
-                Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function master_h(tspan, rho0::Operator, H::AbstractOperator, J;
+                rates=nothing,
+                Jdagger=dagger.(J),
+                fout=nothing,
+                kwargs...)
     check_master(rho0, H, J, Jdagger, rates)
     tmp = copy(rho0)
-    dmaster_(t, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
+    dmaster_(t, rho, drho) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
@@ -29,15 +27,15 @@ H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
 ```
 Further information can be found at [`master`](@ref).
 """
-function master_nh(tspan, rho0::T, Hnh::AbstractOperator{B,B}, J::Vector;
-                rates::DecayRates=nothing,
+function master_nh(tspan, rho0::Operator, Hnh::AbstractOperator, J;
+                rates=nothing,
                 Hnhdagger::AbstractOperator=dagger(Hnh),
-                Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+                Jdagger=dagger.(J),
+                fout=nothing,
+                kwargs...)
     check_master(rho0, Hnh, J, Jdagger, rates)
     tmp = copy(rho0)
-    dmaster_(t, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
+    dmaster_(t, rho, drho) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
@@ -73,21 +71,21 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster.
         be changed.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function master(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector;
-                rates::DecayRates=nothing,
-                Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function master(tspan, rho0::Operator, H::AbstractOperator, J;
+                rates=nothing,
+                Jdagger=dagger.(J),
+                fout=nothing,
+                kwargs...)
     isreducible = check_master(rho0, H, J, Jdagger, rates)
     if !isreducible
         tmp = copy(rho0)
-        dmaster_h_(t, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
+        dmaster_h_(t, rho, drho) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
         return integrate_master(tspan, dmaster_h_, rho0, fout; kwargs...)
     else
         Hnh = nh_hamiltonian(H,J,Jdagger,rates)
         Hnhdagger = dagger(Hnh)
         tmp = copy(rho0)
-        dmaster_nh_(t, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
+        dmaster_nh_(t, rho, drho) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
         return integrate_master(tspan, dmaster_nh_, rho0, fout; kwargs...)
     end
 end
@@ -132,12 +130,12 @@ The given function can either be of the form `f(t, rho) -> (Hnh, Hnhdagger, J, J
 or `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger, rates)` For further information look
 at [`master_dynamic`](@ref).
 """
-function master_nh_dynamic(tspan, rho0::T, f::Function;
-                rates::DecayRates=nothing,
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function master_nh_dynamic(tspan, rho0::Operator, f;
+                rates=nothing,
+                fout=nothing,
+                kwargs...)
     tmp = copy(rho0)
-    dmaster_(t, rho::T, drho::T) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
+    dmaster_(t, rho, drho) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
@@ -166,22 +164,20 @@ operators:
         be changed.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function master_dynamic(tspan, rho0::T, f::Function;
-                rates::DecayRates=nothing,
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function master_dynamic(tspan, rho0::Operator, f;
+                rates=nothing,
+                fout=nothing,
+                kwargs...)
     tmp = copy(rho0)
-    dmaster_(t, rho::T, drho::T) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
+    dmaster_(t, rho, drho) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
 
 # Automatically convert Ket states to density operators
-master(tspan, psi0::Ket{B}, args...; kwargs...) where B<:Basis = master(tspan, dm(psi0), args...; kwargs...)
-master_h(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_h(tspan, dm(psi0), H, J; kwargs...)
-master_nh(tspan, psi0::Ket{B}, Hnh::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_nh(tspan, dm(psi0), Hnh, J; kwargs...)
-master_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_dynamic(tspan, dm(psi0), f; kwargs...)
-master_nh_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_nh_dynamic(tspan, dm(psi0), f; kwargs...)
+for f ∈ [:master,:master_h,:master_nh,:master_dynamic,:master_nh_dynamic]
+    @eval $f(tspan,psi0::Ket,args...;kwargs...) = $f(tspan,dm(psi0),args...;kwargs...)
+end
 
 # Non-hermitian Hamiltonian
 function nh_hamiltonian(H,J,Jdagger,::Nothing)
@@ -207,18 +203,16 @@ function nh_hamiltonian(H,J,Jdagger,rates::AbstractMatrix)
 end
 
 # Recasting needed for the ODE solver is just providing the underlying data
-function recast!(x::T, rho::Operator{B,B,T}) where {B<:Basis,T}
+function recast!(x::T, rho::Operator{B,B,T}) where {B,T}
     rho.data = x
 end
-recast!(rho::Operator{B,B,T}, x::T) where {B<:Basis,T} = nothing
+recast!(rho::Operator{B,B,T}, x::T) where {B,T} = nothing
 
-function integrate_master(tspan, df::Function, rho0,
-                        fout::Union{Nothing, Function}; kwargs...)
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
+function integrate_master(tspan, df, rho0, fout; kwargs...)
     x0 = rho0.data
     state = deepcopy(rho0)
     dstate = deepcopy(rho0)
-    integrate(tspan_, df, x0, state, dstate, fout; kwargs...)
+    integrate(tspan, df, x0, state, dstate, fout; kwargs...)
 end
 
 
@@ -231,9 +225,7 @@ end
 # the type of the given decay rate object which can either be nothing, a vector
 # or a matrix.
 
-function dmaster_h(rho::T, H::AbstractOperator{B,B},
-                    rates::Nothing, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h(rho, H, rates::Nothing, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
     for i=1:length(J)
@@ -248,9 +240,7 @@ function dmaster_h(rho::T, H::AbstractOperator{B,B},
     return drho
 end
 
-function dmaster_h(rho::T, H::AbstractOperator{B,B},
-                    rates::Vector, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h(rho, H, rates::AbstractVector, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
     for i=1:length(J)
@@ -265,9 +255,7 @@ function dmaster_h(rho::T, H::AbstractOperator{B,B},
     return drho
 end
 
-function dmaster_h(rho::T, H::AbstractOperator{B,B},
-                    rates::Matrix, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h(rho, H, rates::AbstractMatrix, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
     for j=1:length(J), i=1:length(J)
@@ -282,9 +270,7 @@ function dmaster_h(rho::T, H::AbstractOperator{B,B},
     return drho
 end
 
-function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
-                    rates::Nothing, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_nh(rho, Hnh, Hnh_dagger, rates::Nothing, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
     for i=1:length(J)
@@ -294,9 +280,7 @@ function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOper
     return drho
 end
 
-function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
-                    rates::Vector, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_nh(rho, Hnh, Hnh_dagger, rates::AbstractVector, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
     for i=1:length(J)
@@ -306,9 +290,7 @@ function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOper
     return drho
 end
 
-function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
-                    rates::Matrix, J::Vector, Jdagger::Vector,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_nh(rho, Hnh, Hnh_dagger, rates::AbstractMatrix, J, Jdagger, drho, tmp)
     QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
     QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
     for j=1:length(J), i=1:length(J)
@@ -323,9 +305,7 @@ function dmaster_liouville(rho,drho,L)
     return drho
 end
 
-function dmaster_h_dynamic(t, rho::T, f::Function,
-                    rates::DecayRates,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
     result = f(t, rho)
     QO_CHECKS[] && @assert 3 <= length(result) <= 4
     if length(result) == 3
@@ -338,9 +318,7 @@ function dmaster_h_dynamic(t, rho::T, f::Function,
     dmaster_h(rho, H, rates_, J, Jdagger, drho, tmp)
 end
 
-function dmaster_nh_dynamic(t, rho::T, f::Function,
-                    rates::DecayRates,
-                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
     result = f(t, rho)
     QO_CHECKS[] && @assert 4 <= length(result) <= 5
     if length(result) == 4
@@ -354,36 +332,35 @@ function dmaster_nh_dynamic(t, rho::T, f::Function,
 end
 
 
-function check_master(rho0::Operator{B,B}, H::AbstractOperator{B,B}, J::Vector, Jdagger::Vector, rates::DecayRates) where B<:Basis
-    # TODO: clean up type checks by dispatch; make type of J known
+function check_master(rho0, H, J, Jdagger, rates)
     isreducible = true # test if all operators are sparse or dense
     if !(isa(H, DenseOpType) || isa(H, SparseOpType))
         isreducible = false
     end
     for j=J
-        @assert isa(j, AbstractOperator{B,B})
+        @assert isa(j, AbstractOperator)
         if !(isa(j, DenseOpType) || isa(j, SparseOpType))
             isreducible = false
         end
         check_samebases(rho0, j)
     end
     for j=Jdagger
-        @assert isa(j, AbstractOperator{B,B})
+        @assert isa(j, AbstractOperator)
         if !(isa(j, DenseOpType) || isa(j, SparseOpType))
             isreducible = false
         end
         check_samebases(rho0, j)
     end
     @assert length(J)==length(Jdagger)
-    if isa(rates,Matrix)
+    if isa(rates,AbstractMatrix)
         @assert size(rates, 1) == size(rates, 2) == length(J)
-    elseif isa(rates,Vector)
+    elseif isa(rates,AbstractVector)
         @assert length(rates) == length(J)
     end
     isreducible
 end
 
 get_type(rho, H, rates::Nothing, J, Jdagger) = promote_type(eltype(rho),eltype(H),eltype.(J)...,eltype.(Jdagger)...)
-get_type(rho, H, rates::Union{Vector,Matrix}, J, Jdagger) = promote_type(eltype(rho),eltype(H),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)
+get_type(rho, H, rates::AbstractVecOrMat, J, Jdagger) = promote_type(eltype(rho),eltype(H),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)
 get_type(rho, Hnh, Hnhdagger, rates::Nothing, J, Jdagger) = promote_type(eltype(rho),eltype(Hnh),eltype(Hnhdagger), eltype.(J)...,eltype.(Jdagger)...)
-get_type(rho, Hnh, Hnhdagger, rates::Union{Vector,Matrix}, J, Jdagger) = promote_type(eltype(rho),eltype(Hnh),eltype(Hnhdagger),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)
+get_type(rho, Hnh, Hnhdagger, rates::AbstractVecOrMat, J, Jdagger) = promote_type(eltype(rho),eltype(Hnh),eltype(Hnhdagger),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)
diff --git a/src/mcwf.jl b/src/mcwf.jl
index 0de7ed2d..e7a6882e 100644
--- a/src/mcwf.jl
+++ b/src/mcwf.jl
@@ -10,15 +10,15 @@ Calculate MCWF trajectory where the Hamiltonian is given in hermitian form.
 
 For more information see: [`mcwf`](@ref)
 """
-function mcwf_h(tspan, psi0::T, H::AbstractOperator{B,B}, J::Vector;
-        seed=rand(UInt), rates::DecayRates=nothing,
-        fout=nothing, Jdagger::Vector=dagger.(J),
-        tmp::T=copy(psi0),
+function mcwf_h(tspan, psi0::Ket, H::AbstractOperator, J;
+        seed=rand(UInt), rates=nothing,
+        fout=nothing, Jdagger=dagger.(J),
+        tmp=copy(psi0),
         display_beforeevent=false, display_afterevent=false,
-        kwargs...) where {B<:Basis,T<:Ket{B}}
+        kwargs...)
     check_mcwf(psi0, H, J, Jdagger, rates)
-    f(t, psi::T, dpsi::T) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates)
-    j(rng, t, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates)
+    f(t, psi, dpsi) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates)
+    j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates)
     integrate_mcwf(f, j, tspan, psi0, seed, fout;
         display_beforeevent=display_beforeevent,
         display_afterevent=display_afterevent,
@@ -36,13 +36,13 @@ H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
 
 For more information see: [`mcwf`](@ref)
 """
-function mcwf_nh(tspan, psi0::T, Hnh::AbstractOperator{B,B}, J::Vector;
+function mcwf_nh(tspan, psi0::Ket, Hnh::AbstractOperator, J;
         seed=rand(UInt), fout=nothing,
         display_beforeevent=false, display_afterevent=false,
-        kwargs...) where {B<:Basis,T<:Ket{B}}
+        kwargs...)
     check_mcwf(psi0, Hnh, J, J, nothing)
-    f(t, psi::T, dpsi::T) = dmcwf_nh(psi, Hnh, dpsi)
-    j(rng, t, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, nothing)
+    f(t, psi, dpsi) = dmcwf_nh(psi, Hnh, dpsi)
+    j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, nothing)
     integrate_mcwf(f, j, tspan, psi0, seed, fout;
         display_beforeevent=display_beforeevent,
         display_afterevent=display_afterevent,
@@ -87,16 +87,16 @@ is returned. These correspond to the times at which a jump occured and the index
 of the jump operators with which the jump occured, respectively.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function mcwf(tspan, psi0::T, H::AbstractOperator{B,B}, J::Vector;
-        seed=rand(UInt), rates::DecayRates=nothing,
-        fout=nothing, Jdagger::Vector=dagger.(J),
+function mcwf(tspan, psi0::Ket, H::AbstractOperator, J;
+        seed=rand(UInt), rates=nothing,
+        fout=nothing, Jdagger=dagger.(J),
         display_beforeevent=false, display_afterevent=false,
-        kwargs...) where {B<:Basis,T<:Ket{B}}
+        kwargs...)
     isreducible = check_mcwf(psi0, H, J, Jdagger, rates)
     if !isreducible
         tmp = copy(psi0)
-        dmcwf_h_(t, psi::T, dpsi::T) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates)
-        j_h(rng, t, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates)
+        dmcwf_h_(t, psi, dpsi) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates)
+        j_h(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates)
         integrate_mcwf(dmcwf_h_, j_h, tspan, psi0, seed,
             fout;
             display_beforeevent=display_beforeevent,
@@ -113,8 +113,8 @@ function mcwf(tspan, psi0::T, H::AbstractOperator{B,B}, J::Vector;
                 Hnh -= complex(float(eltype(H)))(0.5im*rates[i])*Jdagger[i]*J[i]
             end
         end
-        dmcwf_nh_(t, psi::T, dpsi::T) = dmcwf_nh(psi, Hnh, dpsi)
-        j_nh(rng, t, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates)
+        dmcwf_nh_(t, psi, dpsi) = dmcwf_nh(psi, Hnh, dpsi)
+        j_nh(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates)
         integrate_mcwf(dmcwf_nh_, j_nh, tspan, psi0, seed,
             fout;
             display_beforeevent=display_beforeevent,
@@ -152,13 +152,13 @@ is returned. These correspond to the times at which a jump occured and the index
 of the jump operators with which the jump occured, respectively.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function mcwf_dynamic(tspan, psi0::T, f::Function;
-    seed=rand(UInt), rates::DecayRates=nothing,
+function mcwf_dynamic(tspan, psi0::Ket, f;
+    seed=rand(UInt), rates=nothing,
     fout=nothing, display_beforeevent=false, display_afterevent=false,
-    kwargs...) where {T<:Ket}
+    kwargs...)
     tmp = copy(psi0)
-    dmcwf_(t, psi::T, dpsi::T) = dmcwf_h_dynamic(t, psi, f, rates, dpsi, tmp)
-    j_(rng, t, psi::T, psi_new::T) = jump_dynamic(rng, t, psi, f, psi_new, rates)
+    dmcwf_(t, psi, dpsi) = dmcwf_h_dynamic(t, psi, f, rates, dpsi, tmp)
+    j_(rng, t, psi, psi_new) = jump_dynamic(rng, t, psi, f, psi_new, rates)
     integrate_mcwf(dmcwf_, j_, tspan, psi0, seed,
         fout;
         display_beforeevent=display_beforeevent,
@@ -173,12 +173,12 @@ Calculate MCWF trajectory where the dynamic Hamiltonian is given in non-hermitia
 
 For more information see: [`mcwf_dynamic`](@ref)
 """
-function mcwf_nh_dynamic(tspan, psi0::T, f::Function;
-    seed=rand(UInt), rates::DecayRates=nothing,
+function mcwf_nh_dynamic(tspan, psi0::Ket, f;
+    seed=rand(UInt), rates=nothing,
     fout=nothing, display_beforeevent=false, display_afterevent=false,
-    kwargs...) where T<:Ket
-    dmcwf_(t, psi::T, dpsi::T) = dmcwf_nh_dynamic(t, psi, f, dpsi)
-    j_(rng, t, psi::T, psi_new::T) = jump_dynamic(rng, t, psi, f, psi_new, rates)
+    kwargs...)
+    dmcwf_(t, psi, dpsi) = dmcwf_nh_dynamic(t, psi, f, dpsi)
+    j_(rng, t, psi, psi_new) = jump_dynamic(rng, t, psi, f, psi_new, rates)
     integrate_mcwf(dmcwf_, j_, tspan, psi0, seed,
         fout;
         display_beforeevent=display_beforeevent,
@@ -186,8 +186,8 @@ function mcwf_nh_dynamic(tspan, psi0::T, f::Function;
         kwargs...)
 end
 
-function dmcwf_h_dynamic(t, psi::T, f::Function, rates::DecayRates,
-                    dpsi::T, tmp::T) where T<:Ket
+function dmcwf_h_dynamic(t, psi, f, rates,
+                    dpsi, tmp)
     result = f(t, psi)
     QO_CHECKS[] && @assert 3 <= length(result) <= 4
     if length(result) == 3
@@ -200,7 +200,7 @@ function dmcwf_h_dynamic(t, psi::T, f::Function, rates::DecayRates,
     dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates_)
 end
 
-function dmcwf_nh_dynamic(t, psi::T, f::Function, dpsi::T) where T<:Ket
+function dmcwf_nh_dynamic(t, psi, f, dpsi)
     result = f(t, psi)
     QO_CHECKS[] && @assert 3 <= length(result) <= 4
     H, J, Jdagger = result[1:3]
@@ -208,7 +208,7 @@ function dmcwf_nh_dynamic(t, psi::T, f::Function, dpsi::T) where T<:Ket
     dmcwf_nh(psi, H, dpsi)
 end
 
-function jump_dynamic(rng, t, psi::T, f::Function, psi_new::T, rates::DecayRates) where T<:Ket
+function jump_dynamic(rng, t, psi, f, psi_new, rates)
     result = f(t, psi)
     QO_CHECKS[] && @assert 3 <= length(result) <= 4
     J = result[2]
@@ -241,14 +241,14 @@ Integrate a single Monte Carlo wave function trajectory.
         and therefore must not be changed.
 * `kwargs`: Further arguments are passed on to the ode solver.
 """
-function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
-                        psi0::T, seed, fout::Function;
+function integrate_mcwf(dmcwf, jumpfun, tspan,
+                        psi0, seed, fout::Function;
                         display_beforeevent=false, display_afterevent=false,
                         display_jumps=false,
                         save_everystep=false, callback=nothing,
                         saveat=tspan,
                         alg=OrdinaryDiffEq.DP5(),
-                        kwargs...) where T
+                        kwargs...)
 
     tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
     # Display before or after events
@@ -275,7 +275,7 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
         (t,i)->nothing
     end
 
-    function fout_(x::Vector, t, integrator)
+    function fout_(x, t, integrator)
         recast!(x, state)
         fout(t, state)
     end
@@ -291,14 +291,14 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
     cb = jump_callback(jumpfun, seed, scb, save_before!, save_after!, save_t_index, psi0)
     full_cb = OrdinaryDiffEq.CallbackSet(callback,cb,scb)
 
-    function df_(dx::D, x::D, p, t) where D
+    function df_(dx, x, p, t)
         recast!(x, state)
         recast!(dx, dstate)
         dmcwf(t, state, dstate)
         recast!(dstate, dx)
     end
 
-    prob = OrdinaryDiffEq.ODEProblem{true}(df_, as_vector(psi0),(tspan_[1],tspan_[end]))
+    prob = OrdinaryDiffEq.ODEProblem{true}(df_, as_vector(psi0), (tspan_[1],tspan_[end]))
 
     sol = OrdinaryDiffEq.solve(
                 prob,
@@ -316,17 +316,17 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
     end
 end
 
-function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
-                        psi0::T, seed, fout::Nothing;
-                        kwargs...) where T
-    function fout_(t, x::T)
+function integrate_mcwf(dmcwf, jumpfun, tspan,
+                        psi0, seed, fout::Nothing;
+                        kwargs...)
+    function fout_(t, x)
         return normalize(x)
     end
     integrate_mcwf(dmcwf, jumpfun, tspan, psi0, seed, fout_; kwargs...)
 end
 
-function jump_callback(jumpfun::Function, seed, scb, save_before!::Function,
-                        save_after!::Function, save_t_index::Function, psi0::Ket)
+function jump_callback(jumpfun, seed, scb, save_before!,
+                        save_after!, save_t_index, psi0)
 
     tmp = copy(psi0)
     psi_tmp = copy(psi0)
@@ -368,7 +368,7 @@ Default jump function.
 * `J`: List of jump operators.
 * `psi_new`: Result of jump.
 """
-function jump(rng, t, psi::T, J::Vector, psi_new::T, rates::Nothing) where T<:Ket
+function jump(rng, t, psi, J, psi_new, rates::Nothing)
     if length(J)==1
         QuantumOpticsBase.mul!(psi_new,J[1],psi,true,false)
         psi_new.data ./= norm(psi_new)
@@ -387,7 +387,7 @@ function jump(rng, t, psi::T, J::Vector, psi_new::T, rates::Nothing) where T<:Ke
     return i
 end
 
-function jump(rng, t, psi::T, J::Vector, psi_new::T, rates::Vector) where T<:Ket
+function jump(rng, t, psi, J, psi_new, rates::AbstractVector)
     if length(J)==1
         QuantumOpticsBase.mul!(psi_new,J[1],psi,eltype(psi)(sqrt(rates[1])),zero(eltype(psi)))
         psi_new.data ./= norm(psi_new)
@@ -412,8 +412,8 @@ Evaluate non-hermitian Schroedinger equation.
 The non-hermitian Hamiltonian is given in two parts - the hermitian part H and
 the jump operators J.
 """
-function dmcwf_h(psi::T, H::AbstractOperator{B,B},
-                 J::Vector, Jdagger::Vector, dpsi::T, tmp::T, rates::Nothing) where {B<:Basis,T<:Ket{B}}
+function dmcwf_h(psi, H,
+                 J, Jdagger, dpsi, tmp, rates::Nothing)
     QuantumOpticsBase.mul!(dpsi,H,psi,eltype(psi)(-im),zero(eltype(psi)))
     for i=1:length(J)
         QuantumOpticsBase.mul!(tmp,J[i],psi,true,false)
@@ -422,8 +422,8 @@ function dmcwf_h(psi::T, H::AbstractOperator{B,B},
     return dpsi
 end
 
-function dmcwf_h(psi::T, H::AbstractOperator{B,B},
-                 J::Vector, Jdagger::Vector, dpsi::T, tmp::T, rates::Vector) where {B<:Basis,T<:Ket{B}}
+function dmcwf_h(psi, H,
+                 J, Jdagger, dpsi, tmp, rates::AbstractVector)
     QuantumOpticsBase.mul!(dpsi,H,psi,eltype(psi)(-im),zero(eltype(psi)))
     for i=1:length(J)
         QuantumOpticsBase.mul!(tmp,J[i],psi,eltype(psi)(rates[i]),zero(eltype(psi)))
@@ -438,7 +438,7 @@ Evaluate non-hermitian Schroedinger equation.
 
 The given Hamiltonian is already the non-hermitian version.
 """
-function dmcwf_nh(psi::T, Hnh::AbstractOperator{B,B}, dpsi::T) where {B<:Basis,T<:Ket{B}}
+function dmcwf_nh(psi, Hnh, dpsi)
     QuantumOpticsBase.mul!(dpsi,Hnh,psi,eltype(psi)(-im),zero(eltype(psi)))
     return dpsi
 end
@@ -448,20 +448,20 @@ end
 
 Check input of mcwf.
 """
-function check_mcwf(psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector, Jdagger::Vector, rates::DecayRates) where B<:Basis
+function check_mcwf(psi0, H, J, Jdagger, rates)
     # TODO: replace type checks by dispatch; make types of J known
     isreducible = true
     if !(isa(H, DenseOpType) || isa(H, SparseOpType))
         isreducible = false
     end
     for j=J
-        @assert isa(j, AbstractOperator{B,B})
+        @assert isa(j, AbstractOperator)
         if !(isa(j, DenseOpType) || isa(j, SparseOpType))
             isreducible = false
         end
     end
     for j=Jdagger
-        @assert isa(j, AbstractOperator{B,B})
+        @assert isa(j, AbstractOperator)
         if !(isa(j, DenseOpType) || isa(j, SparseOpType))
             isreducible = false
         end
@@ -488,13 +488,13 @@ corresponding set of jump operators is calculated.
 * `rates`: Matrix of decay rates.
 * `J`: Vector of jump operators.
 """
-function diagonaljumps(rates::Matrix, J::Vector{T}) where {B<:Basis,T<:AbstractOperator{B,B}}
+function diagonaljumps(rates::AbstractMatrix, J)
     @assert length(J) == size(rates)[1] == size(rates)[2]
     d, v = eigen(rates)
     d, [sum([v[j, i]*J[j] for j=1:length(d)]) for i=1:length(d)]
 end
 
-function diagonaljumps(rates::Matrix, J::Vector{T}) where {B<:Basis,T<:Union{LazySum{B,B},LazyTensor{B,B},LazyProduct{B,B}}}
+function diagonaljumps(rates::AbstractMatrix, J::Vector{T}) where T<:Union{LazySum,LazyTensor,LazyProduct}
     @assert length(J) == size(rates)[1] == size(rates)[2]
     d, v = eigen(rates)
     d, [LazySum([v[j, i]*J[j] for j=1:length(d)]...) for i=1:length(d)]
diff --git a/src/phasespace.jl b/src/phasespace.jl
index 5bc192ba..1e9b58ce 100644
--- a/src/phasespace.jl
+++ b/src/phasespace.jl
@@ -16,9 +16,9 @@ function can either be evaluated on one point α or on a grid specified by
 the vectors `xvec` and `yvec`. Note that conversion from `x` and `y` to `α` is
 done via the relation ``α = \\frac{1}{\\sqrt{2}}(x + i y)``.
 """
-function qfunc(rho::AbstractOperator{B,B}, alpha::Number) where B<:FockBasis
+function qfunc(rho::AbstractOperator{B,B}, alpha) where B<:FockBasis
     b = basis(rho)
-    _qfunc_operator(rho, convert(ComplexF64, alpha), Ket(b), Ket(b))
+    _qfunc_operator(rho, alpha, Ket(b), Ket(b))
 end
 
 function qfunc(rho::AbstractOperator{B,B}, xvec::AbstractVector, yvec::AbstractVector) where B<:FockBasis
@@ -27,14 +27,14 @@ function qfunc(rho::AbstractOperator{B,B}, xvec::AbstractVector, yvec::AbstractV
     Ny = length(yvec)
     tmp1 = Ket(b)
     tmp2 = Ket(b)
-    result = Matrix{ComplexF64}(undef, Nx, Ny)
+    result = Matrix{eltype(rho)}(undef, Nx, Ny)
     @inbounds for j=1:Ny, i=1:Nx
         result[i, j] = _qfunc_operator(rho, complex(xvec[i], yvec[j])/sqrt(2), tmp1, tmp2)
     end
     result
 end
 
-function qfunc(psi::Ket{B}, alpha::Number) where B<:FockBasis
+function qfunc(psi::Ket{B}, alpha) where B<:FockBasis
     b = basis(psi)
     _conj_alpha = conj(alpha)
     c = one(alpha)
@@ -84,11 +84,11 @@ function qfunc(psi::Ket{B}, xvec::AbstractVector, yvec::AbstractVector) where B<
     return result
 end
 
-function qfunc(state::Union{Ket{B}, AbstractOperator{B,B}}, x::Number, y::Number) where B<:FockBasis
-    qfunc(state, ComplexF64(x, y)/sqrt(2))
+function qfunc(state::Union{Ket{B}, AbstractOperator{B,B}}, x, y) where B<:FockBasis
+    qfunc(state, complex(x, y)/sqrt(2))
 end
 
-function _qfunc_operator(rho::AbstractOperator{B,B}, alpha::Complex, tmp1::Ket, tmp2::Ket) where B<:FockBasis
+function _qfunc_operator(rho, alpha, tmp1, tmp2)
     coherentstate!(tmp1, basis(rho), alpha)
     QuantumOpticsBase.mul!(tmp2,rho,tmp1,true,false)
     a = dot(tmp1.data, tmp2.data)
@@ -106,11 +106,11 @@ function can either be evaluated on one point α or on a grid specified by
 the vectors `xvec` and `yvec`. Note that conversion from `x` and `y` to `α` is
 done via the relation ``α = \\frac{1}{\\sqrt{2}}(x + i y)``.
 """
-function wigner(rho::Operator{B,B}, x::Number, y::Number) where B<:FockBasis
+function wigner(rho::Operator{B,B}, x, y) where B<:FockBasis
     b = basis(rho)
-    N = b.N::Int
+    N = b.N
     N0 = b.offset
-    _2α = complex(convert(Float64, x), convert(Float64, y))*sqrt(2)
+    _2α = complex(x, y)*sqrt(2)
     abs2_2α = abs2(_2α)
     w = complex(0.)
     coefficient = complex(0.)
@@ -125,7 +125,7 @@ end
 
 function wigner(rho::Operator{B,B}, xvec::AbstractVector, yvec::AbstractVector) where B<:FockBasis
     b = basis(rho)
-    N = b.N::Int
+    N = b.N
     N0 = b.offset
     _2α = [complex(x, y)*sqrt(2) for x=xvec, y=yvec]
     abs2_2α = abs2.(_2α)
@@ -287,10 +287,10 @@ systems and can be described by two angles θ and ϕ (although this
 parametrization is not unique), similarly to a qubit on the
 Bloch sphere.
 """
-function coherentspinstate(b::SpinBasis, theta::Real, phi::Real,
-    result = Ket(b, Vector{ComplexF64}(undef, length(b))))
+function coherentspinstate(b::SpinBasis, theta::Real, phi::Real)
+    result = Ket(b)
     data = result.data
-    N = BigInt(length(b)-1)
+    N = length(b)-1
     sinth = sin(0.5theta)
     costh = cos(0.5theta)
     expphi = exp(0.5im*phi)
@@ -316,21 +316,21 @@ resolution `(Ntheta, Nphi)`.
 
 This version calculates the Husimi Q SU(2) function at a position given by θ and ϕ.
 """
-function qfuncsu2(psi::Ket{B}, Ntheta::Int; Nphi::Int=2Ntheta) where B<:SpinBasis
+function qfuncsu2(psi::Ket{B}, Ntheta::Integer; Nphi::Integer=2Ntheta) where B<:SpinBasis
     b = psi.basis
     psi_bra_data = psi.data'
     lb = float(b.spinnumber)
-    result = Array{Float64}(undef, Ntheta,Nphi)
+    result = Array{real(eltype(psi))}(undef, Ntheta,Nphi)
     @inbounds  for i = 0:Ntheta-1, j = 0:Nphi-1
         result[i+1,j+1] = (2*lb+1)/(4pi)*abs2(psi_bra_data*coherentspinstate(b,pi-i*pi/(Ntheta-1),j*2pi/(Nphi-1)-pi).data)
     end
     return result
 end
 
-function qfuncsu2(rho::Operator{B,B}, Ntheta::Int; Nphi::Int=2Ntheta) where B<:SpinBasis
+function qfuncsu2(rho::Operator{B,B}, Ntheta::Integer; Nphi::Integer=2Ntheta) where B<:SpinBasis
     b = basis(rho)
     lb = float(b.spinnumber)
-    result = Array{Float64}(undef, Ntheta,Nphi)
+    result = Array{real(eltype(rho))}(undef, Ntheta,Nphi)
     @inbounds  for i = 0:Ntheta-1, j = 0:Nphi-1
         c = coherentspinstate(b,pi-i*1pi/(Ntheta-1),j*2pi/(Nphi-1)-pi)
         result[i+1,j+1] = abs((2*lb+1)/(4pi)*c.data'*rho.data*c.data)
@@ -374,7 +374,7 @@ function wignersu2(rho::Operator{B,B}, theta::Real, phi::Real) where B<:SpinBasi
     N = length(basis(rho))-1
 
     ### Tensor generation ###
-    BandT = Array{Vector{Float64}}(undef, N,N+1)
+    BandT = Array{Vector{real(eltype(rho))}}(undef, N,N+1)
     BandT[1,1] = collect(range(-N/2, stop=N/2, length=N+1))
     BandT[1,2] = -collect(sqrt.(range(1, stop=N, length=N)).*sqrt.(range((N)/2, stop=1/2, length=N)))
     BandT[2,1] = clebschgordan(1,0,1,0,2,0)*BandT[1,1].*BandT[1,1] -
@@ -419,12 +419,12 @@ function wignersu2(rho::Operator{B,B}, theta::Real, phi::Real) where B<:SpinBasi
     return wignermap*sqrt((N+1)/(4pi))
 end
 
-function wignersu2(rho::Operator{B,B}, Ntheta::Int; Nphi::Int=2Ntheta) where B<:SpinBasis
+function wignersu2(rho::Operator{B,B}, Ntheta::Integer; Nphi::Integer=2Ntheta) where B<:SpinBasis
 
     N = length(basis(rho))-1
 
     ### Tensor generation ###
-    BandT = Array{Vector{Float64}}(undef, N,N+1)
+    BandT = Array{Vector{real(eltype(rho))}}(undef, N,N+1)
     BandT[1,1] = collect(range(-N/2, stop=N/2, length=N+1))
     BandT[1,2] = -collect(sqrt.(range(1, stop=N, length=N)).*sqrt.(range((N)/2, stop=1/2, length=N)))
     BandT[2,1] = clebschgordan(1,0,1,0,2,0)*BandT[1,1].*BandT[1,1] -
@@ -537,14 +537,10 @@ function ylm(l::Integer,m::Integer,theta::Real,phi::Real)
     end
 end
 
-function _calc_ylm_norm(l::Int, m_::Int)
-    # TODO: Clean up Int types
-    if 0 < l+m_ <= 33
-        norm = @. Float64(sqrt(4pi)/sqrt(2*l+1)*sqrt(factorial(Int128(l-m_)))/sqrt(factorial(Int128(l+m_))))
-    elseif 0 < l-m_ <= 33 && l+m_ > 33
-        norm = @. Float64(sqrt(4pi)/sqrt(2*l+1)*sqrt(factorial(Int128(l-m_)))/sqrt(factorial(BigInt(l+m_))))
-    else
-        norm = @. Float64(sqrt(4pi)/sqrt(2*l+1)*sqrt(factorial(BigInt(l-m_)))/sqrt(factorial(BigInt(l+m_))))
+function _calc_ylm_norm(l, m)
+    n = sqrt(4pi/(2l+1))
+    for k=(-m+1):m
+        n /= sqrt(l + k)
     end
-    norm
+    return n
 end
diff --git a/src/schroedinger.jl b/src/schroedinger.jl
index 9aeeb0cc..ef284937 100644
--- a/src/schroedinger.jl
+++ b/src/schroedinger.jl
@@ -14,13 +14,12 @@ Integrate Schroedinger equation to evolve states or compute propagators.
 """
 function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B};
                 fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Union{AbstractOperator{B,B},StateVector{B}}}
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
-    dschroedinger_(t, psi::T, dpsi::T) = dschroedinger(psi, H, dpsi)
+                kwargs...) where {B,T<:Union{AbstractOperator{B,B},StateVector{B}}}
+    dschroedinger_(t, psi, dpsi) = dschroedinger(psi, H, dpsi)
     x0 = psi0.data
     state = copy(psi0)
     dstate = copy(psi0)
-    integrate(tspan_, dschroedinger_, x0, state, dstate, fout; kwargs...)
+    integrate(tspan, dschroedinger_, x0, state, dstate, fout; kwargs...)
 end
 
 
@@ -38,45 +37,44 @@ Integrate time-dependent Schroedinger equation to evolve states or compute propa
         normalized nor permanent! It is still in use by the ode solver and
         therefore must not be changed.
 """
-function schroedinger_dynamic(tspan, psi0::T, f::Function;
+function schroedinger_dynamic(tspan, psi0, f::Function;
                 fout::Union{Function,Nothing}=nothing,
-                kwargs...) where T<:Union{StateVector,AbstractOperator}
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
-    dschroedinger_(t, psi::T, dpsi::T) = dschroedinger_dynamic(t, psi, f, dpsi)
+                kwargs...)
+    dschroedinger_(t, psi, dpsi) = dschroedinger_dynamic(t, psi, f, dpsi)
     x0 = psi0.data
     state = copy(psi0)
     dstate = copy(psi0)
-    integrate(tspan_, dschroedinger_, x0, state, dstate, fout; kwargs...)
+    integrate(tspan, dschroedinger_, x0, state, dstate, fout; kwargs...)
 end
 
 
-recast!(x::D, psi::StateVector{B,D}) where {B<:Basis, D} = (psi.data = x);
-recast!(psi::StateVector{B,D}, x::D) where {B<:Basis, D} = nothing
+recast!(x::D, psi::StateVector{B,D}) where {B, D} = (psi.data = x);
+recast!(psi::StateVector{B,D}, x::D) where {B, D} = nothing
 
 
-function dschroedinger(psi::Union{Ket{B},AbstractOperator{B,B}}, H::AbstractOperator{B,B}, dpsi::Union{Ket{B},AbstractOperator{B,B}}) where B<:Basis
+function dschroedinger(psi, H, dpsi)
     QuantumOpticsBase.mul!(dpsi,H,psi,eltype(psi)(-im),zero(eltype(psi)))
     return dpsi
 end
 
-function dschroedinger(psi::Bra{B}, H::AbstractOperator{B,B}, dpsi::Bra{B}) where B<:Basis
+function dschroedinger(psi::Bra, H, dpsi)
     QuantumOpticsBase.mul!(dpsi,psi,H,eltype(psi)(im),zero(eltype(psi)))
     return dpsi
 end
 
 
-function dschroedinger_dynamic(t, psi0::T, f::Function, dpsi::T) where T<:Union{StateVector,AbstractOperator}
+function dschroedinger_dynamic(t, psi0, f::Function, dpsi)
     H = f(t, psi0)
     dschroedinger(psi0, H, dpsi)
 end
 
 
-function check_schroedinger(psi::Ket, H::AbstractOperator)
+function check_schroedinger(psi::Ket, H)
     check_multiplicable(H, psi)
     check_samebases(H)
 end
 
-function check_schroedinger(psi::Bra, H::AbstractOperator)
+function check_schroedinger(psi::Bra, H)
     check_multiplicable(psi, H)
     check_samebases(H)
 end
diff --git a/src/semiclassical.jl b/src/semiclassical.jl
index ca10a072..a367ff3b 100644
--- a/src/semiclassical.jl
+++ b/src/semiclassical.jl
@@ -16,7 +16,6 @@ using ..timeevolution
 
 
 const QuantumState{B} = Union{Ket{B}, Operator{B,B}}
-const DecayRates = Union{Nothing, Vector, Matrix}
 
 """
 Semi-classical state.
@@ -24,10 +23,10 @@ Semi-classical state.
 It consists of a quantum part, which is either a `Ket` or a `DenseOperator` and
 a classical part that is specified as a complex vector of arbitrary length.
 """
-mutable struct State{B<:Basis,T<:QuantumState{B},C<:Vector}
+mutable struct State{B,T,C}
     quantum::T
     classical::C
-    function State(quantum::T, classical::C) where {B<:Basis,T<:QuantumState{B},C<:Vector}
+    function State(quantum::T, classical::C) where {B,T<:QuantumState{B},C}
         new{B,T,C}(quantum, classical)
     end
 end
@@ -35,8 +34,8 @@ end
 Base.length(state::State) = length(state.quantum) + length(state.classical)
 Base.copy(state::State) = State(copy(state.quantum), copy(state.classical))
 Base.eltype(state::State) = promote_type(eltype(state.quantum),eltype(state.classical))
-normalize!(state::State{B,T}) where {B,T<:Ket} = normalize!(state.quantum)
-normalize(state::T) where {B,K<:Ket,T<:State{B,K}} = State(normalize(state.quantum),copy(state.classical))
+normalize!(state::State) where {B,T} = (normalize!(state.quantum); state)
+normalize(state::State) = State(normalize(state.quantum),copy(state.classical))
 
 function ==(a::State, b::State)
     QuantumOpticsBase.samebases(a.quantum, b.quantum) &&
@@ -47,9 +46,9 @@ end
 
 QuantumOpticsBase.expect(op, state::State) = expect(op, state.quantum)
 QuantumOpticsBase.variance(op, state::State) = variance(op, state.quantum)
-QuantumOpticsBase.ptrace(state::State, indices::Vector{Int}) = State(ptrace(state.quantum, indices), state.classical)
+QuantumOpticsBase.ptrace(state::State, indices) = State(ptrace(state.quantum, indices), state.classical)
 
-QuantumOpticsBase.dm(x::State{B,T}) where {B<:Basis,T<:Ket{B}} = State(dm(x.quantum), x.classical)
+QuantumOpticsBase.dm(x::State) = State(dm(x.quantum), x.classical)
 
 
 """
@@ -71,16 +70,15 @@ Integrate time-dependent Schrödinger equation coupled to a classical system.
         normalized nor permanent!
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function schroedinger_dynamic(tspan, state0::S, fquantum::Function, fclassical::Function;
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
-    dschroedinger_(t, state::S, dstate::S) = dschroedinger_dynamic(t, state, fquantum, fclassical, dstate)
+function schroedinger_dynamic(tspan, state0::State, fquantum, fclassical;
+                fout=nothing,
+                kwargs...)
+    dschroedinger_(t, state, dstate) = dschroedinger_dynamic(t, state, fquantum, fclassical, dstate)
     x0 = Vector{eltype(state0)}(undef, length(state0))
     recast!(state0, x0)
     state = copy(state0)
     dstate = copy(state0)
-    integrate(tspan_, dschroedinger_, x0, state, dstate, fout; kwargs...)
+    integrate(tspan, dschroedinger_, x0, state, dstate, fout; kwargs...)
 end
 
 """
@@ -103,22 +101,19 @@ Integrate time-dependent master equation coupled to a classical system.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
 function master_dynamic(tspan, state0::State{B,T}, fquantum, fclassical;
-                rates::DecayRates=nothing,
-                fout::Union{Function,Nothing}=nothing,
-                tmp::T=copy(state0.quantum),
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
-    function dmaster_(t, state::S, dstate::S) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
-        dmaster_h_dynamic(t, state, fquantum, fclassical, rates, dstate, tmp)
-    end
+                rates=nothing,
+                fout=nothing,
+                tmp=copy(state0.quantum),
+                kwargs...) where {B,T<:Operator}
+    dmaster_(t, state, dstate) = dmaster_h_dynamic(t, state, fquantum, fclassical, rates, dstate, tmp)
     x0 = Vector{eltype(state0)}(undef, length(state0))
     recast!(state0, x0)
     state = copy(state0)
     dstate = copy(state0)
-    integrate(tspan_, dmaster_, x0, state, dstate, fout; kwargs...)
+    integrate(tspan, dmaster_, x0, state, dstate, fout; kwargs...)
 end
 
-function master_dynamic(tspan, state0::State{B,T}, fquantum, fclassical; kwargs...) where {B<:Basis,T<:Ket{B}}
+function master_dynamic(tspan, state0::State{B,T}, fquantum, fclassical; kwargs...) where {B,T<:Ket{B}}
     master_dynamic(tspan, dm(state0), fquantum, fclassical; kwargs...)
 end
 
@@ -155,57 +150,51 @@ sure to take this into account when computing expectation values!
 """
 function mcwf_dynamic(tspan, psi0::State{B,T}, fquantum, fclassical, fjump_classical;
                 seed=rand(UInt),
-                rates::DecayRates=nothing,
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Ket{B}}
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
+                rates=nothing,
+                fout=nothing,
+                kwargs...) where {B,T<:Ket}
     tmp=copy(psi0.quantum)
-    function dmcwf_(t, psi::S, dpsi::S) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
-        dmcwf_h_dynamic(t, psi, fquantum, fclassical, rates, dpsi, tmp)
-    end
+    dmcwf_(t, psi, dpsi) = dmcwf_h_dynamic(t, psi, fquantum, fclassical, rates, dpsi, tmp)
     j_(rng, t, psi, psi_new) = jump_dynamic(rng, t, psi, fquantum, fclassical, fjump_classical, psi_new, rates)
     x0 = Vector{eltype(psi0)}(undef, length(psi0))
     recast!(psi0, x0)
     psi = copy(psi0)
     dpsi = copy(psi0)
-    integrate_mcwf(dmcwf_, j_, tspan_, psi, seed, fout; kwargs...)
+    integrate_mcwf(dmcwf_, j_, tspan, psi, seed, fout; kwargs...)
 end
 
-function recast!(state::State{B,T,C}, x::C) where {B<:Basis,T<:QuantumState{B},C<:Vector}
+function recast!(state::State{B,T,C}, x::C) where {B,T,C}
     N = length(state.quantum)
     copyto!(x, 1, state.quantum.data, 1, N)
     copyto!(x, N+1, state.classical, 1, length(state.classical))
     x
 end
 
-function recast!(x::C, state::State{B,T,C}) where {B<:Basis,T<:QuantumState{B},C<:Vector}
+function recast!(x::C, state::State{B,T,C}) where {B,T,C}
     N = length(state.quantum)
     copyto!(state.quantum.data, 1, x, 1, N)
     copyto!(state.classical, 1, x, N+1, length(state.classical))
 end
 
-function dschroedinger_dynamic(t, state::State{B,T}, fquantum::Function,
-            fclassical::Function, dstate::State{B,T}) where {B<:Basis,T<:Ket{B}}
+function dschroedinger_dynamic(t, state, fquantum, fclassical, dstate)
     fquantum_(t, psi) = fquantum(t, state.quantum, state.classical)
     timeevolution.dschroedinger_dynamic(t, state.quantum, fquantum_, dstate.quantum)
     fclassical(t, state.quantum, state.classical, dstate.classical)
 end
 
-function dmaster_h_dynamic(t, state::State{B,T}, fquantum::Function,
-            fclassical::Function, rates::DecayRates, dstate::State{B,T}, tmp::T) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_h_dynamic(t, state, fquantum, fclassical, rates, dstate, tmp)
     fquantum_(t, rho) = fquantum(t, state.quantum, state.classical)
     timeevolution.dmaster_h_dynamic(t, state.quantum, fquantum_, rates, dstate.quantum, tmp)
     fclassical(t, state.quantum, state.classical, dstate.classical)
 end
 
-function dmcwf_h_dynamic(t, psi::T, fquantum::Function, fclassical::Function, rates::DecayRates,
-                    dpsi::T, tmp::K) where {T,K}
+function dmcwf_h_dynamic(t, psi, fquantum, fclassical, rates, dpsi, tmp)
     fquantum_(t, rho) = fquantum(t, psi.quantum, psi.classical)
     timeevolution.dmcwf_h_dynamic(t, psi.quantum, fquantum_, rates, dpsi.quantum, tmp)
     fclassical(t, psi.quantum, psi.classical, dpsi.classical)
 end
 
-function jump_dynamic(rng, t, psi::T, fquantum::Function, fclassical::Function, fjump_classical::Function, psi_new::T, rates::DecayRates) where T<:State
+function jump_dynamic(rng, t, psi, fquantum, fclassical, fjump_classical, psi_new, rates)
     result = fquantum(t, psi.quantum, psi.classical)
     QO_CHECKS[] && @assert 3 <= length(result) <= 4
     J = result[2]
@@ -220,8 +209,8 @@ function jump_dynamic(rng, t, psi::T, fquantum::Function, fclassical::Function,
     return i
 end
 
-function jump_callback(jumpfun::Function, seed, scb, save_before!::Function,
-                        save_after!::Function, save_t_index::Function, psi0::State)
+function jump_callback(jumpfun, seed, scb, save_before!,
+                        save_after!, save_t_index, psi0::State)
     tmp = copy(psi0)
     psi_tmp = copy(psi0)
     rng = MersenneTwister(convert(UInt, seed))
@@ -248,7 +237,7 @@ function jump_callback(jumpfun::Function, seed, scb, save_before!::Function,
     return OrdinaryDiffEq.ContinuousCallback(djumpnorm,dojump,
                      save_positions = (false,false))
 end
-as_vector(psi::State{B,K}) where {B,K<:Ket} = [psi.quantum.data; psi.classical]
+as_vector(psi::State) = vcat(psi.quantum.data[:], psi.classical)
 
 
 end # module
diff --git a/src/steadystate.jl b/src/steadystate.jl
index 424fd3b9..5d8e2963 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -20,15 +20,11 @@ Calculate steady state using long time master equation evolution.
         density operator `rho` is further used and therefore must not be changed.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function master(H::AbstractOperator{B,B}, J::Vector;
-                rho0::Operator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
+function master(H::AbstractOperator, J;
+                rho0::Operator=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
                 tol=1e-3,
-                rates::Union{Vector, Matrix, Nothing}=nothing,
-                Jdagger::Vector=dagger.(J),
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where B<:Basis
-    t,u = timeevolution.master([0., Inf], rho0, H, J; rates=rates, Jdagger=Jdagger,
-                        fout=fout,
+                kwargs...)
+    t,u = timeevolution.master([0., Inf], rho0, H, J;
                         steady_state = true,
                         tol = tol, kwargs...)
 end
@@ -78,7 +74,7 @@ function liouvillianspectrum(L::SparseSuperOpType; nev::Int = min(10, length(L.b
     return d[indices], ops
 end
 
-liouvillianspectrum(H::AbstractOperator{B,B}, J::Vector; rates::Union{Vector, Matrix}=ones(Float64, length(J)), kwargs...) where B<:Basis = liouvillianspectrum(liouvillian(H, J; rates=rates); kwargs...)
+liouvillianspectrum(H, J; rates=ones(length(J)), kwargs...) = liouvillianspectrum(liouvillian(H, J; rates=rates); kwargs...)
 
 """
     steadystate.eigenvector(L)
@@ -108,4 +104,4 @@ function eigenvector(L::SuperOperator; tol::Real = 1e-9, nev::Int = 2, which::Sy
     return ops[1]/tr(ops[1])
 end
 
-eigenvector(H::AbstractOperator, J::Vector; rates::Union{Vector, Matrix}=ones(Float64, length(J)), kwargs...) = eigenvector(liouvillian(H, J; rates=rates); kwargs...)
+eigenvector(H, J; rates=ones(length(J)), kwargs...) = eigenvector(liouvillian(H, J; rates=rates); kwargs...)
\ No newline at end of file
diff --git a/src/stochastic_base.jl b/src/stochastic_base.jl
index 457ebabe..0e783da3 100644
--- a/src/stochastic_base.jl
+++ b/src/stochastic_base.jl
@@ -1,41 +1,38 @@
 using QuantumOpticsBase
 using QuantumOpticsBase: check_samebases, check_multiplicable
-import ..timeevolution: recast!, QO_CHECKS, pure_inference
+import ..timeevolution: recast!, QO_CHECKS, pure_inference, as_vector
 
 import DiffEqCallbacks, StochasticDiffEq, OrdinaryDiffEq
 
-const DiffArray{T} = Union{AbstractArray{T,1}, AbstractArray{T, 2}}
-
-
 """
-    integrate_stoch(tspan::Vector, df::Function, dg::Vector{Function}, x0::Vector{ComplexF64},
-            state::T, dstate::T, fout::Function; kwargs...)
+    integrate_stoch(tspan, df::Function, dg{Function}, x0{ComplexF64},
+            state, dstate, fout::Function; kwargs...)
 
 Integrate using StochasticDiffEq
 """
-function integrate_stoch(tspan::Vector, df::Function, dg::Function, x0::Vector,
-            state::T, dstate::T, fout::Function, n::Int;
+function integrate_stoch(tspan, df, dg, x0,
+            state, dstate, fout, n;
             save_everystep = false, callback=nothing, saveat=tspan,
-            alg::StochasticDiffEq.StochasticDiffEqAlgorithm=StochasticDiffEq.EM(),
+            alg=StochasticDiffEq.EM(),
             noise_rate_prototype = nothing,
             noise_prototype_classical = nothing,
             noise=nothing,
             ncb=nothing,
-            kwargs...) where T
+            kwargs...)
 
-    function df_(dx::T, x::T, p, t) where T
+    function df_(dx, x, p, t)
         recast!(x, state)
         recast!(dx, dstate)
         df(t, state, dstate)
         recast!(dstate, dx)
     end
 
-    function dg_(dx, x, p, t) where T
+    function dg_(dx, x, p, t)
         recast!(x, state)
         dg(dx, t, state, dstate, n)
     end
 
-    function fout_(x::Vector, t, integrator)
+    function fout_(x, t, integrator)
         recast!(x, state)
         fout(t, state)
     end
@@ -87,9 +84,9 @@ end
 
 Define fout if it was omitted.
 """
-function integrate_stoch(tspan::Vector, df::Function, dg::Function, x0::Vector,
-    state::T, dstate::T, ::Nothing, n::Int; kwargs...) where T
-    function fout(t, state::T)
+function integrate_stoch(tspan, df, dg, x0,
+    state, dstate, ::Nothing, n; kwargs...)
+    function fout(t, state)
         copy(state)
     end
     integrate_stoch(tspan, df, dg, x0, state, dstate, fout, n; kwargs...)
diff --git a/src/stochastic_master.jl b/src/stochastic_master.jl
index 3b96d372..2840ee58 100644
--- a/src/stochastic_master.jl
+++ b/src/stochastic_master.jl
@@ -1,7 +1,5 @@
 import ...timeevolution: dmaster_h, dmaster_nh, dmaster_h_dynamic, check_master
 
-const DecayRates = Union{Vector, Matrix, Nothing}
-
 """
     stochastic.master(tspan, rho0, H, J, C; <keyword arguments>)
 
@@ -32,31 +30,30 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
 function master(tspan, rho0::T, H::AbstractOperator{B,B},
-                J::Vector, C::Vector;
-                rates::DecayRates=nothing,
-                Jdagger::Vector=dagger.(J), Cdagger::Vector=dagger.(C),
-                fout::Union{Function,Nothing}=nothing,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+                J, C;
+                rates=nothing,
+                Jdagger=dagger.(J), Cdagger=dagger.(C),
+                fout=nothing,
+                kwargs...) where {B,T<:Operator{B,B}}
 
     tmp = copy(rho0)
 
     n = length(C)
 
-    dmaster_stoch(dx::DiffArray, t, rho::T, drho::T, n::Int) =
-            dmaster_stochastic(dx, rho, C, Cdagger, drho, n)
+    dmaster_stoch(dx, t, rho, drho, n) = dmaster_stochastic(dx, rho, C, Cdagger, drho, n)
 
     isreducible = check_master(rho0, H, J, Jdagger, rates) && check_master_stoch(rho0, C, Cdagger)
     if !isreducible
-        dmaster_h_determ(t, rho::T, drho::T) =
+        dmaster_h_determ(t, rho, drho) =
             dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
         integrate_master_stoch(tspan, dmaster_h_determ, dmaster_stoch, rho0, fout, n; kwargs...)
     else
         Hnh = copy(H)
-        if isa(rates, Matrix)
+        if isa(rates, AbstractMatrix)
             for i=1:length(J), j=1:length(J)
                 Hnh -= complex(float(eltype(H)))(0.5im*rates[i,j])*Jdagger[i]*J[j]
             end
-        elseif isa(rates, Vector)
+        elseif isa(rates, AbstractVector)
             for i=1:length(J)
                 Hnh -= complex(float(eltype(H)))(0.5im*rates[i])*Jdagger[i]*J[i]
             end
@@ -67,7 +64,7 @@ function master(tspan, rho0::T, H::AbstractOperator{B,B},
         end
         Hnhdagger = dagger(Hnh)
 
-        dmaster_nh_determ(t, rho::T, drho::T) =
+        dmaster_nh_determ(t, rho, drho) =
             dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
         integrate_master_stoch(tspan, dmaster_nh_determ, dmaster_stoch, rho0, fout, n; kwargs...)
     end
@@ -103,11 +100,11 @@ dynamic Hamiltonian and J.
         number if you want to avoid an initial calculation of function outputs!
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function master_dynamic(tspan, rho0::T, fdeterm::Function, fstoch::Function;
-                rates::DecayRates=nothing,
-                fout::Union{Function,Nothing}=nothing,
+function master_dynamic(tspan, rho0::T, fdeterm, fstoch;
+                rates=nothing,
+                fout=nothing,
                 noise_processes::Int=0,
-                kwargs...) where {B<:Basis,T<:Operator{B,B}}
+                kwargs...) where {B,T<:Operator{B,B}}
 
     tmp = copy(rho0)
 
@@ -118,23 +115,20 @@ function master_dynamic(tspan, rho0::T, fdeterm::Function, fstoch::Function;
         n = noise_processes
     end
 
-    dmaster_determ(t, rho::T, drho::T) = dmaster_h_dynamic(t, rho, fdeterm, rates, drho, tmp)
-    dmaster_stoch(dx::DiffArray, t, rho::T, drho::T, n::Int) =
-        dmaster_stoch_dynamic(dx, t, rho, fstoch, drho, n)
+    dmaster_determ(t, rho, drho) = dmaster_h_dynamic(t, rho, fdeterm, rates, drho, tmp)
+    dmaster_stoch(dx, t, rho, drho, n) = dmaster_stoch_dynamic(dx, t, rho, fstoch, drho, n)
     integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n; kwargs...)
 end
 master_dynamic(tspan, psi0::Ket, args...; kwargs...) = master_dynamic(tspan, dm(psi0), args...; kwargs...)
 
 # Derivative functions
-function dmaster_stochastic(dx::Vector, rho::T,
-            C::Vector, Cdagger::Vector, drho::T, ::Int) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_stochastic(dx::AbstractVector, rho, C, Cdagger, drho, n)
     recast!(dx, drho)
     QuantumOpticsBase.mul!(drho,C[1],rho)
     QuantumOpticsBase.mul!(drho,rho,Cdagger[1],true,true)
     drho.data .-= tr(drho)*rho.data
 end
-function dmaster_stochastic(dx::Matrix, rho::T,
-            C::Vector, Cdagger::Vector, drho::T, n::Int) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_stochastic(dx::AbstractMatrix, rho, C, Cdagger, drho, n)
     for i=1:n
         dx_i = @view dx[:, i]
         recast!(dx_i, drho)
@@ -145,8 +139,7 @@ function dmaster_stochastic(dx::Matrix, rho::T,
     end
 end
 
-function dmaster_stoch_dynamic(dx::DiffArray, t, rho::T,
-            f::Function, drho::T, n::Int) where {B<:Basis,T<:Operator{B,B}}
+function dmaster_stoch_dynamic(dx, t, rho, f, drho, n)
     result = f(t, rho)
     QO_CHECKS[] && @assert 2 == length(result)
     C, Cdagger = result
@@ -154,18 +147,18 @@ function dmaster_stoch_dynamic(dx::DiffArray, t, rho::T,
     dmaster_stochastic(dx, rho, C, Cdagger, drho, n)
 end
 
-function integrate_master_stoch(tspan, df::Function, dg::Function,
-                        rho0::Operator, fout::Union{Nothing, Function},
-                        n::Int;
+function integrate_master_stoch(tspan, df, dg,
+                        rho0, fout,
+                        n;
                         kwargs...)
     tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
-    x0 = reshape(rho0.data, length(rho0))
+    x0 = as_vector(rho0)
     state = copy(rho0)
     dstate = copy(rho0)
     integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...)
 end
 
-function check_master_stoch(rho0::Operator{B,B}, C::Vector, Cdagger::Vector) where B<:Basis
+function check_master_stoch(rho0::Operator{B,B}, C, Cdagger) where B
     # TODO: replace type checks by dispatch; make types of C known
     @assert length(C) == length(Cdagger)
     isreducible = true
@@ -184,10 +177,11 @@ function check_master_stoch(rho0::Operator{B,B}, C::Vector, Cdagger::Vector) whe
     isreducible
 end
 
+as_vector(rho::Operator) = reshape(rho.data, length(rho))
 
 # TODO: Speed up by recasting to n-d arrays, remove vector methods
-function recast!(x::Union{Vector, SubArray}, rho::Operator{B,B,T}) where {B<:Basis,T}
+function recast!(x::Union{Vector, SubArray}, rho::Operator{B,B,T}) where {B,T}
     rho.data = reshape(x, size(rho.data))
 end
-recast!(state::Operator{B,B}, x::SubArray) where B<:Basis = (x[:] = state.data)
-recast!(state::Operator{B,B}, x::Vector) where B<:Basis = nothing
+recast!(state::Operator{B,B}, x::SubArray) where B = (x[:] = state.data)
+recast!(state::Operator{B,B}, x::Vector) where B = nothing
diff --git a/src/stochastic_schroedinger.jl b/src/stochastic_schroedinger.jl
index 636c738e..82d2c3e8 100644
--- a/src/stochastic_schroedinger.jl
+++ b/src/stochastic_schroedinger.jl
@@ -20,10 +20,10 @@ Integrate stochastic Schrödinger equation.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
 function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}, Hs::Vector;
-                fout::Union{Function,Nothing}=nothing,
-                normalize_state::Bool=false,
+                fout=nothing,
+                normalize_state=false,
                 calback=nothing,
-                kwargs...) where {B<:Basis,T<:Ket{B}}
+                kwargs...) where {B,T<:Ket{B}}
     tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
 
     n = length(Hs)
@@ -36,12 +36,11 @@ function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}, Hs::Vector;
         @assert isa(h, AbstractOperator{B,B})
     end
 
-    dschroedinger_determ(t, psi::T, dpsi::T) = dschroedinger(psi, H, dpsi)
-    dschroedinger_stoch(dx::DiffArray,
-            t, psi::T, dpsi::T, n::Int) = dschroedinger_stochastic(dx, psi, Hs, dpsi, n)
+    dschroedinger_determ(t, psi, dpsi) = dschroedinger(psi, H, dpsi)
+    dschroedinger_stoch(dx, t, psi, dpsi, n) = dschroedinger_stochastic(dx, psi, Hs, dpsi, n)
 
     if normalize_state
-        norm_func(u::Vector, t, integrator) = normalize!(u)
+        norm_func(u, t, integrator) = normalize!(u)
         ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
                  func_everystep=true,
                  func_start=false)
@@ -53,7 +52,7 @@ function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}, Hs::Vector;
                     ncb=ncb,
                     kwargs...)
 end
-schroedinger(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, Hs::AbstractOperator{B,B}; kwargs...) where B<:Basis = schroedinger(tspan, psi0, H, [Hs]; kwargs...)
+schroedinger(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, Hs::AbstractOperator{B,B}; kwargs...) where B = schroedinger(tspan, psi0, H, [Hs]; kwargs...)
 
 """
     stochastic.schroedinger_dynamic(tspan, state0, fdeterm, fstoch[; fout, ...])
@@ -81,10 +80,10 @@ Integrate stochastic Schrödinger equation with dynamic Hamiltonian.
         each time step taken by the solver.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function schroedinger_dynamic(tspan, psi0::T, fdeterm::Function, fstoch::Function;
-                fout::Union{Function,Nothing}=nothing, noise_processes::Int=0,
-                normalize_state::Bool=false,
-                kwargs...) where T<:Ket
+function schroedinger_dynamic(tspan, psi0::Ket, fdeterm, fstoch;
+                fout=nothing, noise_processes::Int=0,
+                normalize_state=false,
+                kwargs...)
     tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
 
     if noise_processes == 0
@@ -98,13 +97,11 @@ function schroedinger_dynamic(tspan, psi0::T, fdeterm::Function, fstoch::Functio
     x0 = psi0.data
     state = copy(psi0)
 
-    dschroedinger_determ(t, psi::T, dpsi::T) = dschroedinger_dynamic(t, psi, fdeterm, dpsi)
-    dschroedinger_stoch(dx::DiffArray,
-            t, psi::T, dpsi::T, n::Int) =
-        dschroedinger_stochastic(dx, t, psi, fstoch, dpsi, n)
+    dschroedinger_determ(t, psi, dpsi) = dschroedinger_dynamic(t, psi, fdeterm, dpsi)
+    dschroedinger_stoch(dx, t, psi, dpsi, n) = dschroedinger_stochastic(dx, t, psi, fstoch, dpsi, n)
 
     if normalize_state
-        norm_func(u::Vector, t, integrator) = normalize!(u)
+        norm_func(u, t, integrator) = normalize!(u)
         ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
                  func_everystep=true,
                  func_start=false)
@@ -119,13 +116,11 @@ function schroedinger_dynamic(tspan, psi0::T, fdeterm::Function, fstoch::Functio
 end
 
 
-function dschroedinger_stochastic(dx::D, psi::T1, Hs::Vector{T2},
-            dpsi::T1, index::Int) where {D<:Vector,B<:Basis,T1<:Ket{B},T2<:AbstractOperator{B,B}}
+function dschroedinger_stochastic(dx::AbstractVector, psi, Hs, dpsi, index)
     recast!(dx, dpsi)
     dschroedinger(psi, Hs[index], dpsi)
 end
-function dschroedinger_stochastic(dx::Matrix, psi::T1, Hs::Vector{T2},
-            dpsi::T1, n::Int) where {B<:Basis,T1<:Ket{B},T2<:AbstractOperator{B,B}}
+function dschroedinger_stochastic(dx::AbstractMatrix, psi, Hs, dpsi, n)
     for i=1:n
         dx_i = @view dx[:, i]
         recast!(dx_i, dpsi)
@@ -133,8 +128,7 @@ function dschroedinger_stochastic(dx::Matrix, psi::T1, Hs::Vector{T2},
         recast!(dpsi, dx_i)
     end
 end
-function dschroedinger_stochastic(dx::DiffArray,
-            t, psi::T, f::Function, dpsi::T, n::Int) where T<:Ket
+function dschroedinger_stochastic(dx, t, psi, f, dpsi, n)
     ops = f(t, psi)
     if QO_CHECKS[]
         @inbounds for h=ops
diff --git a/src/stochastic_semiclassical.jl b/src/stochastic_semiclassical.jl
index cbffae34..90cae791 100644
--- a/src/stochastic_semiclassical.jl
+++ b/src/stochastic_semiclassical.jl
@@ -37,16 +37,16 @@ Integrate time-dependent Schrödinger equation coupled to a classical system.
         each time step taken by the solver.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function schroedinger_semiclassical(tspan, state0::S, fquantum::Function,
-                fclassical::Function; fstoch_quantum::Union{Nothing, Function}=nothing,
-                fstoch_classical::Union{Nothing, Function}=nothing,
-                fout::Union{Function,Nothing}=nothing,
+function schroedinger_semiclassical(tspan, state0::S, fquantum,
+                fclassical; fstoch_quantum=nothing,
+                fstoch_classical=nothing,
+                fout=nothing,
                 noise_processes::Int=0,
                 noise_prototype_classical=nothing,
                 normalize_state::Bool=false,
                 kwargs...) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
     tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
-    dschroedinger_det(t, state::S, dstate::S) =
+    dschroedinger_det(t, state, dstate) =
             semiclassical.dschroedinger_dynamic(t, state, fquantum, fclassical, dstate)
 
     if isa(fstoch_quantum, Nothing) && isa(fstoch_classical, Nothing)
@@ -76,7 +76,7 @@ function schroedinger_semiclassical(tspan, state0::S, fquantum::Function,
 
     if normalize_state
         len_q = length(state0.quantum)
-        function norm_func(u::Vector, t, integrator)
+        function norm_func(u, t, integrator)
             u .= [normalize!(u[1:len_q]); u[len_q+1:end]]
         end
         ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
@@ -86,7 +86,7 @@ function schroedinger_semiclassical(tspan, state0::S, fquantum::Function,
         ncb = nothing
     end
 
-    dschroedinger_stoch(dx::DiffArray, t, state::S, dstate::S, n::Int) =
+    dschroedinger_stoch(dx, t, state, dstate, n) =
             dschroedinger_stochastic(dx, t, state, fstoch_quantum, fstoch_classical, dstate, n)
 
     integrate_stoch(tspan_, dschroedinger_det, dschroedinger_stoch, x0, state, dstate, fout, n;
@@ -136,11 +136,11 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
 function master_semiclassical(tspan, rho0::S,
-                fquantum::Function, fclassical::Function;
-                fstoch_quantum::Union{Function, Nothing}=nothing,
-                fstoch_classical::Union{Function, Nothing}=nothing,
-                rates::DecayRates=nothing,
-                fout::Union{Function,Nothing}=nothing,
+                fquantum, fclassical;
+                fstoch_quantum=nothing,
+                fstoch_classical=nothing,
+                rates=nothing,
+                fout=nothing,
                 noise_processes::Int=0,
                 noise_prototype_classical=nothing,
                 nonlinear::Bool=true,
@@ -167,11 +167,10 @@ function master_semiclassical(tspan, rho0::S,
         end
     end
 
-    dmaster_determ(t, rho::S, drho::S) =
+    dmaster_determ(t, rho, drho) =
             semiclassical.dmaster_h_dynamic(t, rho, fquantum, fclassical, rates, drho, tmp)
 
-    dmaster_stoch(dx::DiffArray, t, rho::S,
-                    drho::S, n::Int) =
+    dmaster_stoch(dx, t, rho, drho, n) =
         dmaster_stoch_dynamic(dx, t, rho, fstoch_quantum, fstoch_classical, drho, n)
 
     integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n;
@@ -182,18 +181,18 @@ master_semiclassical(tspan, psi0::State{B,T}, args...; kwargs...) where {B<:Basi
         master_semiclassical(tspan, dm(psi0), args...; kwargs...)
 
 # Derivative functions
-function dschroedinger_stochastic(dx::Vector, t,
-        state::S, fstoch_quantum::Function, fstoch_classical::Nothing,
-        dstate::S, ::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
+function dschroedinger_stochastic(dx::AbstractVector, t,
+        state, fstoch_quantum::Function, fstoch_classical::Nothing,
+        dstate, n)
     H = fstoch_quantum(t, state.quantum, state.classical)
     recast!(dx, dstate)
     QO_CHECKS[] && check_schroedinger(state.quantum, H[1])
     dschroedinger(state.quantum, H[1], dstate.quantum)
     recast!(dstate, dx)
 end
-function dschroedinger_stochastic(dx::Matrix,
-        t, state::S, fstoch_quantum::Function,
-        fstoch_classical::Nothing, dstate::S, n::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
+function dschroedinger_stochastic(dx::AbstractMatrix, t,
+        state, fstoch_quantum::Function, fstoch_classical::Nothing,
+        dstate, n)
     H = fstoch_quantum(t, state.quantum, state.classical)
     for i=1:n
         dx_i = @view dx[:, i]
@@ -203,23 +202,23 @@ function dschroedinger_stochastic(dx::Matrix,
         recast!(dstate, dx_i)
     end
 end
-function dschroedinger_stochastic(dx::DiffArray, t,
-            state::S, fstoch_quantum::Nothing, fstoch_classical::Function,
-            dstate::S, ::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
+function dschroedinger_stochastic(dx, t,
+            state, fstoch_quantum::Nothing, fstoch_classical::Function,
+            dstate, n)
     dclassical = @view dx[length(state.quantum)+1:end, :]
     fstoch_classical(t, state.quantum, state.classical, dclassical)
 end
-function dschroedinger_stochastic(dx::Matrix, t, state::S, fstoch_quantum::Function,
-            fstoch_classical::Function, dstate::S, n::Int) where {B<:Basis,T<:Ket{B},S<:State{B,T}}
+function dschroedinger_stochastic(dx, t, state, fstoch_quantum::Function,
+            fstoch_classical::Function, dstate, n)
     dschroedinger_stochastic(dx, t, state, fstoch_quantum, nothing, dstate, n)
 
     dx_i = @view dx[length(state.quantum)+1:end, n+1:end]
     fstoch_classical(t, state.quantum, state.classical, dx_i)
 end
 
-function dmaster_stoch_dynamic(dx::Vector, t,
-            state::S, fstoch_quantum::Function,
-            fstoch_classical::Nothing, dstate::S, ::Int) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
+function dmaster_stoch_dynamic(dx::AbstractVector, t,
+            state, fstoch_quantum::Function,
+            fstoch_classical::Nothing, dstate, n)
     result = fstoch_quantum(t, state.quantum, state.classical)
     QO_CHECKS[] && @assert length(result) == 2
     C, Cdagger = result
@@ -230,9 +229,9 @@ function dmaster_stoch_dynamic(dx::Vector, t,
     dstate.quantum.data .-= tr(dstate.quantum)*state.quantum.data
     recast!(dstate, dx)
 end
-function dmaster_stoch_dynamic(dx::Matrix, t,
-            state::S, fstoch_quantum::Function,
-            fstoch_classical::Nothing, dstate::S, n::Int) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
+function dmaster_stoch_dynamic(dx::AbstractMatrix, t,
+            state, fstoch_quantum::Function,
+            fstoch_classical::Nothing, dstate, n)
     result = fstoch_quantum(t, state.quantum, state.classical)
     QO_CHECKS[] && @assert length(result) == 2
     C, Cdagger = result
@@ -246,33 +245,21 @@ function dmaster_stoch_dynamic(dx::Matrix, t,
         recast!(dstate, dx_i)
     end
 end
-function dmaster_stoch_dynamic(dx::DiffArray, t,
-            state::S, fstoch_quantum::Nothing,
-            fstoch_classical::Function, dstate::S, ::Int) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
+function dmaster_stoch_dynamic(dx, t,
+            state, fstoch_quantum::Nothing,
+            fstoch_classical::Function, dstate, n)
     dclassical = @view dx[length(state.quantum)+1:end, :]
     fstoch_classical(t, state.quantum, state.classical, dclassical)
 end
-function dmaster_stoch_dynamic(dx::Matrix, t,
-            state::S, fstoch_quantum::Function,
-            fstoch_classical::Function, dstate::S, n::Int) where {B<:Basis,T<:Operator{B,B},S<:State{B,T}}
+function dmaster_stoch_dynamic(dx, t,
+            state, fstoch_quantum::Function,
+            fstoch_classical::Function, dstate, n)
     dmaster_stoch_dynamic(dx, t, state, fstoch_quantum, nothing, dstate, n)
 
     dx_i = @view dx[length(state.quantum)+1:end, n+1:end]
     fstoch_classical(t, state.quantum, state.classical, dx_i)
 end
 
-function integrate_master_stoch(tspan, df::Function, dg::Function,
-                        rho0::State{B,T}, fout::Union{Nothing, Function},
-                        n::Int;
-                        kwargs...) where {B<:Basis,T<:Operator{B,B}}
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
-    x0 = Vector{eltype(rho0)}(undef, length(rho0))
-    recast!(rho0, x0)
-    state = copy(rho0)
-    dstate = copy(rho0)
-    integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...)
-end
-
 function recast!(state::State, x::SubArray)
     N = length(state.quantum)
     copyto!(x, 1, state.quantum.data, 1, N)
diff --git a/src/timecorrelations.jl b/src/timecorrelations.jl
index 8c4d6ad5..80567240 100644
--- a/src/timecorrelations.jl
+++ b/src/timecorrelations.jl
@@ -33,26 +33,25 @@ criterion specified in [`steadystate.master`](@ref).
 * `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function correlation(tspan::Vector, rho0::Operator{B,B}, H::AbstractOperator{B,B}, J::Vector,
-                     op1::AbstractOperator{B,B}, op2::AbstractOperator{B,B};
-                     rates::Union{Vector, Matrix, Nothing}=nothing,
-                     Jdagger::Vector=dagger.(J),
-                     kwargs...) where B<:Basis
-    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
+function correlation(tspan, rho0::Operator, H::AbstractOperator, J,
+                     op1, op2;
+                     rates=nothing,
+                     Jdagger=dagger.(J),
+                     kwargs...)
     function fout(t, rho)
         expect(op1, rho)
     end
-    t,u = timeevolution.master(tspan_, op2*rho0, H, J; rates=rates, Jdagger=Jdagger,
+    t,u = timeevolution.master(tspan, op2*rho0, H, J; rates=rates, Jdagger=Jdagger,
                         fout=fout, kwargs...)
     u
 end
 
-function correlation(rho0::Operator{B,B}, H::AbstractOperator{B,B}, J::Vector,
-                     op1::AbstractOperator{B,B}, op2::AbstractOperator{B,B};
+function correlation(rho0::Operator, H::AbstractOperator, J,
+                     op1, op2;
                      tol=1e-4,
-                     rates::Union{Vector, Matrix, Nothing}=nothing,
-                     Jdagger::Vector=dagger.(J),
-                     kwargs...) where B<:Basis
+                     rates=nothing,
+                     Jdagger=dagger.(J),
+                     kwargs...)
     op2rho0 = op2*rho0
     exp1 = expect(op1, op2rho0)
     function fout(t, rho)
@@ -97,11 +96,11 @@ automatically.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
 function spectrum(omega_samplepoints,
-                H::AbstractOperator{B,B}, J::Vector, op::AbstractOperator{B,B};
-                rho0::Operator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
+                H::AbstractOperator, J, op;
+                rho0=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
                 tol=1e-4,
-                rho_ss::Operator{B,B}=steadystate.master(H, J; tol=tol, rho0=rho0)[end][end],
-                kwargs...) where B<:Basis
+                rho_ss=steadystate.master(H, J; tol=tol, rho0=rho0)[end][end],
+                kwargs...)
     domega = minimum(diff(omega_samplepoints))
     dt = 2*pi/abs(omega_samplepoints[end] - omega_samplepoints[1])
     T = 2*pi/domega
@@ -111,11 +110,11 @@ function spectrum(omega_samplepoints,
     return omega_samplepoints, S
 end
 
-function spectrum(H::AbstractOperator{B,B}, J::Vector, op::AbstractOperator{B,B};
-                rho0::Operator{B,B}=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
+function spectrum(H::AbstractOperator, J, op;
+                rho0=tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1))),
                 tol=1e-4,
-                rho_ss::Operator{B,B}=steadystate.master(H, J; tol=tol)[end][end],
-                kwargs...) where B<:Basis
+                rho_ss=steadystate.master(H, J; tol=tol)[end][end],
+                kwargs...)
     tspan, exp_values = correlation(rho_ss, H, J, dagger(op), op, tol=tol, kwargs...)
     dtmin = minimum(diff(tspan))
     T = tspan[end] - tspan[1]
@@ -137,7 +136,7 @@ Calculate spectrum as Fourier transform of a correlation function with a given c
 * `corr`: Correlation function of which the Fourier transform is to be calculated.
 * `normalize_spec`: Specify if spectrum should be normalized to its maximum.
 """
-function correlation2spectrum(tspan, corr::Vector{T}; normalize_spec::Bool=false) where T <: Number
+function correlation2spectrum(tspan, corr; normalize_spec=false)
   n = length(tspan)
   if length(corr) != n
     ArgumentError("tspan and corr must be of same length!")
diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl
index 5885901d..02a1ad48 100644
--- a/src/timeevolution_base.jl
+++ b/src/timeevolution_base.jl
@@ -11,13 +11,13 @@ function recast! end
 
 Integrate using OrdinaryDiffEq
 """
-function integrate(tspan, df::Function, x0::X,
-            state::T, dstate::T, fout::Function;
-            alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(),
+function integrate(tspan, df, x0,
+            state, dstate, fout;
+            alg = OrdinaryDiffEq.DP5(),
             steady_state = false, tol = 1e-3, save_everystep = false, saveat=tspan,
             callback = nothing, kwargs...) where {T,X}
 
-    function df_(dx::T, x::T, p, t) where T
+    function df_(dx, x, p, t)
         recast!(x, state)
         recast!(dx, dstate)
         df(t, state, dstate)
@@ -28,15 +28,16 @@ function integrate(tspan, df::Function, x0::X,
         fout(t, state)
     end
 
-    out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)})
+    tType = float(eltype(tspan))
+    out_type = pure_inference(fout, Tuple{tType,typeof(state)})
 
-    out = DiffEqCallbacks.SavedValues(eltype(tspan),out_type)
+    out = DiffEqCallbacks.SavedValues(tType,out_type)
 
     scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=saveat,
                                          save_everystep=save_everystep,
                                          save_start = false)
 
-    prob = OrdinaryDiffEq.ODEProblem{true}(df_, x0,(tspan[1],tspan[end]))
+    prob = OrdinaryDiffEq.ODEProblem{true}(df_, x0,(convert(tType, tspan[1]),convert(tType, tspan[end])))
 
     if steady_state
         affect! = function (integrator)
@@ -65,9 +66,9 @@ function integrate(tspan, df::Function, x0::X,
     out.t,out.saveval
 end
 
-function integrate(tspan, df::Function, x0::X,
-            state::T, dstate::T, ::Nothing; kwargs...) where {T,X}
-    function fout(t, state::T)
+function integrate(tspan, df, x0,
+            state, dstate, ::Nothing; kwargs...)
+    function fout(t, state)
         copy(state)
     end
     integrate(tspan, df, x0, state, dstate, fout; kwargs...)