Skip to content

preionfrac as a function of z #417

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
SabbahMohammed wants to merge 2 commits into
base: master
Choose a base branch
from
Open

preionfrac as a function of z #417

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d046c33
preionfrac as a function of z
Mar 10, 2026
cc9af28
revert Plotting.jl changes from preionfrac commit
Mar 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 48 additions & 19 deletions src/Nonlinear.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -77,15 +77,15 @@ 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
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
Expand All @@ -94,34 +94,62 @@ 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)
Comment on lines 114 to 115
Copy link

Copilot AI Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This change introduces new behavior where preionfrac may be a z-dependent callable. There are existing tests covering constant preionfrac domain validation (see test/test_ionisation.jl), but there is no coverage for (1) callable preionfrac being accepted and applied at different z values, and (2) callable returns outside [0,1] producing a clear error. Adding tests for these cases would help prevent regressions.

Copilot uses AI. Check for mistakes.
fraction = similar(t)
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
Copy link

Copilot AI Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the callable preionfrac branch, the constructor only emits a warning but does not validate that preionfrac is actually callable with a single z argument. Passing a non-Number non-callable value will fail later at runtime inside getpreionfrac/Plasma* with a MethodError. Consider adding an applicable(preionfrac, 0.0) (or similar) check here and throwing a clear ArgumentError if it isn’t callable.

Suggested change
else # preionfrac is a callable
else # preionfrac is intended to be a callable
if !applicable(preionfrac, 0.0)
throw(ArgumentError("preionfrac must be either a Number in [0, 1] or a callable that accepts a single Float64 argument (z)."))
end

Copilot uses AI. Check for mistakes.
@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)
Copy link

Copilot AI Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

getpreionfrac returns Plas.preionfrac(z) without validating the result. Unlike the constant Number path, this means out-of-range values (<0 or >1) (or non-numeric returns) will silently propagate and can produce unphysical ionisation fractions and downstream computations (e.g. (1 - Plas.fraction[ii]) becoming negative). Consider validating the evaluated value (and throwing DomainError like the constant case) at least once per call site (e.g. inside getpreionfrac).

Suggested change
return Plas.preionfrac(z)
val = Plas.preionfrac(z)
if !(val isa Number)
throw(DomainError(val, "preionfrac(z) must return a numeric value between 0 and 1"))
end
!(0.0 <= val <= 1.0) && throw(DomainError(val, "preionfrac must be between 0 and 1"))
return val

Copilot uses AI. Check for mistakes.
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)
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
22 changes: 11 additions & 11 deletions src/NonlinearRHS.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Comment on lines +121 to 125
Copy link

Copilot AI Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Et_to_Pt! now unconditionally calls each response as resp!(...; z=z). This is a breaking change for any externally-supplied response that previously only implemented (out, E, ρ) (it will now throw a MethodError due to an unexpected keyword). Consider either (a) documenting/enforcing the new response signature at the public API boundary (e.g. wrapping user responses on input so they accept ; z), or (b) adding a fast fallback path that calls without z when the method doesn’t accept that keyword, to preserve backward compatibility.

Copilot uses AI. Check for mistakes.
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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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 -> ω
Expand Down Expand Up @@ -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)
Expand Down
Loading