Skip to content

Commit

Permalink
Refactored dipole plotting
Browse files Browse the repository at this point in the history
  • Loading branch information
@jagot
jagot committed May 13, 2024
1 parent 7fd9c94 commit f8fbb3f
Show file tree
Hide file tree
Showing 2 changed files with 58 additions and 59 deletions.
71 changes: 49 additions & 22 deletions SCIDWrapper/src/dipole_analysis.jl
View file
Original file line number Diff line number Diff line change
@@ -1,19 +1,29 @@
plot_dipole_moment(::Nothing) = nothing

function plot_dipole_moment(results)
t = results.t
function plot_dipole_moment(t, Fv, Av, z, vz=nothing, az=nothing; kwargs...)
au2fs = auconvert(u"fs", 1)
tplot = au2fs*t

pE = plot(tplot, Fv, ylabel=L"$E(t)$ [au]")
pA = plot(tplot, Av, ylabel=L"$A(t)$ [au]")

pz = plot(tplot, z, label=L"$z$")
isnothing(vz) || plot!(pz, tplot, vz, label=L"$v_z$")
isnothing(az) || plot!(pz, tplot, az, label=L"$a_z$")

plot(pE, pA, pz;
xlabel=L"$t$ [fs]", layout=@layout([a;b;c]),
size=(700,800), kwargs...)
end

function plot_dipole_moment(results; kwargs...)
t = results.t

z = real(results.z)
vz = real(results.vz)
az = real(results.az)

pE = plot(tplot, results.Ez, ylabel=L"$E(t)$ [au]")
pA = plot(tplot, results.vp_mag, ylabel=L"$A(t)$ [au]")
pz = plot(tplot, [z vz az], label=[L"$z$" L"$v_z$" L"$a_z$"])
plot(pE, pA, pz, xlabel=L"$t$ [fs]", layout=@layout([a;b;c]),
size=(700,800))
plot_dipole_moment(t, results.Ez, results.vp_mag, z, vz, az; kwargs...)
end

function ωaxis(ωkind, ωunit, Iₚ, ω₀, Uₚ)
Expand All @@ -39,41 +49,58 @@ end
round_quantity(v::Quantity) = ElectricFields.si_round(v)
round_quantity(v) = format("{1:.4f}", v)

plot_dipole_spectrum(::Nothing; kwargs...) = nothing
function plot_dipole_spectrum(results; ωkind=:total, ωunit=u"eV", window=hanning)
ip = results.input_params
δt = ip.dt*ip.detail_frequency
ω₀ = ip.vp_param[1][2][1]
function plot_dipole_spectrum(t,
δt, ω₀, Iₚ, Uₚ, cutoff,
z, vz=nothing, az=nothing;
ωkind=:total, ωunit=u"eV", window=hanning,
kwargs...)

t = results.t
ω = 2π*rfftfreq(length(t), 1/δt)

w = window(length(t))

Z = -ω.^2 .* rfft(w .* real(results.z))
Vz = im*ω .* rfft(w .* real(results.vz))
Az = rfft(w .* real(results.az))
Z = -ω.^2 .* rfft(w .* real(z))
Vz = isnothing(vz) ? nothing : im*ω .* rfft(w .* real(vz))
Az = isnothing(az) ? nothing : rfft(w .* real(az))

# Hack to please Plots.jl that refuses to deal with non-positive
# values for logarithmic axes.
nan_map(v) = v ≤ 0 ? NaN : v

ωf,ωlabel,ωax = ωaxis(ωkind, ωunit, results.Iₚ, ω₀, results.Uₚ)
ωf,ωlabel,ωax = ωaxis(ωkind, ωunit, Iₚ, ω₀, Uₚ)

pz = plot(ωf(ω), [nan_map.(abs.(Z)) nan_map.(abs.(Vz)) nan_map.(abs.(Az))],
pz = plot(ωf(ω), nan_map.(abs.(Z));
xlabel=ωlabel,
ylabel="Dipole spectrum [arb.u.]",
label=[L"$|-\omega^2 Z|$" L"$|-\mathrm{i}\omega V_z|$" L"$|A_z|$"],
label=L"$|-\omega^2 Z|$",
xaxis=ωax, yaxis=:log10,
size=(700,600))
size=(700,600),
kwargs...)

isnothing(vz) || plot!(pz, ωf(ω), nan_map.(abs.(Vz)),
label=L"$|-\mathrm{i}\omega V_z|$")
isnothing(az) || plot!(pz, ωf(ω), nan_map.(abs.(Az)),
label=L"$|A_z|$")

Iₚ = ωf(results.Iₚ)
Iₚ = ωf(Iₚ)
vline!(pz, ([Iₚ]), label="Ionization threshold $(round_quantity(Iₚ))")
cutoff = ωf(results.cutoff)
cutoff = ωf(cutoff)
vline!(pz, ([cutoff]), label="HHG cut-off $(round_quantity(cutoff))")
pz
end

plot_dipole_spectrum(::Nothing; kwargs...) = nothing
function plot_dipole_spectrum(results; kwargs...)
ip = results.input_params
δt = ip.dt*ip.detail_frequency
ω₀ = ip.vp_param[1][2][1]

plot_dipole_spectrum(results.t, δt, ω₀,
results.Iₚ, results.Uₚ, results.cutoff,
results.z, results.vz, results.az;
kwargs...)
end

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)
Expand Down
46 changes: 9 additions & 37 deletions scid-notebook.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,20 +262,10 @@ Since we solved the TDSE for a linearly polarized field, only the

# ╔═╡ 0947d3b4-2c76-4de4-a6cf-f91ae840c4ed
let
ptdse = plot(SCIDWrapper.plot_dipole_moment(results), title="TDSE")

au2fs = auconvert(u"fs", 1)
tplot = au2fs*sfa_t
ptdse = SCIDWrapper.plot_dipole_moment(results, title="TDSE")

Av = vector_potential(F, sfa_t)

pF = plot(tplot, Fv, label=L"F(t)")
pA = plot(tplot, Av, label=L"A(t)")
pd = plot(tplot, -sfa_d, label=L"-d(t)")

psfa = plot(pF, pA, pd, xlabel=L"$t$ [fs]",
layout=@layout([a;b;c]),
title="SFA")
psfa = SCIDWrapper.plot_dipole_moment(sfa_t, Fv, Av, -sfa_d, title="SFA")

plot(ptdse, psfa, size=(900,800))
end
Expand All @@ -298,36 +288,18 @@ end

# ╔═╡ 98dce21a-7df0-4059-857c-2777c6e145e0
let
ptdse = plot(SCIDWrapper.plot_dipole_spectrum(results, ωkind=energy_kind, ωunit=freq_unit,
window=apodizing_window), title="TDSE")
kw = (;ωkind=energy_kind, ωunit=freq_unit, window=apodizing_window)
ptdse = SCIDWrapper.plot_dipole_spectrum(results; kw..., title="TDSE")

xl = xlims(ptdse)

ω₀ = austrip(photon_energy(F))
δt = step(sfa_t)

ω = 2π*rfftfreq(length(sfa_t), 1/step(sfa_t))
w = apodizing_window(length(sfa_t))
Fs = rfft(w .* Fv)
Ds = rfft(w .* sfa_d)

Z = -ω.^2 .* Ds

nan_map(v) = v ≤ 0 ? NaN : v

ωf,ωlabel,ωax = SCIDWrapper.ωaxis(energy_kind, freq_unit, austrip(sfa_Iₚ), ω₀, austrip(sfa_Uₚ))

# pF = plot(ωf(ω), abs.(Fs), label=L"|\hat{F}|(\omega)")
pd = plot(ωf(ω), nan_map.(abs.(Z)), label=L"|-\omega^2\hat{d}(\omega)|")

psfa = plot(pd, # xlabel=ωlabel,
xaxis=ωax, yaxis=:log10,
ylabel="Dipole spectrum [arb.u.]",
title="SFA", xlims=xl)

Iₚ = ωf(austrip(sfa_Iₚ))
vline!(psfa, ([Iₚ]), label="Ionization threshold $(SCIDWrapper.round_quantity(Iₚ))")
cutoff = ωf(austrip(sfa_hhg_cutoff))
vline!(psfa, ([cutoff]), label="HHG cut-off $(SCIDWrapper.round_quantity(cutoff))")
psfa = SCIDWrapper.plot_dipole_spectrum(sfa_t, δt, ω₀,
austrip(sfa_Iₚ), austrip(sfa_Uₚ), austrip(sfa_hhg_cutoff),
-sfa_d;
kw..., xlims=xl, title="SFA")

plot(ptdse, psfa, size=(900,900),
layout=@layout([a;b]))
Expand Down

0 comments on commit f8fbb3f

Please sign in to comment.