Efficient evaluation of BSplineField:s for a vector of times

src/arithmetic.jl

@@ -113,6 +113,8 @@ end
 for fun in [:vector_potential, :field_amplitude, :vector_potential_spectrum]
     @eval $fun(f::SumField, t::Number) =
         $fun(f.a, t) + $fun(f.b, t)
+    @eval $fun(f::SumField, t::AbstractVector) =
+        $fun(f.a, t) + $fun(f.b, t)
 end
 
 polarization(f::SumField) = polarization(f.a)
@@ -226,6 +228,7 @@ Base.parent(f::NegatedField) = f.a
 
 for fun in [:vector_potential, :vector_potential_spectrum]
     @eval $fun(f::NegatedField, t::Number) = -$fun(parent(f), t)
+    @eval $fun(f::NegatedField, t::AbstractVector) = -$fun(parent(f), t)
 end
 
 rotate(f::NegatedField, R) = NegatedField(rotate(f.a, R))
@@ -298,6 +301,8 @@ end
 for fun in [:vector_potential, :field_amplitude, :intensity]
     @eval $fun(f::DelayedField, t::Number) =
         $fun(f.a, t-f.t₀)
+    @eval $fun(f::DelayedField, t::AbstractVector) =
+        $fun(f.a, t .- f.t₀)
 end
 
 vector_potential_spectrum(f::DelayedField, ω) =

src/bspline_field.jl

@@ -39,4 +39,13 @@ linear_combination(V::AbstractVector, C::AbstractMatrix) =
 vector_potential(f::BSplineField, t::Number) =
     linear_combination(f.B[t,:], f.C)
 
+vector_potential(f::BSplineField, t::AbstractVector) =
+    f.B[t,:]*f.C
+
+field_amplitude(f::BSplineField, t::Number) =
+    -linear_combination(derivative(f.B, t, :, 1), f.C)
+
+field_amplitude(f::BSplineField, t::AbstractVector) =
+    -derivative(f.B, t, :, 1)*f.C
+
 export BSplineField

src/bsplines.jl

@@ -146,6 +146,46 @@ end
 getindex(B::BSpline, x, ::Colon) =
     getindex(B, x, 1:size(B,2))
 
+function derivative(B::BSpline, x::Real, j::Integer, m)
+    T = promote_type(eltype(B.t), typeof(x))
+    deBoor(B.t, UnitVector{T}(size(B,2), j),
+           x, find_interval(B.t, x), m)
+end
+
+function derivative(B::BSpline, x::Real, sel::AbstractVector, m)
+    T = promote_type(eltype(B.t), typeof(x))
+    i = find_interval(B.t, x)
+    χ = spzeros(T, length(sel))
+    o = sel[1] - 1
+    for j in sel
+        χ[j-o] = deBoor(B.t, UnitVector{T}(size(B,2), j), x, i, m)
+    end
+    χ
+end
+
+function derivative_basis_function!(χ, B::BSpline, x::AbstractVector, j, m)
+    T = promote_type(eltype(B.t), eltype(x))
+    eⱼ = UnitVector{T}(size(B,2), j)
+    for (is,k) ∈ within_support(x, B.t, j)
+        for i in is
+            χ[i] = deBoor(B.t, eⱼ, x[i], k, m)
+        end
+    end
+end
+
+function derivative(B::BSpline, x::AbstractVector, sel::AbstractVector, m)
+    T = promote_type(eltype(B.t), eltype(x))
+    χ = spzeros(T, length(x), length(sel))
+    o = sel[1] - 1
+    for j in sel
+        derivative_basis_function!(view(χ, :, j-o), B, x, j, m)
+    end
+    χ
+end
+
+derivative(B::BSpline, x, ::Colon, m) =
+    derivative(B, x, 1:size(B,2), m)
+
 struct BSplineView{Bt,Sel}
     B::Bt
     sel::Sel
@@ -174,6 +214,9 @@ end
 getindex(B::BSplineView, x, j) =
     getindex(B.B, x, B.sel[j])
 
+derivative(B::BSplineView, x, j, m) =
+    derivative(B.B, x, B.sel[j], m)
+
 const BSplineOrView = Union{BSpline,BSplineView}
 
 # * Function interpolation

src/dispersed_fields.jl

@@ -98,6 +98,8 @@ Base.parent(f::DispersedField) = f.f
 for fun in [:vector_potential, :field_amplitude]
     @eval $fun(f::DispersedField, t::Number) =
         $fun(f.FB, t)
+    @eval $fun(f::DispersedField, t::AbstractVector) =
+        $fun(f.FB, t)
 end
 
 # Some of these forwards are iffy, i.e. after e.g. a chirp, the

test/bspline_fields.jl

@@ -44,6 +44,17 @@
 
             test_approx_eq(Av, Arec, rtol=1e-4)
             test_approx_eq(Fv, Frec, rtol=1e-4)
+
+            @testset "Scalar evaluation" begin
+                sel = 1:1000
+                correct_tensor(::LinearPolarization, v) = v
+                correct_tensor(::ArbitraryPolarization, v) = transpose(reduce(hcat, v))
+                Arec2 = correct_tensor(polarization(FB), vector_potential.(Ref(FB), t[sel]))
+                Frec2 = correct_tensor(polarization(FB), field_amplitude.(Ref(FB), t[sel]))
+
+                test_approx_eq(selectdim(Arec, 1, sel), Arec2, rtol=1e-14)
+                test_approx_eq(selectdim(Frec, 1, sel), Frec2, rtol=1e-14)
+            end
         end
     end
 end