Added smooth windows to limit time integrals

Author: jagot

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 ElectricFields = "2f84ce32-9bb1-11e8-0d9f-3dce90a4beca"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@@ -18,6 +19,7 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 [compat]
 ElectricFields = "0.1, 0.2"
 FastGaussQuadrature = "0.4.7, 0.5"
+FillArrays = "1"
 HCubature = "1.5"
 ProgressMeter = "1.5,1.6"
 SpecialFunctions = "1.3,2"

src/StrongFieldApproximation.jl

@@ -7,6 +7,7 @@ using FastGaussQuadrature
 
 using LinearAlgebra
 using StaticArrays
+using FillArrays
 
 using ProgressMeter
 using TimerOutputs
@@ -20,6 +21,7 @@ include("momentum_grid.jl")
 include("units.jl")
 include("field.jl")
 include("volkov.jl")
+include("windows.jl")
 
 include("channels.jl")
 include("system.jl")

src/dyson_expansions.jl

@@ -18,7 +18,8 @@ function recombination(system::System, diagram::Diagram, 𝐩, 𝐀, i)
     end
 end
 
-function integrate_diagram(::Type{Amp}, system::System, diagram::Diagram, iref, 𝐤=nothing; memory=typemax(Int), imin=1,
+function integrate_diagram(::Type{Amp}, system::System, diagram::Diagram, iref, 𝐤=nothing;
+                           window=flat_window(), imin=1,
                            to=TimerOutput(), verbosity=1) where Amp
     ions, unique_momenta, momenta, indeterminate_momenta, order = @timeit to "Analyze diagram" analyze_diagram(system, diagram)
 
@@ -39,10 +40,18 @@ function integrate_diagram(::Type{Amp}, system::System, diagram::Diagram, iref,
         ctT = complex(eltype(system.t))
 
         is = vcat(iref, zeros(Int, order))
+
+        memory = length(window)
+        if memory < iref
+            window = vcat(window, zeros(eltype(window), iref))
+        end
     end
 
     @timeit to "Time loop" begin
         for i in TelescopeIterator(max(1,imin):iref-1, order, memory)
+            windw = @timeit "Window" prod(window[(i[1] + 1) .- i])
+            iszero(windw) && continue
+
             is[2:end] .= i
             # is = vcat(iref, i)
             # is = (iref,i...)
@@ -89,7 +98,9 @@ function integrate_diagram(::Type{Amp}, system::System, diagram::Diagram, iref,
                 end
             end
 
-            amplitude += ∂a
+            @timeit to "Accumulate amplitude" begin
+                amplitude += ∂a * windw
+            end
         end
     end
 
@@ -122,7 +133,7 @@ function photoelectron_spectrum(k::AbstractArray{T},
 end
 
 function photoelectron_spectrum(k, args...; kwargs...)
-    system = System(args...; kwargs...)
+    system = System(args...)
     photoelectron_spectrum(k, system, Diagram(system); kwargs...)
 end
 
@@ -140,18 +151,19 @@ function induced_dipole(system::System, diagram::Diagram; verbosity = 1, kwargs.
     DT = eltype(system.𝐀)
     𝐝 = zeros(DT, length(t))
 
-    memory = get(kwargs, :memory, typemax(Int))
+    memory = length(get(kwargs, :window, flat_window()))
 
     @showprogress for (i,t) in enumerate(t)
-        𝐝̃ = integrate_diagram(DT, system, diagram, i; imin=i-memory, verbosity=verbosity-1, kwargs...)
+        𝐝̃ = integrate_diagram(DT, system, diagram, i; imin=i-memory,
+                             verbosity=verbosity-1, kwargs...)
         𝐝[i] = 2real(𝐝̃)
     end
 
     𝐝
 end
 
 function induced_dipole(args...; kwargs...)
-    system = System(args...; kwargs...)
+    system = System(args...)
     diagram = Diagram([(1,0),(1,0)], system)
     induced_dipole(system, diagram; kwargs...)
 end

src/system.jl

@@ -73,8 +73,8 @@ System(ionization_channels::AbstractVector{<:IonizationChannel},
            System(ionization_channels, nothing, NoCouplings,
                   F, fs)
 
-System(Iₚ::Number, args...; kwargs...) =
-    System([IonizationChannel(Iₚ, args...)], args...; kwargs...)
+System(Iₚ::Number, args...) =
+    System([IonizationChannel(Iₚ, args...)], args...)
 
 num_channels(system::System) = length(system.ionization_channels)
 

src/windows.jl

@@ -0,0 +1,40 @@
+flat_window(n=typemax(Int)) = Trues(n)
+
+#=
+
+Smooth turn-off of scalar functions, following Eqs. (18–21) of
+
+- Becke, A. D. (1988). A Multicenter Numerical Integration Scheme for
+  Polyatomic Molecules. The Journal of Chemical Physics, 88(4),
+  2547–2553. http://dx.doi.org/10.1063/1.454033
+
+=#
+
+function becke_smoother(μ, k::Integer)
+    p = μ -> 3μ/2 - μ^3/2
+    f = μ
+    for i = 1:k
+        f = p(f)
+    end
+    (1-f)/2
+end
+
+function becke_smoother(x, a, b, k)
+    if x < a
+        one(x)
+    elseif x > b
+        zero(x)
+    else
+        becke_smoother(2*(x-a)/(b-a)-1, k)
+    end
+end
+
+function smooth_window(::Type{T}, nflat, nsmooth, k) where T
+    w = ones(T, nflat+nsmooth)
+    for i = eachindex(w)
+        w[i] = becke_smoother(i, nflat, nflat+nsmooth, k)
+    end
+    w
+end
+
+export flat_window, smooth_window