Compute strong-field properties for cubic Hermite spline fields

src/arithmetic.jl

@@ -287,6 +287,13 @@ for fun in [:vector_potential, :vector_potential_spectrum]
     @eval $fun(f::ScaledField, t::AbstractVector) = f.η*$fun(parent(f), t)
 end
 
+for fun in [:amplitude, :free_oscillation_amplitude]
+    @eval $fun(f::ScaledField) = f.η*$fun(parent(f))
+end
+for fun in [:intensity, :ponderomotive_potential]
+    @eval $fun(f::ScaledField, args...) = f.η^2*$fun(parent(f), args...)
+end
+
 rotate(f::ScaledField, R) = ScaledField(rotate(parent(f), R), f.η)
 
 function Base.show(io::IO, f::ScaledField)
@@ -296,7 +303,10 @@ end
 
 function Base.show(io::IO, mime::MIME"text/plain", f::ScaledField)
     printfmt(io, "ScaledField: {1:s} × ", f.η)
-    show(io, mime, parent(f))
+    show(io, parent(f))
+    println(io)
+    print(io, "  ")
+    show_strong_field_properties(io, f)
 end
 
 # ** Delayed fields

src/cubic_hermite_spline_field.jl

@@ -62,6 +62,13 @@ Base.show(io::IO, f::CubicHermiteSplineField) =
     printfmt(io, "{1:d}-sample, {2:d}-component Cubic Hermite spline field, t ∈ {3:s}",
              length(f.t), dimensions(f), span(f))
 
+function Base.show(io::IO, ::MIME"text/plain", f::CubicHermiteSplineField)
+    show(io, f)
+    println(io)
+    print(io, "  ")
+    show_strong_field_properties(io, f)
+end
+
 span(f::CubicHermiteSplineField) = first(f.t)..last(f.t)
 
 min_step(t::AbstractRange) = step(t)
@@ -79,6 +86,77 @@ vector_potential(f::CubicHermiteSplineField, t::Number) =
 
 export CubicHermiteSplineField
 
+# * Strong field properties
+
+function _maximize_interpolation(t, v, f, f′; verbosity=0)
+    # This is only a rough approximation, and will most likely only
+    # work reasonably well if there is one isolated pulse in the
+    # trace.
+
+    N = length(t)
+
+    # We begin by finding the maximum point.
+    i = argmax(i -> abs(v[i]), 1:N)
+    tᵢ = t[i]
+    fᵢ = abs(v[i])
+    verbosity > 0 && @info "Initial maximum amplitude" i tᵢ fᵢ
+
+    # Since the electric field is minus the time derivative of the
+    # vector potential, we know that when the vector potential is
+    # extremized, the electric field has a root. The converse is not
+    # necessarily true (only for purely monochromatic fields is that
+    # guaranteed), but we assume it is.
+    tₘₐₓ = find_zero(f′, tᵢ)
+    fₘₐₓ = abs(f(tₘₐₓ))
+
+    verbosity > 0 && @info "Optimized maximum amplitude" tₘₐₓ fₘₐₓ fₘₐₓ ≥ fᵢ
+
+    if fₘₐₓ < fᵢ
+        @warn "Was not able to find the proper maximum"
+        return fᵢ
+    end
+
+    fₘₐₓ
+end
+
+function amplitude(f::CubicHermiteSplineField; kwargs...)
+    # Since the electric field is minus the time derivative of the
+    # vector potential, we know that when the vector potential is
+    # extremized, the electric field has a root. The converse is not
+    # necessarily true (only for purely monochromatic fields is that
+    # guaranteed), but we assume it is.
+    _maximize_interpolation(f.t, f.Aₜ,
+                            Base.Fix1(field_amplitude, f),
+                            (Base.Fix1(vector_potential, f),
+                             t -> -field_amplitude(f, t));
+                            kwargs...)
+end
+
+function ponderomotive_potential(f::CubicHermiteSplineField; kwargs...)
+    # Since the electric field is minus the time derivative of the
+    # vector potential, we know that when the vector potential is
+    # extremized, the electric field has a root. We therefore look for
+    # a root of the electric field, starting from the time we found
+    # above.
+    Aₘₐₓ = _maximize_interpolation(f.t, f.A,
+                                   Base.Fix1(vector_potential, f),
+                                   t -> -field_amplitude(f, t);
+                                   kwargs...)
+
+    Aₘₐₓ^2/4
+end
+
+function free_oscillation_amplitude(f::CubicHermiteSplineField; kwargs...)
+    ∫A = approximate_integral(f.t, f.A)
+    ∫Af = t -> cubic_hermite_interpolation(f.t, ∫A, f.A, t)
+
+    _maximize_interpolation(f.t, ∫A,
+                            ∫Af,
+                            (Base.Fix1(vector_potential, f),
+                             t -> -field_amplitude(f, t));
+                            kwargs...)
+end
+
 # * Cubic Hermite spline interpolation
 
 function cubic_hermite_interpolation(xp::AbstractVector, f::AbstractVector, fₓ::AbstractVector, x::Number)
@@ -101,6 +179,9 @@ function cubic_hermite_interpolation(xp::AbstractVector, f::AbstractVector, fₓ
     h₀₀*f[i] + h₁₀*δx*fₓ[i] + h₀₁ * f[i+1] + h₁₁*δx*fₓ[i+1]
 end
 
+cubic_hermite_interpolation(xp::AbstractVector, f::AbstractVector, fₓ::AbstractVector, x::AbstractVector) =
+    cubic_hermite_interpolation.(Ref(xp), Ref(f), Ref(fₓ), x)
+
 # * FFT derivatives/integrals
 
 # The derivative/integral of a real function is real, but using the

src/strong_field_properties.jl

@@ -75,7 +75,7 @@ free_oscillation_amplitude(F::SumField) =
 
 # * Print info
 
-function show_strong_field_properties(io::IO, F::Union{LinearField,TransverseField})
+function show_strong_field_properties(io::IO, F::AbstractField)
     Uₚ = austrip(ponderomotive_potential(F))
     α = austrip(free_oscillation_amplitude(F))
     printfmt(io, "Uₚ = {1:.4f} Ha = {2:s} => α = {3:.4f} Bohr = {4:s}",