From d046c33915575db5ff0dfd48ca93f0ab916e2b1e Mon Sep 17 00:00:00 2001 From: SabbahMohammed Date: Tue, 10 Mar 2026 15:30:18 +0000 Subject: [PATCH 1/2] preionfrac as a function of z --- src/Nonlinear.jl | 67 ++++++++++++++++++++++++++++++++------------- src/NonlinearRHS.jl | 22 +++++++-------- src/Plotting.jl | 44 ++++++++++++++++++++--------- 3 files changed, 90 insertions(+), 43 deletions(-) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index a90e7198..92c60109 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -23,7 +23,7 @@ end "Kerr response for real field" function Kerr_field(γ3) Kerr = let γ3 = γ3 - function Kerr(out, E, ρ) + function Kerr(out, E, ρ; z=0.0) if size(E,2) == 1 KerrScalar!(out, E, ρ*ε_0*γ3) else @@ -38,7 +38,7 @@ function Kerr_field_nothg(γ3, n) E = Array{Float64}(undef, n) hilbert = Maths.plan_hilbert(E) Kerr = let γ3 = γ3, hilbert = hilbert - function Kerr(out, E, ρ) + function Kerr(out, E, ρ; z=0.0) out .+= ρ*3/4*ε_0*γ3.*abs2.(hilbert(E)).*E end end @@ -62,7 +62,7 @@ end "Kerr response for envelope" function Kerr_env(γ3) Kerr = let γ3 = γ3 - function Kerr(out, E, ρ) + function Kerr(out, E, ρ; z=0.0) if size(E,2) == 1 KerrScalarEnv!(out, E, ρ*ε_0*γ3) else @@ -77,7 +77,7 @@ end function Kerr_env_thg(γ3, ω0, t) C = exp.(2im*ω0.*t) Kerr = let γ3 = γ3, C = C - function Kerr(out, E, ρ) + function Kerr(out, E, ρ; z=0.0) @. out += ρ*ε_0*γ3/4*(3*abs2(E) + C*E^2)*E end end @@ -85,7 +85,7 @@ end "Response type for cumtrapz-based plasma polarisation, adapted from: M. Geissler, G. Tempea, A. Scrinzi, M. Schnürer, F. Krausz, and T. Brabec, Physical Review Letters 83, 2930 (1999)." -struct PlasmaCumtrapz{R, EType, tType} +struct PlasmaCumtrapz{R, EType, tType, PType} ratefunc::R # the ionization rate function ionpot::Float64 # the ionization potential (for calculation of ionization loss) rate::tType # buffer to hold the rate @@ -94,15 +94,22 @@ struct PlasmaCumtrapz{R, EType, tType} J::EType # buffer to hold the plasma current P::EType # buffer to hold the plasma polarisation δt::Float64 # the time step - preionfrac::Float64 # the pre-ionisation fraction + preionfrac::PType # the pre-ionisation fraction (Float64 or callable function of z) end """ - PlasmaCumtrapz(t, E, ratefunc, ionpot) + PlasmaCumtrapz(t, E, ratefunc, ionpot; preionfrac=0.0) Construct the Plasma polarisation response for a field on time grid `t` with example electric field like `E`, an ionization rate callable `ratefunc` and ionization potential `ionpot`. + +The `preionfrac` can be either: +- A `Float64` constant value (default: 0.0) +- A callable (function or interpolation) that takes `z` (position) as input and returns the pre-ionization fraction + +For z-dependent pre-ionization, you can use an interpolation object from Interpolations.jl +or a custom function that maps z to the pre-ionization fraction. """ function PlasmaCumtrapz(t, E, ratefunc, ionpot; preionfrac=0.0) rate = similar(t) @@ -110,18 +117,39 @@ function PlasmaCumtrapz(t, E, ratefunc, ionpot; preionfrac=0.0) phase = similar(E) J = similar(E) P = similar(E) - !(0.0 <= preionfrac <= 1.0) && throw(DomainError(preionfrac, "preionfrac must be between 0 and 1")) - if preionfrac > 0.0 - @warn("Using preionfrac > 0.0 is not a well founded physical model. Use only after careful consideration.") + # Validate preionfrac + if preionfrac isa Number + !(0.0 <= preionfrac <= 1.0) && throw(DomainError(preionfrac, "preionfrac must be between 0 and 1")) + if preionfrac > 0.0 + @warn("Using preionfrac > 0.0 is not a well founded physical model. Use only after careful consideration.") + end + else # preionfrac is a callable + @warn("Using z-dependent preionfrac. Ensure the function returns values between 0 and 1.") end return PlasmaCumtrapz(ratefunc, ionpot, rate, fraction, phase, J, P, t[2]-t[1], preionfrac) end +""" + getpreionfrac(Plas::PlasmaCumtrapz, z) + +Evaluate the pre-ionization fraction at position `z`. +If `preionfrac` is a constant, returns that constant. +If `preionfrac` is a callable, evaluates it at `z`. +""" +function getpreionfrac(Plas::PlasmaCumtrapz, z) + if Plas.preionfrac isa Number + return Plas.preionfrac + else + return Plas.preionfrac(z) + end +end + "The plasma response for a scalar electric field" -function PlasmaScalar!(Plas::PlasmaCumtrapz, E) +function PlasmaScalar!(Plas::PlasmaCumtrapz, E, z=0.0) Plas.ratefunc(Plas.rate, E) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) - @. Plas.fraction = Plas.preionfrac + 1 - exp(-Plas.fraction) + preionfrac_z = getpreionfrac(Plas, z) + @. Plas.fraction = preionfrac_z + 1 - exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt) for ii in eachindex(E) @@ -141,13 +169,14 @@ for the vector field. A similar approach was used in: C Tailliez et al 2020 New J. Phys. 22 103038. """ -function PlasmaVector!(Plas::PlasmaCumtrapz, E) +function PlasmaVector!(Plas::PlasmaCumtrapz, E, z=0.0) Ex = E[:,1] Ey = E[:,2] Em = @. hypot.(Ex, Ey) Plas.ratefunc(Plas.rate, Em) Maths.cumtrapz!(Plas.fraction, Plas.rate, Plas.δt) - @. Plas.fraction = Plas.preionfrac + 1 - exp(-Plas.fraction) + preionfrac_z = getpreionfrac(Plas, z) + @. Plas.fraction = preionfrac_z + 1 - exp(-Plas.fraction) @. Plas.phase = Plas.fraction * e_ratio * E Maths.cumtrapz!(Plas.J, Plas.phase, Plas.δt) for ii in eachindex(Em) @@ -161,17 +190,17 @@ function PlasmaVector!(Plas::PlasmaCumtrapz, E) end "Handle plasma polarisation routing to `PlasmaVector` or `PlasmaScalar`." -function (Plas::PlasmaCumtrapz)(out, Et, ρ) +function (Plas::PlasmaCumtrapz)(out, Et, ρ; z=0.0) if ndims(Et) > 1 if size(Et, 2) == 1 # handle scalar case but within modal simulation - PlasmaScalar!(Plas, reshape(Et, size(Et,1))) + PlasmaScalar!(Plas, reshape(Et, size(Et,1)), z) out .+= ρ .* reshape(Plas.P, size(Et)) else - PlasmaVector!(Plas, Et) # vector case + PlasmaVector!(Plas, Et, z) # vector case out .+= ρ .* Plas.P end else - PlasmaScalar!(Plas, Et) # straight scalar case + PlasmaScalar!(Plas, Et, z) # straight scalar case out .+= ρ .* Plas.P end end @@ -279,7 +308,7 @@ function sqr!(R::RamanPolarEnv, E) end "Calculate Raman polarisation for field/envelope Et" -function (R::RamanPolar)(out, Et, ρ) +function (R::RamanPolar)(out, Et, ρ; z=0.0) # get the field as a 1D Array n = size(Et, 1) if ndims(Et) > 1 diff --git a/src/NonlinearRHS.jl b/src/NonlinearRHS.jl index c020b333..3bebf883 100644 --- a/src/NonlinearRHS.jl +++ b/src/NonlinearRHS.jl @@ -114,29 +114,29 @@ function _cpscb_core(dest, source, N, scale, idcs) end """ - Et_to_Pt!(Pt, Et, responses, density) + Et_to_Pt!(Pt, Et, responses, density; z=0.0) Accumulate responses induced by Et in Pt. """ -function Et_to_Pt!(Pt, Et, responses, density::Number) +function Et_to_Pt!(Pt, Et, responses, density::Number; z=0.0) fill!(Pt, 0) for resp! in responses - resp!(Pt, Et, density) + resp!(Pt, Et, density; z=z) end end -function Et_to_Pt!(Pt, Et, responses, density::AbstractVector) +function Et_to_Pt!(Pt, Et, responses, density::AbstractVector; z=0.0) fill!(Pt, 0) for ii in eachindex(density) for resp! in responses[ii] - resp!(Pt, Et, density[ii]) + resp!(Pt, Et, density[ii]; z=z) end end end -function Et_to_Pt!(Pt, Et, responses, density, idcs) +function Et_to_Pt!(Pt, Et, responses, density, idcs; z=0.0) for i in idcs - Et_to_Pt!(view(Pt, :, i), view(Et, :, i), responses, density) + Et_to_Pt!(view(Pt, :, i), view(Et, :, i), responses, density; z=z) end end @@ -274,7 +274,7 @@ function Erω_to_Prω!(t, x) Modes.to_space!(t.Erω, t.Emω, x, t.ts, z=t.z) to_time!(t.Er, t.Erω, t.Erωo, inv(t.FT)) # get nonlinear pol at r,θ - Et_to_Pt!(t.Pr, t.Er, t.resp, t.density) + Et_to_Pt!(t.Pr, t.Er, t.resp, t.density; z=t.z) @. t.Pr *= t.grid.towin to_freq!(t.Prω, t.Prωo, t.Pr, t.FT) @. t.Prω *= t.grid.ωwin @@ -362,7 +362,7 @@ const nlscale = sqrt(PhysData.ε_0*PhysData.c/2) function (t::TransModeAvg)(nl, Eω, z) to_time!(t.Eto, Eω, t.Eωo, inv(t.FT)) @. t.Eto /= nlscale*sqrt(t.aeff(z)) - Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z)) + Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z); z=z) @. t.Pto *= t.grid.towin to_freq!(nl, t.Pωo, t.Pto, t.FT) t.norm!(nl, z) @@ -466,7 +466,7 @@ place the result in `nl` function (t::TransRadial)(nl, Eω, z) to_time!(t.Eto, Eω, t.Eωo, inv(t.FT)) # transform ω -> t ldiv!(t.Eto, t.QDHT, t.Eto) # transform k -> r - Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z), t.idcs) # add up responses + Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z), t.idcs; z=z) # add up responses @. t.Pto *= t.grid.towin # apodisation mul!(t.Pto, t.QDHT, t.Pto) # transform r -> k to_freq!(nl, t.Pωo, t.Pto, t.FT) # transform t -> ω @@ -597,7 +597,7 @@ function (t::TransFree)(nl, Eωk, z) fill!(t.Eωo, 0) copy_scale!(t.Eωo, Eωk, length(t.grid.ω), t.scale) ldiv!(t.Eto, t.FT, t.Eωo) # transform (ω, ky, kx) -> (t, y, x) - Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z), t.idcs) # add up responses + Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z), t.idcs; z=z) # add up responses @. t.Pto *= t.grid.towin # apodisation mul!(t.Pωo, t.FT, t.Pto) # transform (t, y, x) -> (ω, ky, kx) copy_scale!(nl, t.Pωo, length(t.grid.ω), 1/t.scale) diff --git a/src/Plotting.jl b/src/Plotting.jl index 0fa9763f..6aa6443a 100644 --- a/src/Plotting.jl +++ b/src/Plotting.jl @@ -118,7 +118,7 @@ end Plot all statistics available in `output`. Additional `kwargs` are passed onto `plt.plot()` """ -function stats(output; kwargs...) +function stats(output; savename=nothing, savepath=pwd(), kwargs...) stats = output["stats"] pstats = [] # pulse statistics @@ -160,13 +160,6 @@ function stats(output; kwargs...) haskey(stats, "dz") && push!(fstats, (1e6*stats["dz"], "Stepsize (μm)")) haskey(stats, "core_radius") && push!(fstats, (1e6*stats["core_radius"], "Core radius (μm)")) haskey(stats, "zdw") && push!(fstats, (1e9*stats["zdw"], "ZDW (nm)")) - haskey(stats, "mode_reconstruction_error") && push!( - fstats, (stats["mode_reconstruction_error"], "Mode error")) - haskey(stats, "transverse_points") && push!( - fstats, (stats["transverse_points"], "Transverse grid points")) - haskey(stats, "transverse_integral_error_rel") && push!( - fstats, (stats["transverse_integral_error_rel"], "Transverse integral error (relative)")) - z = stats["z"]*1e2 @@ -187,6 +180,11 @@ function stats(output; kwargs...) end pfig.tight_layout() end + if savename !== nothing + savename1 = savename*"_1.png" + plt.savefig(joinpath(savepath, savename1)) + plt.close(pfig) + end Npl = length(fstats) if Npl > 0 @@ -203,7 +201,12 @@ function stats(output; kwargs...) end ffig.tight_layout() end - [pfig, ffig] + + if savename !== nothing + savename2 = savename*"_2.png" + plt.savefig(joinpath(savepath, savename2)) + plt.close(ffig) + end end """ @@ -240,7 +243,7 @@ the sum of all modes. function prop_2D(output, specaxis=:f; trange=(-50e-15, 50e-15), bandpass=nothing, λrange=(150e-9, 2000e-9), dBmin=-60, - resolution=nothing, modes=nothing, oversampling=4, + resolution=nothing, modes=nothing, oversampling=4, savename=nothing, savepath=pwd(), kwargs...) z = output["z"]*1e2 if specaxis == :λ @@ -267,6 +270,9 @@ function prop_2D(output, specaxis=:f; kwargs...) end fig + if savename !== nothing + plt.savefig(joinpath(savepath, savename)) + end end modeidcs(m::Int, ml) = [m] @@ -379,7 +385,7 @@ Other `kwargs` are passed onto `plt.plot`. function time_1D(output, zslice=maximum(output["z"]); y=:Pt, modes=nothing, oversampling=4, trange=(-50e-15, 50e-15), bandpass=nothing, - FTL=false, propagate=nothing, + FTL=false, propagate=nothing, savename=nothing, savepath=pwd(), kwargs...) t, Et, zactual = getEt(output, zslice, trange=trange, oversampling=oversampling, bandpass=bandpass, @@ -430,6 +436,10 @@ function time_1D(output, zslice=maximum(output["z"]); sfig.set_size_inches(8.5, 5) sfig.tight_layout() sfig + if savename !== nothing + plt.savefig(joinpath(savepath, savename)) + plt.close(sfig) + end end # Automatically find power unit depending on scale of electric field. @@ -461,7 +471,7 @@ Other `kwargs` are passed onto `plt.plot`. """ function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ; modes=nothing, λrange=(150e-9, 1200e-9), - log10=true, log10min=1e-6, resolution=nothing, + log10=true, log10min=1e-6, resolution=nothing, savename=nothing, savepath=pwd(), kwargs...) if specaxis == :λ specx, Iω, zactual = getIω(output, specaxis, zslice, specrange=λrange, resolution=resolution) @@ -504,6 +514,10 @@ function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ; sfig.set_size_inches(8.5, 5) sfig.tight_layout() sfig + if savename !== nothing + plt.savefig(joinpath(savepath, savename)) + plt.close(sfig) + end end dashes = [(0, (10, 1)), @@ -541,7 +555,7 @@ function spectrogram(grid::Grid.AbstractGrid, output, zslice, specaxis=:λ; end function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ; - trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40, + trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40, savename=nothing, savepath=pwd(), kwargs...) ω = Maths.rfftfreq(t)[2:end] tmin, tmax = extrema(trange) @@ -562,6 +576,10 @@ function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ; log && plt.clim(dBmin, 0) plt.colorbar() fig + if savename !== nothing + plt.savefig(joinpath(savepath, savename)) + plt.close(fig) + end end function energy(output; modes=nothing, bandpass=nothing, figsize=(7, 5)) From cc9af2818551834399a7c2cb627733c87fe0d2e9 Mon Sep 17 00:00:00 2001 From: SabbahMohammed Date: Wed, 11 Mar 2026 10:00:27 +0000 Subject: [PATCH 2/2] revert Plotting.jl changes from preionfrac commit --- src/Plotting.jl | 44 +++++++++++++------------------------------- 1 file changed, 13 insertions(+), 31 deletions(-) diff --git a/src/Plotting.jl b/src/Plotting.jl index 6aa6443a..0fa9763f 100644 --- a/src/Plotting.jl +++ b/src/Plotting.jl @@ -118,7 +118,7 @@ end Plot all statistics available in `output`. Additional `kwargs` are passed onto `plt.plot()` """ -function stats(output; savename=nothing, savepath=pwd(), kwargs...) +function stats(output; kwargs...) stats = output["stats"] pstats = [] # pulse statistics @@ -160,6 +160,13 @@ function stats(output; savename=nothing, savepath=pwd(), kwargs...) haskey(stats, "dz") && push!(fstats, (1e6*stats["dz"], "Stepsize (μm)")) haskey(stats, "core_radius") && push!(fstats, (1e6*stats["core_radius"], "Core radius (μm)")) haskey(stats, "zdw") && push!(fstats, (1e9*stats["zdw"], "ZDW (nm)")) + haskey(stats, "mode_reconstruction_error") && push!( + fstats, (stats["mode_reconstruction_error"], "Mode error")) + haskey(stats, "transverse_points") && push!( + fstats, (stats["transverse_points"], "Transverse grid points")) + haskey(stats, "transverse_integral_error_rel") && push!( + fstats, (stats["transverse_integral_error_rel"], "Transverse integral error (relative)")) + z = stats["z"]*1e2 @@ -180,11 +187,6 @@ function stats(output; savename=nothing, savepath=pwd(), kwargs...) end pfig.tight_layout() end - if savename !== nothing - savename1 = savename*"_1.png" - plt.savefig(joinpath(savepath, savename1)) - plt.close(pfig) - end Npl = length(fstats) if Npl > 0 @@ -201,12 +203,7 @@ function stats(output; savename=nothing, savepath=pwd(), kwargs...) end ffig.tight_layout() end - - if savename !== nothing - savename2 = savename*"_2.png" - plt.savefig(joinpath(savepath, savename2)) - plt.close(ffig) - end + [pfig, ffig] end """ @@ -243,7 +240,7 @@ the sum of all modes. function prop_2D(output, specaxis=:f; trange=(-50e-15, 50e-15), bandpass=nothing, λrange=(150e-9, 2000e-9), dBmin=-60, - resolution=nothing, modes=nothing, oversampling=4, savename=nothing, savepath=pwd(), + resolution=nothing, modes=nothing, oversampling=4, kwargs...) z = output["z"]*1e2 if specaxis == :λ @@ -270,9 +267,6 @@ function prop_2D(output, specaxis=:f; kwargs...) end fig - if savename !== nothing - plt.savefig(joinpath(savepath, savename)) - end end modeidcs(m::Int, ml) = [m] @@ -385,7 +379,7 @@ Other `kwargs` are passed onto `plt.plot`. function time_1D(output, zslice=maximum(output["z"]); y=:Pt, modes=nothing, oversampling=4, trange=(-50e-15, 50e-15), bandpass=nothing, - FTL=false, propagate=nothing, savename=nothing, savepath=pwd(), + FTL=false, propagate=nothing, kwargs...) t, Et, zactual = getEt(output, zslice, trange=trange, oversampling=oversampling, bandpass=bandpass, @@ -436,10 +430,6 @@ function time_1D(output, zslice=maximum(output["z"]); sfig.set_size_inches(8.5, 5) sfig.tight_layout() sfig - if savename !== nothing - plt.savefig(joinpath(savepath, savename)) - plt.close(sfig) - end end # Automatically find power unit depending on scale of electric field. @@ -471,7 +461,7 @@ Other `kwargs` are passed onto `plt.plot`. """ function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ; modes=nothing, λrange=(150e-9, 1200e-9), - log10=true, log10min=1e-6, resolution=nothing, savename=nothing, savepath=pwd(), + log10=true, log10min=1e-6, resolution=nothing, kwargs...) if specaxis == :λ specx, Iω, zactual = getIω(output, specaxis, zslice, specrange=λrange, resolution=resolution) @@ -514,10 +504,6 @@ function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ; sfig.set_size_inches(8.5, 5) sfig.tight_layout() sfig - if savename !== nothing - plt.savefig(joinpath(savepath, savename)) - plt.close(sfig) - end end dashes = [(0, (10, 1)), @@ -555,7 +541,7 @@ function spectrogram(grid::Grid.AbstractGrid, output, zslice, specaxis=:λ; end function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ; - trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40, savename=nothing, savepath=pwd(), + trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40, kwargs...) ω = Maths.rfftfreq(t)[2:end] tmin, tmax = extrema(trange) @@ -576,10 +562,6 @@ function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ; log && plt.clim(dBmin, 0) plt.colorbar() fig - if savename !== nothing - plt.savefig(joinpath(savepath, savename)) - plt.close(fig) - end end function energy(output; modes=nothing, bandpass=nothing, figsize=(7, 5))