From 2bd5c3d055dcab48db7698b71d30498b7e29baa7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= Date: Mon, 28 Jun 2021 15:15:12 +0200 Subject: [PATCH] more utilities --- src/scgenerator/data/gas.toml | 91 ++++++++++++++------------- src/scgenerator/io.py | 2 +- src/scgenerator/physics/__init__.py | 97 +++++++++++++++++++++++++++++ src/scgenerator/physics/pulse.py | 72 ++++++++++++++------- 4 files changed, 193 insertions(+), 69 deletions(-) diff --git a/src/scgenerator/data/gas.toml b/src/scgenerator/data/gas.toml index 1e61685..0f4640c 100644 --- a/src/scgenerator/data/gas.toml +++ b/src/scgenerator/data/gas.toml @@ -1,51 +1,7 @@ -[nitrogen] -a = 0.137 -b = 1.709e-5 - -[helium] -a = 0.00346 -b = 2.38e-5 - -[helium_alt] -a = 0.00346 -b = 2.38e-5 - -[hydrogen] -a = 0.02453 -b = 2.651e-5 - -[neon] -a = 0.02135 -b = 1.709e-5 - -[argon] -a = 0.1355 -b = 3.201e-5 - -[argon_alt] -a = 0.1355 -b = 3.201e-5 - -[argon_alt2] -a = 0.1355 -b = 3.201e-5 - -[krypton] -a = 0.2349 -b = 3.978e-5 - -[xenon] -a = 0.425 -b = 5.105e-5 - [air] a = 0.1358 b = 3.64e-5 -[vacuum] -a = 0 -b = 0 - [air.sellmeier] B = [57921050000.0, 1679170000.0] C = [238018500000000.0, 57362000000000.0] @@ -56,7 +12,12 @@ kind = 2 [air.kerr] P0 = 101325 T0 = 293.15 -n2 = 3.01e-23 +n2 = 4.0e-23 +source = "Pigeon, J. J., Tochitsky, S. Y., Welch, E. C., & Joshi, C. (2016). Measurements of the nonlinear refractive index of air, N 2, and O 2 at 10 μm using four-wave mixing. Optics letters, 41(17), 3924-3927." + +[nitrogen] +a = 0.137 +b = 1.709e-5 [nitrogen.sellmeier] B = [32431570000.0] @@ -72,6 +33,10 @@ T0 = 273.15 n2 = 2.2e-23 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904." +[helium] +a = 0.00346 +b = 2.38e-5 + [helium.sellmeier] B = [2.16463842e-5, 2.10561127e-7, 4.7509272e-5] C = [-6.80769781e-16, 5.13251289e-15, 3.18621354e-15] @@ -86,6 +51,10 @@ T0 = 273.15 n2 = 3.1e-25 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904." +[helium_alt] +a = 0.00346 +b = 2.38e-5 + [helium_alt.sellmeier] B = [14755297000.0] C = [426297400000000.0] @@ -99,6 +68,10 @@ P0 = 30400.0 T0 = 273.15 n2 = 3.1e-25 +[hydrogen] +a = 0.02453 +b = 2.651e-5 + [hydrogen.sellmeier] B = [0.0148956, 0.0049037] C = [1.807e-10, 9.2e-11] @@ -113,6 +86,10 @@ T0 = 273.15 n2 = 6.36e-24 source = "Shelton, D. P., & Rice, J. E. (1994). Measurements and calculations of the hyperpolarizabilities of atoms and small molecules in the gas phase. Chemical Reviews, 94(1), 3-29" +[neon] +a = 0.02135 +b = 1.709e-5 + [neon.sellmeier] B = [1281450000.0, 22048600000.0] C = [184661000000000.0, 376840000000000.0] @@ -126,6 +103,10 @@ T0 = 273.15 n2 = 8.7e-25 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904." +[argon] +a = 0.1355 +b = 3.201e-5 + [argon.sellmeier] B = [2501410000.0, 500283000.0, 52234300000.0] C = [91012000000000.0, 87892000000000.0, 214020000000000.0] @@ -140,6 +121,10 @@ T0 = 273.15 n2 = 9.7e-24 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904." +[argon_alt] +a = 0.1355 +b = 3.201e-5 + [argon_alt.sellmeier] B = [0.0002033229, 0.0003445831] C = [2.0612e-16, 8.066e-15] @@ -154,6 +139,10 @@ T0 = 273.15 n2 = 9.7e-24 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904." +[argon_alt2] +a = 0.1355 +b = 3.201e-5 + [argon_alt2.sellmeier] B = [0.030182943] C = [144000000000000.0] @@ -169,6 +158,10 @@ T0 = 273.15 n2 = 9.7e-24 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904." +[krypton] +a = 0.2349 +b = 3.978e-5 + [krypton.sellmeier] B = [2536370000.0, 2736490000.0, 62080200000.0] C = [65474200000000.0, 73698000000000.0, 181080000000000.0] @@ -182,6 +175,10 @@ T0 = 273.15 n2 = 2.2e-23 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904." +[xenon] +a = 0.425 +b = 5.105e-5 + [xenon.sellmeier] B = [3228690000.0, 3553930000.0, 60676400000.0] C = [46301000000000.0, 59578000000000.0, 112740000000000.0] @@ -195,6 +192,10 @@ T0 = 273.15 n2 = 5.8e-23 source = "Wahlstrand, J. K., Cheng, Y. H., & Milchberg, H. M. (2012). High field optical nonlinearity and the Kramers-Kronig relations. Physical review letters, 109(11), 113904." +[vacuum] +a = 0 +b = 0 + [vacuum.sellmeier] B = [] C = [] diff --git a/src/scgenerator/io.py b/src/scgenerator/io.py index 1018f16..83f4b71 100644 --- a/src/scgenerator/io.py +++ b/src/scgenerator/io.py @@ -162,7 +162,7 @@ def load_config(path: os.PathLike) -> BareConfig: return BareConfig(**config) -def load_material_dico(name): +def load_material_dico(name: str) -> dict[str, Any]: """loads a material dictionary Parameters ---------- diff --git a/src/scgenerator/physics/__init__.py b/src/scgenerator/physics/__init__.py index e69de29..392b603 100644 --- a/src/scgenerator/physics/__init__.py +++ b/src/scgenerator/physics/__init__.py @@ -0,0 +1,97 @@ +""" +This file contains functions that don't properly fit into the 'fiber' or 'pulse' category. +They're also not necessarily 'pure' function as they do some io and stuff. +""" + +from typing import TypeVar + +import numpy as np +from scipy.optimize import minimize_scalar + +from ..io import load_material_dico +from .. import math +from . import fiber, materials, units, pulse + +T = TypeVar("T") + + +def group_delay_to_gdd(wavelength: np.ndarray, group_delay: np.ndarray) -> np.ndarray: + w = units.m.inv(wavelength) + gdd = np.gradient(group_delay, w) + return gdd + + +def material_dispersion( + wavelengths, + material: str, + pressure=None, + temperature=None, + ideal=False, +): + """returns the dispersion profile (beta_2) of a bulk material. + + Parameters + ---------- + wavelengths : ndarray, shape (n, ) + wavelengths over which to calculate the dispersion + material : str + material name in lower case + temperature : float, optional + Temperature of the material + pressure : float, optional + constant pressure + ideal : bool, optional + whether to use the ideal gas law instead of the van der Waals equation, by default False + + Returns + ------- + out : ndarray, shape (n, ) + beta2 as function of wavelength + """ + + w = units.m(wavelengths) + material_dico = load_material_dico(material) + if ideal: + n_gas_2 = materials.sellmeier(wavelengths, material_dico, pressure, temperature) + 1 + else: + N_1 = materials.number_density_van_der_waals( + pressure=pressure, temperature=temperature, material_dico=material_dico + ) + N_0 = materials.number_density_van_der_waals(material_dico=material_dico) + n_gas_2 = materials.sellmeier(wavelengths, material_dico) * N_1 / N_0 + 1 + + return fiber.beta2(w, np.sqrt(n_gas_2)) + + +def find_optimal_depth( + spectrum: T, w_c: np.ndarray, w0: float, material: str = "silica", max_z: float = 1.0 +) -> tuple[T, pulse.OptimizeResult]: + """finds the optimal silica depth to compress a pulse + + Parameters + ---------- + spectrum : np.ndarray or Spectrum, shape (n, ) + spectrum from which to remove 2nd order dispersion + w_c : np.ndarray, shape (n, ) + corresponding centered angular frequencies (w-w0) + w0 : float + pump central angular frequency + material : str + material name, by default 'silica' + max_z : float + maximum propagation distance in m + + Returns + ------- + float + distance in m + """ + silica_disp = material_dispersion(units.To.m(w_c + w0), material) + + propagate = lambda z: spectrum * np.exp(-0.5j * silica_disp * w_c ** 2 * z) + + def score(z): + return 1 / np.max(math.abs2(np.fft.ifft(propagate(z)))) + + opti = minimize_scalar(score, bracket=(0, max_z)) + return opti.x diff --git a/src/scgenerator/physics/pulse.py b/src/scgenerator/physics/pulse.py index 207d973..0b903a1 100644 --- a/src/scgenerator/physics/pulse.py +++ b/src/scgenerator/physics/pulse.py @@ -13,6 +13,7 @@ import itertools import os from pathlib import Path from typing import Literal, Tuple, TypeVar +from collections import namedtuple import matplotlib.pyplot as plt import numpy as np @@ -48,6 +49,11 @@ P0T0_to_E0_fac = dict( ) """relates the total energy (amplitue^2) to the t0 parameter of the amplitude and the peak intensity (peak_amplitude^2)""" +PulseProperties = namedtuple( + "PulseProperties", + "quality mean_coherence fwhm_noise mean_fwhm peak_rin energy_rin timing_jitter", +) + def initial_field(t, shape, t0, peak_power): """returns the initial field @@ -475,12 +481,14 @@ def g12(values): # Create all the possible pairs of values n = len(values) field_pairs = itertools.combinations(values, 2) + mean_spec = np.mean(abs2(values), axis=0) + mask = mean_spec > 1e-15 * mean_spec.max() corr = np.zeros_like(values[0]) for pair in field_pairs: - corr += pair[0].conj() * pair[1] - g12_arr = corr / (n * (n - 1) / 2 * np.mean(abs2(values), axis=0)) + corr[mask] += pair[0][mask].conj() * pair[1][mask] + corr[mask] = corr[mask] / (n * (n - 1) / 2 * mean_spec[mask]) - return np.abs(g12_arr) + return np.abs(corr) def avg_g12(values): @@ -842,39 +850,45 @@ def _detailed_find_lobe_limits( ) -def measure_properties(spectra, t, compress=True, debug=""): +def measure_properties(spectra, t, compress=True, return_limits=False, debug="") -> PulseProperties: """measure the quality factor, the fwhm variation, the peak power variation, Parameters ---------- - spectra : 2D array - set of n spectra in sc-ordering that differ only by noise - t : 1D array - time axis of the simulation - compress : bool, optional - whether to perform pulse compression. Default value is True, but this - should be set to False to measure the initial pulse as output by gaussian_pulse - or sech_pulse because compressing it would result in glitches and wrong measurements + spectra : np.ndarray, shape (n, nt) + set of n spectra on a grid of nt angular frequency points + t : np.ndarray, shape (nt, ) + time axis of the simulation + compress : bool, optional + whether to perform pulse compression. Default value is True, but this + should be set to False to measure the initial pulse as output by gaussian_pulse + or sech_pulse because compressing it would result in glitches and wrong measurements + return_limits : bool, optional + return the time values of the limits Returns ---------- - qf : float + PulseProperties object: + quality : float quality factor of the pulse ensemble - mean_g12 : float + mean_gmean_coherence12 : float mean coherence of the spectra ensemble - fwhm_var : float - relative noise in temporal width of the compressed pulse - fwhm_abs : float - width of the mean compressed pulse - int_var : flaot - relative noise in the compressed pulse peak intensity - t_jitter : float + fwhm_noise : float + relative noise in temporal width of the (compressed) pulse + mean_fwhm : float + width of the mean (compressed) pulse + peak_rin : float + relative noise in the (compressed) pulse peak intensity + energy_rin : float + relative noise in the (compressed) pulse total_energy + timing_jitter : float standard deviantion in absolute temporal peak position + all_limits : list[tuple[np.ndarray, np.ndarray]], only if return_limits = True + list of tuples of the form ([left_lobe_lim, right_lobe_lim, lobe_pos], [left_hm, right_hm]) """ if compress: fields = abs2(compress_pulse(spectra)) else: - print("Skipping compression") fields = abs2(ifft(spectra)) field = np.mean(fields, axis=0) @@ -897,16 +911,28 @@ def measure_properties(spectra, t, compress=True, debug=""): P0 = [] fwhm = [] t_offset = [] + energies = [] + all_lims: list[tuple[np.ndarray, np.ndarray]] = [] for f in fields: lobe_lim, fwhm_lim, _, big_spline = find_lobe_limits(t, f, debug) + all_lims.append((lobe_lim, fwhm_lim)) P0.append(big_spline(lobe_lim[2])) fwhm.append(length(fwhm_lim)) t_offset.append(lobe_lim[2]) + energies.append(np.trapz(fields, t)) fwhm_var = np.std(fwhm) / np.mean(fwhm) int_var = np.std(P0) / np.mean(P0) + en_var = np.std(energies) / np.mean(energies) t_jitter = np.std(t_offset) - return qf, mean_g12, fwhm_var, fwhm_abs, int_var, t_jitter + if isinstance(mean_g12, np.ndarray) and mean_g12.ndim == 0: + mean_g12 = mean_g12[()] + pp = PulseProperties(qf, mean_g12, fwhm_var, fwhm_abs, int_var, en_var, t_jitter) + + if return_limits: + return pp, all_lims + else: + return pp def measure_field(t: np.ndarray, field: np.ndarray) -> Tuple[float, float, float]: