diff --git a/src/Compressor.jl b/src/Compressor.jl index d49f351..60b8124 100644 --- a/src/Compressor.jl +++ b/src/Compressor.jl @@ -4,31 +4,52 @@ import HiSol.Solitons: T0P0 import HiSol.HCF: Aeff0, Leff, transmission, intensity_modeavg, α import HiSol.Focusing: max_flength, diverged_beam import HiSol.Data: n2_gas -import Luna.PhysData: pressure +import HiSol.GasFlow: tube_PVflow_choked, bar, slpm +import Luna.PhysData: pressure, density import Luna.Tools: τfw_to_τ0 import Logging: @debug, @info import Printf: @sprintf import Roots: find_zero import PyPlot: plt +function fix_pressure(pr, a, flength, gas, max_pressure, max_flow, ngradient) + pr = min(pr, max_pressure) + flow = 0.0 + (flength <= 0) && return 0.0, 0.0 + if ngradient > 0 + flow = ngradient*slpm(tube_PVflow_choked(a, flength/ngradient, gas, pr*bar)) + if flow > max_flow + pr = find_zero(pr) do prg + flow = ngradient*slpm(tube_PVflow_choked(a, flength/ngradient, + gas, prg*bar)) + max_flow - flow + end + end + end + flow = ngradient*slpm(tube_PVflow_choked(a, flength/ngradient, gas, pr*bar)) + pr, flow +end + function optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; thickness=1e-3, material=:SiO2, Bmax=0.2, LIDT=2000, S_fluence=5, S_sf=1.5, S_ion=10, - entrance_window=true, exit_window=true) + entrance_window=true, exit_window=true, ngradient=0, + max_pressure=60.0, max_flow=1.0, polarisation=:linear) factor = τfwhm_in/τfwhm_out ρcrit = critical_density(gas, λ0, τfwhm_in, energy) ρ = ρcrit/S_sf - pr = pressure(gas, ρ) - n2 = n2_gas(gas, pr) + prm = pressure(gas, ρ) @debug(@sprintf("ρcrit: %.4e m⁻³", ρcrit)) @debug(@sprintf("ρ: %.4e m⁻³", ρ)) - @debug(@sprintf("pressure: %.3f bar", pr)) - @debug(@sprintf("n₂: %.4e m²/W", n2)) + @debug(@sprintf("crit pressure: %.3f bar", pr)) + @debug(@sprintf("crit n₂: %.4e m²/W", n2)) _, P0 = T0P0(τfwhm_in, energy) φm = nonlinear_phase(factor) + φm *= (ngradient > 0 ? 3/2 : 1.0) + φm *= (polarisation == :circular ? 3/2 : 1.0) # φm = γP0Leff γLeff = φm/P0 @@ -53,6 +74,8 @@ function optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; aguess *= 0.99 continue end + pr, flow = fix_pressure(prm, aguess, flength, gas, max_pressure, max_flow, ngradient) + n2 = n2_gas(gas, pr) γLeff_this = n2*k0/(A0*aguess^2)*Leff(flength, aguess, λ0) @debug( @sprintf( @@ -64,9 +87,11 @@ function optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; enough = γLeff_this > γLeff if aguess < 10e-6 error("""Could not find core radius: not enough nonlinearity. - Decrease the compression factor or increase the maximum length.""") + Decrease the compression factor or increase the maximum length, + maximum pressure or maximum flow.""") end - if P0/(A0*aguess^2) > Isupp/S_ion + I_this = P0/(A0*aguess^2) * (polarisation == :circular ? 2/3 : 1.0) + if I_this > Isupp/S_ion error("""Could not find core radius: intensity too high. Decrease the energy or increase the maximum length.""") end @@ -78,6 +103,8 @@ function optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; global flength = max_flength(a, λ0, energy, τfwhm_in, maxlength; thickness, material, Bmax, LIDT, S_fluence, entrance_window, exit_window) + global pr, flow = fix_pressure(prm, a, flength, gas, max_pressure, max_flow, ngradient) + n2 = n2_gas(gas, pr) γLeff_this = n2*k0/(A0*a^2)*Leff(flength, a, λ0) γLeff_this - γLeff end @@ -85,30 +112,31 @@ function optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; t = transmission(flength, aopt, λ0) @debug(@sprintf("Results: a = %.3f μm, %.2f m long, transmission %.4f", 1e6aopt, flength, t)) - return aopt, flength, pr, t + return aopt, flength, pr, flow, t end function params_maxlength(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; thickness=1e-3, material=:SiO2, Bmax=0.2, S_sf=1.5, - entrance_window=true, exit_window=true, LIDT=2000, S_fluence=5) - - ρcrit = critical_density(gas, λ0, τfwhm_in, energy) - ρ = ρcrit/S_sf - pr = pressure(gas, ρ) - n2 = n2_gas(gas, pr) - + entrance_window=true, exit_window=true, LIDT=2000, S_fluence=5, + ngradient=0, max_pressure=60.0, max_flow=1.0, polarisation=:linear) _, P0 = T0P0(τfwhm_in, energy) - broadfac_req = τfwhm_in/τfwhm_out φnl_req = nonlinear_phase(broadfac_req) - + φnl_req *= (ngradient > 0 ? 3/2 : 1.0) + φnl_req *= (polarisation == :circular ? 3/2 : 1.0) k0 = 2π/λ0 A0 = Aeff0() function params(a) maxflength = max_flength(a, λ0, energy, τfwhm_in, maxlength; - thickness, material, Bmax, LIDT, S_fluence, - entrance_window, exit_window) + thickness, material, Bmax, LIDT, S_fluence, + entrance_window, exit_window) + ρcrit = critical_density(gas, λ0, τfwhm_in, energy) + ρ = ρcrit/S_sf + prm = pressure(gas, ρ) + pr, flow = fix_pressure(prm, a, maxflength, gas, max_pressure, max_flow, ngradient) + n2 = n2_gas(gas, pr) + γthis = n2*k0/(A0*a^2) Leff_req = φnl_req/P0/γthis maxLeff = Leff(maxflength, a, λ0) @@ -121,26 +149,27 @@ function params_maxlength(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; γLeff = γthis * Leff(flength, a, λ0) t = transmission.(flength, a, λ0) - intensity = P0/(A0*a^2) + intensity = P0/(A0*a^2) * (polarisation == :circular ? 2/3 : 1.0) φnl = P0*γLeff window_distance = (maxlength-maxflength)/2 beamsize_w0 = diverged_beam(a, λ0, window_distance) (;φnl=P0*γLeff, transmission=t, intensity, flength, broadening_factor=broadening_factor(φnl), pressure=pr, n2=n2, - P0=P0, window_distance, beamsize_w0) + P0=P0, window_distance, beamsize_w0, flow) end end function plot_optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; thickness=1e-3, material=:SiO2, Bmax=0.2, S_sf=1.5, S_ion=10, entrance_window=true, exit_window=true, LIDT=2000, S_fluence=5, - amin=25e-6, amax=500e-6, Na=512) + amin=25e-6, amax=500e-6, Na=512, + ngradient=0, max_pressure=60.0, max_flow=1.0, polarisation=:linear) factor = τfwhm_in/τfwhm_out f = params_maxlength(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; thickness, material, Bmax, S_sf, - entrance_window, exit_window, LIDT, S_fluence) - + entrance_window, exit_window, LIDT, S_fluence, ngradient, + max_pressure, max_flow, polarisation) Isupp = barrier_suppression_intensity(gas) A0 = Aeff0() @@ -152,21 +181,23 @@ function plot_optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; t = getindex.(params, :transmission) intensity = getindex.(params, :intensity) P0 = params[1].P0 - pr = params[1].pressure + pr = getindex.(params, :pressure) + flow = getindex.(params, :flow) t[flength .== 0] .= NaN try - global aopt, flopt, _, topt = optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; - thickness, material, Bmax, S_sf, S_ion, LIDT, S_fluence, entrance_window, exit_window) - global intopt = P0/(A0*aopt^2) + global aopt, flopt, propt, flowopt, topt = optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; + thickness, material, Bmax, S_sf, S_ion, LIDT, S_fluence, entrance_window, exit_window, + ngradient, max_pressure, max_flow, polarisation) + global intopt = P0/(A0*aopt^2) * (polarisation == :circular ? 2/3 : 1.0) catch e global aopt = missing end fig = plt.figure() - fig.set_size_inches(15, 4) - plt.subplot(1, 4, 1) + fig.set_size_inches(15, 7) + plt.subplot(2, 4, 1) plt.plot(1e6a, broadfac) if ~ismissing(aopt) plt.plot(1e6aopt, factor, "o"; color="k", label=@sprintf("%.1f μm", 1e6aopt)) @@ -179,7 +210,7 @@ function plot_optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; plt.ylabel("Broadening factor") - plt.subplot(1, 4, 2) + plt.subplot(2, 4, 2) plt.plot(1e6a, flength) if ~ismissing(aopt) plt.plot(1e6aopt, flopt, "o"; color="k", label=@sprintf("%.3f m", flopt)) @@ -190,7 +221,7 @@ function plot_optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; plt.xlabel("Core radius (μm)") plt.ylabel("Fibre length (m)") - plt.subplot(1, 4, 3) + plt.subplot(2, 4, 3) plt.plot(1e6a, 100*t) if ~ismissing(aopt) plt.plot(1e6aopt, 100*topt, "o"; color="k", label=@sprintf("%.1f %%", 100*topt)) @@ -201,7 +232,7 @@ function plot_optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; plt.xlabel("Core radius (μm)") plt.ylabel("Transmission (%)") - plt.subplot(1, 4, 4) + plt.subplot(2, 4, 4) plt.plot(1e6a, intensity*1e-4) plt.axhline(Isupp*1e-4/S_ion; linestyle="--", color="r", label="Limit") if ~ismissing(aopt) @@ -214,9 +245,33 @@ function plot_optimise(τfwhm_in, τfwhm_out, gas, λ0, energy, maxlength; plt.xlabel("Core radius (μm)") plt.ylabel("Intensity (W/cm\$^{-2}\$)") + plt.subplot(2, 4, 5) + plt.plot(1e6a, flow) + plt.axhline(max_flow; linestyle="--", color="r", label="Limit") + if ~ismissing(aopt) + plt.plot(1e6aopt, flowopt, "o"; color="k", label=@sprintf("%.3f slpm", flowopt)) + plt.legend() + end + plt.xlim(1e6.*extrema(a)) + plt.ylim(ymin=0) + plt.xlabel("Core radius (μm)") + plt.ylabel("Gas flow (slpm)") + + plt.subplot(2, 4, 6) + plt.plot(1e6a, pr) + plt.axhline(max_pressure; linestyle="--", color="r", label="Limit") + if ~ismissing(aopt) + plt.plot(1e6aopt, propt, "o"; color="k", label=@sprintf("%.3f bar", propt)) + plt.legend() + end + plt.xlim(1e6.*extrema(a)) + plt.ylim(ymin=0) + plt.xlabel("Core radius (μm)") + plt.ylabel("Pressure (bar)") + fig.tight_layout() plt.subplots_adjust(top=0.9) - plt.suptitle(@sprintf("Input peak power: %.2f GW | Pressure: %.2f bar", P0*1e-9, pr)) + plt.suptitle(@sprintf("Input peak power: %.2f GW", P0*1e-9)) return fig, f end diff --git a/src/HiSol.jl b/src/HiSol.jl index 932775d..49bb38c 100644 --- a/src/HiSol.jl +++ b/src/HiSol.jl @@ -5,8 +5,8 @@ include("HCF.jl") include("Solitons.jl") include("Limits.jl") include("Focusing.jl") -include("Compressor.jl") include("GasFlow.jl") +include("Compressor.jl") include("Design.jl") end