From cce1458a11ac97dda6b85757c85e72e17846a63f Mon Sep 17 00:00:00 2001 From: cbrahms Date: Tue, 28 Sep 2021 20:40:18 +0100 Subject: [PATCH 1/4] cycle average for ADK --- src/Ionisation.jl | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/Ionisation.jl b/src/Ionisation.jl index 1e483933..3d7361ec 100644 --- a/src/Ionisation.jl +++ b/src/Ionisation.jl @@ -17,7 +17,7 @@ Return a closure `ionrate!(out, E)` which calculates the ADK ionisation rate for field `E` and places the result in `out`. If `threshold` is true, use [`ADK_threshold`](@ref) to avoid calculation below floating-point precision. """ -function ionrate_fun!_ADK(ionpot::Float64, threshold=true) +function ionrate_fun!_ADK(ionpot::Float64, threshold=true; cycle_average=false) nstar = sqrt(0.5/(ionpot/au_energy)) cn_sq = 2^(2*nstar)/(nstar*gamma(nstar+1)*gamma(nstar)) ω_p = ionpot/ħ @@ -29,12 +29,24 @@ function ionrate_fun!_ADK(ionpot::Float64, threshold=true) thr = 0 end + Ip_au = ionpot / au_energy + F0_au = (2Ip_au)^(3/2) + F0 = F0_au*au_Efield + avfac = sqrt.(3/(π*F0)) + + ionrate! = let nstar=nstar, cn_sq=cn_sq, ω_p=ω_p, ω_t_prefac=ω_t_prefac, thr=thr function ir(E) if abs(E) >= thr - (ω_p*cn_sq* - (4*ω_p/(ω_t_prefac*abs(E)))^(2*nstar-1) - *exp(-4/3*ω_p/(ω_t_prefac*abs(E)))) + if cycle_average + (avfac*sqrt(abs(E))*ω_p*cn_sq* + (4*ω_p/(ω_t_prefac*abs(E)))^(2*nstar-1) + *exp(-4/3*ω_p/(ω_t_prefac*abs(E)))) + else + (ω_p*cn_sq* + (4*ω_p/(ω_t_prefac*abs(E)))^(2*nstar-1) + *exp(-4/3*ω_p/(ω_t_prefac*abs(E)))) + end else zero(E) end @@ -119,9 +131,9 @@ function ionrate_fun!_PPTcached(material::Symbol, λ0; kwargs...) end function ionrate_fun!_PPTcached(ionpot::Float64, λ0, Z, l; - sum_tol=1e-4, rcycle=true, N=2^16, Emax=nothing, + sum_tol=1e-4, cycle_average=false, N=2^16, Emax=nothing, cachedir=joinpath(homedir(), ".luna", "pptcache")) - h = hash((ionpot, λ0, Z, l, sum_tol, rcycle, N, Emax)) + h = hash((ionpot, λ0, Z, l, sum_tol, cycle_average, N, Emax)) fname = string(h, base=16)*".h5" fpath = joinpath(cachedir, fname) lockpath = joinpath(cachedir, "pptlock") @@ -134,7 +146,7 @@ function ionrate_fun!_PPTcached(ionpot::Float64, λ0, Z, l; return rate else E, rate = makePPTcache(ionpot::Float64, λ0, Z, l; - sum_tol=sum_tol, rcycle=rcycle, N=N, Emax=Emax) + sum_tol=sum_tol, cycle_average, N=N, Emax=Emax) @info "Saving PPT rate cache for $(ionpot/electron) eV, $(λ0*1e9) nm in $cachedir" pidlock = mkpidlock(lockpath) if isfile(fpath) # makePPTcache takes a while - has another process saved first? @@ -160,7 +172,7 @@ function loadPPTaccel(fpath) end function makePPTcache(ionpot::Float64, λ0, Z, l; - sum_tol=1e-4, rcycle=true, N=2^16, Emax=nothing) + sum_tol=1e-4, cycle_average=false, N=2^16, Emax=nothing) Emax = isnothing(Emax) ? 2*barrier_suppression(ionpot, Z) : Emax # ω0 = 2π*c/λ0 @@ -169,7 +181,7 @@ function makePPTcache(ionpot::Float64, λ0, Z, l; E = collect(range(Emin, stop=Emax, length=N)); @info "Pre-calculating PPT rate for $(ionpot/electron) eV, $(λ0*1e9) nm" - rate = ionrate_PPT(ionpot, λ0, Z, l, E; sum_tol=sum_tol, rcycle=rcycle); + rate = ionrate_PPT(ionpot, λ0, Z, l, E; sum_tol=sum_tol, cycle_average); @info "PPT pre-calcuation done" return E, rate end @@ -235,13 +247,13 @@ function ionrate_fun!_PPT(args...) end """ - ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, rcycle=true) + ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, cycle_average=false) Create closure to calculate PPT ionisation rate. # Keyword arguments - `sum_tol::Number`: Relative tolerance used to truncate the infinite sum. -- `rcycle::Bool`: If `true` (default), remove the cycle-average factor. +- `cycle_average::Bool`: If `false` (default), calculate the cycle-averaged rate # References [1] Ilkov, F. A., Decker, J. E. & Chin, S. L. @@ -254,7 +266,7 @@ Ultrashort filaments of light in weakly ionized, optically transparent media. Rep. Prog. Phys. 70, 1633–1713 (2007) (Appendix A) """ -function ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, rcycle=true) +function ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, cycle_average=false) Ip_au = ionpot / au_energy ns = Z/sqrt(2Ip_au) ls = ns-1 @@ -294,7 +306,7 @@ function ionrate_fun_PPT(ionpot::Float64, λ0, Z, l; sum_tol=1e-4, rcycle=true) lret *= exp(-2v*(asinh(γ) - γ*sqrt(1+γ2)/(1+2γ2))) lret *= Ip_au * γ2/(1+γ2) # Remove cycle average factor, see eq. (2) of [1] - if rcycle + if !cycle_average lret *= sqrt(π*E0_au/(3E_au)) end k = ceil(v) From 94b8f27b8bc027ff5d5170642bde4a0c0c87062a Mon Sep 17 00:00:00 2001 From: cbrahms Date: Tue, 28 Sep 2021 21:02:27 +0100 Subject: [PATCH 2/4] add test --- test/test_ionisation.jl | 33 ++++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/test/test_ionisation.jl b/test/test_ionisation.jl index a28e4ace..e0ca9f6c 100644 --- a/test/test_ionisation.jl +++ b/test/test_ionisation.jl @@ -1,5 +1,6 @@ import Test: @test, @testset -import Luna: Ionisation, Maths, PhysData +import Luna: Ionisation, Maths, PhysData, Tools +import NumericalIntegration: integrate, SimpsonEven @test Ionisation.ionrate_ADK(:He, 1e10) ≈ 1.2416371415312408e-18 @test Ionisation.ionrate_ADK(:He, 2e10) ≈ 1.0772390893742478 @@ -44,6 +45,36 @@ outneg = similar(out) ratefun!(outneg, -E) @test out == outneg +## +@testset "cycle-averaging $gas" for gas in (:He, :Ar, :Xe) + + bsi = Tools.field_to_intensity( + Ionisation.barrier_suppression( + PhysData.ionisation_potential(gas), + 1)) + isy = 10 .^ collect(16:0.01:log10(bsi)) + E0 = Tools.intensity_to_field.(isy) + + Ip_au = PhysData.ionisation_potential(gas; unit=:atomic) + F0_au = (2Ip_au)^(3/2) + F0 = F0_au*PhysData.au_Efield + + t = collect(range(0, 2π; length=2^12)) + Et = sin.(t) + + adk_avg = zero(E0) + rf! = Ionisation.ionrate_fun!_ADK(gas) + out = zero(Et) + for (idx, E0i) in enumerate(E0) + rf!(out, E0i*Et) + adk_avg[idx] = 1/2π * integrate(t, out, SimpsonEven()) + end + + adk = Ionisation.ionrate_ADK.(gas, E0) + adk_avg_kw = Ionisation.ionrate_ADK.(gas, E0; cycle_average=true) + @test all(isapprox.(adk_avg, adk_avg_kw; rtol=0.05)) +end +## # this_folder = dirname(@__FILE__) # import CSV # import PyPlot: plt, pygui From 96826cc30c3d3ab556a2783754729f05a77e3c2c Mon Sep 17 00:00:00 2001 From: cbrahms Date: Tue, 28 Sep 2021 21:03:41 +0100 Subject: [PATCH 3/4] docstring --- src/Ionisation.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Ionisation.jl b/src/Ionisation.jl index 3d7361ec..d94d25b9 100644 --- a/src/Ionisation.jl +++ b/src/Ionisation.jl @@ -15,7 +15,8 @@ import Luna: @hlock Return a closure `ionrate!(out, E)` which calculates the ADK ionisation rate for the electric field `E` and places the result in `out`. If `threshold` is true, use [`ADK_threshold`](@ref) -to avoid calculation below floating-point precision. +to avoid calculation below floating-point precision. If `cycle_average` is `true`, calculate +the cycle-averaged ADK ionisation rate instead. """ function ionrate_fun!_ADK(ionpot::Float64, threshold=true; cycle_average=false) nstar = sqrt(0.5/(ionpot/au_energy)) From aa3964628e03f73f8f6a41d1e63a37007fd6cfb5 Mon Sep 17 00:00:00 2001 From: cbrahms Date: Thu, 7 Oct 2021 10:02:37 +0100 Subject: [PATCH 4/4] review edits --- src/Ionisation.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/Ionisation.jl b/src/Ionisation.jl index d94d25b9..e4d2cb78 100644 --- a/src/Ionisation.jl +++ b/src/Ionisation.jl @@ -30,6 +30,9 @@ function ionrate_fun!_ADK(ionpot::Float64, threshold=true; cycle_average=false) thr = 0 end + # Zenghu Chang: Fundamentals of Attosecond Optics (2011) p. 184 + # Section 4.2.3.1 Cycle-Averaged Rate + # ̄w_ADK(Fₐ) = √(3/π) √(Fₐ/F₀) w_ADK(Fₐ) where Fₐ is the field amplitude Ip_au = ionpot / au_energy F0_au = (2Ip_au)^(3/2) F0 = F0_au*au_Efield @@ -39,17 +42,15 @@ function ionrate_fun!_ADK(ionpot::Float64, threshold=true; cycle_average=false) ionrate! = let nstar=nstar, cn_sq=cn_sq, ω_p=ω_p, ω_t_prefac=ω_t_prefac, thr=thr function ir(E) if abs(E) >= thr - if cycle_average - (avfac*sqrt(abs(E))*ω_p*cn_sq* - (4*ω_p/(ω_t_prefac*abs(E)))^(2*nstar-1) - *exp(-4/3*ω_p/(ω_t_prefac*abs(E)))) - else - (ω_p*cn_sq* + r = (ω_p*cn_sq* (4*ω_p/(ω_t_prefac*abs(E)))^(2*nstar-1) *exp(-4/3*ω_p/(ω_t_prefac*abs(E)))) + if cycle_average + r *= avfac*sqrt(abs(E)) end + return r else - zero(E) + return zero(E) end end function ionrate!(out, E)