Simplified ramps, down-ramp, ensure continuity (#16)

Project.toml

@@ -1,7 +1,7 @@
 name = "ElectricFields"
 uuid = "2f84ce32-9bb1-11e8-0d9f-3dce90a4beca"
 authors = ["Stefanos Carlström <[email protected]>"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

src/field_types.jl

@@ -681,14 +681,17 @@ Constant field of
              f.E₀, au2si_round(f.E₀, u"V/m"))
 end
 
-field_amplitude(f::ConstantField, t::Number) =
-    f.E₀*(0 ≤ t && t ≤ f.tmax)
-
 intensity(f::ConstantField, t::Number) = field_amplitude(f, t)^2
 
 function vector_potential(f::ConstantField{T}, t::Number) where T
-    t = clamp(t, zero(T), f.tmax)
-    -f.E₀*t
+    z = zero(T)
+    if t < z
+        z
+    elseif t < f.tmax
+        -f.E₀*t
+    else
+        -f.E₀*f.tmax
+    end
 end
 
 intensity(f::ConstantField) = f.E₀^2
@@ -712,57 +715,65 @@ field_types[:constant] = ConstantField
 # * Ramps
 
 @doc raw"""
-    Ramp(tmax, E₀, f, g, name)
+    Ramp(tmax, E₀, f, name)
 
 The field amplitude of a ramp is defined as
 
 ```math
 F(t) = \begin{cases}
-E_0f(\tau), & 0 \leq t_{\textrm{max}}, \\
+E_0f'(\tau), & 0 \leq t_{\textrm{max}}, \\
 0, \textrm{else},
 \end{cases}
 \implies
 A(t) = \begin{cases}
--E_0t_{\textrm{max}}g(\tau), & 0 \leq t \leq t_{\textrm{max}},\\
+-E_0t_{\textrm{max}}f(\tau), & 0 \leq t \leq t_{\textrm{max}},\\
 A(t_{\textrm{max}}), & t_{\textrm{max}} < t,\\
 0, & \textrm{else},
 \end{cases}
 \quad
 \tau = \frac{t}{t_{\textrm{max}}}.
 ```
 
-To define a new ramp, one thus only needs to define a pair of
-functions on the unit interval, ``f(\tau)`` that rises from ``0`` to
-``1``, and its integral ``g(\tau)``. Similar to
-[`ConstantField`](@ref), `Ramp` is a _non-propagating_ field, but is
-realizable in e.g. a capacitor.
+To define a new ramp, one thus only needs to define one function on
+the unit interval, ``f(\tau)`` whose derivative ``f'(\tau)`` rises
+from ``0`` to ``1``. Similar to [`ConstantField`](@ref), `Ramp` is a
+_non-propagating_ field, but is realizable in e.g. a capacitor.
 
 Three kinds of ramps are predefined, `:linear_ramp`,
 `:parabolic_ramp`, `:sin²_ramp` (with the alias `:sin2_ramp`).
 
-# Example
+# Examples
 
 ```jldoctest
 julia> @field(F) do
            E₀ = 1.0
            tmax = 4.0u"fs"
            kind = :sin²_ramp
        end
-sin² ramp of
+sin² up-ramp of
+  - 165.3655 jiffies = 4.0000 fs duration, and
+  - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
+
+julia> @field(F) do
+           E₀ = 1.0
+           tmax = 4.0u"fs"
+           kind = :linear_ramp
+           ramp = :down
+       end
+Linear down-ramp of
   - 165.3655 jiffies = 4.0000 fs duration, and
   - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
 ```
 """
-struct Ramp{T,EFunc,AFunc} <: AbstractField
+struct Ramp{T,AFunc} <: AbstractField
     tmax::T
     E₀::T
-    f::EFunc
-    g::AFunc
+    f::AFunc
     name::String
     params::Dict{Symbol, Any}
 end
 
-function Ramp(f, g, name, field_params)
+function Ramp(f, name, field_params)
     test_field_parameters(field_params, [:I₀, :E₀])
     test_field_parameters(field_params, [:tmax])
 
@@ -778,41 +789,50 @@ function Ramp(f, g, name, field_params)
 
     @unpack tmax, E₀ = field_params
 
-    Ramp(austrip(tmax), austrip(E₀), f, g, name, field_params)
+    r = get(field_params, :ramp, :up)
+    a,suffix = if r == :up
+        f, "up"
+    elseif r == :down
+        A = f(1)
+        τ -> A - f(1-τ), "down"
+    else
+        throw(ArgumentError("Unknown :ramp kind $(r)"))
+    end
+
+    Ramp(austrip(tmax), austrip(E₀), a, name*" "*suffix*"-", field_params)
 end
 
 function Base.show(io::IO, f::Ramp)
     printfmt(io, """
-$(f.name) ramp of
+$(f.name)ramp of
   - {1:.4f} jiffies = {2:s} duration, and
   - E₀ = {3:.4e} au = {4:s}""",
              f.tmax, au2si_round(f.tmax, u"s"),
              f.E₀, au2si_round(f.E₀, u"V/m"))
 end
 
 ElectricFields.field_types[:linear_ramp] =
-    p -> Ramp(identity, t -> t^2/2, "Linear", p)
+    p -> Ramp(t -> t^2/2, "Linear", p)
 
 ElectricFields.field_types[:parabolic_ramp] =
-    p -> Ramp(t -> 2t - t^2, t -> t^2 - t^3/3, "Parabolic", p)
+    p -> Ramp(t -> t^2 - t^3/3, "Parabolic", p)
 
 ElectricFields.field_types[:sin2_ramp] =
     ElectricFields.field_types[:sin²_ramp] =
-    p -> Ramp(t -> sinpi(t/2)^2, t -> t/2 - sinpi(t)/2π, "sin²", p)
-
-function field_amplitude(f::Ramp, t::Number)
-    if 0 ≤ t ≤ f.tmax
-        f.E₀*f.f(t/f.tmax)
-    else
-        zero(f.E₀)
-    end
-end
+    p -> Ramp(t -> t/2 - sinpi(t)/2π, "sin²", p)
 
 intensity(f::Ramp, t::Number) = field_amplitude(f, t)^2
 
 function vector_potential(f::Ramp{T}, t::Number) where T
-    t = clamp(t, zero(T), f.tmax)/f.tmax
-    -f.E₀*f.tmax*f.g(t)
+    z = zero(T)
+    if t < z
+        z
+    elseif t < f.tmax
+        τ = clamp(t, zero(T), f.tmax)/f.tmax
+        -f.E₀*f.tmax*f.f(τ)
+    else
+        -f.E₀*f.tmax*f.f(one(T))
+    end
 end
 
 polarization(::Ramp) = LinearPolarization()

test/field_types.jl

@@ -174,31 +174,61 @@
             kind = :linear_ramp
         end
 
+       @field(F2) do
+            I₀ = 4.0
+            tmax = 4.0
+            kind = :linear_ramp
+            ramp = :down
+        end
+
         withenv("UNITFUL_FANCY_EXPONENTS" => true) do
             @test string(F) == """
-                               Linear ramp of
+                               Linear up-ramp of
                                  - 4.0000 jiffies = 96.7554 as duration, and
                                  - E₀ = 2.0000e+00 au = 1.0284 TV m⁻¹"""
+            @test string(F2) == """
+                                Linear down-ramp of
+                                  - 4.0000 jiffies = 96.7554 as duration, and
+                                  - E₀ = 2.0000e+00 au = 1.0284 TV m⁻¹"""
         end
         withenv("UNITFUL_FANCY_EXPONENTS" => false) do
             @test string(F) == """
-                               Linear ramp of
+                               Linear up-ramp of
                                  - 4.0000 jiffies = 96.7554 as duration, and
                                  - E₀ = 2.0000e+00 au = 1.0284 TV m^-1"""
+
+            @test string(F2) == """
+                                Linear down-ramp of
+                                  - 4.0000 jiffies = 96.7554 as duration, and
+                                  - E₀ = 2.0000e+00 au = 1.0284 TV m^-1"""
         end
 
         @test F isa ElectricFields.Ramp
 
         @test field_amplitude(F, -1) == 0
+        @test field_amplitude(F, 0) ≈ 0 rtol=1e-7
         @test field_amplitude(F, 2) ≈ 1 rtol=1e-7
-        @test field_amplitude(F, 4) ≈ 2 rtol=1e-7
+        @test field_amplitude(F, 4-2eps()) ≈ 2 rtol=1e-7
         @test field_amplitude(F, 5) == 0
 
+        @test field_amplitude(F2, -1) == 0
+        @test field_amplitude(F2, 0) ≈ 2 rtol=1e-7
+        @test field_amplitude(F2, 2) ≈ 1 rtol=1e-7
+        @test field_amplitude(F2, 4-2eps()) ≈ 0 atol=1e-10
+        @test field_amplitude(F2, 5) == 0
+
         @test intensity(F, -1) == 0
+        @test intensity(F, 0) ≈ 0 rtol=1e-7
         @test intensity(F, 2) ≈ 1 rtol=1e-7
-        @test intensity(F, 4) ≈ 4 rtol=1e-7
+        @test intensity(F, 4-2eps()) ≈ 4 rtol=1e-7
         @test intensity(F, 5) == 0
 
+        @test intensity(F2, -1) == 0
+        @test intensity(F2, 0) ≈ 4 rtol=1e-7
+        @test intensity(F2, 2) ≈ 1 rtol=1e-7
+        @test intensity(F2, 4-2eps()) ≈ 0 atol=1e-10
+        @test intensity(F2, 5) == 0
+
         @test field_amplitude(F, 0, 4.0) ≈ 4.0 rtol=1e-7
 
         @test polarization(F) == LinearPolarization()
@@ -212,15 +242,22 @@
 
         @test dimensions(F) == 1
 
-        @testset "Ramp kind = $(kind)" for kind in (:linear_ramp,
-                                                    :parabolic_ramp,
-                                                    :sin²_ramp)
+        @testset "Ramp kind = $(kind)" for (kind,half) in ((:linear_ramp, 0.5),
+                                                           (:parabolic_ramp, 0.75),
+                                                           (:sin²_ramp, 0.5))
             @field(R) do
                 E₀ = 1.0
                 tmax = 4.0u"fs"
                 kind = kind
             end
 
+            @field(R2) do
+                E₀ = 1.0
+                tmax = 4.0u"fs"
+                kind = kind
+                ramp = :down
+            end
+
             @field(C) do
                 E₀ = 1.0
                 tmax = 4.0u"fs"
@@ -234,6 +271,9 @@
             @test iszero(field_amplitude(F, s.left-1))
             @test iszero(field_amplitude(F, s.right+1))
 
+            @test field_amplitude(R, austrip(2u"fs")) ≈ half
+            @test field_amplitude(R2, austrip(2u"fs")) ≈ half
+
             @test vector_potential(F, s.right+1) ≈ vector_potential(F, s.right) rtol=1e-7
 
             t = range(0, stop=1.1span(F).right, length=1000)