From b58e0f40aecc74ab5ed275a800da74919e0bb448 Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 22 Apr 2026 13:57:20 +0100 Subject: [PATCH 01/13] Adding data for Ne, Ar, and Kr from Bideau-Mehu. --- src/PhysData.jl | 47 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/src/PhysData.jl b/src/PhysData.jl index f05372f4..bb750c9d 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -47,15 +47,18 @@ const m_u = ustrip(CODATA2014.m_u) "Atomic unit of electric polarisability" const au_polarisability = electron^2*ustrip(CODATA2014.a_0)^2/au_energy -const gas = (:Air, :He, :HeJ, :HeB, :Ne, :Ar, :ArB, :Kr, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) +const gas = (:Air, :He, :HeJ, :HeB, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrBideauMehu, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) const gas_str = Dict( :He => "He", :HeB => "He", :HeJ => "He", :Ar => "Ar", :ArB => "Ar", + :ArBideauMehu => "Ar", :Ne => "Neon", + :NeBideauMehu => "Neon", :Kr => "Krypton", + :KrBideauMehu => "Krypton", :Xe => "Xenon", :Air => "Air", :N2 => "Nitrogen", @@ -138,6 +141,26 @@ function γ_QuanfuHe(A, B, C, dens) return μm -> complex((1 + 1e-8*(A + B/(C - (1e4/μm)^2))))/dens end +""" + γ_BideauMehu2(B1, C1, B2, C2, dens) + +2-term Sellmeier expression for Ne from Journal of Quantitative Spectroscopy and Radiative Transfer, volume 25, issue 5, May 1981, pages 395-402 +https://doi.org/10.1016/0022-4073(81)90057-1 +""" +function γ_BideauMehu2(B1, C1, B2, C2, dens) + return μm -> (((B1/(C1-1/μm^2) + B2/(C2-1/μm^2))+1)^2 - 1)/dens +end + +""" + γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) + +3-term Sellmeier expression for Ar and Kr from Journal of Quantitative Spectroscopy and Radiative Transfer, volume 25, issue 5, May 1981, pages 395-402 +https://doi.org/10.1016/0022-4073(81)90057-1 +""" +function γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) + return μm -> (((B1/(C1-1/μm^2) + B2/(C2-1/μm^2) + B3/(C3-1/μm^2))+1)^2 - 1)/dens +end + """ sellmeier_gas(material::Symbol) @@ -166,6 +189,12 @@ function sellmeier_gas(material::Symbol) B2 = 4018.63e-8 C2 = 5.728e-3 return γ_Börzsönyi(B1/dens, C1, B2/dens, C2) + elseif material == :NeBideauMehu + B1 = 0.00128145 + C1 = 184.6661 + B2 = 0.0220486 + C2 = 376.84 + return γ_BideauMehu2(B1, C1, B2, C2, dens) elseif material == :Ar B1 = 0.00032323117217767093 C1 = 0.0045416501944977915 @@ -180,12 +209,28 @@ function sellmeier_gas(material::Symbol) B2 = 34458.31e-8 C2 = 8.066e-3 return γ_Börzsönyi(B1/dens, C1, B2/dens, C2) + elseif material == :ArBideauMehu + B1 = 2.50141e-3 + C1 = 91.012 + B2 = 5.00283e-4 + C2 = 87.892 + B3 = 5.22343e-2 + C3 = 214.02 + return γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) elseif material == :Kr B1 = 26102.88e-8 C1 = 2.01e-6 B2 = 56946.82e-8 C2 = 10.043e-3 return γ_Börzsönyi(B1/dens, C1, B2/dens, C2) + elseif material == :KrBideauMehu + B1 = 0.00253637 + C1 = 65.4742 + B2 = 0.00273649 + C2 = 73.698 + B3 = 0.0620802 + C3 = 181.08 + return γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) elseif material == :Xe B1 = 103701.61e-8 C1 = 12750e-6 From 6bb88e728d6e497e979ffd69754936c437c98d2b Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 29 Apr 2026 09:43:52 +0100 Subject: [PATCH 02/13] Adding a RDW phase-match finding function including nonlinear and ionisation terms. --- src/PhysData.jl | 11 ++- src/Tools.jl | 196 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 206 insertions(+), 1 deletion(-) diff --git a/src/PhysData.jl b/src/PhysData.jl index bb750c9d..024f4e29 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -47,7 +47,7 @@ const m_u = ustrip(CODATA2014.m_u) "Atomic unit of electric polarisability" const au_polarisability = electron^2*ustrip(CODATA2014.a_0)^2/au_energy -const gas = (:Air, :He, :HeJ, :HeB, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrBideauMehu, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) +const gas = (:Air, :He, :HeJ, :HeB, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrB, :KrBideauMehu, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) const gas_str = Dict( :He => "He", :HeB => "He", @@ -58,6 +58,7 @@ const gas_str = Dict( :Ne => "Neon", :NeBideauMehu => "Neon", :Kr => "Krypton", + :KrB => "Krypton", :KrBideauMehu => "Krypton", :Xe => "Xenon", :Air => "Air", @@ -218,6 +219,14 @@ function sellmeier_gas(material::Symbol) C3 = 214.02 return γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) elseif material == :Kr + B1 = 7.156207185088921e-6 + C1 = 0.016126479267946525 + B2 = 0.00023020614937540427 + C2 = 0.01348822574907029 + B3 = 0.000604323109684193 + C3 = 0.004204161467601559 + return γ_JCT(B1/dens, C1, B2/dens, C2, B3/dens, C3) + elseif material == :KrB B1 = 26102.88e-8 C1 = 2.01e-6 B2 = 56946.82e-8 diff --git a/src/Tools.jl b/src/Tools.jl index 04c2e9eb..bd020e17 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -235,6 +235,202 @@ function pressureZDW(a::Number, gas::Symbol, λzd; Pmax=100, clad=:SiO2, kwargs. end end +""" + Calculates the critical electron density for a given wavelength, which is the density at which the plasma frequency equals the frequency of the light + Ref: https://doi.org/10.1016/j.physrep.2006.12.005, top of page 56 + Arguments: + λ: the wavelength of light +""" +function critical_electron_density(λ) + ω = PhysData.wlfreq(λ) + return (PhysData.ε_0*PhysData.m_e*ω^2)/PhysData.electron^2 +end + + +""" + Calculates the energy required for a certain soliton order + Arguments: + a: radius of the capillary core + gas: type of gas in the capillary + pressure: pressure of the gas in the capillary + τFWHM: pulse duration (FWHM) + λp: central wavelength of the pump (soliton) pulse + N: soliton order + energy_lims: limits for the root finding algorithm that looks for the energy; by default these are set to 1 nJ to 1 J +""" +function energyN(a, gas, pressure, τFWHM, λp, N; energy_lims=(1e-9, 1.0)) + + EforN = find_zero(energy_lims) do energy + params = Tools.capillary_params(energy, τFWHM, λp, a, gas; P=pressure) + params.N - N + end + + return EforN + +end + +""" + Helper function to calculate the linear contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission +""" +function Δβlin(a, gas, pressure, λp; kwargs...) + + mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) + ωsol = PhysData.wlfreq(λp) + β1 = Modes.dispersion(mode, 1, ωsol) + β0 = Modes.β(mode, ωsol) + Δβlin(ω) = Modes.β(mode, ω) - β1*(ω - ωsol) - β0 + + return Δβlin + +end + +""" + Helper function to calculate the nonlinear contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission +""" +function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; use_PpeakSelfCompressed=false) + + energy = energyN(a, gas, pressure, τFWHM, λp, soliton_order) + params = Tools.capillary_params(energy, τFWHM, λp, a, gas; P=pressure) + + if use_PpeakSelfCompressed + Ppeak_sc = params.P0*(4.6*soliton_order) + else + Ppeak_sc = params.P0 + end + + ωsol = PhysData.wlfreq(λp) + return params.γ*Ppeak_sc*(1/ωsol) # γ*Pp*(1/ωsol) +end + +""" + Helper function to calculate the ionisation contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission +""" +function Δβion(a, gas, pressure, λp, ionisation_fraction; kwargs...) + neutral_density = PhysData.density(gas, pressure) + electron_density = ionisation_fraction*neutral_density + critical_density = critical_electron_density(λp) + ωsol = PhysData.wlfreq(λp) + mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) + nlin = real.(Capillary.neff(mode, ωsol)) # taking the real part here, otherwise the ionisation term becomes complex + return (1/(2*PhysData.c*nlin))*(electron_density/critical_density)*(ωsol^2) +end + +""" + A helper type to store parameters for phase-matching +""" +struct PhaseMatching + mode + ωsol + β0 + β1 + ΔβnonlinCoeff + ΔβionCoeff +end + +""" + A helper function to calculate the phase-matching curve from pre-calculated parameters before passing to the root-finding algorithm; this speeds up the calculation a lot compared to when just passing Δβ as a full function, that has to be evaluated at every iteration of the root-finding +""" +function (pm::PhaseMatching)(ω) + + return Modes.β(pm.mode, ω) - pm.β1*(ω - pm.ωsol) - pm.β0 - pm.ΔβnonlinCoeff*ω + pm.ΔβionCoeff*(1/ω) + +end + +""" + A helper function to pre-compute parameters for phase-matching +""" +function make_PhaseMatching(a, gas, pressure, λp; + include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, include_Δβion=false, ionisation_fraction=nothing, use_PpeakSelfCompressed=false, kwargs...) + + mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) + ωsol = PhysData.wlfreq(λp) + β1 = Modes.dispersion(mode, 1, ωsol) + β0 = Modes.β(mode, ωsol) + + if include_Δβnonlin # perform the assertions only if the nonlinear contribution is included, otherwise these parameters are not needed and can be left as nothing + @assert !isnothing(soliton_order) "If include_Δβnonlin is true, soliton_order must be given" + @assert soliton_order > 0 "Soliton order must be positive" + @assert !isnothing(τFWHM) "If include_Δβnonlin is true, τFWHM must be given" + end + + if include_Δβion # perform the assertions only if the ionisation contribution is included, otherwise these parameters are not needed and can be left as nothing + @assert !isnothing(ionisation_fraction) "If include_Δβion is true, ionisation_fraction must be given" + @assert ionisation_fraction >= 0 && ionisation_fraction <= 1 "Ionisation fraction must be between 0 and 1" + end + + ΔβnonlinCoeff = include_Δβnonlin ? Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; use_PpeakSelfCompressed) : 0.0 + ΔβionCoeff = include_Δβion ? Δβion(a, gas, pressure, λp, ionisation_fraction; kwargs...) : 0.0 + + return PhaseMatching(mode, ωsol, β0, β1, ΔβnonlinCoeff, ΔβionCoeff) +end + +""" + Function for calculating the phase-matched resonant dispersive wave (RDW) wavelength in a gas-filled capillary, taking into account the nonlinear and ionisation contributions to phase-matching. + Ref: PRL 115, 033901 (2015), https://doi.org/10.1103/PhysRevLett.115.033901 + Arguments: + a: radius of the capillary core + gas: type of gas in the capillary + pressure: pressure of the gas in the capillary + λp: central wavelength of the pump (soliton) pulse + kwargs: additional keyword arguments for the mode (e.g. m, n, ϕ, kind, clad; used in calling MarcatiliMode() line 95 in Luna.Capillary) + soliton_order: used to calculate the energy/peak power + τFWHM: pulse duration (FWHM) + ionisation_fraction: electron fraction + λlims: limits for the root finding algorithm that looks for the phase-matching; by default these are set to the minimum wavelength just above the first resonance in the Sellmeier equation for the given gas and the maximum for the pump wavelength minus 1 nm +""" +function λRDWfull(a, gas, pressure, λp; + include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, include_Δβion=false, ionisation_fraction=nothing, use_PpeakSelfCompressed=false, + λlims=nothing, kwargs...) + + if isnothing(λlims) + λUVlim = Dict( + :He => 90e-9, + :HeB => 90e-9, + :HeJ => 90e-9, + :Ne => 100e-9, + :NeBideauMehu => 100e-9, + :Ar => 110e-9, + :ArB => 110e-9, + :ArBideauMehu => 110e-9, + :Kr => 130e-9, + :KrB => 130e-9, + :KrBideauMehu => 130e-9 + ) + λlims = (λUVlim[gas], λp-10e-9) # pump wavelength minus 10 nm for safety + end + + phase_matching = make_PhaseMatching(a, gas, pressure, λp; include_Δβnonlin=include_Δβnonlin, soliton_order=soliton_order, τFWHM=τFWHM, include_Δβion=include_Δβion, ionisation_fraction=ionisation_fraction, use_PpeakSelfCompressed=use_PpeakSelfCompressed, kwargs...) + + # if ionization effects are included, shrink λlims, because the dispersion curve gets lifted in the IR and a new phase-matched point appears near the pump + if include_Δβion + + # first find the linear phase-matching wavelength + λlin = Tools.λRDW(a, gas, pressure, λp; λlims=λlims) # use the OG Luna function to calculate the linear phase-matched RDW + + # find the root close to the pump wavelength + ωmax = NaN + try + @assert !ismissing(λlin) "λlin should not be NaN here, the first root-finding has failed" + ωmax = find_zero(phase_matching, extrema(PhysData.wlfreq.((λlin+50e-9, λp-10e-9)))) # linear RDW + 50 nm for safety + λlims = (minimum(λlims), PhysData.wlfreq(ωmax) - 10e-9) + catch err + @warn "A phase-matching wavelength couldn't be found for Δβ..." exception=(err, catch_backtrace()) + end + + end + + # find RDW phase-matched wavelength with the full Δβ and potentially shrunk λlims + ωRDW = NaN + try + ωRDW = find_zero(phase_matching, extrema(PhysData.wlfreq.(λlims))) + catch err + @warn "A phase-matching wavelength couldn't be found for the given parameters..." exception=(err, catch_backtrace()) + end + + return PhysData.wlfreq(ωRDW) + +end + function aperture_filter(a, dist, radius) w0 = 0.64*a function filter(λ) From ca094620567b917a90f97a508ff05c5c2d5b3ccd Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 29 Apr 2026 12:14:14 +0100 Subject: [PATCH 03/13] Adding the new gas symbols to the lookup lists for gas nonlineaity and ionisation parameters --- src/PhysData.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/PhysData.jl b/src/PhysData.jl index 024f4e29..5cc2c62d 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -642,7 +642,7 @@ References: function γ3_gas(material::Symbol; source=nothing) # TODO: More Bishop/Shelton; Wahlstrand updated values. if source === nothing - if material in (:He, :HeB, :HeJ, :Ne, :Ar, :ArB, :Kr, :Xe, :N2) + if material in (:He, :HeB, :HeJ, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrB, :KrBideauMehu,:Xe, :N2) source = :Lehmeier elseif material in (:H2, :CH4, :SF6, :D2) source = :Shelton @@ -659,11 +659,11 @@ function γ3_gas(material::Symbol; source=nothing) # Table 1 in [3] if material in (:He, :HeB, :HeJ) fac = 1 - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) fac = 1.8 - elseif material in (:Ar, :ArB) + elseif material in (:Ar, :ArB, :ArBideauMehu) fac = 23.5 - elseif material == :Kr + elseif material in (:Kr, :KrB, :KrBideauMehu) fac = 64.0 elseif material == :Xe fac = 188.2 @@ -785,11 +785,11 @@ Possible units are `:SI`, `:atomic` and `:eV`. function ionisation_potential(material; unit=:SI) if material in (:He, :HeB, :HeJ) Ip = 0.9036 - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) Ip = 0.7925 - elseif material in (:Ar, :ArB) + elseif material in (:Ar, :ArB, :ArBideauMehu) Ip = 0.5792 - elseif material == :Kr + elseif material in (:Kr, :KrB, :KrBideauMehu) Ip = 0.5142 elseif material == :Xe Ip = 0.4458 @@ -831,11 +831,11 @@ Return the quantum numbers of the `material` for use in the PPT ionisation rate. """ function quantum_numbers(material) # Returns n, l, ion Z - if material in (:Ar, :ArB) + if material in (:Ar, :ArB, :ArBideauMehu) return 3, 1, 1 - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) return 2, 1, 1; - elseif material == :Kr + elseif material in (:Kr, :KrB, :KrBideauMehu) return 4, 1, 1 elseif material == :Xe return 5, 1, 1 @@ -876,9 +876,9 @@ function polarisability(material, ion=false; unit=:SI) end if material in (:He, :HeB, :HeJ) return (ion ? 0.2811 : 1.3207)*factor - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) return (ion ? 1.2417 : 2.376)*factor - elseif material in (:Ar, :ArB) + elseif material in (:Ar, :ArB, :ArBideauMehu) return (ion ? 6.807 : 10.762)*factor else return missing @@ -915,11 +915,11 @@ Ammosov, M. V., Delone, N. B. & Krainov, V. P. Tunnel Ionization Of Complex Atom function Cnl_ADK(material) if material in (:He, :HeB, :HeJ) return 1.99 - elseif material == :Ne + elseif material in (:Ne, :NeBideauMehu) return 1.31 - elseif material in (:Ar, :ArB) + elseif material in (:Ar, :ArB, :ArBideauMehu) return 1.9 - elseif material == :Kr + elseif material in (:Kr, :KrB, :KrBideauMehu) return 2.17 elseif material == :Xe return 2.27 From 294af2900e6869342bed13966f9167070b53602f Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 29 Apr 2026 14:31:47 +0100 Subject: [PATCH 04/13] =?UTF-8?q?Krypton=20B=C3=B6rzs=C3=B6nyi=20is=20:Kr?= =?UTF-8?q?=20and=20the=20Kingston+B=C3=B6rzs=C3=B6nyi=20extension=20is=20?= =?UTF-8?q?:KrJCT?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/PhysData.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/PhysData.jl b/src/PhysData.jl index 5cc2c62d..ab770d5d 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -47,7 +47,7 @@ const m_u = ustrip(CODATA2014.m_u) "Atomic unit of electric polarisability" const au_polarisability = electron^2*ustrip(CODATA2014.a_0)^2/au_energy -const gas = (:Air, :He, :HeJ, :HeB, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrB, :KrBideauMehu, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) +const gas = (:Air, :He, :HeJ, :HeB, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrJCT, :KrBideauMehu, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) const gas_str = Dict( :He => "He", :HeB => "He", @@ -58,7 +58,7 @@ const gas_str = Dict( :Ne => "Neon", :NeBideauMehu => "Neon", :Kr => "Krypton", - :KrB => "Krypton", + :KrJCT => "Krypton", :KrBideauMehu => "Krypton", :Xe => "Xenon", :Air => "Air", @@ -218,7 +218,7 @@ function sellmeier_gas(material::Symbol) B3 = 5.22343e-2 C3 = 214.02 return γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) - elseif material == :Kr + elseif material == :KrJCT B1 = 7.156207185088921e-6 C1 = 0.016126479267946525 B2 = 0.00023020614937540427 @@ -226,7 +226,7 @@ function sellmeier_gas(material::Symbol) B3 = 0.000604323109684193 C3 = 0.004204161467601559 return γ_JCT(B1/dens, C1, B2/dens, C2, B3/dens, C3) - elseif material == :KrB + elseif material == :Kr B1 = 26102.88e-8 C1 = 2.01e-6 B2 = 56946.82e-8 @@ -642,7 +642,7 @@ References: function γ3_gas(material::Symbol; source=nothing) # TODO: More Bishop/Shelton; Wahlstrand updated values. if source === nothing - if material in (:He, :HeB, :HeJ, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrB, :KrBideauMehu,:Xe, :N2) + if material in (:He, :HeB, :HeJ, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrJCT, :KrBideauMehu,:Xe, :N2) source = :Lehmeier elseif material in (:H2, :CH4, :SF6, :D2) source = :Shelton @@ -663,7 +663,7 @@ function γ3_gas(material::Symbol; source=nothing) fac = 1.8 elseif material in (:Ar, :ArB, :ArBideauMehu) fac = 23.5 - elseif material in (:Kr, :KrB, :KrBideauMehu) + elseif material in (:Kr, :KrJCT, :KrBideauMehu) fac = 64.0 elseif material == :Xe fac = 188.2 @@ -789,7 +789,7 @@ function ionisation_potential(material; unit=:SI) Ip = 0.7925 elseif material in (:Ar, :ArB, :ArBideauMehu) Ip = 0.5792 - elseif material in (:Kr, :KrB, :KrBideauMehu) + elseif material in (:Kr, :KrJCT, :KrBideauMehu) Ip = 0.5142 elseif material == :Xe Ip = 0.4458 @@ -835,7 +835,7 @@ function quantum_numbers(material) return 3, 1, 1 elseif material in (:Ne, :NeBideauMehu) return 2, 1, 1; - elseif material in (:Kr, :KrB, :KrBideauMehu) + elseif material in (:Kr, :KrJCT, :KrBideauMehu) return 4, 1, 1 elseif material == :Xe return 5, 1, 1 @@ -919,7 +919,7 @@ function Cnl_ADK(material) return 1.31 elseif material in (:Ar, :ArB, :ArBideauMehu) return 1.9 - elseif material in (:Kr, :KrB, :KrBideauMehu) + elseif material in (:Kr, :KrJCT, :KrBideauMehu) return 2.17 elseif material == :Xe return 2.27 From 4c8fa67f1155e231129810cc3cd3de6d86148689 Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 29 Apr 2026 14:45:57 +0100 Subject: [PATCH 05/13] =?UTF-8?q?Going=20back:=20all=20B=C3=B6rzs=C3=B6nyi?= =?UTF-8?q?=20is=20explicitly=20B=20and=20the=20new=20Kingston+B=C3=B6rzs?= =?UTF-8?q?=C3=B6nyi=20is=20the=20default?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/PhysData.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/PhysData.jl b/src/PhysData.jl index ab770d5d..5cc2c62d 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -47,7 +47,7 @@ const m_u = ustrip(CODATA2014.m_u) "Atomic unit of electric polarisability" const au_polarisability = electron^2*ustrip(CODATA2014.a_0)^2/au_energy -const gas = (:Air, :He, :HeJ, :HeB, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrJCT, :KrBideauMehu, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) +const gas = (:Air, :He, :HeJ, :HeB, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrB, :KrBideauMehu, :Xe, :N2, :H2, :O2, :CH4, :SF6, :N2O, :D2) const gas_str = Dict( :He => "He", :HeB => "He", @@ -58,7 +58,7 @@ const gas_str = Dict( :Ne => "Neon", :NeBideauMehu => "Neon", :Kr => "Krypton", - :KrJCT => "Krypton", + :KrB => "Krypton", :KrBideauMehu => "Krypton", :Xe => "Xenon", :Air => "Air", @@ -218,7 +218,7 @@ function sellmeier_gas(material::Symbol) B3 = 5.22343e-2 C3 = 214.02 return γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) - elseif material == :KrJCT + elseif material == :Kr B1 = 7.156207185088921e-6 C1 = 0.016126479267946525 B2 = 0.00023020614937540427 @@ -226,7 +226,7 @@ function sellmeier_gas(material::Symbol) B3 = 0.000604323109684193 C3 = 0.004204161467601559 return γ_JCT(B1/dens, C1, B2/dens, C2, B3/dens, C3) - elseif material == :Kr + elseif material == :KrB B1 = 26102.88e-8 C1 = 2.01e-6 B2 = 56946.82e-8 @@ -642,7 +642,7 @@ References: function γ3_gas(material::Symbol; source=nothing) # TODO: More Bishop/Shelton; Wahlstrand updated values. if source === nothing - if material in (:He, :HeB, :HeJ, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrJCT, :KrBideauMehu,:Xe, :N2) + if material in (:He, :HeB, :HeJ, :Ne, :NeBideauMehu, :Ar, :ArB, :ArBideauMehu, :Kr, :KrB, :KrBideauMehu,:Xe, :N2) source = :Lehmeier elseif material in (:H2, :CH4, :SF6, :D2) source = :Shelton @@ -663,7 +663,7 @@ function γ3_gas(material::Symbol; source=nothing) fac = 1.8 elseif material in (:Ar, :ArB, :ArBideauMehu) fac = 23.5 - elseif material in (:Kr, :KrJCT, :KrBideauMehu) + elseif material in (:Kr, :KrB, :KrBideauMehu) fac = 64.0 elseif material == :Xe fac = 188.2 @@ -789,7 +789,7 @@ function ionisation_potential(material; unit=:SI) Ip = 0.7925 elseif material in (:Ar, :ArB, :ArBideauMehu) Ip = 0.5792 - elseif material in (:Kr, :KrJCT, :KrBideauMehu) + elseif material in (:Kr, :KrB, :KrBideauMehu) Ip = 0.5142 elseif material == :Xe Ip = 0.4458 @@ -835,7 +835,7 @@ function quantum_numbers(material) return 3, 1, 1 elseif material in (:Ne, :NeBideauMehu) return 2, 1, 1; - elseif material in (:Kr, :KrJCT, :KrBideauMehu) + elseif material in (:Kr, :KrB, :KrBideauMehu) return 4, 1, 1 elseif material == :Xe return 5, 1, 1 @@ -919,7 +919,7 @@ function Cnl_ADK(material) return 1.31 elseif material in (:Ar, :ArB, :ArBideauMehu) return 1.9 - elseif material in (:Kr, :KrJCT, :KrBideauMehu) + elseif material in (:Kr, :KrB, :KrBideauMehu) return 2.17 elseif material == :Xe return 2.27 From 0f24db1ace4b6af0524a7c889858ee78efce7372 Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 29 Apr 2026 14:52:28 +0100 Subject: [PATCH 06/13] Tiny change to run the tests again (switched back from draft after the last commit and they won't run) --- src/PhysData.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/PhysData.jl b/src/PhysData.jl index 5cc2c62d..90a2a2fa 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -145,8 +145,7 @@ end """ γ_BideauMehu2(B1, C1, B2, C2, dens) -2-term Sellmeier expression for Ne from Journal of Quantitative Spectroscopy and Radiative Transfer, volume 25, issue 5, May 1981, pages 395-402 -https://doi.org/10.1016/0022-4073(81)90057-1 +2-term Sellmeier expression for Ne from https://doi.org/10.1016/0022-4073(81)90057-1 """ function γ_BideauMehu2(B1, C1, B2, C2, dens) return μm -> (((B1/(C1-1/μm^2) + B2/(C2-1/μm^2))+1)^2 - 1)/dens @@ -155,8 +154,7 @@ end """ γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) -3-term Sellmeier expression for Ar and Kr from Journal of Quantitative Spectroscopy and Radiative Transfer, volume 25, issue 5, May 1981, pages 395-402 -https://doi.org/10.1016/0022-4073(81)90057-1 +3-term Sellmeier expression for Ar and Kr from https://doi.org/10.1016/0022-4073(81)90057-1 """ function γ_BideauMehu3(B1, C1, B2, C2, B3, C3, dens) return μm -> (((B1/(C1-1/μm^2) + B2/(C2-1/μm^2) + B3/(C3-1/μm^2))+1)^2 - 1)/dens From d14b111a82493a810b219d504d0f71720d9baa25 Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 29 Apr 2026 16:30:22 +0100 Subject: [PATCH 07/13] =?UTF-8?q?Failing=20tests=20because=20krypton=20B?= =?UTF-8?q?=C3=B6rzs=C3=B6nyi=20is=20now=20:KrB=20instead=20of=20:Kr?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test/test_tools.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_tools.jl b/test/test_tools.jl index 85e5aaf3..ec9db806 100644 --- a/test/test_tools.jl +++ b/test/test_tools.jl @@ -10,7 +10,7 @@ p = Tools.params(300e-6, 10e-15, 800e-9, m, :HeB, P=0.4) @test isapprox(p.zdw, 378.8e-9, rtol=1e-2) @test isapprox(p.P0/p.Pcr, 0.0398, rtol=2e-2) # test backup zdw method -p = Tools.capillary_params(6e-9, 20e-15, 800e-9, 14e-6, :Kr, P=15.0) +p = Tools.capillary_params(6e-9, 20e-15, 800e-9, 14e-6, :KrB, P=15.0) @test isapprox(p.zdw, 7.693023014958748e-7, rtol=1e-7) end From 2cb52ab49b223634aa15fefe2a9b9f21606e7b5c Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 13 May 2026 11:02:47 +0100 Subject: [PATCH 08/13] Added functions for Nmax ion and self-focusing --- src/PhysData.jl | 11 +++++ src/Tools.jl | 110 +++++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 110 insertions(+), 11 deletions(-) diff --git a/src/PhysData.jl b/src/PhysData.jl index 90a2a2fa..12497f95 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -941,6 +941,17 @@ function lookup_glass(material::Symbol) return spl end + +function Ith(material; ionLevel=1) + + Ip = ionisation_potential(material; unit=:eV) + Ith = 3.8e9*(Ip^4)/(ionLevel^2) + + return Ith*1e4 # convert from W/cm^2 to W/m^2 + +end + + """ lookup_metal(material::Symbol) diff --git a/src/Tools.jl b/src/Tools.jl index bd020e17..2f51048a 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -235,6 +235,8 @@ function pressureZDW(a::Number, gas::Symbol, λzd; Pmax=100, clad=:SiO2, kwargs. end end + + """ Calculates the critical electron density for a given wavelength, which is the density at which the plasma frequency equals the frequency of the light Ref: https://doi.org/10.1016/j.physrep.2006.12.005, top of page 56 @@ -284,22 +286,66 @@ function Δβlin(a, gas, pressure, λp; kwargs...) end +""" + Calculate the estimated compression factor and quality factor of the soliton self-compression + Schade et al. "Scaling rules for high quality soliton self-compression in hollow-core fibers" https://doi.org/10.1364/OE.426307 + Arguments: + N: soliton order + Returns compression factor and compression quality (Fc, Q) +""" +function compression_estimation(a, gas, pressure, λp, τFWHM, N; include_ξ=false) + + A(N) = (1/(2*N)+1.7/(N^2)) + Q(N) = 3.7/(N + 2.2) + + ωsol = PhysData.wlfreq(λp) + mode = Capillary.MarcatiliMode(a, gas, pressure) + β2 = Modes.dispersion(mode, 2, ωsol) + β3 = Modes.dispersion(mode, 3, ωsol) + ξ = β3/(τFWHM*abs(β2)) + + Fc = 3/(A(N)) + + if include_ξ + Fc *= (1-N*ξ) + end + + return Fc, Q(N) +end + +function Ppeak(τFWHM, energy; pulse_shape=:sech) + + if pulse_shape == :sech + return 0.88*(energy/τFWHM) + elseif pulse_shape == :gauss + return 0.94*(energy/τFWHM) + else + error("pulse_shape must be one of: :sech, :gauss") + end + +end + """ Helper function to calculate the nonlinear contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission """ -function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; use_PpeakSelfCompressed=false) +function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; input_pulse_shape=:gauss, kwargs...) energy = energyN(a, gas, pressure, τFWHM, λp, soliton_order) params = Tools.capillary_params(energy, τFWHM, λp, a, gas; P=pressure) - - if use_PpeakSelfCompressed - Ppeak_sc = params.P0*(4.6*soliton_order) + ωsol = PhysData.wlfreq(λp) + + if input_pulse_shape == :sech + Pp = Ppeak(τFWHM, energy, pulse_shape=:sech) + elseif input_pulse_shape == :gauss + Pp = 0.936*Ppeak(τFWHM, energy, pulse_shape=:gauss) # 0.88/0.94 else - Ppeak_sc = params.P0 + error("input_pulse_shape must be one of: :sech, :gauss") end - ωsol = PhysData.wlfreq(λp) - return params.γ*Ppeak_sc*(1/ωsol) # γ*Pp*(1/ωsol) + soliton_factor = ((2*soliton_order-1)/soliton_order)^2 # ((2N-1)/N)^2 + + return params.γ*(soliton_factor*Pp)*(1/ωsol) + end """ @@ -340,7 +386,8 @@ end A helper function to pre-compute parameters for phase-matching """ function make_PhaseMatching(a, gas, pressure, λp; - include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, include_Δβion=false, ionisation_fraction=nothing, use_PpeakSelfCompressed=false, kwargs...) + include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, input_pulse_shape=:gauss, + include_Δβion=false, ionisation_fraction=nothing, kwargs...) mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) ωsol = PhysData.wlfreq(λp) @@ -358,7 +405,7 @@ function make_PhaseMatching(a, gas, pressure, λp; @assert ionisation_fraction >= 0 && ionisation_fraction <= 1 "Ionisation fraction must be between 0 and 1" end - ΔβnonlinCoeff = include_Δβnonlin ? Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; use_PpeakSelfCompressed) : 0.0 + ΔβnonlinCoeff = include_Δβnonlin ? Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; input_pulse_shape=input_pulse_shape) : 0.0 ΔβionCoeff = include_Δβion ? Δβion(a, gas, pressure, λp, ionisation_fraction; kwargs...) : 0.0 return PhaseMatching(mode, ωsol, β0, β1, ΔβnonlinCoeff, ΔβionCoeff) @@ -379,7 +426,8 @@ end λlims: limits for the root finding algorithm that looks for the phase-matching; by default these are set to the minimum wavelength just above the first resonance in the Sellmeier equation for the given gas and the maximum for the pump wavelength minus 1 nm """ function λRDWfull(a, gas, pressure, λp; - include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, include_Δβion=false, ionisation_fraction=nothing, use_PpeakSelfCompressed=false, + include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, input_pulse_shape=:gauss, + include_Δβion=false, ionisation_fraction=nothing, λlims=nothing, kwargs...) if isnothing(λlims) @@ -399,7 +447,9 @@ function λRDWfull(a, gas, pressure, λp; λlims = (λUVlim[gas], λp-10e-9) # pump wavelength minus 10 nm for safety end - phase_matching = make_PhaseMatching(a, gas, pressure, λp; include_Δβnonlin=include_Δβnonlin, soliton_order=soliton_order, τFWHM=τFWHM, include_Δβion=include_Δβion, ionisation_fraction=ionisation_fraction, use_PpeakSelfCompressed=use_PpeakSelfCompressed, kwargs...) + phase_matching = make_PhaseMatching(a, gas, pressure, λp; + include_Δβnonlin=include_Δβnonlin, soliton_order=soliton_order, τFWHM=τFWHM, input_pulse_shape=input_pulse_shape, + include_Δβion=include_Δβion, ionisation_fraction=ionisation_fraction, kwargs...) # if ionization effects are included, shrink λlims, because the dispersion curve gets lifted in the IR and a new phase-matched point appears near the pump if include_Δβion @@ -431,6 +481,44 @@ function λRDWfull(a, gas, pressure, λp; end +function fλ(λ, material, pressure; temperature=PhysData.roomtemp) + + χ1function = PhysData.χ1_fun(material, pressure, temperature) + + return Maths.derivative(χ1function, λ, 2) + +end + + +function δ(λ, zdw, material, pressure; kind=:HE, n=1, m=1) + + unm = Capillary.get_unm(n, m, kind) + + return ((unm^2*λ^3)/(8*π^3*PhysData.c^2))*((fλ(λ, material, pressure))/(fλ(zdw, material, pressure)) - 1) + +end + +function NmaxSelfFocusing(λp, zdw, gas, τFWHM, a; safetyFactor=10, pulseShape=:gauss, kind=:HE, n=1, m=1) + + pressure = pressureZDW(a, gas, zdw) + τ0 = τfw_to_τ0(τFWHM, pulseShape) + + return sqrt((τ0^2*λp)/(safetyFactor*abs(δ(λp, zdw, gas, pressure; kind=kind, n=n, m=m)))) + +end + +function NmaxIonisation(λp, zdw, gas, τFWHM, a; safetyFactor=10, pulseShape=:gauss, kind=:HE, n=1, m=1) + + pressure = pressureZDW(a, gas, zdw) + τ0 = τfw_to_τ0(τFWHM, pulseShape) + ωpump = PhysData.wlfreq(λp) + _, _, n2 = getN0n0n2(ωpump, gas; P=pressure) + unm = Capillary.get_unm(n, m, kind) + + return sqrt((τ0^2*n2*PhysData.PhysData.Ith(gas)*unm^2)/(safetyFactor*π*λp*abs(δ(λp, zdw, gas, pressure; kind=kind, n=n, m=m))*fλ(zdw, gas, pressure))) + +end + function aperture_filter(a, dist, radius) w0 = 0.64*a function filter(λ) From 39c9b4b9f01a486935d04fc81b7d1926485ede52 Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 13 May 2026 15:52:49 +0100 Subject: [PATCH 09/13] Add loss up to Lfiss in the Ppeak estimation for the nonlinear phase-matching term --- src/Tools.jl | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/Tools.jl b/src/Tools.jl index 2f51048a..6b64560f 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -328,7 +328,7 @@ end """ Helper function to calculate the nonlinear contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission """ -function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; input_pulse_shape=:gauss, kwargs...) +function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=false, input_pulse_shape=:gauss, kwargs...) energy = energyN(a, gas, pressure, τFWHM, λp, soliton_order) params = Tools.capillary_params(energy, τFWHM, λp, a, gas; P=pressure) @@ -344,7 +344,15 @@ function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; input_pulse_sh soliton_factor = ((2*soliton_order-1)/soliton_order)^2 # ((2N-1)/N)^2 - return params.γ*(soliton_factor*Pp)*(1/ωsol) + attenuation = 1.0 + + if includeLoss + Lfiss = params.Lfiss + mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) + attenuation = α(mode, PhysData.wlfreq(λp); z=Lfiss) + end + + return params.γ*(soliton_factor*Pp*attenuation)*(1/ωsol) end @@ -386,7 +394,7 @@ end A helper function to pre-compute parameters for phase-matching """ function make_PhaseMatching(a, gas, pressure, λp; - include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, input_pulse_shape=:gauss, + include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, includeLoss=false, input_pulse_shape=:gauss, include_Δβion=false, ionisation_fraction=nothing, kwargs...) mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) @@ -405,7 +413,7 @@ function make_PhaseMatching(a, gas, pressure, λp; @assert ionisation_fraction >= 0 && ionisation_fraction <= 1 "Ionisation fraction must be between 0 and 1" end - ΔβnonlinCoeff = include_Δβnonlin ? Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; input_pulse_shape=input_pulse_shape) : 0.0 + ΔβnonlinCoeff = include_Δβnonlin ? Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=includeLoss,input_pulse_shape=input_pulse_shape) : 0.0 ΔβionCoeff = include_Δβion ? Δβion(a, gas, pressure, λp, ionisation_fraction; kwargs...) : 0.0 return PhaseMatching(mode, ωsol, β0, β1, ΔβnonlinCoeff, ΔβionCoeff) @@ -426,7 +434,7 @@ end λlims: limits for the root finding algorithm that looks for the phase-matching; by default these are set to the minimum wavelength just above the first resonance in the Sellmeier equation for the given gas and the maximum for the pump wavelength minus 1 nm """ function λRDWfull(a, gas, pressure, λp; - include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, input_pulse_shape=:gauss, + include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, includeLoss=false, input_pulse_shape=:gauss, include_Δβion=false, ionisation_fraction=nothing, λlims=nothing, kwargs...) @@ -448,7 +456,7 @@ function λRDWfull(a, gas, pressure, λp; end phase_matching = make_PhaseMatching(a, gas, pressure, λp; - include_Δβnonlin=include_Δβnonlin, soliton_order=soliton_order, τFWHM=τFWHM, input_pulse_shape=input_pulse_shape, + include_Δβnonlin=include_Δβnonlin, soliton_order=soliton_order, τFWHM=τFWHM, includeLoss=includeLoss, input_pulse_shape=input_pulse_shape, include_Δβion=include_Δβion, ionisation_fraction=ionisation_fraction, kwargs...) # if ionization effects are included, shrink λlims, because the dispersion curve gets lifted in the IR and a new phase-matched point appears near the pump From ae59211bae9a2bd376bcb5229ba15e397c431c0e Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 13 May 2026 15:53:42 +0100 Subject: [PATCH 10/13] Missed to import the attenuation constant from the right module :( --- src/Tools.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Tools.jl b/src/Tools.jl index 6b64560f..f38d548a 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -349,7 +349,7 @@ function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=fa if includeLoss Lfiss = params.Lfiss mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) - attenuation = α(mode, PhysData.wlfreq(λp); z=Lfiss) + attenuation = Modes.α(mode, PhysData.wlfreq(λp); z=Lfiss) end return params.γ*(soliton_factor*Pp*attenuation)*(1/ωsol) From e810742a12fa0de23b3ac5e56a7a220e17bd185f Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 13 May 2026 16:13:47 +0100 Subject: [PATCH 11/13] Removing the kwargs... fromthe phase-matching functions --- src/Tools.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/Tools.jl b/src/Tools.jl index f38d548a..ccb094f9 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -274,9 +274,9 @@ end """ Helper function to calculate the linear contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission """ -function Δβlin(a, gas, pressure, λp; kwargs...) +function Δβlin(a, gas, pressure, λp) - mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) + mode = Capillary.MarcatiliMode(a, gas, pressure) ωsol = PhysData.wlfreq(λp) β1 = Modes.dispersion(mode, 1, ωsol) β0 = Modes.β(mode, ωsol) @@ -328,7 +328,7 @@ end """ Helper function to calculate the nonlinear contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission """ -function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=false, input_pulse_shape=:gauss, kwargs...) +function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=false, input_pulse_shape=:gauss) energy = energyN(a, gas, pressure, τFWHM, λp, soliton_order) params = Tools.capillary_params(energy, τFWHM, λp, a, gas; P=pressure) @@ -348,7 +348,7 @@ function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=fa if includeLoss Lfiss = params.Lfiss - mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) + mode = Capillary.MarcatiliMode(a, gas, pressure) attenuation = Modes.α(mode, PhysData.wlfreq(λp); z=Lfiss) end @@ -359,12 +359,12 @@ end """ Helper function to calculate the ionisation contribution to the phase-mismatch Δβ, which is used in the phase-matching condition for RDW emission """ -function Δβion(a, gas, pressure, λp, ionisation_fraction; kwargs...) +function Δβion(a, gas, pressure, λp, ionisation_fraction) neutral_density = PhysData.density(gas, pressure) electron_density = ionisation_fraction*neutral_density critical_density = critical_electron_density(λp) ωsol = PhysData.wlfreq(λp) - mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) + mode = Capillary.MarcatiliMode(a, gas, pressure) nlin = real.(Capillary.neff(mode, ωsol)) # taking the real part here, otherwise the ionisation term becomes complex return (1/(2*PhysData.c*nlin))*(electron_density/critical_density)*(ωsol^2) end @@ -395,9 +395,9 @@ end """ function make_PhaseMatching(a, gas, pressure, λp; include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, includeLoss=false, input_pulse_shape=:gauss, - include_Δβion=false, ionisation_fraction=nothing, kwargs...) + include_Δβion=false, ionisation_fraction=nothing) - mode = Capillary.MarcatiliMode(a, gas, pressure; kwargs...) + mode = Capillary.MarcatiliMode(a, gas, pressure) ωsol = PhysData.wlfreq(λp) β1 = Modes.dispersion(mode, 1, ωsol) β0 = Modes.β(mode, ωsol) @@ -414,7 +414,7 @@ function make_PhaseMatching(a, gas, pressure, λp; end ΔβnonlinCoeff = include_Δβnonlin ? Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=includeLoss,input_pulse_shape=input_pulse_shape) : 0.0 - ΔβionCoeff = include_Δβion ? Δβion(a, gas, pressure, λp, ionisation_fraction; kwargs...) : 0.0 + ΔβionCoeff = include_Δβion ? Δβion(a, gas, pressure, λp, ionisation_fraction) : 0.0 return PhaseMatching(mode, ωsol, β0, β1, ΔβnonlinCoeff, ΔβionCoeff) end @@ -436,7 +436,7 @@ end function λRDWfull(a, gas, pressure, λp; include_Δβnonlin=false, soliton_order=nothing, τFWHM=nothing, includeLoss=false, input_pulse_shape=:gauss, include_Δβion=false, ionisation_fraction=nothing, - λlims=nothing, kwargs...) + λlims=nothing) if isnothing(λlims) λUVlim = Dict( @@ -457,7 +457,7 @@ function λRDWfull(a, gas, pressure, λp; phase_matching = make_PhaseMatching(a, gas, pressure, λp; include_Δβnonlin=include_Δβnonlin, soliton_order=soliton_order, τFWHM=τFWHM, includeLoss=includeLoss, input_pulse_shape=input_pulse_shape, - include_Δβion=include_Δβion, ionisation_fraction=ionisation_fraction, kwargs...) + include_Δβion=include_Δβion, ionisation_fraction=ionisation_fraction) # if ionization effects are included, shrink λlims, because the dispersion curve gets lifted in the IR and a new phase-matched point appears near the pump if include_Δβion From 5106a613aac95866ac87e9431c2350bfaf1afad2 Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 13 May 2026 16:24:28 +0100 Subject: [PATCH 12/13] Transmission is (1 - attenuation) --- src/Tools.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Tools.jl b/src/Tools.jl index ccb094f9..1c8f4ece 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -349,7 +349,7 @@ function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=fa if includeLoss Lfiss = params.Lfiss mode = Capillary.MarcatiliMode(a, gas, pressure) - attenuation = Modes.α(mode, PhysData.wlfreq(λp); z=Lfiss) + attenuation = 1.0 - Modes.α(mode, PhysData.wlfreq(λp); z=Lfiss) end return params.γ*(soliton_factor*Pp*attenuation)*(1/ωsol) From 399f756ed2872817d3a50058497066adbcd1fb71 Mon Sep 17 00:00:00 2001 From: Teodora Grigorova Date: Wed, 13 May 2026 16:40:02 +0100 Subject: [PATCH 13/13] Slightly changed how to include the loss --- src/Tools.jl | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/Tools.jl b/src/Tools.jl index 1c8f4ece..9fa07afc 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -334,6 +334,12 @@ function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=fa params = Tools.capillary_params(energy, τFWHM, λp, a, gas; P=pressure) ωsol = PhysData.wlfreq(λp) + if includeLoss # if include loss is true, decrease the energy for which the peak power is calculated + Lfiss = params.Lfiss + mode = Capillary.MarcatiliMode(a, gas, pressure) + energy = (1.0 - Modes.α(mode, PhysData.wlfreq(λp); z=Lfiss))*energy + end + if input_pulse_shape == :sech Pp = Ppeak(τFWHM, energy, pulse_shape=:sech) elseif input_pulse_shape == :gauss @@ -344,15 +350,7 @@ function Δβnonlin(a, gas, pressure, τFWHM, λp, soliton_order; includeLoss=fa soliton_factor = ((2*soliton_order-1)/soliton_order)^2 # ((2N-1)/N)^2 - attenuation = 1.0 - - if includeLoss - Lfiss = params.Lfiss - mode = Capillary.MarcatiliMode(a, gas, pressure) - attenuation = 1.0 - Modes.α(mode, PhysData.wlfreq(λp); z=Lfiss) - end - - return params.γ*(soliton_factor*Pp*attenuation)*(1/ωsol) + return params.γ*(soliton_factor*Pp)*(1/ωsol) end