diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 11a59d8..17f7558 100644 --- a/.github/workflows/ci.yml +++ b/.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)' diff --git a/Project.toml b/Project.toml index f0e207c..1a15cb4 100644 --- a/Project.toml +++ b/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" diff --git a/docs/Project.toml b/docs/Project.toml index c1d08e7..04b25c3 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index bdd3840..2d675d7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,7 +6,6 @@ isdefined(Main, :NOPLOTS) && NOPLOTS || include("plots.jl") makedocs(; modules=[ElectricFields], authors="Stefanos Carlström 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(; diff --git a/docs/plots.jl b/docs/plots.jl index 9e3397f..e6898ab 100644 --- a/docs/plots.jl +++ b/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() diff --git a/docs/src/dispersion.md b/docs/src/dispersion.md new file mode 100644 index 0000000..12fa406 --- /dev/null +++ b/docs/src/dispersion.md @@ -0,0 +1,264 @@ +# Dispersion + +Dispersion is most easily calculated in the frequency domain, where it +simply amounts to multiplication of the spectral amplitudes by a +transfer function: +```math +\hat{\vec{F}}_2(\omega) = +H(\omega) +\hat{\vec{F}}_1(\omega) +``` +where ``H(\omega)`` is the transfer function, ``\hat{\vec{F}}_1(\omega)`` is +the field before dispersion, and ``\hat{\vec{F}}_2(\omega)`` after. Since +``\hat{\vec{F}}(\omega)\equiv-\im\omega\hat{\vec{A}}(\omega)``, the same relation +holds for the vector potential: +```math +\hat{\vec{A}}_2(\omega) = +H(\omega) +\hat{\vec{A}}_1(\omega). +``` + +Dispersion is implemented in ElectricFields.jl via the transform pair +[`rfft`](@ref)/[`irfft`](@ref): + +```math +\vec{A}_1\{t\} +\overset{\texttt{rfft}}{\to} +\hat{\vec{A}}_1\{\omega\} +\to +\hat{\vec{A}}_2\{\omega\} = +H\{\omega\} +\hat{\vec{A}}_1\{\omega\} +\overset{\texttt{irfft}}{\to} +\vec{A}_2\{t\} +``` +where the notation ``\{t\}`` and ``\{\omega\}`` signifies uniformly +spaced vectors, since the FFT is a numeric algorithm. Since +[`field_amplitude`](@ref) and [`vector_potential`](@ref) allow the +evaluation of a field at arbitrary time points (and integrals over +arbitrary time intervals), we fit ``\vec{A}_2\{t\}`` to a +[`BSplineField`](@ref). The next complication is that the +transfer function ``H`` will in general introduce a temporal spread, +such that the compact support of the dispersed field is much wider +than that of the original one. We therefore successively increase the +time span (using [`ElectricFields.find_time_span`](@ref)) we evaluate +``\vec{A}_2\{t\}`` on until it has converged. Before we fit this +function to a [`BSplineField`](@ref), we truncate the time span again to a window +where ``|\vec{A}_2\{t\}|>\epsilon``. + +```@docs +DispersedField +``` + +## Simple dispersive elements + +Any purely dispersive medium (i.e. no loss or gain) can be written as +```math +H(\omega) = \exp[-\im\phi(\omega)] +``` +where the phase can be Taylor-expanded as +```math +\phi(\omega) = +\phi_0 + +a(\omega-\omega_0) + +b(\omega-\omega_0)^2 + +..., +``` +``\phi_0`` is a [`PhaseShift`](@ref), ``a=2\pi t_d`` amounts to a time +[`delay`](@ref) by ``t_d`` (by the Fourier shift theorem), and ``b`` +introduced a [`Chirp`](@ref). + +```julia-repl +julia> @field(F) do + λ = 800u"nm" + I₀ = 1.0 + τ = 3u"fs" + σoff = 4.0 + σmax = 6.0 + env = :trunc_gauss + ϕ = π + end +Linearly polarized field with + - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² => + - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹ + - A₀ = 17.5580 au + – a Fixed carrier @ λ = 800.0000 nm (T = 2.6685 fs, ω = 0.0570 Ha = 1.5498 eV, f = 374.7406 THz); CEP = 1.00π + – and a Truncated Gaussian envelope of duration 124.0241 jiffies = 3.0000 fs (intensity FWHM; turn-off from 5.0959 fs to 7.6439 fs) + – and a bandwidth of 0.0224 Ha = 608.3170 meV ⟺ 147.0904 THz ⟺ 5933.9307 Bohr = 314.0101 nm + – Uₚ = 77.0706 Ha = 2.0972 keV => α = 308.2823 Bohr = 16.3136 nm + +julia> Fc = chirp(F, austrip(5u"fs^2"), verbosity=4) +┌ Info: Finding large enough time span to encompass dispersed field +│ f = +│ Linearly polarized field with +│ - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² => +│ - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹ +│ - A₀ = 17.5580 au +│ – a Fixed carrier @ λ = 800.0000 nm (T = 2.6685 fs, ω = 0.0570 Ha = 1.5498 eV, f = 374.7406 THz); CEP = 1.00π +│ – and a Truncated Gaussian envelope of duration 124.0241 jiffies = 3.0000 fs (intensity FWHM; turn-off from 5.0959 fs to 7.6439 fs) +│ – and a bandwidth of 0.0224 Ha = 608.3170 meV ⟺ 147.0904 THz ⟺ 5933.9307 Bohr = 314.0101 nm +│ – Uₚ = 77.0706 Ha = 2.0972 keV => α = 308.2823 Bohr = 16.3136 nm +│ de = Chirp(b = 8545.5457 = 5.0000 fs², ω₀ = 0.0570 = 1.5498 eV) +│ max_iter = 7 +│ ξ = 2.0 +└ tol = 0.0005 +---------------------------------------------------------------------------------------------------- +1 +t′ = -632.0183332958276:1.1049271561115868:633.1232604519392 +R = 2.095947834009663 +---------------------------------------------------------------------------------------------------- +2 +t′ = -1264.0366665916551:1.1049271561115868:1265.1415937477668 +R = 0.16138039790711892 +---------------------------------------------------------------------------------------------------- +3 +t′ = -2528.0733331833103:1.1049271561115868:2529.178260339422 +R = 0.0013929399886955529 +---------------------------------------------------------------------------------------------------- +4 +t′ = -5056.146666366621:1.1049271561115868:5057.251593522733 +R = 0.00012089884366390845 +┌ Info: Truncated to time interval -896.38155595379 .. 888.6605640985183 +│ a = 1480 +│ b = 3098 +│ cutoff = 0.0014901161193847656 +└ abs_cutoff = 0.014115901307391281 +┌ Info: Generated B-spline +│ num_knots = 324 +└ B = BSpline basis with typename(ElectricFields.LinearKnotSet)(Float64) of order k = 3 (parabolic) on -896.38155595379 .. 888.6605640985183 (324 intervals) +DispersedField: +Linearly polarized field with + - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² => + - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹ + - A₀ = 17.5580 au + – a Fixed carrier @ λ = 800.0000 nm (T = 2.6685 fs, ω = 0.0570 Ha = 1.5498 eV, f = 374.7406 THz); CEP = 1.00π + – and a Truncated Gaussian envelope of duration 124.0241 jiffies = 3.0000 fs (intensity FWHM; turn-off from 5.0959 fs to 7.6439 fs) + – and a bandwidth of 0.0224 Ha = 608.3170 meV ⟺ 147.0904 THz ⟺ 5933.9307 Bohr = 314.0101 nm + – Uₚ = 77.0706 Ha = 2.0972 keV => α = 308.2823 Bohr = 16.3136 nm + – dispersed through Chirp(b = 8545.5457 = 5.0000 fs², ω₀ = 0.0570 = 1.5498 eV) +``` + +![Chirped field](figures/chirped_field.svg) + +In black, we see the original field, in red the chirped field, and the +bottom panel shows a [Gabor +transform](https://en.wikipedia.org/wiki/Gabor_transform) of the +chirped field, along with horizontal line at the carrier frequency +``\omega_0``, and a diagonal line at the expected instantaneous +frequency ``\omega=\omega_0 + \frac{1}{2}\frac{b}{\gamma^2 + b^2}t``, +where ``\gamma = \frac{\tau^2}{8\ln2}``. + +```@docs +ElectricFields.DispersiveElement +PhaseShift +phase_shift(f::ElectricFields.AbstractField, ϕ; kwargs...) +Chirp +chirp +ElectricFields.CascadedDispersiveElement +ElectricFields.find_time_span +``` + +## Dispersive media + +In addition to the simple, analytic dispersive elements listed above, +[`PhaseShift`](@ref) and [`Chirp`](@ref), we also support calculating +the dispersion resulting from pulse propagation through a dispersive +[`ElectricFields.Medium`](@ref), which models a physical medium +phenomenologically. Typically, one employs a [`Sellmeier`](@ref) +equation to describe the variation of the refractive index with the +wavelength of the incident light, and the compute the dispersion +according to +```math +H(\omega) = +\exp[-\im k(\omega)d], +``` +where ``d`` is the thickness of the medium, +```math +k(\omega) = n(\omega)k_0 +``` +is the wavevector in the medium with refractive index +``n(\omega)``. For convenience, we typically subtract the linear +component at the carrier energy ``\omega_0`` according to +```math +\tilde{k}(\omega) = +k(\omega) - +\omega +\left.\frac{\partial k}{\partial\omega}\right|_{\omega_0} +``` +to keep the pulse centred in the frame of reference. + +For an [`IsotropicMedium`](@ref), the square of the refractive index, +``n^2(\lambda)`` does not depend on the direction of polarization. For +an anistropic medium, a [`Crystal`](@ref), different polarization axes +may have different refractive indices. This is modelled using the +[index ellipsoid](https://en.wikipedia.org/wiki/Index_ellipsoid)[^demo], +which obeys +```math +\frac{x^2}{n_a^2} + \frac{y^2}{n_b^2} + \frac{z^2}{n_c^2} = 1. +``` +By choosing the orientation of the crystal appropriately (or +equivalently rotate the field), we can achieve different dispersion +for the different components of the field: + +```julia-repl +julia> @field(F) do + λ = 800u"nm" + I₀ = 1.0 + τ = 3u"fs" + σoff = 4.0 + σmax = 6.0 + env = :trunc_gauss + ϕ = π + end; + +julia> F = rotate(F, ElectricFields.compute_rotation((π/3, [0.4,1,0]))); + +julia> de = Crystal(KTP, 12u"μm", ω₀=photon_energy(F)); + +julia> Fc = DispersedField(F, de) +DispersedField: +Transversely polarized field with + - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² => + - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹ + - A₀ = 17.5580 au + – a LinearTransverseCarrier: Fixed carrier @ λ = 800.0000 nm (T = 2.6685 fs, ω = 0.0570 Ha = 1.5498 eV, f = 374.7406 THz); CEP = 1.00π + – a Truncated Gaussian envelope of duration 124.0241 jiffies = 3.0000 fs (intensity FWHM; turn-off from 5.0959 fs to 7.6439 fs) + – and a rotation of 0.33π about [0.371, 0.928, 0.000] + – and a bandwidth of 0.0224 Ha = 608.3170 meV ⟺ 147.0904 THz ⟺ 5933.9307 Bohr = 314.0101 nm + – Uₚ = 77.0706 Ha = 2.0972 keV => α = 308.2823 Bohr = 16.3136 nm + – dispersed through Crystal(226767.13 Bohr = 12.00 μm of Sellmeier{Float64, Float64, Quantity{Float64, 𝐋², Unitful.FreeUnits{(μm²,), 𝐋², nothing}}, Tuple{Quantity{Float64, 𝐋⁻², Unitful.FreeUnits{(μm⁻²,), 𝐋⁻², nothing}}}}[Sellmeier(1.10468, [0.89342], [2], [0.04438 μm²], [-0.01036 μm⁻²], [2]), Sellmeier(1.14559, [0.87629], [2], [0.0485 μm²], [-0.01173 μm⁻²], [2]), Sellmeier(0.9446, [1.3617], [2], [0.047 μm²], [-0.01491 μm⁻²], [2])], R = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], ∂k∂ω₀ = 0.0134) +``` + +![Dispersed field](figures/dispersed_field.svg) + +[^demo]: [Interactive demonstration of the index ellipsoid.](https://micro.magnet.fsu.edu/primer/java/polarizedlight/ellipsoid/index.html) + +```@docs +ElectricFields.Medium +IsotropicMedium +IsotropicMedium(material, d; ω₀=nothing) +Crystal +Crystal(material, d, R=I; ω₀=nothing) +``` + +### Sellmeier equations + +```@docs +ElectricFields.Sellmeier +BK7 +SiO₂ +Calcite +Quartz +KTP +``` + +## B-splines + +The present B-spline implementation is cannibalized from +[CompactBases.jl](https://github.com/JuliaApproximation/CompactBases.jl), +to avoid dragging in all the dependencies. + +```@docs +BSplineField +ElectricFields.BSpline +``` diff --git a/docs/src/field_types.md b/docs/src/field_types.md index 31c7519..05a03cb 100644 --- a/docs/src/field_types.md +++ b/docs/src/field_types.md @@ -15,6 +15,7 @@ ElectricFields.SumField ElectricFields.WrappedField ElectricFields.NegatedField ElectricFields.DelayedField +delay ElectricFields.PaddedField ElectricFields.WindowedField ElectricFields.phase_shift diff --git a/docs/src/interface.md b/docs/src/interface.md index 7ea6c4a..9f60ed5 100644 --- a/docs/src/interface.md +++ b/docs/src/interface.md @@ -30,4 +30,8 @@ ElectricFields.nfft ElectricFields.fft_vector_potential ElectricFields.nfft_vector_potential ElectricFields.fftω +ElectricFields.AbstractFFTs.rfft +ElectricFields.AbstractFFTs.irfft +ElectricFields.rfft_vector_potential +ElectricFields.rfftω ``` diff --git a/src/ElectricFields.jl b/src/ElectricFields.jl index e76d7dd..46c3250 100644 --- a/src/ElectricFields.jl +++ b/src/ElectricFields.jl @@ -2,6 +2,8 @@ module ElectricFields using LinearAlgebra using StaticArrays +using SparseArrays +using Statistics using Unitful using UnitfulAtomic @@ -36,8 +38,14 @@ include("arithmetic.jl") include("field_dsl.jl") -# include("sellmeier.jl") -# include("dispersed_fields.jl") +include("knot_sets.jl") +include("bsplines.jl") +include("bspline_field.jl") + +include("dispersed_fields.jl") +include("dispersive_elements.jl") +include("sellmeier.jl") +include("dispersive_media.jl") include("strong_field_properties.jl") diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 0186407..6815def 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -113,6 +113,8 @@ end for fun in [:vector_potential, :field_amplitude, :vector_potential_spectrum] @eval $fun(f::SumField, t::Number) = $fun(f.a, t) + $fun(f.b, t) + @eval $fun(f::SumField, t::AbstractVector) = + $fun(f.a, t) + $fun(f.b, t) end polarization(f::SumField) = polarization(f.a) @@ -226,6 +228,7 @@ Base.parent(f::NegatedField) = f.a for fun in [:vector_potential, :vector_potential_spectrum] @eval $fun(f::NegatedField, t::Number) = -$fun(parent(f), t) + @eval $fun(f::NegatedField, t::AbstractVector) = -$fun(parent(f), t) end rotate(f::NegatedField, R) = NegatedField(rotate(f.a, R)) @@ -298,6 +301,8 @@ end for fun in [:vector_potential, :field_amplitude, :intensity] @eval $fun(f::DelayedField, t::Number) = $fun(f.a, t-f.t₀) + @eval $fun(f::DelayedField, t::AbstractVector) = + $fun(f.a, t .- f.t₀) end vector_potential_spectrum(f::DelayedField, ω) = @@ -322,9 +327,26 @@ rotate(f::DelayedField, R) = DelayedField(rotate(f.a, R), f.t₀) # Convention for delayed fields: a field delayed by a /positive/ # time, comes /later/, i.e. we write \(f(t-\delta t)\). +""" + delay(a, t₀) + +Delay the field `a` by the time `t₀`. +""" delay(a::AbstractField, t₀::Union{<:Real,<:Unitful.Time}) = DelayedField(a, austrip(t₀)) + +""" + delay(a, ϕ) + +Delay the field `a` by the phase `ϕ` with respect to its +[`period`](@ref). +""" delay(a::AbstractField, ϕ::Quantity{<:Real, NoDims}) = delay(a, (ϕ/(2π*u"rad"))*period(a)) +""" + delay(a) + +Return the delay of the field `a`. +""" delay(a::DelayedField) = a.t₀ delay(a::AbstractField) = 0 @@ -637,7 +659,7 @@ end # **** Kaiser @doc raw""" - Kaser(α) + Kaiser(α) ```math W(x) = \frac{I_0\left[\pi\alpha\sqrt{1 - (2x)^2}\right]}{I_0(\pi\alpha)}, @@ -663,7 +685,7 @@ function window_derivative(w::Kaiser, x) abs(2x) > 1 && return zero(x) a = π*w.α f = √(1 - (2x)^2) - pf = 2a*x/(besseli(0,a)*f) + pf = 4a*x/(besseli(0,a)*f) -pf*besseli(1, a*f) end diff --git a/src/bspline_field.jl b/src/bspline_field.jl new file mode 100644 index 0000000..1857718 --- /dev/null +++ b/src/bspline_field.jl @@ -0,0 +1,51 @@ +""" + BSplineField(B, C) + +Represents an electric field using an expansion over a B-spline basis +`B`, with `C` being the array (vector for 1d fields, 3-column matrix +for 3d fields) expansion coefficients for the vector potential. +""" +struct BSplineField{Bt,Ct<:AbstractArray} <: AbstractField + B::Bt + C::Ct +end + +function BSplineField(B::BSplineOrView, t::AbstractVector, A::AbstractVecOrMat) + BSplineField(B, interpolate(B, t, A)) +end + +Base.show(io::IO, f::BSplineField) = + write(io, "B-spline field") + +function Base.show(io::IO, ::MIME"text/plain", f::BSplineField) + show(io, f) + println(io, " expanded over") + write(io, " ") + show(io, f.B) +end + +dimensions(::BSplineField{<:Any,<:AbstractVector}) = 1 +dimensions(::BSplineField{<:Any,<:AbstractMatrix}) = 3 + +polarization(::BSplineField{<:Any,<:AbstractVector}) = LinearPolarization() +polarization(::BSplineField{<:Any,<:AbstractMatrix}) = ArbitraryPolarization() + +linear_combination(V::AbstractVector, C::AbstractVector) = + dot(V, C) + +linear_combination(V::AbstractVector, C::AbstractMatrix) = + SVector(dot(V, view(C, :, 1)), dot(V, view(C, :, 2)), dot(V, view(C, :, 3))) + +vector_potential(f::BSplineField, t::Number) = + linear_combination(f.B[t,:], f.C) + +vector_potential(f::BSplineField, t::AbstractVector) = + f.B[t,:]*f.C + +field_amplitude(f::BSplineField, t::Number) = + -linear_combination(derivative(f.B, t, :, 1), f.C) + +field_amplitude(f::BSplineField, t::AbstractVector) = + -derivative(f.B, t, :, 1)*f.C + +export BSplineField diff --git a/src/bsplines.jl b/src/bsplines.jl new file mode 100644 index 0000000..bb0eb30 --- /dev/null +++ b/src/bsplines.jl @@ -0,0 +1,236 @@ +# * Unit Vectors +""" + UnitVector{T}(N, k) + +Helper vector type of length `N` where the `k`th element is `one(T)` +and all the others `zero(T)`. +""" +struct UnitVector{T} <: AbstractVector{T} + N::Int + k::Int +end + +Base.size(e::UnitVector) = (e.N,) +Base.getindex(e::UnitVector{T}, i::Int) where T = i == e.k ? one(T) : zero(T) +Base.show(io::IO, e::UnitVector{T}) where T = write(io, "ê{$T}($(e.k)|$(e.N))") + +# * B-Splines + +# This is a cannibalized version of the code found in CompactBases.jl, +# to avoid dragging in all dependencies. + +""" + BSpline(t, x, w, B, S) + +Basis structure for the B-splines generated by the knot set `t`. +""" +struct BSpline{KnotSet<:AbstractKnotSet} + t::KnotSet +end + +# * Properties + +function firstindex(B::BSpline, i) + @assert i == 2 + 1 +end +function lastindex(B::BSpline, i) + @assert i == 2 + numfunctions(B.t) +end +function size(B::BSpline, i) + @assert i == 2 + numfunctions(B.t) +end + +knotset(B::BSpline) = B.t +order(B::BSpline) = order(knotset(B)) + +function show(io::IO, B::BSpline) + write(io, "BSpline basis with $(B.t)") +end + +# # * Basis functions + +""" + deBoor(t, c, x[, i[, m=0]]) + +Evaluate the spline given by the knot set `t` and the set of control +points `c` at `x` using de Boor's algorithm. `i` is the index of the +knot interval containing `x`. If `m≠0`, calculate the `m`th derivative +at `x` instead. +""" +function deBoor(t::AbstractKnotSet, c::AbstractVector, x::Number, + i=find_interval(t, x), m=0) + T = promote_type(eltype(t), eltype(c), typeof(x)) + isnothing(i) && return zero(T) + k = order(t) + nc = length(c) + k == 1 && return (nc < i || m > 0) ? zero(T) : c[i] + + α = [r > 0 && r ≤ nc ? T(c[r]) : zero(T) + for r ∈ i-k+1:i] + nt = length(t) + nf = numfunctions(t) + for j = 1:k-1 + for r = i:-1:max(i-k+j,1) + jjj = r+k-j + (jjj > nt || jjj < 1) && continue + r′ = r - i + k + r ≠ 1 && r′ == 1 && continue + + a = t[r+k-j] + b = t[r] + + α[r′] = if j ≤ m + (k-j)*if r == 1 + -α[r′] + elseif r == nf + j + α[r′-1] + else + α[r′-1] - α[r′] + end + else + if r == 1 + (b-x)*α[r′] + elseif r == nf + j + (x-a)*α[r′-1] + else + ((x-a)*α[r′-1] + (b-x)*α[r′]) + end + end + α[r′] /= (b-a) + end + end + + α[end] +end + +function getindex(B::BSpline, x::Real, j::Integer) + T = promote_type(eltype(B.t), typeof(x)) + deBoor(B.t, UnitVector{T}(size(B,2), j), + x, find_interval(B.t, x)) +end + +function getindex(B::BSpline, x::Real, sel::AbstractVector) + T = promote_type(eltype(B.t), typeof(x)) + i = find_interval(B.t, x) + χ = spzeros(T, length(sel)) + o = sel[1] - 1 + for j in sel + χ[j-o] = deBoor(B.t, UnitVector{T}(size(B,2), j), x, i) + end + χ +end + +function basis_function!(χ, B::BSpline, x::AbstractVector, j) + T = promote_type(eltype(B.t), eltype(x)) + eⱼ = UnitVector{T}(size(B,2), j) + for (is,k) ∈ within_support(x, B.t, j) + for i in is + χ[i] = deBoor(B.t, eⱼ, x[i], k) + end + end +end + +function getindex(B::BSpline, x::AbstractVector, sel::AbstractVector) + T = promote_type(eltype(B.t), eltype(x)) + χ = spzeros(T, length(x), length(sel)) + o = sel[1] - 1 + for j in sel + basis_function!(view(χ, :, j-o), B, x, j) + end + χ +end + +getindex(B::BSpline, x, ::Colon) = + getindex(B, x, 1:size(B,2)) + +function derivative(B::BSpline, x::Real, j::Integer, m) + T = promote_type(eltype(B.t), typeof(x)) + deBoor(B.t, UnitVector{T}(size(B,2), j), + x, find_interval(B.t, x), m) +end + +function derivative(B::BSpline, x::Real, sel::AbstractVector, m) + T = promote_type(eltype(B.t), typeof(x)) + i = find_interval(B.t, x) + χ = spzeros(T, length(sel)) + o = sel[1] - 1 + for j in sel + χ[j-o] = deBoor(B.t, UnitVector{T}(size(B,2), j), x, i, m) + end + χ +end + +function derivative_basis_function!(χ, B::BSpline, x::AbstractVector, j, m) + T = promote_type(eltype(B.t), eltype(x)) + eⱼ = UnitVector{T}(size(B,2), j) + for (is,k) ∈ within_support(x, B.t, j) + for i in is + χ[i] = deBoor(B.t, eⱼ, x[i], k, m) + end + end +end + +function derivative(B::BSpline, x::AbstractVector, sel::AbstractVector, m) + T = promote_type(eltype(B.t), eltype(x)) + χ = spzeros(T, length(x), length(sel)) + o = sel[1] - 1 + for j in sel + derivative_basis_function!(view(χ, :, j-o), B, x, j, m) + end + χ +end + +derivative(B::BSpline, x, ::Colon, m) = + derivative(B, x, 1:size(B,2), m) + +struct BSplineView{Bt,Sel} + B::Bt + sel::Sel +end + +knotset(B::BSplineView) = knotset(B.B) + +function size(B::BSplineView, i) + @assert i == 2 + length(B.sel) +end + +function Base.show(io::IO, B::BSplineView) + write(io, "Restriction to basis functions ") + show(io, B.sel) + write(io, " of ") + show(io, B.B) +end + +function getindex(B::BSpline, ::Colon, sel::AbstractVector) + all(∈(1:size(B,2)), sel) || + throw(ArgumentError("Invalid selection")) + BSplineView(B, sel) +end + +getindex(B::BSplineView, x, j) = + getindex(B.B, x, B.sel[j]) + +derivative(B::BSplineView, x, j, m) = + derivative(B.B, x, B.sel[j], m) + +const BSplineOrView = Union{BSpline,BSplineView} + +# * Function interpolation + +interpolate(B::BSplineOrView, x::AbstractVector, y::AbstractVector) = + B[x,:] \ y + +function interpolate(B::BSplineOrView, x::AbstractVector, y::AbstractMatrix) + T = promote_type(eltype(knotset(B)), eltype(x), eltype(y)) + n = size(y,2) + C = zeros(T, size(B,2), n) + V = B[x,:] + for j = 1:n + C[:,j] = V \ y[:,j] + end + C +end diff --git a/src/dispersed_fields.jl b/src/dispersed_fields.jl index 41bf736..3c98dfc 100644 --- a/src/dispersed_fields.jl +++ b/src/dispersed_fields.jl @@ -1,40 +1,118 @@ -# * Dispersed fields -# We calculate the effect of dispersion described by the Sellmeier -# equations in the frequency domain, to which we transform via an -# FFT. For this to be possible, we only allow the evaluation of the -# field for a specified sampling frequency, i.e. we don't provide an -# implementation for evaluating a =DispersedField= at an arbitrary -# time point. - -struct DispersedField <: AbstractField - a::AbstractField - m::Medium - d::Unitful.Length +""" + DispersedField(f, de, B, s) + +Represents the field `f` dispersed through the +[`DispersiveElement`](@ref) `de`. Since the dispersion is most +efficiently calculated in the frequency domain, this is done via +[`rfft`](@ref)/[`irfft`](@ref), and the dispersed field is +interpolated as the [`BSplineField`](@ref) `B`. The temporal [`span`](@ref) +of the dispersed field is given by `s`. +""" +struct DispersedField{Ft,Dt,Bt,S} <: AbstractField + f::Ft + de::Dt + FB::Bt + s::S end -function show(io::IO, f::DispersedField) - show(io, f.a) - write(io, "\n – dispersed through $(f.d) of $(f.m)") +""" + DispersedField(f, de; spline_order=7, Bfs=20/period(F)) + +Convenience constructor for [`DispersedField`](@ref) that +automatically finds the appropriate time span for `f` after it has +been dispersed through `de`, and then fits a [`BSpline`](@ref) to the +resultant vector potential, such that [`vector_potential`](@ref) and +[`field_amplitude`](@ref) can be evaluated at arbitrary times. The +number of knots can be influenced by the \"sampling frequency\" `Bfs`. +""" +function DispersedField(f::AbstractField, de; + spline_order=3, Bfs=20/period(parent(f)), + cutoff=1e5*√(eps()), + verbosity=0, kwargs...) + s = find_time_span(f, de; verbosity=verbosity-2, kwargs...) + t = timeaxis(f, s=s) + s₀ = span(f) + + ω = rfftω(t) + H = frequency_response(de, ω) + + As = rfft_vector_potential(f, t) + Arec = irfft(H.*As, t) + + # Truncate time span such that |A| > cutoff*max|A|, but still + # cover the whole time span of the undispersed field. + if cutoff > 0 + Amag = vec(sum(abs, Arec, dims=2)) + Amax = maximum(Amag) + abs_cutoff = cutoff*Amax + + sel = findall(∈(s₀), t) + a = min(findfirst(>(abs_cutoff), Amag), first(sel)) + b = max(findlast(>(abs_cutoff), Amag), last(sel)) + + s = t[a] .. t[b] + t = t[a:b] + Arec = selectdim(Arec, 1, a:b) + + verbosity > 1 && @info "Truncated to time interval $(s)" a b cutoff abs_cutoff + end + + num_knots = ceil(Int, austrip(Bfs)*(s.right-s.left)) + B = BSpline(LinearKnotSet(spline_order, s.left, s.right, num_knots)) + verbosity > 1 && @info "Generated B-spline" num_knots B + + FB = BSplineField(B, t, Arec) + + DispersedField(f, de, FB, s) end -disperse(a::AbstractField, m::Medium, d::Unitful.Length) = - DispersedField(a, m, d) +@doc raw""" + DispersedField(f::DispersedField, de) -(f::DispersedField)(t::Union{Unitful.Time,AbstractVector{<:Unitful.Time}}) = - error("Dispersed fields can only be evaluated by specifying a sampling frequency") +Stacking another [`DispersiveElement`](@ref) `de` onto an already +dispersed field `f` is accomlished by multiplying the transfer +functions together, with the earliest applied as the right-most factor: +```math +H(\omega) = H_n(\omega) ... H_3(\omega)H_2(\omega)H_1(\omega). +``` +""" +DispersedField(f::DispersedField, de; kwargs...) = + DispersedField(parent(f), de*f.de; kwargs...) -function (f::DispersedField)(fs::Unitful.Frequency = default_sampling_frequency(f.a)) - F̂ = fft(f.a, fs) - ifft(F̂.*dispersion(f.m, f.d, fftfreq(f.a, fs), frequency(f.a)))*(u"V"*u"m") +function show(io::IO, f::DispersedField) + write(io, "DispersedField(") + show(io, f.f) + write(io, ", ") + show(io, f.de) + write(io, ")") +end +function show(io::IO, mime::MIME"text/plain", f::DispersedField) + println(io, "DispersedField:") + show(io, mime, f.f) + write(io, "\n – dispersed through ") + show(io, mime, f.de) end -max_frequency(f::DispersedField) = max_frequency(f.a) +Base.parent(f::DispersedField) = f.f + +for fun in [:vector_potential, :field_amplitude] + @eval $fun(f::DispersedField, t::Number) = + $fun(f.FB, t) + @eval $fun(f::DispersedField, t::AbstractVector) = + $fun(f.FB, t) +end + +# Some of these forwards are iffy, i.e. after e.g. a chirp, the +# carrier is most certainly not a FixedCarrier anymore, and the +# duration is not the Fourier-limited one. +for fun in [:params, :carrier, :envelope, :polarization, + :wavelength, :period, :frequency, :max_frequency, + :wavenumber, :fundamental, :photon_energy, + :intensity, :amplitude, :duration, :continuity, + :dimensions, :rotation_matrix] + @eval $fun(f::DispersedField) = $fun(parent(f)) +end -# This is, strictly speaking, not true, since dispersing a pulse will -# in general stretch it in the time domain. It is up to the user to -# ensure that the time window is large enough to contain even the -# stretched pulse, by setting the σmax (or friends) to a large enough -# value. -span(f::DispersedField) = span(f.a) +span(f::DispersedField) = f.s -export disperse +export DispersedField diff --git a/src/dispersive_elements.jl b/src/dispersive_elements.jl new file mode 100644 index 0000000..894c018 --- /dev/null +++ b/src/dispersive_elements.jl @@ -0,0 +1,288 @@ +# * (An)isotropy + +abstract type Tropy end +struct Isotropic <: Tropy end +struct Anisotropic <: Tropy end + +Base.:(*)(::Isotropic, ::Isotropic) = Isotropic() +Base.:(*)(::Tropy, ::Tropy) = Anisotropic() + +# * Dispersive elements + +""" + DispersiveElement + +Base type for all dispersive elements +""" +abstract type DispersiveElement end +tropy(::DispersiveElement) = Anisotropic() +Base.Broadcast.broadcastable(de::DispersiveElement) = Ref(de) + +frequency_response(::Isotropic, de::DispersiveElement, ω::AbstractVector) = + frequency_response.(de, ω) + +function frequency_response(::Anisotropic, de::DispersiveElement, ω::AbstractVector) + T = eltype(frequency_response(de, first(ω))) + H = zeros(T, length(ω), 3) + for (i,ω) in enumerate(ω) + H[i,:] .= frequency_response(de, ω) + end + H +end + +frequency_response(de::DispersiveElement, ω::AbstractVector) = + frequency_response(tropy(de), de, ω) + +# ** Phase shift + +@doc raw""" + PhaseShift(ϕ) + +Represents a phase shift according to +```math +H(\omega) = +\exp(-\im\phi). +``` + +# Example + +```julia +julia> PhaseShift(6) +PhaseShift(ϕ = 6.0000 rad) +``` +""" +struct PhaseShift{T} <: DispersiveElement + ϕ::T +end + +Base.show(io::IO, e::PhaseShift) = + printfmt(io, "PhaseShift(ϕ = {1:0.4g} rad)", e.ϕ) + +tropy(::PhaseShift) = Isotropic() + +frequency_response(e::PhaseShift, ::Number) = + exp(-im*e.ϕ) + +""" + phase_shift(f, ϕ) + +Returns the field resulting from applying a [`PhaseShift`](@ref) to +the field `f`. +""" +phase_shift(f::AbstractField, ϕ; kwargs...) = + DispersedField(f, PhaseShift(ϕ); kwargs...) + +# ** Chirp + +@doc raw""" + Chirp(b, ω₀) + +Represents chirp according to +```math +H(\omega) = +\exp[-\im b(ω-ω_0)^2]. +``` + +# Example + +```julia +julia> Chirp(austrip(5u"fs^2"), austrip(1.5u"eV")) +Chirp(b = 8545.5457 = 5.0000 fs², ω₀ = 0.0551 = 1.5000 eV) +``` +""" +struct Chirp{T,U} <: DispersiveElement + b::T + ω₀::U +end + +function Base.show(io::IO, c::Chirp) + bu = u"fs^2" + ω₀u = u"eV" + printfmt(io, "Chirp(b = {1:0.4g} = {2:0.4g} {3:s}, ω₀ = {4:0.4g} = {5:0.4g} {6:s})", + c.b, ustrip(auconvert(bu, c.b)), bu, + c.ω₀, ustrip(auconvert(ω₀u, c.ω₀)), ω₀u) +end + +tropy(::Chirp) = Isotropic() + +frequency_response(e::Chirp, ω::Number) = + exp(-im*e.b*(ω-e.ω₀)^2) + +""" + chirp(f, b, ω₀=photon_energy(f)) + +Returns the field resulting from applying a [`Chirp`](@ref) to the +field `f`. + +# Example + +```julia +julia> @field(F) do + I₀ = 1.0 + T = 2.0u"fs" + σ = 3.0 + Tmax = 3.0 + end +Linearly polarized field with + - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² => + - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹ + - A₀ = 13.1594 au + – a Fixed carrier @ λ = 599.5849 nm (T = 2.0000 fs, ω = 0.0760 Ha = 2.0678 eV, f = 500.0000 THz) + – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±82.68σ) + – and a bandwidth of 0.3925 Ha = 10.6797 eV ⟺ 2.5823 PHz ⟺ 58518.2144 Bohr = 3.0967 μm + – Uₚ = 43.2922 Ha = 1.1780 keV => α = 173.1690 Bohr = 9.1637 nm + +julia> chirp(F, austrip(1u"fs^2")) +DispersedField: +Linearly polarized field with + - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² => + - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹ + - A₀ = 13.1594 au + – a Fixed carrier @ λ = 599.5849 nm (T = 2.0000 fs, ω = 0.0760 Ha = 2.0678 eV, f = 500.0000 THz) + – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±82.68σ) + – and a bandwidth of 0.3925 Ha = 10.6797 eV ⟺ 2.5823 PHz ⟺ 58518.2144 Bohr = 3.0967 μm + – Uₚ = 43.2922 Ha = 1.1780 keV => α = 173.1690 Bohr = 9.1637 nm + – dispersed through Chirp(b = 1709.1091 = 1.0000 fs², ω₀ = 0.0760 = 2.0678 eV) +``` +""" +chirp(f::AbstractField, b, ω₀=photon_energy(f); kwargs...) = + DispersedField(f, Chirp(b, ω₀); kwargs...) + +# ** Cascade of dispersive elements + +""" + CascadedDispersiveElement(elements) + +Represents the combination of multiple [`DispersiveElement`](@ref)s. + +# Example + +```julia +julia> PhaseShift(6)*Chirp(austrip(5u"fs^2"), austrip(1.5u"eV")) +[PhaseShift(ϕ = 6.0000 rad) ∗ Chirp(b = 8545.5457 = 5.0000 fs², ω₀ = 0.0551 = 1.5000 eV)] +``` +""" +struct CascadedDispersiveElement{Elements} <: DispersiveElement + elements::Elements +end + +function Base.show(io::IO, ce::CascadedDispersiveElement) + if isempty(ce.elements) + write(io, "CascadedDispersiveElement(())") + return + end + write(io, "[") + n = length(ce.elements) + for (i,e) in enumerate(ce.elements) + show(io, e) + i == n || write(io, " ∗ ") + end + write(io, "]") +end + +tropy(ce::CascadedDispersiveElement) = + foldr((a,t) -> tropy(a) * t, ce.elements, init=Isotropic()) + +function frequency_response(ce::CascadedDispersiveElement, ω::Number) + v = 1 + # We "have to" loop over the list in reverse, since the + # multiplication rules below ensure the right-most factor in a + # multiplication ends up as the last element of the list. Since + # all frequency responses are scalars, they commute however. + for e in reverse(ce.elements) + v *= frequency_response(e, ω) + end + v +end + +Base.:(*)(a::DispersiveElement, b::DispersiveElement) = + CascadedDispersiveElement((a,b)) + +Base.:(*)(a::CascadedDispersiveElement, b::DispersiveElement) = + CascadedDispersiveElement((a.elements...,b)) + +Base.:(*)(a::DispersiveElement, b::CascadedDispersiveElement) = + CascadedDispersiveElement((a,b.elements...)) + +Base.:(*)(a::CascadedDispersiveElement, b::CascadedDispersiveElement) = + CascadedDispersiveElement((a.elements...,b.elements...)) + +# * Find time span + +function join_ranges(r1::StepRangeLen, r2::StepRangeLen) + r1.ref == r2.ref || throw(ArgumentError("Ranges do not have the same origin")) + step(r1) == step(r2) || throw(ArgumentError("Ranges have different steps")) + + n1 = length(r1) + n2 = length(r2) + + StepRangeLen(r2.ref, r2.step, n1+n2, n1) +end + +""" +find_overlapping_range(r1::StepRangeLen, r2::StepRangeLen) + +Find the index range of `r2` which overlap with all of `r1`. +""" +function find_overlapping_range(r1::StepRangeLen, r2::StepRangeLen) + r1.ref == r2.ref || throw(ArgumentError("Ranges do not have the same origin")) + step(r1) == step(r2) || throw(ArgumentError("Ranges have different steps")) + + n1 = length(r1) + n2 = length(r2) + + sel = (1:n1) .+ (r2.offset-r1.offset) + @assert r1 == r2[sel] + sel +end + +""" + find_time_span(f, de[, fs]; max_iter=10, ξ=2.0, tol) + +Find a suitable time span ``[a,b]`` such that when `f` is dispersed +through the [`DispersiveElement`](@ref) `de`, the compact support of +the resulting field is contained in the time interval. This is done by +successively multiplying the time span on which `F` is evaluated `ξ` +before the `RFFT`, until the `IRFFT` has converged. +""" +function find_time_span(f, de, args...; max_iter=7, ξ = 2.0, tol=5e-4, verbosity=0) + verbosity > 0 && @info "Finding large enough time span to encompass dispersed field" f de max_iter ξ tol + a,b = endpoints(span(f)) + a == b && throw(ArgumentError("Cannot successively double an infinitesimal interval $(span(f))")) + δt = step(timeaxis(f, args...)) + + Δ = (b-a)/2 + t₀ = (a+b)/2 + + tf = l -> join_ranges(reverse(range(0, stop=-l, step=-δt)), range(0, stop=l, step=δt)) .+ t₀ + + t = nothing + F = nothing + + R = Inf + for i = 0:max_iter + t′ = tf(Δ*ξ^i) + if verbosity > 1 && i > 0 + println(repeat("-", 100)) + println(i) + @show t′ + end + H′ = frequency_response(de, rfftω(t′)) + F′ = irfft(H′.*rfft(f, t′), t′) + + if i > 0 + sel = find_overlapping_range(t, t′) + R = norm(F - selectdim(F′, 1, sel)) + verbosity > 1 && @show R + R < tol && break + end + + t = t′ + F = F′ + end + R > tol && @warn "Could not find large enough time span in $(max_iter) iterations" + t[1] .. t[end] +end + +# * Exports +export PhaseShift, Chirp, chirp diff --git a/src/dispersive_media.jl b/src/dispersive_media.jl new file mode 100644 index 0000000..b95dc90 --- /dev/null +++ b/src/dispersive_media.jl @@ -0,0 +1,125 @@ +""" + Medium + +Base type for all dispersive media +""" +abstract type Medium <: DispersiveElement end + +""" + IsotropicMedium(material, d[, ∂k∂ω₀ = 0]) + +Describes the dispersion through an isotropic medium of thickness `d` +(such as an atomic gas), where the refractive index, given by +`material`, is the same in all directions. + +The optional `∂k∂ω₀` can be used to subtract the linear slope from the +dispersion relation, i.e. remove the time shift, such that a pulse +with central angular frequency `ω₀` stays centred in the frame of +reference. +""" +struct IsotropicMedium{Material,D<:Length,U} <: Medium + material::Material + d::D + ∂k∂ω₀::U +end + +""" + IsotropicMedium(material, d; ω₀) + +Convenience constructor for [`IsotropicMedium`](@ref); if a central +angular frequency `ω₀` is provided, the linear slope of the dispersion +of the `material` is computed at `ω₀`. This slope is subsequently +subtracted from the dispersion. +""" +IsotropicMedium(material, d; ω₀=nothing) = + IsotropicMedium(material, d, isnothing(ω₀) ? 0 : dispersion_slope(material, ω₀)) + +tropy(::IsotropicMedium) = Isotropic() + +function Base.show(io::IO, m::IsotropicMedium) + dau = austrip(m.d) + du = u"μm" + printfmt(io, "IsotropicMedium({1:0.2g} Bohr = {2:0.2g} {3:s} of ", + dau, ustrip(auconvert(du, dau)), du) + show(io, m.material) + printfmt(io, ", ∂k∂ω₀ = {1:0.4g})", m.∂k∂ω₀) +end + +function frequency_response(m::IsotropicMedium, ω::Number) + f = ω/2π + n = real(m.material(auconvert(u"Hz", f))) + + # Wavevector in vacuum + k₀ = ω/austrip(1u"c") + # Wavevector in the medium, subtracting the linear component to + # centre the pulse in the frame of reference. + k = n*k₀ - ω*m.∂k∂ω₀ + + d = austrip(m.d) + + exp(-im*k*d) +end + +""" + Crystal(material, d, R[, ∂k∂ω₀ = 0]) + +Describes the dispersion through a general crystal of thickness `d` +and orientation given by the rotation `R`, with different refractive +indices, given by the components of `material`, along the different +crystal axes. + +In the special case where `material` has two components, we have a +_uniaxial_ crystal, which is birefringent, and `material[1]` is +referred to as _ordinary_ the refractive index, and `material[2]` as +the _extraordinary_ refractive index. + +The optional `∂k∂ω₀` can be used to subtract the linear slope from the +dispersion relation, i.e. remove the time shift, such that a pulse +with central angular frequency `ω₀` stays centred in the frame of +reference. +""" +struct Crystal{Material,D<:Length,Rotation<:AbstractMatrix,U} <: Medium + material::Material + d::D + R::Rotation + ∂k∂ω₀::U +end + +""" + Crystal(material, d, R=I; ω₀) + +Convenience constructor for [`Crystal`](@ref); if a central +angular frequency `ω₀` is provided, the linear slope of the dispersion +of the `material` is computed at `ω₀`. This slope is subsequently +subtracted from the dispersion. +""" +Crystal(material, d, R=I; ω₀=nothing) = + Crystal(material, d, rotation_matrix(R), + isnothing(ω₀) ? 0 : dispersion_slope(material, ω₀)) + +function Base.show(io::IO, m::Crystal) + dau = austrip(m.d) + du = u"μm" + printfmt(io, "Crystal({1:0.2g} Bohr = {2:0.2g} {3:s} of ", + dau, ustrip(auconvert(du, dau)), du) + show(io, m.material) + printfmt(io, ", R = {1:s}, ∂k∂ω₀ = {2:0.4g})", m.R, m.∂k∂ω₀) +end + +function frequency_response(m::Crystal, ω::Number) + f = ω/2π + n⁻² = inv.(n²_3d(m.material, auconvert(u"Hz", f))) + n = .√(maybe_complex.(inv.(m.R*n⁻²))) + + # Wavevector in vacuum + k₀ = ω/austrip(1u"c") + # Wavevector in the medium, subtracting the linear component to + # centre the pulse in the frame of reference. + k = n*k₀ .- ω*m.∂k∂ω₀ + + d = austrip(m.d) + + exp.(-im*k*d) +end + +export IsotropicMedium, Crystal diff --git a/src/envelopes.jl b/src/envelopes.jl index 55a1ab0..5e9f1ab 100644 --- a/src/envelopes.jl +++ b/src/envelopes.jl @@ -328,29 +328,42 @@ Required parameters: Beware that this envelope can introduce artifacts at the ends of the pulse, such that the electric field is non-vanishing, depending of e.g. the phase of the carrier. + +The shape of the ramp can be influenced by chaining ``g[f(t)]``, where +``g(x)`` should map smoothly to ``[0,1]`` for ``x\in[0,1]``. Currently +implemented ramps are +- `ramp_kind = :linear` ``\implies g(x) = x`` (default) +- `ramp_kind = :sin² | ramp_kind = :sin2` ``\implies g(x) = \sin^2(\pi x/2)`` """ -struct TrapezoidalEnvelope{T} <: AbstractEnvelope +struct TrapezoidalEnvelope{T,G} <: AbstractEnvelope ramp_up::T flat::T ramp_down::T period::T + g::G + ramp_kind::Symbol - function TrapezoidalEnvelope(ramp_up::T, flat::T, ramp_down::T, period::T) where T + function TrapezoidalEnvelope(ramp_up::T, flat::T, ramp_down::T, period::T, + g::G, ramp_kind) where {T,G<:Function} ramp_up >= 0 || error("Negative up-ramp not supported") flat >= 0 || error("Negative flat region not supported") ramp_down >= 0 || error("Negative down-ramp not supported") ramp_up + flat + ramp_down > 0 || error("Pulse length must be non-zero") - new{T}(ramp_up, flat, ramp_down, period) + new{T,G}(ramp_up, flat, ramp_down, period, g, ramp_kind) end end + +TrapezoidalEnvelope(ramp_up, flat, ramp_down, period) = + TrapezoidalEnvelope(ramp_up, flat, ramp_down, period, identity, :linear) + envelope_types[:trapezoidal] = TrapezoidalEnvelope # Common alias envelope_types[:tophat] = TrapezoidalEnvelope function (env::TrapezoidalEnvelope{T})(t) where T t = real(t)/env.period - if t < 0 + f = if t < 0 zero(T) elseif t < env.ramp_up t/env.ramp_up @@ -361,11 +374,12 @@ function (env::TrapezoidalEnvelope{T})(t) where T else zero(T) end + env.g(f) end show(io::IO, env::TrapezoidalEnvelope) = - printfmt(io, "/{1:d}‾{2:d}‾{3:d}\\ cycles trapezoidal envelope", - env.ramp_up, env.flat, env.ramp_down) + printfmt(io, "/{1:d}‾{2:d}‾{3:d}\\ cycles trapezoidal envelope, with {4:s} ramps", + env.ramp_up, env.flat, env.ramp_down, env.ramp_kind) function TrapezoidalEnvelope(field_params::Dict{Symbol,Any}, args...) test_field_parameters(field_params, [:T]) # Period time required to relate ramps/flat to cycles @@ -378,10 +392,22 @@ function TrapezoidalEnvelope(field_params::Dict{Symbol,Any}, args...) ramp_up = ramp ramp_down = ramp end + if !ramp_kind + ramp_kind = :linear + end + end + + @unpack ramp_up, flat, ramp_down, T, ramp_kind = field_params + + g = if ramp_kind == :linear + identity + elseif ramp_kind == :sin² || ramp_kind == :sin2 + x -> sinpi(x/2)^2 + else + throw(ArgumentError("Unknown ramp kind $(ramp_kind)")) end - @unpack ramp_up, flat, ramp_down, T = field_params - TrapezoidalEnvelope(ramp_up, flat, ramp_down, austrip(T)) + TrapezoidalEnvelope(ramp_up, flat, ramp_down, austrip(T), g, ramp_kind) end continuity(::TrapezoidalEnvelope) = 0 diff --git a/src/field_types.jl b/src/field_types.jl index d1cef26..a3ff8ba 100644 --- a/src/field_types.jl +++ b/src/field_types.jl @@ -481,7 +481,7 @@ Linearly polarized field with - E₀ = 2.2361e-01 au = 114.9832 GV m^-1 - A₀ = 0.2236 au – a Fixed carrier @ λ = 45.5634 nm (T = 151.9830 as, ω = 1.0000 Ha = 27.2114 eV, f = 6.5797 PHz) - – and a /1‾3‾1\\ cycles trapezoidal envelope + – and a /1‾3‾1\\ cycles trapezoidal envelope, with linear ramps – and a bandwidth of Inf Ha = Inf eV ⟺ Inf Hz ⟺ Inf Bohr = Inf m – Uₚ = 0.0125 Ha = 340.1423 meV => α = 0.2236 Bohr = 11.8328 pm ``` diff --git a/src/knot_sets.jl b/src/knot_sets.jl new file mode 100644 index 0000000..415d4dc --- /dev/null +++ b/src/knot_sets.jl @@ -0,0 +1,375 @@ +import Base: first, last, length, size, + getindex, firstindex, lastindex, eachindex, iterate, + eltype, similar, + show + +# * AbstractKnotSet +""" + AbstractKnotSet{k,ml,mr,T} + +Abstract base for B-spline knot sets. `T` is the `eltype` of the knot +set, `k` is the order of the piecewise polynomials (`order = degree + +1`) and `ml` and `mr` are the knot multiplicities of the left and +right endpoints, respectively. """ +abstract type AbstractKnotSet{k,ml,mr,T} <: AbstractVector{T} end + +""" + assert_multiplicities(k,ml,mr,t) + +Assert that the multiplicities at the endpoints, `ml` and `mr`, +respectively, are consistent with the order `k`. Also check that the +amount of knots in the knot set `t` are enough to support the +requested order `k`. +""" +function assert_multiplicities(k,ml,mr,t) + k > 0 || throw(ArgumentError("Polynomial order must be a positive integer, got k = $k")) + ml ∈ 1..k || throw(ArgumentError("Left multiplicity has to be in range 1..$(k), got $ml")) + mr ∈ 1..k || throw(ArgumentError("Right multiplicity has to be in range 1..$(k), got $mr")) + n = length(t) + nreq = max(2, k + 3 - ml - mr) + n ≥ nreq || throw(ArgumentError("Polynomial order k = $(k) and left,right multiplicities $(ml),$(mr) require a knot sequence of at least $(nreq) points")) +end + +""" + order(t) + +Returns the order `k` of the knot set `t`. +""" +order(t::AbstractKnotSet{k}) where k = k + +""" + numinterval(t) + +Returns the number of intervals generated by the knot set `t`. + +# Examples + +```jldoctest +julia> ElectricFields.numintervals(ElectricFields.LinearKnotSet(3, 0, 1, 2)) +2 +``` +""" +numintervals(t::AbstractKnotSet) = length(t.t)-1 + +""" + numfunctions(t) + +Returns the number of basis functions generated by knot set `t`. + +# Examples + +```jldoctest +julia> ElectricFields.numfunctions(ElectricFields.LinearKnotSet(3, 0, 1, 2)) +4 + +julia> ElectricFields.numfunctions(ElectricFields.LinearKnotSet(5, 0, 1, 2)) +6 +``` +""" +numfunctions(t::AbstractKnotSet) = length(t) - order(t) + +""" + leftmultiplicity(t) + +Return the multiplicity of the left endpoint. +""" +leftmultiplicity(::AbstractKnotSet{k,ml,mr}) where {k,ml,mr} = ml + +""" + rightmultiplicity(t) + +Return the multiplicity of the right endpoint. +""" +rightmultiplicity(::AbstractKnotSet{k,ml,mr}) where {k,ml,mr} = mr + +""" + first(t) + +Return the first knot of `t`. +""" +first(t::AbstractKnotSet) = first(t.t) + +""" + last(t) + +Return the last knot of `t`. +""" +last(t::AbstractKnotSet) = last(t.t) + +""" + length(t) + +Return the number of knots of `t`. + +# Examples + +```jldoctest +julia> length(ElectricFields.LinearKnotSet(3, 0, 1, 3)) +8 + +julia> length(ElectricFields.LinearKnotSet(3, 0, 1, 3, 1, 1)) +4 +``` +""" +length(t::AbstractKnotSet) = length(t.t) + leftmultiplicity(t) + rightmultiplicity(t) - 2 +# size(t) is not very useful for knot sets, but necessary to fulfill +# the interface expected from AbstractVectors. +size(t::AbstractKnotSet) = (length(t),) + +""" + getindex(t, i) + +Return the `i`th knot of the knot set `t`, accounting for the +multiplicities of the endpoints. + +# Examples + +```jldoctest +julia> ElectricFields.LinearKnotSet(3, 0, 1, 3)[2] +0.0 + +julia> ElectricFields.LinearKnotSet(3, 0, 1, 3, 1, 1)[2] +0.3333333333333333 +``` + +""" +function getindex(t::AbstractKnotSet{k,ml,mr}, i::Integer) where {k,ml,mr} + i < 1 && throw(BoundsError(t, i)) + + ni = numintervals(t) + if i < ml + first(t) + elseif i < ml + ni + t.t[i-ml+1] + elseif i < ml + ni + mr + last(t) + else + throw(BoundsError(t, i)) + end +end +lastindex(t::AbstractKnotSet) = length(t) +eachindex(t::AbstractKnotSet) = 1:length(t) +# function iterate(t::AbstractKnotSet, (element,i)=(t[1],0)) +# i ≥ length(t) && return nothing +# element, (t[i+1], i+1) +# end + +eltype(t::AbstractKnotSet{k,ml,mr,T}) where {k,ml,mr,T} = T +similar(t::AbstractKnotSet{k,ml,mr,T}) where {k,ml,mr,T} = Vector{T}(undef, length(t)) + +function show_order_multiplicities(io::IO, t::AbstractKnotSet{k,ml,mr}) where {k,ml,mr} + write(io, "order k = $(k)") + k < 7 && write(io, " (", + ["constant", "linear", "parabolic", "cubic", "quartic", "quintic"][k], + ")") + ml == k && mr == k && return + write(io, " (left, right multiplicities: $(ml), $(mr))") +end + +function show(io::IO, t::AbstractKnotSet) + # https://stackoverflow.com/a/76363569 + write(io, "$(nameof(typeof(t)))($(eltype(t))) of ") + show_order_multiplicities(io, t) + write(io, " on $(first(t)..last(t)) ($(numintervals(t)) intervals)") +end + +# TODO: Implement via the `iterate` interface? +""" + nonempty_intervals(t) + +Return the indices of all intervals of the knot set `t` that are +non-empty. + +# Examples + +```jldoctest +julia> ElectricFields.nonempty_intervals(ElectricFields.ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3)) +4-element $(Vector{Int}): + 1 + 3 + 4 + 5 +``` + +""" +nonempty_intervals(t::AbstractKnotSet) = + [j for j in 1:length(t)-1 + if t[j] ≠ t[j+1]] + +""" + find_interval(t, x[, i=ml]) + +Find the interval in the knot set `t` that includes `x`, starting from +interval `i` (which by default is the first non-zero interval of the +knot set). The search complexity is linear, but by storing the result +and using it as starting point for the next call to `find_interval`, +the knot set need only be traversed once. +""" +function find_interval(t::AbstractKnotSet{k,ml,mr,T}, x, i=ml) where {T,k,ml,mr} + (x < first(t) || x > last(t) || i > length(t) || x < t[i]) && return nothing + x == last(t) && return length(t) - mr + for r ∈ i:length(t) + t[r] > x && return r-1 + end + @assert false +end + +const RightContinuous{T} = Interval{:closed,:open,T} + +""" + within_support(x, t, j) + +Return the indices of the elements of `x` that lie withing the compact +support of the `j`th basis function (enumerated `1..n`), given the +knot set `t`. For each index of `x` that is covered, the index `k` of +the interval within which `x[i]` falls is also returned. +""" +function within_support(x::AbstractVector, t::AbstractKnotSet, j::Integer) + isempty(x) && return 1:0 + k = order(t) + # The last interval includes the right endpoint as well, whereas + # the other intervals are only right-continuous. + IntervalKind(i) = i == numintervals(t) ? ClosedInterval : RightContinuous + ml = leftmultiplicity(t) + supports = [(findall(in(IntervalKind(i-ml+1)(t[i], t[i+1])), x), i) + for i = j:j+k-1] + filter(s -> !isempty(s[1]), supports) +end + +# * Specific knot sets +# ** Arbitrary knot set +""" + ArbitraryKnotSet{k,ml,mr}(t) + +An arbitrary knot set of order `k` and left and right multiplicities +of `ml` and `mr`, respectively. The knot set is specified by the +non-decreasing vector `t`; the same knot may appear multiple times, +which influences the continuity of the B-splines at that location. +""" +struct ArbitraryKnotSet{k,ml,mr,T,V<:AbstractVector{T}} <: AbstractKnotSet{k,ml,mr,T} + t::V + ArbitraryKnotSet{k,ml,mr}(t::V) where {k,ml,mr,T,V<:AbstractVector{T}} = + new{k,ml,mr,T,V}(sort(t)) +end + +""" + ArbitraryKnotset(k, t[, ml=k, mr=k]) + +Construct an order-`k` arbitrary knot set, with the locations of the +knots given by non-decreasing vector `t`. +""" +function ArbitraryKnotSet(k::Integer, t::AbstractVector, + ml::Integer=k, mr::Integer=k) + assert_multiplicities(k, ml, mr, t) + ArbitraryKnotSet{k,ml,mr}(t) +end + +# ** Linear knot set +""" + LinearKnotSet{k,ml,mr}(t) + +A knot set of order `k` and left and right multiplicities of `ml` and +`mr`, respectively, whose knots are uniformly distributed according to +the range `t`. The interior basis functions are thus Cᵏ⁻²-continuous. +""" +struct LinearKnotSet{k,ml,mr,T,R<:AbstractRange{T}} <: AbstractKnotSet{k,ml,mr,T} + t::R + LinearKnotSet{k,ml,mr}(t::R) where {k,ml,mr,T,R<:AbstractRange{T}} = + new{k,ml,mr,T,R}(t) +end + +""" + LinearKnotSet(k, a, b, N[, ml=k, mr=k]) + +Construct an order-`k` linear knot set spanning from `a` to `b`, with +`N` intervals. + +# Examples + +```jldoctest +julia> ElectricFields.LinearKnotSet(2, 0, 1, 3) +6-element ElectricFields.LinearKnotSet{2, 2, 2, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}: + 0.0 + 0.0 + 0.3333333333333333 + 0.6666666666666666 + 1.0 + 1.0 +``` +""" +function LinearKnotSet(k::Integer, a, b, N::Integer, + ml::Integer=k, mr::Integer=k) + t = range(a, stop=b, length=N+1) + assert_multiplicities(k, ml, mr, t) + LinearKnotSet{k,ml,mr}(t) +end + +# ** Exponential knot set +""" + ExpKnotSet{k,ml,mr}(exponents, base, t, include0) + +A knot set of order `k` and left and right multiplicities of `ml` and +`mr`, respectively, whose knots are exponentially distributed +according to `t = base .^ exponents`, optionally including 0 as the +left endpoint.""" +struct ExpKnotSet{k,ml,mr,T,R<:AbstractRange{T},TV<:AbstractVector{T}} <: AbstractKnotSet{k,ml,mr,T} + exponents::R + base::T + t::TV + include0::Bool + ExpKnotSet{k,ml,mr}(exponents::R, base::T, t::TV, include0) where {k,ml,mr,T,R,TV} = + new{k,ml,mr,T,R,TV}(exponents, base, t, include0) +end + +""" + ExpKnotSet(k, a, b, N[, ml=k, mr=k, base=10, include0=true]) + +Construct an order-`k` knot spanning from `base^a` to `base^b` in `N` +intervals, optionally including 0 as the left endpoint. + +# Examples + +```jldoctest +julia> ElectricFields.ExpKnotSet(2, -4, 2, 7) +10-element ElectricFields.ExpKnotSet{2, 2, 2, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}}: + 0.0 + 0.0 + 0.0001 + 0.001 + 0.010000000000000002 + 0.1 + 1.0 + 10.0 + 100.0 + 100.0 +``` +""" +function ExpKnotSet(k::Integer, a::T, b::T, N::Integer, + ml::Integer=k, mr::Integer=k; + base::T=T(10), include0::Bool=true) where T + exponents = range(a, stop=b, length=include0 ? N : N+1) + t = base .^ exponents + assert_multiplicities(k, ml, mr, t) + ExpKnotSet{k,ml,mr}(exponents, eltype(t)(base), include0 ? vcat(0,t) : t, include0) +end + +function show(io::IO, t::ExpKnotSet) + write(io, "$(nameof(typeof(t)))($(eltype(t))) of ") + show_order_multiplicities(io, t) + write(io, " on ") + t.include0 && write(io, "0,") + write(io, "$(t.base^first(t.exponents))..$(t.base^last(t.exponents)) ($(numintervals(t)) intervals)") +end + +# function arcsin_knot_set(k::Integer, a::Integer, b::Integer, N::Integer) +# N2 = N/2 +# arcsin.(range(-N2, stop=N2, lengthN+1)/N2)*(b-a)/π+(a+b)/2 +# end + +# arcsin_half_knot_set(a, b, N) = arcsin(linspace(0,N,N+1)/N)*(b-a)*2.0/π+a + +# function exp_linear_knot_set(a, b, N) +# ap = a != 0 ? a : 1e-1 +# @assert t[2][1]+t[2][2] == N +# [logspace(log10(ap),log10(t[1]), t[2][1]); linspace(t[1],b,t[2][2]+1)[2:end]] +# end diff --git a/src/rotations.jl b/src/rotations.jl index 5bb383b..7a72e0b 100644 --- a/src/rotations.jl +++ b/src/rotations.jl @@ -37,4 +37,7 @@ function rotation_axis(R::AbstractMatrix) normalize(real(ee.vectors[:,i])) end +rotation_matrix(R::AbstractMatrix) = compute_rotation(R) +rotation_matrix(I::UniformScaling) = SMatrix{3,3}(compute_rotation(I)) + export rotation_angle, rotation_axis, rotation_matrix diff --git a/src/sellmeier.jl b/src/sellmeier.jl index fb85685..fff1e39 100644 --- a/src/sellmeier.jl +++ b/src/sellmeier.jl @@ -3,93 +3,249 @@ using Unitful import Unitful: Length, Area, Frequency import Base: show -using Printf +@derived_dimension InverseArea Unitful.𝐋^(-2) @doc raw""" - Medium(B, C) + Sellmeier(A, B, p, C, D, q) The [Sellmeier equations](https://en.wikipedia.org/wiki/Sellmeier_equation) are used -to describe dispersion in glass using a series of resonances: +to describe dispersion in glass using a series of resonances, here +generalized to support variations found in the literature: ```math n^2(\lambda) = -1 + \sum_i \frac {B_i\lambda^2}{\lambda^2-C_i}. +1 + A + \sum_i \frac {B_i\lambda^{p_i}}{\lambda^2-C_i} + \sum_j D_j \lambda^{q_j}. ``` """ -struct Medium{T<:Real,A<:Area{T}} - B::Vector{T} - C::Vector{A} +struct Sellmeier{T<:Real,Bt,Ct<:Area{T},Dt} + A::T + B::Vector{Bt} + p::Vector{Int} + C::Vector{Ct} + D::Dt + q::Vector{Int} end -function show(io::IO, m::Medium) - write(io, "Medium(") +@doc raw""" + Sellmeier(B, C) + +Convenience constructor for [`Sellmeier`](@ref) for the canonical form of +the Sellmeier equation: + +```math +n^2(\lambda) = +1 + \sum_i \frac {B_i\lambda^2}{\lambda^2-C_i}. +``` +""" +Sellmeier(B::Vector{T}, C::Vector{Ct}) where {T<:Real,Ct<:Area{T}} = + Sellmeier(zero(T), B, fill(2, length(B)), C, zeros(T,0), zeros(Int, 0)) + +function show(io::IO, m::Sellmeier) + write(io, "Sellmeier(") + show(io, m.A) + write(io, ", ") show(io, m.B) - write(io, @sprintf(", [%s])", - join(string.(m.C), ", "))) + write(io, ", ") + show(io, m.p) + printfmt(io, ", [{1:s}]", + join(string.(m.C), ", ")) + printfmt(io, ", [{1:s}]", + join(string.(m.D), ", ")) + write(io, ", ") + show(io, m.q) + write(io, ")") end +Base.Broadcast.broadcastable(m::Sellmeier) = Ref(m) + +const AnisotropicMaterial = SVector{<:Any,<:Sellmeier} + # ** BK7 -""" +@doc raw""" BK7 [Borosilicate glass](https://en.wikipedia.org/wiki/Borosilicate_glass) is commonly used in optical lenses + +```math +n^2(\lambda) = +1 + +\frac{1.03961212\lambda^2}{\lambda^2 - 6.00069867\times10^{-3}\;\textrm{μm}^2} + +\frac{0.231792344\lambda^2}{\lambda^2 - 2.00179144\times10^{-2}\;\textrm{μm}^2} + +\frac{1.01046945}{\lambda^2 - 1.03560653\times10^{2}\;\textrm{μm}^2}. +``` """ -BK7 = Medium([1.03961212, 0.231792344, 1.01046945], - [6.00069867e-3u"μm"^2, 2.00179144e-2u"μm"^2, 1.03560653e2u"μm"^2]) +const BK7 = Sellmeier([1.03961212, 0.231792344, 1.01046945], + [6.00069867e-3u"μm"^2, 2.00179144e-2u"μm"^2, 1.03560653e2u"μm"^2]) # ** SiO₂ -""" +@doc raw""" SiO₂ [Fused silica](https://en.wikipedia.org/wiki/Fused_quartz#Optical_properties) + +```math +n^2(\lambda) = +1 + +\frac{0.6961663\lambda^2}{\lambda^2-(0.0684043\;\textrm{μm})^2} + +\frac{0.4079426\lambda^2}{\lambda^2-(0.1162414\;\textrm{μm})^2} + +\frac{0.8974794\lambda^2}{\lambda^2-(9.896161\;\textrm{μm})^2}. +``` """ -SiO₂ = Medium([0.6961663, 0.4079426, 0.8974794], - [0.0684043u"μm", 0.1162414u"μm", 9.896161u"μm"].^2) +const SiO₂ = Sellmeier([0.6961663, 0.4079426, 0.8974794], + [0.0684043u"μm", 0.1162414u"μm", 9.896161u"μm"].^2) -# * Refractive index +# ** Calcite -function n²(m::Medium, λ::Length) - if !isfinite(λ) - 1.0 - else - 1.0 + sum([Bᵢ*λ^2/(λ^2-Cᵢ) - for (Bᵢ,Cᵢ) ∈ zip(m.B, m.C)]) |> NoUnits - end -end +@doc raw""" + Calcite -maybe_complex(x) = x<0 ? complex(x) : x -refractive_index(m::Medium, λ::Length) = √(maybe_complex(n²(m, λ))) -refractive_index(m::Medium, f::Frequency) = refractive_index(m, u"c"/f) +Calcite is a uniaxial, birefringent crystal. -(m::Medium)(x) = refractive_index(m, x) +Numbers taken from . -# * Dispersion +```math +\begin{aligned} +n_{\mathrm{o}}^2(\lambda) +&= +2.69705 + +\frac{0.0192064}{\lambda^2-0.01820\;\textrm{μm}^2} - +0.0151624\lambda^2, \\ +n_{\mathrm{e}}^2(\lambda) +&= +2.18438 + +\frac{0.0087309}{\lambda^2-0.01018\;\textrm{μm}^2} - +0.0024411\lambda^2. +\end{aligned} +``` +""" +const Calcite = SVector(Sellmeier(1.69705, [0.0192064u"μm^2"], [0], [0.01820u"μm^2"], (-0.0151624u"μm^-2",), [2]), + Sellmeier(1.18438, [0.0087309u"μm^2"], [0], [0.01018u"μm^2"], (-0.0024411u"μm^-2",), [2])) + +# ** Quartz +@doc raw""" + Quartz + +Quartz is a uniaxial, birefrigent crystal. + +Numbers taken from . + +```math +\begin{aligned} +n_{\mathrm{o}}^2(\lambda) +&=2.3573 - +0.01170\;\textrm{μm}^{-2}\lambda^2 + +\frac{0.01054\;\textrm{μm}^{2}}{\lambda^2} + +\frac{1.3414\times10^{-4}\;\textrm{μm}^{4}}{\lambda^4} - +\frac{4.4537\times10^{-7}\;\textrm{μm}^{6}}{\lambda^6} + +\frac{5.9236\times10^{-8}\;\textrm{μm}^{8}}{\lambda^8}, \\ +n_{\mathrm{e}}^2(\lambda) +&=2.3849 - +0.01259\;\textrm{μm}^{-2}\lambda^2 + +\frac{0.01079\;\textrm{μm}^{2}}{\lambda^2} + +\frac{1.6518\times10^{-4}\;\textrm{μm}^{4}}{\lambda^4} - +\frac{1.9474\times10^{-6}\;\textrm{μm}^{6}}{\lambda^6} + +\frac{9.3648\times10^{-8}\;\textrm{μm}^{8}}{\lambda^8}. +\end{aligned} +``` """ - dispersion(m, d, f[, f₀=0u"Hz]) +const Quartz = SVector(Sellmeier(1.3573, Float64[], Int[], fill(0.0u"μm^2", 0), + (-0.01170u"μm^-2", 0.01054u"μm^2", 1.3414e-4u"μm^4", -4.4537e-7u"μm^6", 5.9236e-8u"μm^8"), [2,-2,-4,-6,-8] ), + Sellmeier(1.3849, Float64[], Int[], fill(0.0u"μm^2", 0), + (-0.01259u"μm^-2", 0.01079u"μm^2", 1.6518e-4u"μm^4", -1.9474e-6u"μm^6", 9.3648e-8u"μm^8"), [2,-2,-4,-6,-8])) + +# ** KTP + +# """ +# KTP + +# KTP is a biaxial crystal. -Calculate dispersion through a medium `m` of length `d`. Optionally, -remove central frequency k-vector, to keep pulse temporally centred. +# Numbers taken from . +# """ +# const KTP = SVector(Sellmeier(2.0065, [0.03901u"μm^2"], [0], [0.04251u"μm^2"], (-0.01327u"μm^-2",), [2]), +# Sellmeier(2.0333, [0.04154u"μm^2"], [0], [0.04547u"μm^2"], (-0.01408u"μm^-2",), [2]), +# Sellmeier(2.0065, [0.05694u"μm^2"], [0], [0.05658u"μm^2"], (-0.01682u"μm^-2",), [2])) + +@doc raw""" + KTP + +KTP is a biaxial crystal. + +Numbers taken from . + +```math +\begin{aligned} +n_x^2(\lambda) +&= 2.10468 + +\frac{0.89342\lambda^2}{\lambda^2-0.04438\;\textrm{μm}^2} - +0.01036\;\textrm{μm}^{-2}\lambda^2, \\ +n_y^2(\lambda) +&= 2.14559 + +\frac{0.87629\lambda^2}{\lambda^2-0.0485\;\textrm{μm}^2} - +0.01173\;\textrm{μm}^{-2}\lambda^2, \\ +n_z^2(\lambda) +&= 1.9446 + +\frac{1.3617\lambda^2}{\lambda^2-0.047\;\textrm{μm}^2} - +0.01491\;\textrm{μm}^{-2}\lambda^2. \\ +\end{aligned} +``` """ -function dispersion(m::Medium, d::Length, f::AbstractVector{F}, f₀::Frequency = 0u"Hz") where {F<:Frequency} - n = real(m.(f)) +const KTP = SVector(Sellmeier(1.10468, [0.89342], [2], [0.04438u"μm^2"], (-0.01036u"μm^-2",), [2]), + Sellmeier(1.14559, [0.87629], [2], [0.0485u"μm^2"], (-0.01173u"μm^-2",), [2]), + Sellmeier(0.9446, [1.3617], [2], [0.047u"μm^2"], (-0.01491u"μm^-2",), [2])) - # Wavevector in vacuum - k₀ = 2π*f./u"c" - # Wavevector in the medium - k = n.*k₀ +# * Refractive index - if !iszero(f₀) - # Slope of dispersion relation at central frequency - ∂k∂f = (2π/u"c")*ForwardDiff.derivative(f -> real(m(f*u"Hz"))*f, f₀/u"Hz" .|> NoUnits) - k -= ∂k∂f*f - end +function n²(m::Sellmeier, λ::Length) + O = one(ustrip(λ)) + v = O + m.A + isinf(λ) && return v - exp.(-im*k*d) + v += sum([NoUnits(Bᵢ*λ^pᵢ/(λ^2*(1+im*eps())-Cᵢ)) + for (Bᵢ,pᵢ,Cᵢ) ∈ zip(m.B, m.p, m.C)], init=zero(O)) + v += sum([NoUnits(Dⱼ*λ^qⱼ) + for (Dⱼ,qⱼ) ∈ zip(m.D,m.q)], init=zero(O)) + + v +end +n²(m::Sellmeier, f::Frequency) = n²(m, u"c"/f) + +n²(m::AbstractVector{<:Sellmeier}, v) = n².(m, v) + +function n²_3d(m::Sellmeier, v) + nn = n²(m, v) + SVector(nn, nn, nn) +end + +function n²_3d(m::SVector{2}, v) + nₒ²,nₑ² = n².(m, v) + SVector(nₒ², nₒ², nₑ²) +end + +n²_3d(m::SVector{3}, v) = n²(m, v) + +maybe_complex(x) = x<0 ? complex(x) : x +maybe_complex(x::Complex) = x +refractive_index(m::Sellmeier, v) = √(maybe_complex(n²(m, v))) + +(m::Sellmeier)(x) = refractive_index(m, x) + +# * Dispersion + +function dispersion_slope(m::Sellmeier, ω₀) + f₀ = ustrip(auconvert(u"Hz", ω₀/2π)) + inv(austrip(1u"c"))*ForwardDiff.derivative(f -> real(m(f*u"Hz"))*f, f₀) end + +# We can only remove a common slope, since the different axes _should_ +# be delayed with respect to each other. +dispersion_slope(m::AnisotropicMaterial, ω₀) = + mean(dispersion_slope.(m, ω₀)) + # * Exports -export Medium, BK7, SiO₂, refractive_index, dispersion +export Sellmeier, BK7, SiO₂, Calcite, Quartz, KTP, refractive_index, dispersion, dispersion_slope diff --git a/src/spectra.jl b/src/spectra.jl index 6aa9e05..3406338 100644 --- a/src/spectra.jl +++ b/src/spectra.jl @@ -120,9 +120,50 @@ temporal grid `t`. nfft_vector_potential(f::AbstractField, t::AbstractRange) = _nfft(fft_vector_potential(f, t), t) + +rfftω(args...) = 2π*rfftfreq(args...) +""" + rfftω(t::AbstractRange) + +Return the angular frequency grid corresponding to uniform sampling in +time with the tempolar grid `t`, in the case of real-valued signals. +""" +rfftω(t::AbstractRange) = rfftω(length(t), 1/step(t)) + +""" + rfft(f::AbstractField, t::AbstractRange) + +Compute the RFFT of the [`field_amplitude`](@ref) of the field `f` +sampled on the uniform temporal grid `t`, in the case of real-valued +signals. +""" +FFTW.rfft(f::AbstractField, t::AbstractRange) = + rfft(field_amplitude(f, t), 1) + +""" + rfft_vector_potential(f::AbstractField, t::AbstractRange) + +Compute the RFFT of the [`vector_potential`](@ref) of the field `f` +sampled on the uniform temporal grid `t`, in the case of real-valued +signals. +""" +rfft_vector_potential(f::AbstractField, t::AbstractRange) = + rfft(vector_potential(f, t), 1) + +""" + irfft(F, t) + +Compute the IRFFT of the spectral amplitudes `F`, with the end-result +resolved on the time axis `t`, in the case of real-valued signals. +""" +FFTW.irfft(F::AbstractVecOrMat, t::AbstractRange) = + irfft(F, length(t), 1) + # * Exports export spectrum, fftω, fftshift, fft, fft_vector_potential, - nfft, nfft_vector_potential + nfft, nfft_vector_potential, + rfftω, + rfft, rfft_vector_potential, irfft diff --git a/src/time_axis.jl b/src/time_axis.jl index 3e94059..dbb2073 100644 --- a/src/time_axis.jl +++ b/src/time_axis.jl @@ -22,9 +22,9 @@ as default, ``f_s=100f_{\textrm{max}}``. This also makes plots nicer. If `fs<:Integer`, the sampling frequency is computed as `fs/T`, i.e. `fs` steps per cycle. """ -steps(f::AbstractField, fs::Real=default_sampling_frequency(f)) = - ceil(Int, fs*abs(-(endpoints(span(f))...))) -steps(f::AbstractField, ndt::Int) = steps(f, ndt/austrip(period(f))) +steps(f::AbstractField, fs::Real=default_sampling_frequency(f); s=span(f)) = + ceil(Int, fs*abs(-(endpoints(s)...))) +steps(f::AbstractField, ndt::Int; kwargs...) = steps(f, ndt/austrip(period(f)); kwargs...) """ timeaxis(f[, fs]) @@ -33,9 +33,8 @@ Construct a time axis for the field `f`, covering the [`span`](@ref) of the envelope in the number of [`steps`](@ref) given by the sample frequency `fs`. """ -function timeaxis(f::AbstractField, fs::Number=default_sampling_frequency(f)) - num_steps = steps(f, fs) - s = span(f) +function timeaxis(f::AbstractField, fs::Number=default_sampling_frequency(f); s=span(f)) + num_steps = steps(f, fs; s=s) if num_steps > 1 || s.left == s.right range(s, length=num_steps) else @@ -47,6 +46,7 @@ function timeaxis(f::AbstractField, fs::Number=default_sampling_frequency(f)) end end + """ timeaxis(f, (N,dt)) diff --git a/test/apodized_fields.jl b/test/apodized_fields.jl index 82ca27a..4a0d01e 100644 --- a/test/apodized_fields.jl +++ b/test/apodized_fields.jl @@ -21,6 +21,59 @@ mid_sel2 = findall(∈(7.5 .. 14.0), t) post_sel = findall(>(14.0), t) + @testset "Cosine sum window macro" begin + rl!(ex) = Base.remove_linenums!(ex) + csw(a) = rl!(ElectricFields.cosine_sum_window(:TestWindow, a, "Test Window")) + + function compare_terms(a::Number, b::Number; kwargs...) + @test isapprox(a, b; kwargs...) + end + + function compare_terms(a, b) + if a.head == :call + @assert b.head == :call + f = first(a.args) + @assert f == first(b.args) + @assert length(a.args) == length(b.args) + if f == :(*) || f == :(+) + for j = 2:length(a.args) + compare_terms(a.args[j], b.args[j]) + end + else + @test a == b + end + else + @test a == b + end + end + + function test_window(name, a, exp_doc, exp_expr, exp_der) + q = csw(a) + @testset "$(name)" begin + doc = q.args[1].args[3] + @test doc.head == :string + @test doc.args[2] == :TestWindow + @test doc.args[4] == "Test Window" + @test doc.args[6] == exp_doc + + val = q.args[3].args[2].args[2] + compare_terms(val, exp_expr) + + der = q.args[4].args[2].args[2] + compare_terms(der, exp_der) + end + end + test_window("Zero", (), "0", :(zero(x)), :(zero(x))) + test_window("Half", (0.5), "0.5", :(0.5), :(zero(x))) + test_window("One", (1), "1", :(1), :(zero(x))) + test_window("Two", (1,1), "1 + 1\\cos(2 \\pi x)", + :(1 + 1 * cospi(2x)), + :(-6.283185307179586 * sinpi(2x))) + test_window("Three", (1,1,1), "1 + 1\\cos(2 \\pi x) + 1\\cos(4 \\pi x)", + :(1 + 1 * cospi(2x) + 1 * cospi(4x)), + :(-6.283185307179586 * sinpi(2x) + -12.566370614359172 * sinpi(4x))) + end + ElectricFields.@cosine_sum_window Zero () "Zero" ElectricFields.@cosine_sum_window One (1.0,) "One" @@ -32,7 +85,8 @@ ElectricFields.BlackmanHarris(), ElectricFields.Kaiser(3), ElectricFields.Kaiser(2), ElectricFields.Rect(), Zero(), One()) - w = ElectricFields.window_value.(window, x) + wx = Base.Fix1(ElectricFields.window_value, window) + w = wx.(x) ∂w = ElectricFields.window_derivative.(window, x) @test all(iszero, w[pre_sel]) @@ -50,6 +104,9 @@ # the derivative, just that the behaviour is reasonable. @test all(≥(0), ∂w[mid_sel1]) @test all(≤(0), ∂w[mid_sel2]) + # Here we however compare the analytic derivatives with + # AD. + test_approx_eq(ForwardDiff.derivative.(wx, x[2:end-1]), ∂w[2:end-1]) end @test all(iszero, w[post_sel]) @@ -82,7 +139,7 @@ # For the same reason, the intensity is not exactly the original # intensity (flat for this particular field) times the window # value squared, hence the high tolerance. - @test abs2.(w[mid_sel]) ≈ Iwv[mid_sel] rtol=5e-2 + test_approx_eq(abs2.(w[mid_sel]), Iwv[mid_sel], rtol=0.2) @test all(iszero, Fwv[post_sel]) @test all(iszero, Awv[pre_sel]) @@ -97,7 +154,7 @@ - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹ - A₀ = 1.0000 au – a Fixed carrier @ λ = 45.5634 nm (T = 151.9830 as, ω = 1.0000 Ha = 27.2114 eV, f = 6.5797 PHz) - – and a /0‾3‾0\\ cycles trapezoidal envelope + – and a /0‾3‾0\\ cycles trapezoidal envelope, with linear ramps – 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""" end diff --git a/test/arithmetic.jl b/test/arithmetic.jl index 6e99a1d..1d9369c 100644 --- a/test/arithmetic.jl +++ b/test/arithmetic.jl @@ -310,6 +310,9 @@ end @test field_amplitude(phase_shift(B, 2.0), 0.4) == field_amplitude(phase_shift(A, 2.0), 0.4) @test field_amplitude(phase_shift(B, 2.0), 3.0) == zero(field_amplitude(phase_shift(A, 2.0), 0.4)) + @test field_amplitude(B, -5, -4) ≈ field_amplitude(A, austrip(-0.1u"fs"), -4) + @test field_amplitude(B2, -5, -4) ≈ field_amplitude(A2, austrip(-0.1u"fs"), -4) + R = [1 0 0; 0 1 1; 0 -1 1] rB = test_rotated_field(B, R) @test rB isa ElectricFields.WindowedField diff --git a/test/bspline_fields.jl b/test/bspline_fields.jl new file mode 100644 index 0000000..de93629 --- /dev/null +++ b/test/bspline_fields.jl @@ -0,0 +1,60 @@ +@testset "B-spline fields" begin + @testset "$(n)d fields" for n in (1,3) + @field(F) do + ω = 1.0 + I₀ = 1.0 + σ = 7.0 + σoff = 4.0 + σmax = 6.0 + env = :trunc_gauss + end + if n == 3 + F = rotate(F, ElectricFields.compute_rotation((π/3, [0.4,1,0]))) + end + + t = timeaxis(F) + Av = vector_potential(F, t) + Fv = field_amplitude(F, t) + + spline_order = 7 + num_steps = length(t) + num_knots = 100 + + BB = ElectricFields.BSpline(ElectricFields.LinearKnotSet(spline_order, t[1], t[end], num_knots)) + @testset "B-spline $(case)" for case in ("full range", "restricted") + B = if case == "restricted" + BB[:,begin+1:end-1] + else + BB + end + + FB = BSplineField(B, t, Av) + + @test dimensions(FB) == n + + @test pretty_print_object(FB) == + if case == "restricted" + "B-spline field expanded over\n Restriction to basis functions 2:105 of BSpline basis with LinearKnotSet(Float64) of order k = 7 on -42.0 .. 42.0 (100 intervals)" + else + "B-spline field expanded over\n BSpline basis with LinearKnotSet(Float64) of order k = 7 on -42.0 .. 42.0 (100 intervals)" + end + + Arec = vector_potential(FB, t) + Frec = field_amplitude(FB, t) + + test_approx_eq(Av, Arec, rtol=1e-4) + test_approx_eq(Fv, Frec, rtol=1e-4) + + @testset "Scalar evaluation" begin + sel = 1:1000 + correct_tensor(::LinearPolarization, v) = v + correct_tensor(::ArbitraryPolarization, v) = transpose(reduce(hcat, v)) + Arec2 = correct_tensor(polarization(FB), vector_potential.(Ref(FB), t[sel])) + Frec2 = correct_tensor(polarization(FB), field_amplitude.(Ref(FB), t[sel])) + + test_approx_eq(selectdim(Arec, 1, sel), Arec2, rtol=1e-14) + test_approx_eq(selectdim(Frec, 1, sel), Frec2, rtol=1e-14) + end + end + end +end diff --git a/test/bsplines.jl b/test/bsplines.jl new file mode 100644 index 0000000..3eb39e5 --- /dev/null +++ b/test/bsplines.jl @@ -0,0 +1,365 @@ +using IntervalSets +import ElectricFields: RightContinuous, find_interval, within_support, +LinearKnotSet, ArbitraryKnotSet, ExpKnotSet, +numfunctions, BSpline, nonempty_intervals, numintervals, +knotset, order +using SparseArrays + +within_interval(x::AbstractRange, interval::Interval) = findall(in(interval), x) + +""" +within_interval_linear(x, interval) + +Return the range of indices of elements of the sorted vector/range `x` +that fall within `interval`. This loops through all elements of `x` +(linear complexity), so it should not be used for large number of +elements. +""" +function within_interval_linear(x, interval) + sel = findall(e -> e ∈ interval, x) + isempty(sel) ? (1:0) : (minimum(sel):maximum(sel)) +end + +function test_within_interval(x, interval, expected=nothing; reversed=true) + if reversed + x = reverse(x) + end + N = length(x) + + linear_expected = within_interval_linear(x, interval) + if isnothing(expected) + expected = linear_expected + else + if reversed + expected = (N-expected[end]+1):(N-expected[1]+1) + end + # This is mostly a sanity check to test if the expectations + # are actually correct. + expected ≠ linear_expected && + throw(ArgumentError("Invalid expectations: $(expected), actual: $(linear_expected)")) + end + + result = :(within_interval($x, $interval)) + expr = :($result == $expected || isempty($result) && isempty($expected)) + if !(@eval $expr) + println("Looking for elements of $x ∈ $interval, got $(@eval $result), expected $expected") + length(x) < 30 && println(" x = ", collect(enumerate(x)), "\n") + end + @eval @test $expr +end + +@testset "B-splines" begin + @testset "Unit vectors" begin + e₆ = ElectricFields.UnitVector{Bool}(10, 6) + @test size(e₆) == (10,) + @test e₆ == vcat(zeros(Bool, 5), one(Bool), zeros(Bool, 4)) + @test string(e₆) == "ê{Bool}(6|10)" + end + + @testset "Knot sets" begin + @testset "Simple tests" begin + @test ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3) ≈ [0.0, 1, 1, 3, 4, 6, 6, 6] + @test LinearKnotSet(3, 0, 1, 2) ≈ [0,0,0,0.5,1,1,1] + @test LinearKnotSet(2, 0, 1, 2, 1, 1) ≈ [0,0.5,1] + @test ExpKnotSet(2, -4, 2, 7) ≈ [0,0,0.0001,0.001,0.01,0.1,1,10,100,100] + @test string(ExpKnotSet(2, -4, 2, 7)) == "ExpKnotSet(Float64) of order k = 2 (linear) on 0,0.0001..100.0 (7 intervals)" + end + + @testset "k = $k" for k ∈ 1:6 + @testset "Full multiplicity" begin + t = LinearKnotSet(k, 0, 1, 1) + @test length(t) == 2k + @test collect(t) == vcat(fill(0, k), fill(1, k)) + @test length(nonempty_intervals(t)) == 1 + end + @testset "Simple endpoints" begin + t = LinearKnotSet(k, 0, 1, k, 1, 1) + @test length(t) == k+1 + @test collect(t) == range(0, stop=1, length=k+1) + @test length(nonempty_intervals(t)) == k + end + end + @test length(nonempty_intervals(ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3))) == 4 + + @testset "Invalid order/multiplicities" begin + @testset "Invalid order" begin + @test_throws ArgumentError LinearKnotSet(0, 0, 1, 4) + @test_throws ArgumentError LinearKnotSet(-1, 0, 1, 4) + end + @testset "Invalid multiplicities" begin + @test_throws ArgumentError LinearKnotSet(1, 0, 1, 4, 0, 1) + @test_throws ArgumentError LinearKnotSet(1, 0, 1, 4, 1, 0) + @test_throws ArgumentError LinearKnotSet(1, 0, 1, 4, 1, 2) + end + @testset "Not enough interior points" begin + @test_throws ArgumentError LinearKnotSet(2, 0, 1, 1, 1, 1) + end + end + + @testset "Properties" begin + @testset "#intervals = $N" for N = 1:4 + @testset "k = $k" for k = 1:N + @testset "Simple multiplicities" begin + t = LinearKnotSet(k, 0.0, 1.0, N, 1, 1) + @test numintervals(t) == N + @test numfunctions(t) == N - k + 1 + end + end + end + end + + @testset "Compact support" begin + a,b = 0,1 + x = range(a, stop=b, length=21) + @testset "Interval coverage" begin + @testset "Two intervals" begin + @testset "Reversed: $reversed" for reversed in [false, true] + test_within_interval(x, 0..0.5, 1:11, reversed=reversed) + end + @testset "L=$L" for L=[:closed,:open] + @testset "R=$R" for R=[:closed,:open] + @testset "Reversed: $reversed" for reversed in [false, true] + test_within_interval(x, Interval{L,R}(0,0.5), reversed=reversed) + test_within_interval(x, Interval{L,R}(0.25,0.5), reversed=reversed) + end + end + end + end + @testset "Three intervals" begin + @testset "Reversed: $reversed" for reversed in [false, true] + test_within_interval(x, RightContinuous(0,1/3), 1:7, reversed=reversed) + test_within_interval(x, RightContinuous(1/3,2/3), 8:14, reversed=reversed) + test_within_interval(x, 2/3..1, 15:21, reversed=reversed) + end + end + @testset "Open interval" begin + @testset "Reversed: $reversed" for reversed in [false, true] + test_within_interval(x, OpenInterval(0.2,0.4), 6:8, reversed=reversed) + end + end + @testset "Random intervals" begin + @testset "L=$L" for L=[:closed,:open] + @testset "R=$R" for R=[:closed,:open] + @testset "Reversed: $reversed" for reversed in [false, true] + for i = 1:20 + interval = Interval{L,R}(minmax(rand(),rand())...) + test_within_interval(x, interval, reversed=reversed) + end + end + end + end + end + end + + @testset "Partially covered intervals" begin + k = 3 + t = LinearKnotSet(k, 0, 1, 2) + @testset "$name, x = $x" for (name,x) in [ + ("Outside left",range(-1,stop=-0.5,length=10)), + ("Touching left",range(-1,stop=0,length=10)), + ("Touching left-ϵ",range(-1,stop=0-eps(),length=10)), + ("Touching left+ϵ",range(-1,stop=0+eps(),length=10)), + + ("Outside right",range(1.5,stop=2,length=10)), + ("Touching right",range(1,stop=2,length=10)), + ("Touching right-ϵ",range(1-eps(),stop=2,length=10)), + ("Touching right+ϵ",range(1+eps(),stop=2,length=10)), + + ("Other right",range(0.5,stop=1,length=10)), + ("Other right-ϵ",range(0.5-eps(0.5),stop=1,length=10)), + ("Other right+ϵ",range(0.5+eps(0.5),stop=1,length=10)), + + ("Complete", range(0,stop=1,length=10)), + ("Complete-ϵ", range(eps(),stop=1-eps(),length=10)), + ("Complete+ϵ", range(-eps(),stop=1+eps(),length=10)), + + ("Left partial", range(-0.5,stop=0.6,length=10)), + ("Left", range(-0.5,stop=1.0,length=10)), + ("Right partial", range(0.5,stop=1.6,length=10)), + ("Right", range(0,stop=1.6,length=10))] + @testset "L=$L" for L=[:closed,:open] + @testset "R=$R" for R=[:closed,:open] + @testset "Reversed: $reversed" for reversed in [false, true] + for i in nonempty_intervals(t) + interval = Interval{L,R}(t[i], t[i+1]) + test_within_interval(x, interval, reversed=reversed) + end + end + end + end + end + end + end + + @testset "Finding intervals" begin + t = ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3) + R = BSpline(t) + # Test without specifying initial interval + @test find_interval(t, 0.5) == 1 + for (i,xs) in enumerate([[0, 0.5, 1.0-eps()], + [], + [1.0, 2.0, 3.0-eps(3.0)], + [3, 3.5, 4.0-eps(4.0)], + [4, 5, 6.0-eps(6.0),6]]) + for j in [1,i] # Test starting from the beginning, and from the interval itself + for x in xs + result = :(find_interval($t, $x, $j)) + res = @eval $result + expr = :($result == $i) + (@eval $expr) || + println("Expected to find $x ∈ interval #$i, got $(isnothing(res) ? "nothing" : res)") + @eval @test $expr + # @test find_interval(t, x) == + find_interval(t, x) + within_support(x:x, t, i)[1][2] + end + end + end + @test isnothing(find_interval(t, 0.5, 2)) + end + + @testset "Support of Heavyside splines" begin + a,b = 0,1 + x = range(a, stop=b, length=21) + t = LinearKnotSet(1, a, b, 4) + supports = [within_support(x, t, j) + for j = 1:numfunctions(t)] + @test length(supports) == 4 + # Each basis function should cover one interval only (since order = 1). + @test all(length.(supports) .== 1) + # Test that all elements of x are covered by the basis + # functions, and that none of the basis functions overlap. + @test first.(first.(supports)) == [1:5, 6:10, 11:15, 16:21] + end + + @testset "Support of linear splines" begin + a,b = 0,1 + x = range(a, stop=b, length=21) + @testset "Simple multiplicity" begin + t = LinearKnotSet(2, a, b, 2, 1, 1) + supports = [within_support(x, t, j) + for j = 1:numfunctions(t)] + @test length(supports) == 1 + @test supports[1] == [(1:10,1), (11:21,2)] + end + @testset "Full multiplicity" begin + @testset "#intervals = 2" begin + t = LinearKnotSet(2, a, b, 2) + supports = [within_support(x, t, j) + for j = 1:numfunctions(t)] + @test length(supports) == 3 + @test supports[1] == [(1:10,2)] + @test supports[2] == [(1:10,2), (11:21,3)] + @test supports[3] == [(11:21,3)] + end + @testset "#intervals = 3" begin + t = LinearKnotSet(2, a, b, 3) + supports = [within_support(x, t, j) + for j = 1:numfunctions(t)] + @test length(supports) == 4 + @test supports[1] == [(1:7,2)] + @test supports[2] == [(1:7,2), (8:14,3)] + @test supports[3] == [(8:14,3), (15:21,4)] + @test supports[4] == [(15:21,4)] + end + end + end + end + + @testset "Splines" begin + @testset "Discontinuous knot set" begin + t = ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3) + B = BSpline(t) + + @test firstindex(B,2) == 1 + @test lastindex(B,2) == 5 + @test size(B,2) == 5 + + @test knotset(B) == t + @test order(B) == 3 + + @test all(B[0.0, :] .== 0) + @test all(B[0.0:0.5:1.0, :][1, :] .== 0) + + @test B[0.0, 1] == 0.0 + @test B[0.1, 1] ≈ 0.01 + end + + @testset "Evaluate B-splines" begin + k = 3 + t = LinearKnotSet(k, 0, 1, 2) + B = BSpline(t) + @testset "Eval on subintervals" begin + x₁ = range(-1,stop=-0.5,length=10) + χ₁ = B[x₁, :] + @test norm(χ₁) == 0 + + function testbasis(x) + χ = B[x, :] + B̃ = spzeros(Float64, length(x), 4) + B̃[:,1] = (x .>= 0) .* (x .< 0.5) .* ((2x .- 1).^2) + B̃[:,2] = (x .>= 0) .* (x .< 0.5) .* (2/3*(1 .- (3x .- 1).^2)) + + (x .>= 0.5) .* (x .< 1) .* (2*(x .- 1).^2) + B̃[:,3] = (x .>= 0) .* (x .< 0.5) .* (2*x.^2) + + (x .>= 0.5) .* (x .< 1) .* (2/3*(1 .- (3x .- 2).^2)) + B̃[:,4] = (x .>= 0.5) .* (x .<= 1) .* ((2x .- 1).^2) + for j = 1:4 + test_approx_eq(χ[:,j],B̃[:,j], atol=1e-15, rtol=1e-15) + end + end + testbasis(range(0,stop=1,length=50)) + testbasis(range(-0.5,stop=0.6,length=40)) + testbasis(range(0.5,stop=1.6,length=40)) + end + end + + @testset "Restricted basis" begin + t = LinearKnotSet(3, 0.0, 1.0, 6) + B = BSpline(t) + B̃ = B[:,2:end-1] + @test size(B̃,2) == 6 + χ̃ = B̃[range(0, stop=1, length=10), :] + @test size(χ̃,2) == 6 + end + + @testset "Evaluate spline" begin + @testset "Linear knot set" begin + k = 4 + t = LinearKnotSet(k, 0, 1, 5) + B = BSpline(t) + x = range(0, stop=1, length=301) + χ = B[x, :] + + i = 1:size(B,2) + c = sin.(i) + + s′ = χ*c + @test all(isfinite, s′) + @test all(!isnan, s′) + end + @testset "Discontinuous knot set" begin + t = ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3) + B = BSpline(t) + + x = range(0, stop=6, length=301) + χ = B[x, :] + @test all(isfinite, χ) + @test all(!isnan, χ) + end + end + + @testset "Function interpolation" begin + B = BSpline(LinearKnotSet(7, 0, 1, 10)) + x = range(0, stop=1, length=100) + c = B[x,:] \ sin.(2π*x) + x̃ = range(0, stop=1, length=301) + χ = B[x̃,:] + test_approx_eq(χ*c, sin.(2π*x̃), rtol=2e-7) + + B̃ = B[:,2:end-1] + c̃ = B̃[x,:] \ sin.(2π*x) + χ̃ = B̃[x̃,:] + test_approx_eq(χ̃*c̃, sin.(2π*x̃), rtol=2e-7) + end + end +end diff --git a/test/dispersed_fields.jl b/test/dispersed_fields.jl new file mode 100644 index 0000000..93b3531 --- /dev/null +++ b/test/dispersed_fields.jl @@ -0,0 +1,177 @@ +import ElectricFields: CascadedDispersiveElement, frequency_response + +@testset "Dispersed fields" begin + @testset "Chirped fields" begin + function chirped_gaussian(t::Number, ω₀, ϕ₀, τ, η; verbose = false) + γ = τ^2/(8log(2)) + g = inv(γ^2 + η^2) + + ϕη = atan(-η,γ)/2 + A = √(γ*√g)*exp(im*ϕη) + + aη = γ*g/4 + bη = η*g/2 + + verbose && @info "Pulse params" ω₀ γ η ϕ₀ ϕη A aη bη + + ca = exp(im*(ω₀*t + ϕ₀ + bη/2*t^2)) + env = exp(-aη*t^2) + + ca′ = im*(ω₀ + bη*t)*exp(im*(ω₀*t + ϕ₀ + bη/2*t^2)) + env′ = (-2aη*t)*exp(-aη*t^2) + 1/ω₀*imag(A*ca*env), -1/ω₀*imag(A*(ca′*env + ca*env′)) + end + + function chirped_gaussian(t::AbstractVector, args...; kwargs...) + T = eltype(t) + if !isempty(t) + A,F = chirped_gaussian(first(t), args...; kwargs...) + T = promote_type(T, typeof(A), typeof(F)) + end + V = zeros(T, length(t), 2) + for (i,t) in enumerate(t) + V[i,:] .= chirped_gaussian(t, args...; kwargs...) + end + V[:,1],V[:,2] + end + + @testset "$(name) pulse" for (name,τ,ηs,rtol) in ( + ("Short", 3u"fs", 5.0u"fs^2"*[1,-1], 4e-2), + ("Long", 30u"fs", 15.0u"fs^2"*[1,-1], 5e-3)) + @testset "η = $(η)" for η in ηs + τ = austrip(τ) + η = austrip(η) + + @field(F) do + λ = 800u"nm" + I₀ = 1.0 + τ = τ + σoff = 4.0 + σmax = 6.0 + env = :trunc_gauss + end + + ω₀ = photon_energy(F) + Fc = chirp(F, η, ω₀) + + withenv("UNITFUL_FANCY_EXPONENTS" => true) do + s = η > 0 ? "" : "-" + v1,v2 = s .* (name == "Short" ? ("8545.5457", "5.0000") : ("25636.6372","15.0000")) + @test string(Fc.de) == "Chirp(b = $(v1) = $(v2) fs², ω₀ = 0.0570 = 1.5498 eV)" + end + + torig = timeaxis(F) + t = timeaxis(Fc) + + @test t[1] ≲ torig[1] rtol=5e-3 + @test t[end] ≳ torig[end] rtol=5e-3 + + Fv = field_amplitude(Fc, t) + Av = vector_potential(Fc, t) + + F₀ = amplitude(F) + Avref,Fvref = F₀ .* chirped_gaussian(t, ω₀, 0.0, τ, η) + + test_approx_eq(Fv, Fvref, rtol=rtol) + test_approx_eq(Av, Avref, rtol=rtol) + end + end + end + + @testset "Cascaded dispersive elements" begin + a = PhaseShift(π) + b = Chirp(austrip(5u"fs^2"), 1.0) + c = CascadedDispersiveElement(()) + + @test string(c) == "CascadedDispersiveElement(())" + + d = a*b + e = d*b + f = a*d + g = d*d + + @test frequency_response(c, 0.4) == 1 + + @test d == CascadedDispersiveElement((a,b)) + @test frequency_response(d, 0.4) ≈ frequency_response(a, 0.4)*frequency_response(b, 0.4) + @test string(d) == "[$(a) ∗ $(b)]" + + @test e == CascadedDispersiveElement((a,b,b)) + @test frequency_response(e, 0.4) ≈ frequency_response(a, 0.4)*frequency_response(b, 0.4)^2 + @test string(e) == "[$(a) ∗ $(b) ∗ $(b)]" + + @test f == CascadedDispersiveElement((a,a,b)) + @test frequency_response(f, 0.4) ≈ frequency_response(a, 0.4)^2*frequency_response(b, 0.4) + @test string(f) == "[$(a) ∗ $(a) ∗ $(b)]" + + @test g == CascadedDispersiveElement((a,b,a,b)) + @test frequency_response(g, 0.4) ≈ frequency_response(a, 0.4)^2*frequency_response(b, 0.4)^2 + @test string(g) == "[$(a) ∗ $(b) ∗ $(a) ∗ $(b)]" + end + + @testset "Isotropic material" begin + τ = austrip(3u"fs") + + @field(F) do + λ = 800u"nm" + I₀ = 1.0 + τ = τ + σoff = 4.0 + σmax = 12.0 + env = :trunc_gauss + ϕ = π + end + + ω₀ = photon_energy(F) + de = IsotropicMedium(BK7, 220u"μm", ω₀=ω₀) + Fc = DispersedField(F, de, verbosity=4) + Fc2 = phase_shift(Fc, π) + + torig = timeaxis(F) + t = timeaxis(Fc) + + @test t[1] ≲ torig[1] rtol=5e-3 + @test t[end] ≳ torig[end] rtol=5e-3 + + @test field_amplitude(Fc2, t) ≈ -field_amplitude(Fc, t) + + @test string(Fc) == "DispersedField($(F), $(de))" + @test string(Fc2) == "DispersedField($(F), $(Fc2.de))" + + @test pretty_print_object(Fc) == "DispersedField:\n$(pretty_print_object(F))\n – dispersed through $(pretty_print_object(Fc.de))" + @test pretty_print_object(Fc2) == "DispersedField:\n$(pretty_print_object(F))\n – dispersed through $(pretty_print_object(Fc2.de))" + end + + @testset "Crystal($(material))" for material in (BK7,Quartz,KTP) + τ = austrip(3u"fs") + + @field(F) do + λ = 800u"nm" + I₀ = 1.0 + τ = τ + σoff = 4.0 + σmax = 12.0 + env = :trunc_gauss + ϕ = π + end + + ω₀ = photon_energy(F) + de = Crystal(material, 12u"μm", ω₀=ω₀) + Fc = DispersedField(F, de, verbosity=4) + Fc2 = phase_shift(Fc, π) + + torig = timeaxis(F) + t = timeaxis(Fc) + + @test t[1] ≲ torig[1] rtol=5e-3 + @test t[end] ≳ torig[end] rtol=5e-3 + + @test field_amplitude(Fc2, t) ≈ -field_amplitude(Fc, t) + + @test string(Fc) == "DispersedField($(F), $(de))" + @test string(Fc2) == "DispersedField($(F), $(Fc2.de))" + + @test pretty_print_object(Fc) == "DispersedField:\n$(pretty_print_object(F))\n – dispersed through $(pretty_print_object(Fc.de))" + @test pretty_print_object(Fc2) == "DispersedField:\n$(pretty_print_object(F))\n – dispersed through $(pretty_print_object(Fc2.de))" + end +end diff --git a/test/envelopes.jl b/test/envelopes.jl index df36726..278f9b8 100644 --- a/test/envelopes.jl +++ b/test/envelopes.jl @@ -24,7 +24,7 @@ @test continuity(env) == 0 @test span(env) == 0..12 - @test string(env) == "/1‾3‾2\\ cycles trapezoidal envelope" + @test string(env) == "/1‾3‾2\\ cycles trapezoidal envelope, with linear ramps" @field(F2) do I₀ = 1.0 @@ -32,13 +32,29 @@ ramp = 1.0 flat = 3.0 env = :trapezoidal + ramp_kind = :sin² end env2 = envelope(F2) + @test env2(-1.0) == 0.0 + @test env2(0.0) == 0.0 + @test env2(2/3) ≈ 0.25 + @test env2(1.0) ≈ 0.5 + @test env2(4/3) ≈ 0.75 + @test env2(2.5) == 1.0 + @test env2(7.9) == 1.0 + @test env2(8.0+2/3) ≈ 0.75 + @test env2(9.0) ≈ 0.5 + @test env2(8.0+4/3) ≈ 0.25 + @test env2(10.0) == 0.0 + @test env2(11.0) == 0.0 + @test env2.ramp_up == env2.ramp_down == 1.0 @test span(env2) == 0..10 + @test string(env2) == "/1‾3‾1\\ cycles trapezoidal envelope, with sin² ramps" + @test_throws ErrorException ElectricFields.TrapezoidalEnvelope(-1.0, 1.0, 1.0, 1.0) @test_throws ErrorException ElectricFields.TrapezoidalEnvelope(1.0, -1.0, 1.0, 1.0) @test_throws ErrorException ElectricFields.TrapezoidalEnvelope(1.0, 1.0, -1.0, 1.0) diff --git a/test/field_types.jl b/test/field_types.jl index c56f563..ba4fde8 100644 --- a/test/field_types.jl +++ b/test/field_types.jl @@ -60,7 +60,7 @@ - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹ - A₀ = 0.3183 au – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV, f = 20.6707 PHz) - – and a /1‾3‾2\\ cycles trapezoidal envelope + – and a /1‾3‾2\\ cycles trapezoidal envelope, with linear ramps – and a bandwidth of Inf Ha = Inf eV ⟺ Inf Hz ⟺ Inf Bohr = Inf m – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm""" end @@ -71,7 +71,7 @@ - E₀ = 1.0000e+00 au = 514.2207 GV m^-1 - A₀ = 0.3183 au – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV, f = 20.6707 PHz) - – and a /1‾3‾2\\ cycles trapezoidal envelope + – and a /1‾3‾2\\ cycles trapezoidal envelope, with linear ramps – and a bandwidth of Inf Ha = Inf eV ⟺ Inf Hz ⟺ Inf Bohr = Inf m – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm""" end @@ -183,6 +183,7 @@ @test F isa ElectricFields.ConstantField + @test vector_potential(F, -1) == 0 @test field_amplitude(F, -1) == 0 @test field_amplitude(F, 2) ≈ 2 rtol=1e-7 @test field_amplitude(F, 5) == 0 @@ -253,6 +254,8 @@ @test F isa ElectricFields.Ramp + @test vector_potential(F, -1) == 0 + @test field_amplitude(F, -1) == 0 @test field_amplitude(F, 0) ≈ 0 rtol=1e-7 @test field_amplitude(F, 2) ≈ 1 rtol=1e-7 @@ -379,6 +382,13 @@ tCpD = transverse_field(CpD) @test tCpD isa ElectricFields.LinearTransverseField + rtCpD = rotate(tCpD, [0 -1 0 + 1 0 0 + 0 0 1]) + + withenv("UNITFUL_FANCY_EXPONENTS" => true) do + @test pretty_print_object(rtCpD) == "Linearly polarized field in transverse plane constructed from\n┌ Constant field of\n│ - 124.0241 jiffies = 3.0000 fs duration, and\n│ - E₀ = 1.0000e-01 au = 51.4221 GV m⁻¹\n│ – delayed by -124.0241 jiffies = -3.0000 fs\n⊕\n│ sin² down-ramp of\n│ - 124.0241 jiffies = 3.0000 fs duration, and\n└ - E₀ = 1.0000e-01 au = 51.4221 GV m⁻¹\n – and a rotation of 0.50π about [0.000, 0.000, 1.000]" + end end @testset "Rotation of fields" begin diff --git a/test/rotations.jl b/test/rotations.jl index 9840aa0..05dc90d 100644 --- a/test/rotations.jl +++ b/test/rotations.jl @@ -1,5 +1,4 @@ using StaticArrays -using LinearAlgebra import ElectricFields: compute_rotation, rotation_angle, rotation_axis @testset "Rotations" begin @@ -53,4 +52,17 @@ import ElectricFields: compute_rotation, rotation_angle, rotation_axis @test rotation_axis(R) ≈ normalize(axis) end end + + @testset "Rotation matrix" begin + @test rotation_matrix([1.0 1.0 0.0 + 0.0 1.0 1.0 + 0.0 0.0 1.0]) ≈ SMatrix{3,3}(I) rtol=1e-14 + @test rotation_matrix(2*[1.0 0.0 0.0 + 0.0 1.0 0.0 + 0.0 0.0 1.0]) ≈ SMatrix{3,3}(I) rtol=1e-14 + @test rotation_matrix([1.0 1.0 1.0 + 0.0 1.0 1.0 + 0.0 0.0 1.0]) ≈ SMatrix{3,3}(I) rtol=1e-14 + @test rotation_matrix(π*I) ≈ SMatrix{3,3}(I) rtol=1e-14 + end end diff --git a/test/runtests.jl b/test/runtests.jl index e1389bf..5335bc0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,8 +2,10 @@ using ElectricFields using Test using Unitful using UnitfulAtomic +using LinearAlgebra using IntervalSets using FFTW +using ForwardDiff using PrettyTables function test_approx_eq(a, b; on_fail::Union{Nothing,Function}=nothing, isbroken=false, kwargs...) @@ -30,6 +32,15 @@ function test_approx_eq(a, b; on_fail::Union{Nothing,Function}=nothing, isbroken end end +≲(a,b; kwargs...) = a < b || isapprox(a, b; kwargs...) +≳(a,b; kwargs...) = a > b || isapprox(a, b; kwargs...) + +function pretty_print_object(obj) + buf = IOBuffer() + show(buf, MIME"text/plain"(), obj) + String(take!(buf)) +end + @testset "ElectricFields.jl" begin include("namespace_macro.jl") include("units.jl") @@ -83,5 +94,9 @@ end include("spectra.jl") - # include("sellmeier.jl") + include("bsplines.jl") + include("bspline_fields.jl") + + include("sellmeier.jl") + include("dispersed_fields.jl") end diff --git a/test/sellmeier.jl b/test/sellmeier.jl index 8aeee16..6abc90e 100644 --- a/test/sellmeier.jl +++ b/test/sellmeier.jl @@ -2,36 +2,29 @@ @testset "String representation" begin Bstr = string(BK7.B) Cstr = join(string.(BK7.C), ", ") - @test string(BK7) == "Medium($(Bstr), [$(Cstr)])" + @test string(BK7) == "Sellmeier(0.0, $(Bstr), [2, 2, 2], [$(Cstr)], [], $(Int)[])" end @test BK7(0.5876u"μm") ≈ 1.5168 atol=1e-3 @test SiO₂(0.5876u"μm") ≈ 1.4585 atol=1e-3 - @testset "Dispersion" begin - # Dispersing a pulse through a positive amount of glass should - # /delay/ the pulse, i.e. its maximum should arrive /later/, - # and vice versa for a negative amount of glass (achievable - # through precompensation, common in pulse compressors). + Cn² = ElectricFields.n².(Ref(Calcite), [1064,800,632.8,532,400,355,266]*u"nm") + Cno = .√(first.(Cn²)) + Cne = .√(last.(Cn²)) + test_approx_eq(Cno, [1.6423,1.6487,1.6557,1.6629,1.6823,1.6951,1.7497], rtol=5e-5) + test_approx_eq(Cne, [1.4797,1.4821,1.4852,1.4886,1.4974,1.5032,1.5259], rtol=5e-5) - λ = 500.0u"nm" - f₀ = u"c"/λ |> u"THz" - ω₀ = 2π*f₀ - τ = 6.2u"fs" # Pulse duration, intensity FWHM - γ = τ^2/8log(2) + Qn² = ElectricFields.n².(Ref(Quartz), [1064,800,632.8,532,400,355,266]*u"nm") + Qno = .√(first.(Qn²)) + Qne = .√(last.(Qn²)) + test_approx_eq(Qno, [1.5341,1.5384,1.5427,1.5469,1.5577,1.5646,1.5916], rtol=5e-5) + test_approx_eq(Qne, [1.5428,1.5473,1.5517,1.5561,1.5673,1.5744,1.6024], rtol=5e-5) - f = range(0,stop=30,length=2000)*f₀ - ω = 2π*f - - Ê = exp.(-(ω .- ω₀).^2*γ) - Ê′ = Ê.*dispersion(BK7, 6u"μm", f) - Ê′′ = Ê.*dispersion(BK7, -6u"μm", f) - Ê′′′ = Ê.*dispersion(BK7, -6u"μm", f, f₀) - - time_domain_envelope(spectrum) = abs.(fftshift(ifft(spectrum)*√(length(spectrum)))) - - @test argmax(time_domain_envelope(Ê′)) > argmax(time_domain_envelope(Ê)) - @test argmax(time_domain_envelope(Ê′′)) < argmax(time_domain_envelope(Ê)) - @test argmax(time_domain_envelope(Ê′′′)) == argmax(time_domain_envelope(Ê)) - end + Kn² = ElectricFields.n².(Ref(KTP), [1064,532]*u"nm") + Knx = real(.√(ElectricFields.maybe_complex.(first.(Kn²)))) + Kny = real(.√(ElectricFields.maybe_complex.((e -> e[2]).(Kn²)))) + Knz = real(.√(ElectricFields.maybe_complex.(last.(Kn²)))) + test_approx_eq(Knx, [1.7377,1.7780], rtol=5e-3) + test_approx_eq(Kny, [1.7453,1.7886], rtol=5e-3) + test_approx_eq(Knz, [1.8297,1.8887], rtol=5e-3) end diff --git a/test/spectra.jl b/test/spectra.jl index 8f65641..6ee032f 100644 --- a/test/spectra.jl +++ b/test/spectra.jl @@ -1,5 +1,5 @@ @testset "Field spectra" begin - @testset "Gaussian, $(kind)" for kind = [:linear, :transverse, :elliptical] + @testset "Gaussian, $(kind)" for kind = [:linear, :linear_transverse, :transverse, :elliptical] if kind == :elliptical @field(F) do I₀ = 1.0 @@ -9,12 +9,16 @@ ξ = 0.67 end else + k = (kind == :linear_transverse ? :linear : kind) @field(F) do I₀ = 1.0 T = 1.0 τ = 2.0 σmax = 6.0 - kind = kind + kind = k + end + if kind == :linear_transverse + F = transverse_field(F) end end @@ -26,7 +30,7 @@ F̂v = exp.(im*ω*t[1]) .* fftshift(nfft(F, t), 1) Âv = exp.(im*ω*t[1]) .* fftshift(nfft_vector_potential(F, t), 1) - sel = kind == :linear ? 1 : (kind == :transverse ? 3 : Colon()) + sel = kind == :linear ? 1 : (kind ∈ (:linear_transverse, :transverse) ? 3 : Colon()) @test norm(F̂v[:,sel]) ≈ 1.4251872372067347 @test norm(Âv[:,sel]) ≈ 0.22581870287199948