Analytic computation of spectra of some fields (#22)

* Analytic computation of spectra of some fields

* Analytic spectra for arithmetic fields

* Fixed normalization of convolution

* Implemented (normalized) FFTs of fields

* WIP Analytic spectrum for ConstantField

* Spectra for transverse fields

* Better test of FWHM in view of unit conversion inaccuracy (#24)

* Added test of analytic and numeric field spectra

* Bump version to 0.2.0 since new Gaussian exponents can lead to different results

* Increased test coverage

* Removed unused methods

* Fixed sign of spectrum

* Some more tests

* Fix for PrettyTables < 2

* Update compat for Documenter

* Drop Julia < 1.6

* Fix doctests

* Some docstrings

* Preview PR docs

* Added spectrum example plot

* Missing plot and imports

* Missing import

* Increase test coverage

* Fixed file name

* Improved plot

.github/workflows/ci.yml

@@ -10,7 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
+          - '1.6'
           - '1'
         os:
           - ubuntu-latest
@@ -74,8 +74,6 @@ jobs:
             Pkg.develop(PackageSpec(path=pwd()))
             Pkg.instantiate()'
       - name: Run Doctests
-        # See https://github.com/actions/toolkit/issues/399
-        continue-on-error: true
         run: |
           julia --project=docs -e '
             using Documenter: DocMeta, doctest

Project.toml

@@ -1,9 +1,10 @@
 name = "ElectricFields"
 uuid = "2f84ce32-9bb1-11e8-0d9f-3dce90a4beca"
 authors = ["Stefanos Carlström <[email protected]>"]
-version = "0.1.6"
+version = "0.2.0"
 
 [deps]
+AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -20,17 +21,19 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 FFTW = "1.3"
 Formatting = "0.4"
 ForwardDiff = "0.10"
-IntervalSets = "0.5.3,0.7"
+IntervalSets = "0.7"
 Optim = "1.3"
 Parameters = "0.12"
+PrettyTables = "2"
 Roots = "2"
 StaticArrays = "1.1"
 Unitful = "1.7"
 UnitfulAtomic = "1.0"
-julia = "1.5"
+julia = "1.6"
 
 [extras]
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Test", "PrettyTables"]

docs/Project.toml

@@ -1,9 +1,10 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
-Documenter = "0.26"
+Documenter = "0.27"
 Unitful = "1.7"
 UnitfulAtomic = "1.0"

docs/make.jl

@@ -41,4 +41,5 @@ makedocs(;
 
 deploydocs(;
     repo="github.com/jagot/ElectricFields.jl",
+    push_preview = true,
 )

docs/plots.jl

@@ -1,7 +1,9 @@
 using ElectricFields
 using Unitful
+using FFTW
 
 using PyPlot
+using Jagot
 using Jagot.plotting
 plot_style("ggplot")
 
@@ -98,6 +100,71 @@ function index_polarized_example()
     end
 end
 
+function index_spectrum_example()
+    @field(F) do
+        I₀ = 1.0
+        T = 1.0
+        τ = 2.0
+        σmax = 6.0
+    end
+
+    t = timeaxis(F)
+    Fv = field_amplitude(F, t);
+    Av = vector_potential(F, t);
+
+    ω = fftshift(fftω(t));
+    # We need to undo the phase, since the FFT does not care that
+    # pulse is centred around zero.
+    F̂v = exp.(im*ω*t[1]) .* fftshift(nfft(F, t), 1);
+    Âv = exp.(im*ω*t[1]) .* fftshift(nfft_vector_potential(F, t), 1);
+
+    F̂v_exact = field_amplitude_spectrum(F, ω);
+    Âv_exact = vector_potential_spectrum(F, ω);
+
+    sel = ind(ω, -20):ind(ω, 20)
+
+    savedfigure("index_spectrum_example", figsize=(8,10)) do
+        csubplot(311) do
+            plot(t, Fv, label=L"F(t)")
+            plot(t, Av, label=L"A(t)")
+            axes_labels_opposite(:x)
+            xlabel(L"t/T")
+            legend(loc=1)
+        end
+        csubplot(323, nox=true) do
+            semilogy(ω[sel], abs2.(F̂v[sel,:]), label=L"$|F(\omega)|^2$, FFT")
+            semilogy(ω[sel], abs2.(Âv[sel,:]), label=L"$|A(\omega)|^2$, FFT")
+            ax = axis()
+            semilogy(ω[sel], abs2.(F̂v_exact[sel,:]), "--", label=L"$|F(\omega)|^2$, exact")
+            semilogy(ω[sel], abs2.(Âv_exact[sel,:]), "--", label=L"$|A(\omega)|^2$, exact")
+            axis(ax)
+            legend(loc=3)
+        end
+        csubplot(325) do
+            semilogy(ω[sel], abs.(abs2.(F̂v[sel,:])-abs2.(F̂v_exact[sel,:])), label=L"||\hat{F}_{\mathrm{FFT}}|^2-|\hat{F}_{\mathrm{exact}}|^2|")
+            semilogy(ω[sel], abs.(abs2.(Âv[sel,:])-abs2.(Âv_exact[sel,:])), label=L"||\hat{A}_{\mathrm{FFT}}|^2-|\hat{A}_{\mathrm{exact}}|^2|")
+            xlabel(L"$\omega$ [rad/jiffies]")
+            legend(loc=3)
+        end
+        csubplot(324, nox=true) do
+            plot(ω[sel], (angle.(F̂v[sel,:])), label=L"$\arg\{F(\omega)\}$, FFT")
+            plot(ω[sel], (angle.(Âv[sel,:])), label=L"$\arg\{A(\omega)\}$, FFT")
+            plot(ω[sel], (angle.(F̂v_exact[sel,:])), "--", label=L"$\arg\{F(\omega)\}$, exact")
+            plot(ω[sel], (angle.(Âv_exact[sel,:])), "--", label=L"$\arg\{A(\omega)\}$, exact")
+            axes_labels_opposite(:y)
+            legend(loc=3)
+            π_labels(:y)
+        end
+        csubplot(326) do
+            semilogy(ω[sel], abs2.(F̂v[sel,:] - F̂v_exact[sel,:]), label=L"|\hat{F}_{\mathrm{FFT}}-\hat{F}_{\mathrm{exact}}|")
+            semilogy(ω[sel], abs2.(Âv[sel,:] - Âv_exact[sel,:]), label=L"|\hat{A}_{\mathrm{FFT}}-\hat{A}_{\mathrm{exact}}|")
+            axes_labels_opposite(:y)
+            xlabel(L"$\omega$ [rad/jiffies]")
+            legend(loc=3)
+        end
+    end
+end
+
 macro echo(expr)
     println(expr)
     :(@time $expr)
@@ -107,3 +174,4 @@ end
 mkpath("docs/src/figures")
 @echo index_example()
 @echo index_polarized_example()
+@echo index_spectrum_example()

docs/src/index.md

@@ -63,6 +63,8 @@ Fourier transform, will receive better support in a future release.
 
 ## Examples
 
+### Linear polarization
+
 The simplest pulse is created thus:
 
 ```jldoctest index_examples
@@ -75,27 +77,28 @@ julia> @field(IR) do
        σmax = 6.0
        end
 Linearly polarized field with
-  - I₀ = 2.8495e-03 au = 1.0e14 W cm⁻² =>
-    - E₀ = 5.3380e-02 au = 27.4492 GV m⁻¹
+  - I₀ = 2.8495e-03 au = 1.0e14 W cm^-2 =>
+    - E₀ = 5.3380e-02 au = 27.4492 GV m^-1
     - A₀ = 0.9372 au
-  – a Fixed carrier @ λ = 800.0000 nm (T = 2.6685 fs, ω = 0.0570 Ha = 1.5498 eV)
+  – a Fixed carrier @ λ = 800.0000 nm (T = 2.6685 fs, ω = 0.0570 Ha = 1.5498 eV, f = 374.7406 THz)
   – and a Gaussian envelope of duration 6.2000 fs (intensity FWHM; ±6.08σ)
+  – and a bandwidth of 0.0108 Ha = 294.3469 meV ⟺ 71.1728 THz ⟺ 2871.2568 Bohr = 151.9404 nm
   – Uₚ = 0.2196 Ha = 5.9759 eV => α = 16.4562 Bohr = 870.8242 pm
 
 julia> vector_potential(IR, 4.0)
-0.2116055371709056
+0.21160647961322301
 
 julia> field_amplitude(IR, 4.0), field_envelope(IR, 4.0)
-(-0.05194633272360931, 0.05336201848846937)
+(-0.05194703530281701, 0.05336224984361336)
 
 julia> instantaneous_intensity(IR, 4.0), intensity(IR, 4.0)
-(0.0026984214834319233, 0.002847505017163747)
+(0.0026984944767521166, 0.002847529708372214)
 
 julia> span(IR)
--661.9198939608041..661.9198939608041
+-661.9198939608042..661.9198939608042
 
 julia> timeaxis(IR)
--661.9198939608041:1.1041199232040102:661.9198939608041
+-661.9198939608042:1.1041199232040104:661.9198939608042
 ```
 
 ![Simple example](figures/index_example.svg)
@@ -110,14 +113,17 @@ julia> @field(XUV) do
        σmax = 4
        end
 Linearly polarized field with
-  - I₀ = 4.0000e-02 au = 1.40377808e15 W cm⁻² =>
-    - E₀ = 2.0000e-01 au = 102.8441 GV m⁻¹
+  - I₀ = 4.0000e-02 au = 1.40377808e15 W cm^-2 =>
+    - E₀ = 2.0000e-01 au = 102.8441 GV m^-1
     - A₀ = 0.2000 au
-  – a Fixed carrier @ λ = 45.5634 nm (T = 151.9830 as, ω = 1.0000 Ha = 27.2114 eV)
+  – a Fixed carrier @ λ = 45.5634 nm (T = 151.9830 as, ω = 1.0000 Ha = 27.2114 eV, f = 6.5797 PHz)
   – and a Gaussian envelope of duration 3.6283 fs (intensity FWHM; ±4.04σ)
+  – and a bandwidth of 0.0185 Ha = 502.9732 meV ⟺ 121.6184 THz ⟺ 15.9151 Bohr = 842.1896 pm
   – Uₚ = 0.0100 Ha = 272.1138 meV => α = 0.2000 Bohr = 10.5835 pm
 ```
 
+### Arbitrary polarization, composite fields
+
 Other [Envelopes](@ref) and polarizations, as well as some simple arithmetic is possible:
 ```jldoctest index_examples
 julia> @field(A) do
@@ -128,11 +134,12 @@ julia> @field(A) do
            ξ = 1.0
        end
 Transversely polarized field with
-  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
-    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
+  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
+    - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
     - A₀ = 1.0000 au
-  – a Elliptical carrier with ξ = 1.00 (RCP) @ λ = 45.5634 nm (T = 151.9830 as, ω = 1.0000 Ha = 27.2114 eV)
+  – a Elliptical carrier with ξ = 1.00 (RCP) @ λ = 45.5634 nm (T = 151.9830 as, ω = 1.0000 Ha = 27.2114 eV, f = 6.5797 PHz)
   – and a 6.00 cycles cos² envelope
+  – and a bandwidth of Inf Ha = Inf eV ⟺ Inf Hz ⟺ Inf Bohr = Inf m
   – Uₚ = 0.2500 Ha = 6.8028 eV => α = 1.0000 Bohr = 52.9177 pm
 
 julia> @field(B) do
@@ -143,41 +150,88 @@ julia> @field(B) do
            ξ = -1.0
        end
 Transversely polarized field with
-  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
-    - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
+  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
+    - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
     - A₀ = 0.5000 au
-  – a Elliptical carrier with ξ = -1.00 (LCP) @ λ = 22.7817 nm (T = 75.9915 as, ω = 2.0000 Ha = 54.4228 eV)
+  – a Elliptical carrier with ξ = -1.00 (LCP) @ λ = 22.7817 nm (T = 75.9915 as, ω = 2.0000 Ha = 54.4228 eV, f = 13.1594 PHz)
   – and a 6.00 cycles cos² envelope
+  – and a bandwidth of Inf Ha = Inf eV ⟺ Inf Hz ⟺ Inf Bohr = Inf m
   – Uₚ = 0.0625 Ha = 1.7007 eV => α = 0.2500 Bohr = 13.2294 pm
 
 julia> F = A + delay(B, 3/2π)
 ┌ Transversely polarized field with
-│   - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
-│     - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
+│   - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
+│     - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
 │     - A₀ = 1.0000 au
-│   – a Elliptical carrier with ξ = 1.00 (RCP) @ λ = 45.5634 nm (T = 151.9830 as, ω = 1.0000 Ha = 27.2114 eV)
+│   – a Elliptical carrier with ξ = 1.00 (RCP) @ λ = 45.5634 nm (T = 151.9830 as, ω = 1.0000 Ha = 27.2114 eV, f = 6.5797 PHz)
 │   – and a 6.00 cycles cos² envelope
+│   – and a bandwidth of Inf Ha = Inf eV ⟺ Inf Hz ⟺ Inf Bohr = Inf m
 │   – Uₚ = 0.2500 Ha = 6.8028 eV => α = 1.0000 Bohr = 52.9177 pm
 ⊕
 │ Transversely polarized field with
-│   - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
-│     - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
+│   - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
+│     - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
 │     - A₀ = 0.5000 au
-│   – a Elliptical carrier with ξ = -1.00 (LCP) @ λ = 22.7817 nm (T = 75.9915 as, ω = 2.0000 Ha = 54.4228 eV)
+│   – a Elliptical carrier with ξ = -1.00 (LCP) @ λ = 22.7817 nm (T = 75.9915 as, ω = 2.0000 Ha = 54.4228 eV, f = 13.1594 PHz)
 │   – and a 6.00 cycles cos² envelope
+│   – and a bandwidth of Inf Ha = Inf eV ⟺ Inf Hz ⟺ Inf Bohr = Inf m
 │   – Uₚ = 0.0625 Ha = 1.7007 eV => α = 0.2500 Bohr = 13.2294 pm
 └   – delayed by 0.4775 jiffies = 11.5493 as
 
 
 julia> field_amplitude(F, 4.0)
-3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
- -0.8793235934912678
+3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
+ -0.8793235934912674
  -0.0
   0.06802883592577502
 ```
 
 ![Simple polarized example](figures/index_polarized_example.svg)
 
+### Spectra
+
+We can compute the spectrum of a field using either the Fast Fourier
+Transform [`fft`](@ref), [`nfft`](@ref), or for some fields
+analytically [`field_amplitude_spectrum`](@ref):
+
+```jldoctest
+julia> @field(F) do
+           I₀ = 1.0
+           T = 1.0
+           τ = 2.0
+           σmax = 6.0
+       end
+Linearly polarized field with
+  - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
+    - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
+    - A₀ = 0.1592 au
+  – a Fixed carrier @ λ = 7.2516 nm (T = 24.1888 as, ω = 6.2832 Ha = 170.9742 eV, f = 41.3414 PHz)
+  – and a Gaussian envelope of duration 48.3777 as (intensity FWHM; ±7.06σ)
+  – and a bandwidth of 1.3863 Ha = 37.7230 eV ⟺ 9.1214 PHz ⟺ 30.2350 Bohr = 1.6000 nm
+  – Uₚ = 0.0063 Ha = 172.3181 meV => α = 0.0253 Bohr = 1.3404 pm
+
+julia> t = timeaxis(F)
+-6.0:0.010008340283569641:6.0
+
+julia> Fv = field_amplitude(F, t);
+
+julia> Av = vector_potential(F, t);
+
+julia> ω = fftshift(fftω(t));
+
+julia> # We need to undo the phase, since the FFT does not care that
+       # pulse is centred around zero.
+       F̂v = exp.(im*ω*t[1]) .* fftshift(nfft(F, t), 1);
+
+julia> Âv = exp.(im*ω*t[1]) .* fftshift(nfft_vector_potential(F, t), 1);
+
+julia> F̂v_exact = field_amplitude_spectrum(F, ω);
+
+julia> Âv_exact = vector_potential_spectrum(F, ω);
+```
+
+![Simple spectrum example](figures/index_spectrum_example.svg)
+
 ## Reference
 
 ```@index

src/ElectricFields.jl

@@ -6,23 +6,27 @@ using StaticArrays
 using Unitful
 using UnitfulAtomic
 
-# using FFTW
+using AbstractFFTs
+using FFTW
 using ForwardDiff
 using Optim
 using Roots
 
 using IntervalSets
-import IntervalSets: duration
 
 using Parameters
 
 import Base: show
 using Formatting
 
+abstract type AbstractField end
+
 include("units.jl")
 include("rotations.jl")
 include("relation_dsl.jl")
 
+include("spectra.jl")
+
 include("field_types.jl")
 include("time_axis.jl")
 include("carriers.jl")