Glass materials in radial.jl

[glamacc]

Hi, I am trying to use radial.jl to simulate the nonlinear propagation (third order) through a SiO2 window. In particular I still haven't figured out how to modify "responses" and "linop" for the glass case. It would be awesome if you could point me in the right direction since I am very new to Jiulia. By the way, thank you for this very nice code!

[chrisbrahms]

Hi! Here is some example code to propagate through a window. This uses envelope propagation (without THG) to make it faster. If you want to include THG, replace EnvGrid with RealGrid and Kerr_env with Kerr_field.

Hope this helps! But feel free to follow up if something doesn't work.

using Luna
import Luna: Hankel
import PyPlot: plt

λ0 = 800e-9 # central wavelength

thickness = 1e-3 # window thickness
material = :SiO2 # window material

R = 15e-3 # pupil radius for Hankel transform
N = 2^10 # number of samples for Hankel transform

w0 = 3e-3 # 1/e² radius of the beam
energy = 3e-3 # pulse energy
τfwhm = 30e-15 # pulse duration FWHM

grid = Grid.EnvGrid(thickness, λ0, (400e-9, 6e-6), 400e-15) # time grid
q = Hankel.QDHT(R, N, dim=2) # space grid

χ3 = PhysData.χ3(material) # χ3 of silica
responses = (Nonlinear.Kerr_env(χ3),) # nonlinear responses 

nfun = PhysData.ref_index_fun(material) # function returning ref index of the window
nfunreal(λ) = real(nfun(λ)) # use only the real part--ignore (negligible) absorption

linop = LinearOps.make_const_linop(grid, q, nfunreal) # linear operator

normfun = NonlinearRHS.const_norm_radial(grid, q, nfunreal)
densityfun = z -> 1 # densityfun is required but not used for a window, so just set it to 1


inputs = Fields.GaussGaussField(;λ0, τfwhm, energy, w0)
Eω, transform, FT = Luna.setup(grid, q, densityfun, normfun, responses, inputs)

# the below line propagates *linearly* backwards through the window,
# such that the pulse without nonlinearity comes out of the window compressed
Fields.prop_material!(Eω, grid, material, -thickness, λ0)


output = Output.MemoryOutput(0, grid.zmax, 51) # 51 z slices in the output

Luna.run(Eω, grid, linop, transform, FT, output) # run the simulation

##
Eωk = output["Eω"][:, :, end] # ω-k-domain field at the window output
Eωr0 = Hankel.onaxis(Eωk, q) # frequency-domain field on the optical axis
Eωr = q \ Eωk # inv Hankel transform
idx = argmin(abs.(q.r .- w0)) # index into q.r closest to w0
Eω_w0 = Eωr[:, idx] # frequency-domain field one 1/e² radius away from the axis

λ, Iλr0 = Processing.getIω(Processing.getEω(grid, Eωr0)..., :λ)
_, Iλ_w0 = Processing.getIω(Processing.getEω(grid, Eω_w0)..., :λ)

plt.figure()
plt.semilogy(1e9λ, Maths.normbymax(Iλr0); label="On axis")
plt.semilogy(1e9λ, Maths.normbymax(Iλ_w0); label="At 1x 1/e² radius")
ymax = max(maximum(Iλ_w0), maximum(Iλr0))
plt.xlim(600, 1000)
plt.ylim(1e-4, 5)
plt.xlabel("Wavelength (nm)")
plt.ylabel("Spectral energy density (norm.)")
plt.legend()