Merge branch 'master' into version

README.md

@@ -5,7 +5,7 @@
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://lupo-lab.com/Luna.jl)
 
 
-Luna.jl is a flexible platform for the simulation of nonlinear optical dynamics—both in waveguides (such as optical fibres) and free-space geometries—using the unidirectional pulse propagation equation (UPPE). Some of the key features of Luna:
+Luna.jl is a flexible platform for the simulation of nonlinear optical dynamics—both in waveguides (such as optical fibres) and free-space geometries—using the unidirectional pulse propagation equation (UPPE) and its approximate forms, such as the commonly used generalised nonlinear Schrödinger equation (GNLSE). Some of the key features of Luna:
 - A variety of propagation geometries treated in a unified way:
     - Single-mode (mode-averaged) propagation in waveguides
     - Multi-mode propagation in waveguides with arbitrary (including non-symmetric) mode-shapes, full polarisation resolution, and intermodal coupling for arbitrary nonlinear polarisation terms
@@ -16,22 +16,22 @@ Luna.jl is a flexible platform for the simulation of nonlinear optical dynamics
 - A range of linear and nonlinear optical effects:
     - Modal dispersion and loss in waveguides
     - Optical Kerr effect (third-order nonlinearity)
-    - Raman scattering in molecular gases
+    - Raman scattering in molecular gases or glasses
     - Strong-field photoionisation and plasma dynamics
 - A built-in interface for the running and processing of multi-dimensional [parameter scans](#running-parameter-scans) in serial or parallel
 - A standard library of [plotting](#plotting-results) and [processing](#output-processing) functions, including calculation of spectrograms and beam properties
 
 Luna is designed to be extensible: adding e.g. a new type of waveguide or a new nonlinear effect is straightforward, even without editing the main source code.
 
-Luna was originally developed for modelling ultrafast pulse propagation in gas-filled hollow capillary fibres and hollow-core photonic crystal fibres. Therefore such simulations have particularly good support.
+Luna was originally developed for modelling ultrafast pulse propagation in gas-filled hollow capillary fibres and hollow-core photonic crystal fibres. Therefore, such simulations have particularly good support. It is also excellent for modelling propagation in solid-core fibres.
 
 Luna is written in the [Julia programming language](https://julialang.org/), chosen for its unique combination of readability, ease of use, and speed. If you want to use Luna but are new to Julia, see [the relevant section of this README](#new-to-julia).
 
 There are two ways of using Luna:
-1. A very simple high-level interface for the most heavily optimised application of Luna—propagation in gas-filled hollow capillary fibres and hollow-core photonic crystal fibres—consisting of the function [`prop_capillary`](#quickstart) and some helper functions to create input pulses.
+1. A very simple high-level interface for the most heavily optimised applications of Luna: propagation in gas-filled hollow capillary fibres and hollow-core photonic crystal fibres (consisting of the function [`prop_capillary`](#quickstart) and some helper functions to create input pulses); or propagation of simple GNLSE simulations (consisting of the function [`prop_gnlse`](#gnlse)).
 2. A low-level interface which allows for full control and customisation of the simulation parameters, the use of custom waveguide modes and gas fills (including gas mixtures), and free-space propagation simulations.
 
-For a short introduction on how to use the simple interface, see the [Quickstart](#quickstart) section below. More information, including on the internals of Luna, can be found in the [Documentation](http://lupo-lab.com/Luna.jl).
+For a short introduction on how to use the simple interface, see the [Quickstart](#quickstart) or [GNLSE](#gnlse) sections below. More information, including on the internals of Luna, can be found in the [Documentation](http://lupo-lab.com/Luna.jl).
 
 ## Installation
 Luna requires Julia v1.7, which can be obtained from [here](https://julialang.org/downloads/). Once Julia is installed, open a new Julia terminal and install Luna:
@@ -126,6 +126,24 @@ The `Processing` module contains many useful functions for more detailed process
 - FWHM widths in frequency (`Processing.fwhm_f`) and time (`Processing.fwhm_t`) as well as time-bandwidth product (`Processing.time_bandwidth`)
 - g₁₂ coherence between multiple fields (`Processing.coherence`)
 
+## GNLSE propagation
+To run a simple simulation of nonlinear pulse propagation in an optical fibre using the generalised nonlinear Schrödinger equation (GNLSE), you can use `prop_gnlse`. As an example, we can model supercontinuum generation in a solid-core photonic crystal fibre for parameters corresponding to the simulations in Fig. 3 of Dudley et. al, RMP 78 1135 (2006).
+```julia
+julia> using Luna
+julia> γ = 0.11
+julia> flength = 15e-2
+julia> βs = [0.0, 0.0, -1.1830e-26, 8.1038e-41, -9.5205e-56,  2.0737e-70, -5.3943e-85,  1.3486e-99, -2.5495e-114,  3.0524e-129, -1.7140e-144]
+julia> output = prop_gnlse(γ, flength, βs; λ0=835e-9, τfwhm=50e-15, power=10e3, pulseshape=:sech, λlims=(400e-9, 2400e-9), trange=12.5e-12)
+```
+After this has run, you can visualise the output, with e.g.
+```julia
+julia> Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
+PyPlot.Figure(PyObject <Figure size 2400x800 with 4 Axes>)
+```
+This should show a plot like this:
+![GNLSE propagation example](assets/readme_gnlse_scg.png)
+
+
 ## Examples
 The [examples folder](examples/) contains complete simulation examples for a variety of scenarios, both for the [simple interface](examples/simple_interface/) and the [low-level interface](examples/low_level_interface). Some of the simple interface examples require the `PyPlot` package to be present, and many of the low-level examples require other packages as well--you can install these by simply typing `] add PyPlot` at the Julia REPL or the equivalent for other packages.
 

examples/low_level_interface/gnlse/simplescg_modeAvg.jl

@@ -0,0 +1,40 @@
+# supercontinuum from simple GNLSE parameters
+# Fig.3 of Dudley et. al, RMP 78 1135 (2006)
+
+using Luna
+
+βs =  [0.0, 0.0, -1.1830e-26, 8.1038e-41, -9.5205e-56,  2.0737e-70, -5.3943e-85,  1.3486e-99, -2.5495e-114,  3.0524e-129, -1.7140e-144]
+γ = 0.11
+flength = 15e-2
+fr = 0.18
+τfwhm = 50e-15
+λ0 = 835e-9
+energy = 568e-12
+
+grid = Grid.RealGrid(flength, λ0, (400e-9, 1400e-9), 10e-12)
+
+m = SimpleFibre.SimpleMode(PhysData.wlfreq(λ0), βs)
+aeff = z -> 1.0
+densityfun = z -> 1.0
+
+linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)
+
+k0 = 2π/λ0
+n2 = γ/k0*aeff(0.0)
+χ3 = 4/3 * n2 * (PhysData.ε_0*PhysData.c)
+responses = (Nonlinear.Kerr_field((1 - fr)*χ3),
+             Nonlinear.RamanPolarField(grid.to, Raman.raman_response(grid.to, :SiO2, fr*χ3*PhysData.ε_0)))
+
+inputs = (Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy), Fields.ShotNoise())
+norm! = NonlinearRHS.norm_mode_average_gnlse(grid, aeff)
+Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff, norm! = norm!)
+
+output = Output.MemoryOutput(0, grid.zmax, 201)
+Luna.run(Eω, grid, linop, transform, FT, output)
+
+##
+Plotting.pygui(true)
+#Plotting.stats(output)
+#Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
+#Plotting.time_1D(output, range(0.0, 1.0, length=5).*flength, trange=(-1e-12, 5e-12))
+Plotting.spec_1D(output, range(0.0, 1.0, length=5).*flength, λrange=(400e-9, 1300e-9))

examples/low_level_interface/gnlse/simplescg_modeAvg_env.jl

@@ -0,0 +1,40 @@
+# supercontinuum from simple GNLSE parameters
+# Fig.3 of Dudley et. al, RMP 78 1135 (2006)
+
+using Luna
+
+βs =  [0.0, 0.0, -1.1830e-26, 8.1038e-41, -9.5205e-56,  2.0737e-70, -5.3943e-85,  1.3486e-99, -2.5495e-114,  3.0524e-129, -1.7140e-144]
+γ = 0.11
+flength = 15e-2
+fr = 0.18
+τfwhm = 50e-15
+λ0 = 835e-9
+energy = 568e-12
+
+grid = Grid.EnvGrid(flength, λ0, (400e-9, 1400e-9), 10e-12)
+
+m = SimpleFibre.SimpleMode(PhysData.wlfreq(λ0), βs)
+aeff = z -> 1.0
+densityfun = z -> 1.0
+
+linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)
+
+k0 = 2π/λ0
+n2 = γ/k0*aeff(0.0)
+χ3 = 4/3 * n2 * (PhysData.ε_0*PhysData.c)
+responses = (Nonlinear.Kerr_env((1 - fr)*χ3),
+             Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, :SiO2, fr*χ3*PhysData.ε_0)))
+
+inputs = (Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy), Fields.ShotNoise())
+norm! = NonlinearRHS.norm_mode_average_gnlse(grid, aeff)
+Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff, norm! = norm!)
+
+output = Output.MemoryOutput(0, grid.zmax, 201)
+Luna.run(Eω, grid, linop, transform, FT, output)
+
+##
+Plotting.pygui(true)
+#Plotting.stats(output)
+Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
+#Plotting.time_1D(output, [0.0, 2.5, 5.0], trange=(-5e-12, 5e-12))
+Plotting.spec_1D(output, range(0.0, 1.0, length=5).*flength, λrange=(400e-9, 1300e-9))

examples/low_level_interface/stepindex/step_modeAvg_env.jl

@@ -0,0 +1,39 @@
+# propagation in step index fibre
+
+using Luna
+
+# single mode fibre at 1030 nm
+a = 5e-6
+NA = 0.08
+flength = 2.0
+fr = 0.18
+τfwhm = 1e-12
+λ0 = 1030e-9
+energy = 10e-9
+
+grid = Grid.EnvGrid(flength, λ0, (980e-9, 1200e-9), 10e-12)
+
+m = StepIndexFibre.StepIndexMode(a, NA, accellims=(900e-9, 1200e-9, 100))
+aeff = let aeffc=Modes.Aeff(m, z=0.0)
+    z -> aeffc
+end
+densityfun = z -> 1.0
+
+linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)
+
+responses = (Nonlinear.Kerr_env((1 - fr)*PhysData.χ3(:SiO2)),
+             Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, :SiO2, fr*PhysData.ε_0*PhysData.χ3(:SiO2))))
+
+inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)
+Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
+
+statsfun = Stats.default(grid, Eω, m, linop, transform)
+output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
+Luna.run(Eω, grid, linop, transform, FT, output)
+
+##
+Plotting.pygui(true)
+#Plotting.stats(output)
+Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(980e-9, 1200e-9), trange=(-2e-12, 2e-12))
+#Plotting.time_1D(output, [0.0, 2.5, 5.0], trange=(-5e-12, 5e-12))
+Plotting.spec_1D(output, [0.0, 2.5, 5.0], λrange=(980e-9, 1080e-9))

examples/low_level_interface/stepindex/stepscg_modeAvg.jl

@@ -0,0 +1,38 @@
+# supercontinuum in strand of silica in air
+
+using Luna
+
+# single mode fibre at 1030 nm
+a = 1.25e-6
+flength = 15e-2
+fr = 0.18
+τfwhm = 50e-15
+λ0 = 835e-9
+energy = 568e-12
+
+grid = Grid.RealGrid(flength, λ0, (400e-9, 1400e-9), 10e-12)
+
+m = StepIndexFibre.StepIndexMode(a, accellims=(400e-9, 1400e-9, 100))
+aeff = let aeffc=Modes.Aeff(m, z=0.0)
+    z -> aeffc
+end
+densityfun = z -> 1.0
+
+linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)
+
+responses = (Nonlinear.Kerr_field((1 - fr)*PhysData.χ3(:SiO2)),
+             Nonlinear.RamanPolarField(grid.to, Raman.raman_response(grid.to, :SiO2, fr*PhysData.ε_0*PhysData.χ3(:SiO2))))
+
+inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)
+Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
+
+statsfun = Stats.default(grid, Eω, m, linop, transform)
+output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
+Luna.run(Eω, grid, linop, transform, FT, output)
+
+##
+Plotting.pygui(true)
+#Plotting.stats(output)
+#Plotting.prop_2D(output)
+#Plotting.time_1D(output, [0.0, 2.5, 5.0], trange=(-5e-12, 5e-12))
+Plotting.spec_1D(output, range(0.0, 1.0, length=5).*flength, λrange=(400e-9, 1300e-9))

examples/low_level_interface/stepindex/stepscg_modeAvg_env.jl

@@ -0,0 +1,38 @@
+# supercontinuum in strand of silica in air
+
+using Luna
+
+# single mode fibre at 1030 nm
+a = 1.25e-6
+flength = 15e-2
+fr = 0.18
+τfwhm = 50e-15
+λ0 = 835e-9
+energy = 568e-12
+
+grid = Grid.EnvGrid(flength, λ0, (400e-9, 1400e-9), 10e-12)
+
+m = StepIndexFibre.StepIndexMode(a, accellims=(400e-9, 1400e-9, 100))
+aeff = let aeffc = Modes.Aeff(m, z=0)
+    z -> aeffc
+end
+densityfun = z -> 1.0
+
+linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)
+
+responses = (Nonlinear.Kerr_env((1 - fr)*PhysData.χ3(:SiO2)),
+             Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, :SiO2, fr*PhysData.ε_0*PhysData.χ3(:SiO2))))
+
+inputs = (Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy), Fields.ShotNoise())
+Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
+
+
+output = Output.MemoryOutput(0, grid.zmax, 201)
+Luna.run(Eω, grid, linop, transform, FT, output)
+
+##
+Plotting.pygui(true)
+#Plotting.stats(output)
+Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
+#Plotting.time_1D(output, [0.0, 2.5, 5.0], trange=(-5e-12, 5e-12))
+Plotting.spec_1D(output, range(0.0, 1.0, length=5).*flength, λrange=(400e-9, 1300e-9))

examples/low_level_interface/stepindex/stepsimplecmp.jl

@@ -0,0 +1,86 @@
+using Luna
+using Luna.PhysData: Polynomials
+import PyPlot: plt
+
+a = 1.25e-6
+flength = 15e-2
+fr = 0.18
+τfwhm = 50e-15
+λ0 = 835e-9
+energy = 568e-12
+grid = Grid.EnvGrid(flength, λ0, (400e-9, 1400e-9), 10e-12)
+
+# supercontinuum in a strand of silica in air
+m = StepIndexFibre.StepIndexMode(a, accellims=(400e-9, 1400e-9, 100))
+aeff = let aeffc = Modes.Aeff(m, z=0)
+    z -> aeffc
+end
+densityfun = z -> 1.0
+linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)
+responses = (Nonlinear.Kerr_env((1 - fr)*PhysData.χ3(:SiO2)),
+             Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, :SiO2, fr*PhysData.ε_0*PhysData.χ3(:SiO2))))
+inputs = (Fields.SechField(λ0=λ0, τfwhm=τfwhm, energy=energy), Fields.ShotNoise())
+Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
+outputm = Output.MemoryOutput(0, grid.zmax, 201)
+Luna.run(Eω, grid, linop, transform, FT, outputm)
+
+Plotting.prop_2D(outputm, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
+
+##
+# Approximate step-index mode with β coefficients
+λmin = 400e-9
+λmax = 1400e-9
+order = 11
+ωs = range(PhysData.wlfreq(λmax), PhysData.wlfreq(λmin), length=512)
+ω0 = PhysData.wlfreq(λ0)
+βs = Modes.β.(m, ωs)
+# help the fit by scaling frequency down to around unity
+p = Polynomials.fit(1e-15*(ωs .- ω0), βs, order)
+# each term n=0...order is scaled by 1e15ⁿ by scaling of ω - undo that
+# also need to multiply by n! to match the form of β expansion, which is:
+# β(ω) ≈ ∑₀ⁿ βₙ(ω-ω₀)/n! with βₙ = ∂β/∂ω at ω0
+βcoeffs = p.coeffs .* 1e-15 .^ (collect(0:order)) .* factorial.(0:order)
+
+s = SimpleFibre.SimpleMode(ω0, βcoeffs)
+plt.figure()
+plt.plot(ωs .- ω0, Modes.β_ret.(m, ωs; λ0), label="Full")
+plt.plot(ωs .- ω0, Modes.β_ret(s, ωs; λ0), "--", label="Polynomial fit")
+plt.legend()
+plt.xlabel("ω-ω0 (rad/s)")
+plt.ylabel("β (1/m)")
+
+N0, n0, n2 = Tools.getN0n0n2(ω0, :SiO2)
+γ = Tools.getγ(ω0, m, n2)
+
+##
+aeff = z -> 1.0
+linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, s, λ0)
+k0 = 2π/λ0
+n2 = γ/k0*aeff(0.0)
+n0 = real(PhysData.ref_index(:SiO2, λ0))
+χ3 = 4/3 * n2 * (PhysData.ε_0*PhysData.c) * n0 * n0 / Modes.neff(m, ω0)
+responses = (Nonlinear.Kerr_env((1 - fr)*χ3),
+             Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, :SiO2, fr*χ3*PhysData.ε_0)))
+norm! = NonlinearRHS.norm_mode_average_gnlse(grid, aeff)
+Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff, norm! = norm!)
+outputs = Output.MemoryOutput(0, grid.zmax, 201)
+Luna.run(Eω, grid, linop, transform, FT, outputs)
+
+Plotting.prop_2D(outputs, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
+
+#isapprox(Processing.getIω(outputm, :λ, flength)[2], Processing.getIω(outputs, :λ, flength)[2], rtol=1.1e-1)
+
+##
+Plotting.prop_2D(outputm)
+Plotting.prop_2D(outputs)
+##
+λ, Iλm = Processing.getIω(outputm, :λ, flength)
+_, Iλs = Processing.getIω(outputs, :λ, flength)
+plt.figure()
+plt.semilogy(1e9λ, Iλm, label="step-index mode")
+plt.semilogy(1e9λ, Iλs, "--",  label="GNLSE approximation")
+plt.legend()
+plt.xlabel("Wavelength (nm)")
+plt.ylabel("SED (a.u.)")
+
+# close but not exact. This is because we cannot fully cancel the frequency dependence of neff.

examples/simple_interface/gnlse_scg.jl

@@ -0,0 +1,19 @@
+# supercontinuum from simple GNLSE parameters
+# Fig.3 of Dudley et. al, RMP 78 1135 (2006)
+
+using Luna
+
+γ = 0.11
+flength = 15e-2
+βs =  [0.0, 0.0, -1.1830e-26, 8.1038e-41, -9.5205e-56,  2.0737e-70, -5.3943e-85,  1.3486e-99, -2.5495e-114,  3.0524e-129, -1.7140e-144]
+
+τfwhm = 50e-15
+λ0 = 835e-9
+power = 10000.0
+
+output = prop_gnlse(γ, flength, βs; λ0, τfwhm, power, pulseshape=:sech, λlims=(400e-9, 2400e-9), trange=12.5e-12, ramanmodel=:sdo, τ1=12.2e-15, τ2=32e-15)
+
+##
+Plotting.pygui(true)
+Plotting.prop_2D(output, :λ, dBmin=-40.0,  λrange=(400e-9, 1300e-9), trange=(-1e-12, 5e-12))
+Plotting.spec_1D(output, range(0.0, 1.0, length=5).*flength, λrange=(400e-9, 1300e-9))

examples/simple_interface/gnlse_sol.jl

@@ -0,0 +1,24 @@
+# N = 5 soliton
+
+using Luna
+
+γ = 0.1
+β2 = -1e-26
+N = 4.0
+τ0 = 280e-15
+τfwhm = (2*log(1 + sqrt(2)))*τ0
+fr = 0.18
+P0 = N^2*abs(β2)/((1-fr)*γ*τ0^2)
+flength = π*τ0^2/abs(β2)
+βs =  [0.0, 0.0, β2]
+
+λ0 = 835e-9
+λlims = [450e-9, 8000e-9]
+trange = 4e-12
+
+output = prop_gnlse(γ, flength, βs; λ0, τfwhm, power=P0, pulseshape=:sech, λlims, trange,
+                    raman=false, shock=false, fr, shotnoise=false)
+
+##
+Plotting.pygui(true)
+Plotting.prop_2D(output, :ω, dBmin=-100.0,  λrange=(720e-9,1000e-9), trange=(-300e-15, 300e-15), oversampling=1)