diff --git a/src/bases/timebasis.jl b/src/bases/timebasis.jl
index 5b087ece..f1350faf 100644
--- a/src/bases/timebasis.jl
+++ b/src/bases/timebasis.jl
@@ -204,13 +204,13 @@ function assemblydata(tbf::TimeBasisFunction)
     els = [ simplex(point((i-1)*Δt),point(i*Δt)) for i in 1:num_cells ]
 
     for k in 1 : numfunctions(tbf)
-        tk = (k-1) * Δt
+        tk = k * Δt
         for i in 1 : numintervals(tbf)
             # Focus on interval [(i-2)Δt,(i-1)Δt]
             p = tbf.polys[i]
             q = substitute(p,t-tk)
 
-            c = k + i - 2
+            c = k + i - 1
             1 <= c <= num_cells || continue
             for d = 0 : degree(q)
                 r = d + 1
@@ -301,10 +301,10 @@ function assemblydata(tbf::TimeBasisDelta)
     num_funcs = zeros(Int, num_cells, num_refs)
     data = fill((0,z), max_num_funcs, num_refs, num_cells)
 
-    els = [ simplex(point((i-0)*Δt),point((i+1)*Δt)) for i in 1:num_cells ]
+    els = [ simplex(point((i-1)*Δt),point(i*Δt)) for i in 1:num_cells ]
 
-    for k in 1 : numfunctions(tbf)-1
-        data[1,1,k] = (k+1,w)
+    for k in 1 : numfunctions(tbf)
+        data[1,1,k] = (k,w)
     end
 
     return els, AssemblyData(data)
diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl
index fc869ca6..e6d284c4 100644
--- a/src/maxwell/mwops.jl
+++ b/src/maxwell/mwops.jl
@@ -14,6 +14,16 @@ struct MWSingleLayer3D{T,U} <: MaxwellOperator3D{T,U}
   β::U
 end
 
+struct MWStaticExtractedHyperSingular{T,U} <: MaxwellOperator3D{T,U}
+  gamma::T
+  β::U
+end
+
+# Maxwell hypersingular operator, with its static kernel removed, can be assembled using basis functions other than div-conforming ones
+MWStaticExtractedHyperSingular(gamma) = MWStaticExtractedHyperSingular(gamma, -1/(gamma))
+
+defaultquadstrat(op::MWStaticExtractedHyperSingular, tfs::RTRefSpace, bfs::RTRefSpace) = DoubleNumSauterQstrat(6,7,5,5,4,3)
+
 gamma(op::MWSingleLayer3D{Val{0}, U}) where {U} = zero(U)
 
 scalartype(op::MWSingleLayer3D{T,U}) where {T,U} = promote_type(T,U)
@@ -137,6 +147,44 @@ function (igd::Integrand{<:MWSingleLayer3DReg})(x,y,f,g)
     end
 end
 
+# MWStaticExtractedHyperSingular operator, discretized and tested by RT. Integration-by-parts is performed
+function (igd::Integrand{<:MWStaticExtractedHyperSingular, <:RTRefSpace, <:RTRefSpace})(x,y,f,g)
+    β = igd.operator.β
+    γ = igd.operator.gamma
+
+    r = cartesian(x) - cartesian(y)
+    R = norm(r)
+    iR = 1/R
+    γR = γ*R
+
+    green = expm1(-γR) * (i4pi * iR) + γ * i4pi
+    βG = β * green
+
+    _integrands(f,g) do fi,gj
+        βG * dot(fi.divergence, gj.divergence)
+    end
+end
+
+# MWStaticExtractedHyperSingular operator, discretized by RT and tested by other functions. Integration-by-parts is not performed. Instead, grad of regularized Green's function is used
+function (igd::Integrand{<:MWStaticExtractedHyperSingular, <:RefSpace, <:RTRefSpace})(x,y,f,g)
+    β = igd.operator.β
+    γ = igd.operator.gamma
+
+    r = cartesian(x) - cartesian(y)
+    R = norm(r)
+    iR = 1/R
+    γR = γ*R
+
+    gradgreen = -(expm1(-γR) * (1 + γR) + γR) * (i4pi * iR^3) * r
+
+    # Minus sign is because integration-by-part is not performed
+    βgG = -β * gradgreen
+
+    _integrands(f,g) do fi,gj
+        dot(fi.value, βgG*gj.divergence)
+    end
+end
+
 
 function (igd::Integrand{<:MWDoubleLayer3D})(x,y,f,g)
 
diff --git a/src/utils/specialfns.jl b/src/utils/specialfns.jl
index 18cda970..fc4d40b8 100644
--- a/src/utils/specialfns.jl
+++ b/src/utils/specialfns.jl
@@ -28,16 +28,23 @@ end
 function fouriertransform(a::Array, dt, t0, dim=1)
     n = size(a,dim)
     dω = 2π / (n*dt)
-    b = fftshift(fft(a, dim), dim) * dt / sqrt(2π)
     ω0 = -dω * div(n,2)
+    ω = ω0 .+ dω .* (0:n-1)
+
+    b = fftshift(fft(a, dim), dim) * dt / sqrt(2π)
+    b .*= reshape(exp.(-im .* ω .* t0), ntuple(j -> j == dim ? n : 1, ndims(a))) 
+
     b, dω, ω0
 end
 
 function inversefouriertransform(a::Array, dω, ω0, dim=1)
     n = size(a,dim)
     dt = 2π/ (n*dω)
-    b = ifft(a,dim) * sqrt(2π) / dt
+    ω = ω0 .+ dω .* (0:n-1)
     t0 = -dt * div(n,2)
+
+    b = ifft(ifftshift(a .* reshape(exp.(im .* ω .* t0), ntuple(j -> j == dim ? n : 1, ndims(a))), dim), dim) * sqrt(2π) / dt
+
     b, dt, t0
 end
 
diff --git a/test/test_assemble_tdexcitation.jl b/test/test_assemble_tdexcitation.jl
new file mode 100644
index 00000000..80e657a6
--- /dev/null
+++ b/test/test_assemble_tdexcitation.jl
@@ -0,0 +1,83 @@
+@testitem "testing TDFunctional" begin
+
+    using CompScienceMeshes
+    
+    Γ = meshcuboid(1.0, 1.0, 1.0, 2.0; generator=:gmsh)
+    X = raviartthomas(Γ)
+
+    struct ConstFunctional{T} <: BEAST.Functional{T}
+        constant::T
+    end
+
+    function(f::ConstFunctional)(r)
+        f.constant * [1.0, 0.0, 0.0]
+    end
+
+    struct FuncXGaussian{T} <: BEAST.TDFunctional{T}
+        functional::ConstFunctional{T}
+        gaussian::BEAST.Gaussian{T}
+    end
+
+    function(f::FuncXGaussian)(r,t)
+        r = cartesian(r)
+        t = cartesian(t)[1]
+
+        f.functional(r) * f.gaussian(t)     
+    end
+
+    Δt, Nt = 0.5, 20
+    duration = 4 * Δt                                
+    delay = 10 * Δt                                        
+    
+    func = ConstFunctional(1.0)
+    gaussian = creategaussian(duration, delay)
+    constxgaussian = FuncXGaussian(func, gaussian)
+
+
+    ### Testing with Dirac delta 
+    δ = timebasisdelta(Δt, Nt)
+    A = assemble(constxgaussian, X⊗δ)
+    # Analytically temporal testing
+    Ats = assemble(func, X) * gaussian.(Δt*[1:1:Nt;])'
+
+    @test maximum(abs.(A - Ats)) < 1e-12
+
+
+    ### Testing with pulse functions
+    p = timebasiscxd0(Δt, Nt) 	                			                      
+    B = assemble(constxgaussian, X⊗p)
+    igaussian = integrate(gaussian)
+    Bts = assemble(func, X) * (igaussian.(Δt*[1:1:Nt;]) - igaussian.(Δt*[0:1:Nt-1;]))'
+    
+    @test maximum(abs.(B - Bts)) < 1e-12
+
+    
+    ### Testing with hat functions
+    import SpecialFunctions: erf 
+    function integratewithh(g::BEAST.Gaussian, ts, Δt)
+        A = g.scaling
+        t0 = g.delay
+        w = g.width
+        
+        r = zeros(length(ts))
+
+        for i in ts
+            a = (i-1)*Δt
+            x = i*Δt
+            b = (i+1)*Δt
+
+            ua = 4*(a-t0)/w
+            ux = 4*(x-t0)/w
+            ub = 4*(b-t0)/w
+
+            r[i] = A/(2*Δt) * ((t0 - a)*(erf(ux) - erf(ua)) + (b - t0)*(erf(ub) - erf(ux))) + A*w/(8*√(π)*Δt)*(exp(-ua^2) + exp(-ub^2) - 2*exp(-ux^2))
+        end
+        return r
+    end
+
+    h = timebasisc0d1(Δt, Nt) 	 
+    C = assemble(constxgaussian, X⊗h)
+    Cts = assemble(func, X) * integratewithh(gaussian, [1:1:Nt;], Δt)'
+
+    @test maximum(abs.(C - Cts)) < 1e-12
+end
diff --git a/test/test_tdrhs_scaling.jl b/test/test_tdrhs_scaling.jl
index 5f7ea258..211c1e69 100644
--- a/test/test_tdrhs_scaling.jl
+++ b/test/test_tdrhs_scaling.jl
@@ -44,12 +44,12 @@ derive(gaussian2).(taxis2)[120]/sol2
 @test all(derive(gaussian1).(taxis1)/sol1 .≈ derive(gaussian2).(taxis2)/sol2)
 
 timeels1, timead1 = BEAST.assemblydata(U1)
-@test timeels1[100][1][1] ≈ Δt1*100
-@test timeels1[100][2][1] ≈ Δt1*101
+@test timeels1[100][1][1] ≈ Δt1*99
+@test timeels1[100][2][1] ≈ Δt1*100
 
 timeels2, timead2 = BEAST.assemblydata(U2)
-@test timeels2[100][1][1] ≈ Δt2*100
-@test timeels2[100][2][1] ≈ Δt2*101
+@test timeels2[100][1][1] ≈ Δt2*99
+@test timeels2[100][2][1] ≈ Δt2*100
 
 b1[120]/sol1
 b2[120]/sol2