Dispersed fields (#32)

* Implemented B-spline field

* Generalized Sellmeier equations, multiaxial media

* Implemented dispersed fields

Implementation of anisotropic crystals (through the index) ellipsoid
still remains.

* Fix broken timeaxis

* Dispersion tests; truncate time interval to minimize B-spline expansion

* Fix documentation

* Dispersion plots for documentation

* Don't crash on arbitrarily polarized fields

* Julia 1.6 compatible multiplication

* Fix documentation build

* Fix doctests

* Implemented dispersion for multiaxial crystals

* Increase test coverage

* Fixed test

* Increased test coverage

* Better test of chirped pulses

* Increased test coverage

* Improved documentation for dispersed fields

* Efficient evaluation of BSplineField:s for a vector of times

* Fixed plot label position

* Speling mistak

* Fixed factor in derivative of Kaiser windows; test using AD

* Allow changing ramps of TrapezoidalEnvelope

.github/workflows/ci.yml

@@ -75,6 +75,6 @@ jobs:
           julia --project=docs -e '
             using Documenter: DocMeta, doctest
             using ElectricFields
-            DocMeta.setdocmeta!(ElectricFields, :DocTestSetup, :(using ElectricFields, Unitful); recursive=true)
+            DocMeta.setdocmeta!(ElectricFields, :DocTestSetup, :(using ElectricFields, Unitful, UnitfulAtomic); recursive=true)
             doctest(ElectricFields)'
 

Project.toml

@@ -13,8 +13,10 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 

docs/Project.toml

@@ -3,6 +3,7 @@ CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 

docs/make.jl

@@ -6,7 +6,6 @@ isdefined(Main, :NOPLOTS) && NOPLOTS || include("plots.jl")
 makedocs(;
     modules=[ElectricFields],
     authors="Stefanos Carlström <[email protected]> and contributors",
-    repo="https://github.com/jagot/ElectricFields.jl/blob/{commit}{path}#L{line}",
     sitename="ElectricFields.jl",
     format=Documenter.HTML(;
         prettyurls=get(ENV, "CI", "false") == "true",
@@ -35,9 +34,11 @@ makedocs(;
         "Envelopes" => "envelopes.md",
         "Carriers" => "carriers.md",
         "Field properties" => "properties.md",
+        "Dispersion" => "dispersion.md",
         "Reference" => "reference.md",
     ],
-    doctest=false
+    doctest=false,
+    checkdocs=:exports
 )
 
 deploydocs(;

docs/plots.jl

@@ -1,12 +1,23 @@
 using ElectricFields
 using Unitful
+using UnitfulAtomic
 using FFTW
+using Statistics
 
 using PythonPlot
 using Jagot
 using Jagot.plotting
 plot_style("ggplot")
 
+function wind_transf(t, x, w)
+    f = rfftfreq(length(t),1.0/(t[2]-t[1]))
+    ecat(a,b) = cat(a,b,dims=ndims(x)+1)
+    Cₓ = reduce(ecat, x.*circshift(w,i) for i = 1:length(t))
+    Wₓ = rfft(Cₓ, 1)
+    f,Wₓ
+end
+gabor(t, x, σ) = wind_transf(t, x, fftshift(exp.(-(t .- mean(t)).^2/2σ^2)))
+
 function savedocfig(name,dir="figures")
     fig = gcf()
     filename = joinpath(@__DIR__, "src", dir, "$(name).svg")
@@ -190,7 +201,7 @@ function apodized_field()
 
     w = ElectricFields.window_value.(Fw.window, 1.0, 14.0, t)
 
-   savedfigure("apodized_field", figsize=(7,8)) do
+    savedfigure("apodized_field", figsize=(7,8)) do
         csubplot(211, nox=true) do
             plot(tplot, Fv)
             plot(tplot, Fwv)
@@ -237,6 +248,130 @@ function apodizing_windows()
     end
 end
 
+function chirped_field()
+    τ = austrip(3u"fs")
+    η = austrip(5.0u"fs^2")
+
+    @field(F) do
+        λ = 800u"nm"
+        I₀ = 1.0
+        τ = τ
+        σoff = 4.0
+        σmax = 6.0
+        env = :trunc_gauss
+        ϕ = π
+    end
+    Fc = chirp(F, η)
+
+    γ = τ^2/(8log(2))
+    ω₀ = photon_energy(F)
+
+    t = range(span(Fc), length=1500)
+    tplot = ustrip(auconvert.(u"fs", 1))*t
+
+    Av = vector_potential(F, t)
+    Fv = field_amplitude(F, t)
+
+    Frec = @time field_amplitude(Fc, t)
+    Arec = @time vector_potential(Fc, t)
+
+    freq,G = gabor(t, Frec, 1austrip(period(F)))
+    fsel = ind(freq, 0):ind(freq, 2ω₀/2π)
+
+    savedfigure("chirped_field", figsize=(7,8)) do
+        csubplot(311) do
+            plot(tplot, Av, "k")
+            plot(tplot, Arec)
+            xlabel(L"$t$ [fs]")
+            axes_labels_opposite(:x)
+            ylabel(L"A(t)")
+        end
+        csubplot(312,nox=true) do
+            plot(tplot, Fv, "k")
+            plot(tplot, Frec)
+            ylabel(L"F(t)")
+        end
+        csubplot(313) do
+            plot_map(tplot, freq[fsel], abs2.(G[fsel,:]))
+            hline(ω₀/2π, color="white")
+            yl = ylim()
+            plot(tplot, (ω₀ .+ 1/2*η/(γ^2 + η^2)*t)/2π, color="white")
+            ylim(yl)
+            xlabel(L"$t$ [fs]")
+            ylabel(L"\omega/2\pi")
+        end
+    end
+end
+
+function dispersed_field()
+    τ = austrip(3u"fs")
+
+    @field(F) do
+        λ = 800u"nm"
+        I₀ = 1.0
+        τ = τ
+        σoff = 4.0
+        σmax = 6.0
+        env = :trunc_gauss
+        ϕ = π
+    end
+    F = rotate(F, ElectricFields.compute_rotation((π/3, [0.4,1,0])))
+
+    ω₀ = photon_energy(F)
+
+    de = Crystal(KTP, 12u"μm", ω₀=ω₀)
+    Fc = DispersedField(F, de, spline_order=3, verbosity=4)
+
+    t = timeaxis(Fc)
+    tplot = ustrip(auconvert.(u"fs", 1))*t
+
+    Fv = field_amplitude(F, t)
+    Av = vector_potential(F, t)
+
+    Frec = @time field_amplitude(Fc, t)
+    Arec = @time vector_potential(Fc, t)
+
+    savedfigure("dispersed_field", figsize=(7,8)) do
+        csubplot(3,2,(1,1), nox=true) do
+            plot(tplot, Fv[:,1], "k")
+            plot(tplot, Frec[:,1])
+            ylabel(L"F_x(t)")
+        end
+        csubplot(3,2,(1,2), nox=true) do
+            plot(tplot, Av[:,1], "k")
+            plot(tplot, Arec[:,1])
+            ylabel(L"A_x(t)")
+            axes_labels_opposite(:y)
+        end
+
+        csubplot(3,2,(2,1), nox=true) do
+            plot(tplot, Fv[:,2], "k")
+            plot(tplot, Frec[:,2])
+            ylabel(L"F_y(t)")
+        end
+        csubplot(3,2,(2,2), nox=true) do
+            plot(tplot, Av[:,2], "k")
+            plot(tplot, Arec[:,2])
+            ylabel(L"A_y(t)")
+            axes_labels_opposite(:y)
+        end
+
+        csubplot(3,2,(3,1)) do
+            plot(tplot, Fv[:,3], "k")
+            plot(tplot, Frec[:,3])
+            ylabel(L"F_z(t)")
+            xlabel(L"$t$ [fs]")
+        end
+        csubplot(3,2,(3,2)) do
+            plot(tplot, Av[:,3], "k")
+            plot(tplot, Arec[:,3])
+            ylabel(L"A_z(t)")
+            axes_labels_opposite(:y)
+            xlabel(L"$t$ [fs]")
+        end
+    end
+end
+
 macro echo(expr)
     println(expr)
     :(@time $expr)
@@ -250,3 +385,5 @@ mkpath(fig_dir)
 @echo index_spectrum_example()
 @echo apodized_field()
 @echo apodizing_windows()
+@echo chirped_field()
+@echo dispersed_field()