LupoLab / Luna.jl

Nonlinear optical pulse propagator
MIT License
65 stars 27 forks source link

How to correctly propagate in vacuum speed of light reference frame #370

Open CaptainKubi opened 5 days ago

CaptainKubi commented 5 days ago

For some of our simulations, we would need the pulses propagated in the reference frame of the vacuum speed of light instead of the medium group velocity to incorporate the result later in another simulation.

We are doing free space simulations with radial symmetry using z-dependent media.

I found a section in the code here: https://github.com/LupoLab/Luna.jl/blob/98aa6d8bad97408270e6cd58b62be2a720c2e2df/src/LinearOps.jl#L167C1-L183C4

"""
    make_linop(grid, q::QDHT, nfun)

Make z-dependent linear operator for radial free-space propagation. `nfun(ω; z)` should
return the refractive index as a function of frequency `ω` and (kwarg) propagation
distance `z`.
"""
function make_linop(grid::Grid.RealGrid, q::Hankel.QDHT, nfun)
    kr2 = q.k.^2
    k2 = zero(grid.ω)
    nfunλ(z) = λ -> nfun(wlfreq(λ), z=z)
    function linop!(out, z)
        β1 = PhysData.dispersion_func(1, nfunλ(z))(grid.referenceλ)
        k2[grid.sidx] .= (nfun.(grid.ω[grid.sidx]; z=z) .* grid.ω[grid.sidx]./PhysData.c).^2
        _fill_linop_r!(out, grid, β1, k2, kr2, q.N)
    end
end

The parameter β1 seems to be the 1 / group velocity in the medium at position z for the reference wavelength λ.

I assume that if I change this function to always use β1 = 1/PhysData.c, the propagation should happen in the reference frame of the vacuum speed of light.

Would this be the correct approach?

Is there some other point in the code that influences the reference frame?

In our first tests, we had some inconsistencies in the form of small shifts, that lead us to believe that we did not implement the different reference frame correctly.

If you could give us a hint how to implement this reference frame correctly, we would be very thankful!

chrisbrahms commented 2 days ago

Hi!

Yes, make_linop is where the frame velocity is defined, and setting β1 = 1/c should result in a frame moving at the speed of light. I don't think it shows up anywhere else. My approach to add this option would probably be an extra keyword argument to make_linop which can override/replace the calculated value when given (pull request welcome if you implement this!).

Regarding the unexpected shifts, you could put up a minimum working example here which reproduces the problem and I could have a look. I've never needed to do this, so haven't had cause to look into it.