From 9d036ae608fcad40e5eb5d84b989aa372b939bab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= Date: Wed, 30 Jun 2021 09:37:49 +0200 Subject: [PATCH] started working on plasma --- README.md | 4 +- src/scgenerator/__init__.py | 2 +- src/scgenerator/physics/__init__.py | 45 +++++++++++++++----- src/scgenerator/physics/fiber.py | 2 +- src/scgenerator/physics/materials.py | 61 +++++++++++++++++++++++++++- src/scgenerator/physics/pulse.py | 58 ++++++++++++++++++++------ src/scgenerator/physics/simulate.py | 5 ++- src/scgenerator/plotting.py | 1 + 8 files changed, 149 insertions(+), 29 deletions(-) diff --git a/README.md b/README.md index 2a52923..d24b027 100644 --- a/README.md +++ b/README.md @@ -118,13 +118,13 @@ capillary_nested : int, optional gamma: float, optional unless beta is directly provided - nonlinear parameter in m^2 / W. Will overwrite any computed gamma parameter. + nonlinear parameter in m^-1*W^-1. Will overwrite any computed gamma parameter. effective_mode_diameter : float, optional effective mode field diameter in m n2 : float, optional - non linear refractive index + non linear refractive index in m^2/W A_eff : float, optional effective mode field area diff --git a/src/scgenerator/__init__.py b/src/scgenerator/__init__.py index e8e721e..097aeb5 100644 --- a/src/scgenerator/__init__.py +++ b/src/scgenerator/__init__.py @@ -10,6 +10,6 @@ from .math import abs2, argclosest, span from .physics import fiber, materials, pulse, simulate, units from .physics.simulate import RK4IP, new_simulation, resume_simulations from .plotting import plot_avg, plot_results_1D, plot_results_2D, plot_spectrogram -from .spectra import Pulse +from .spectra import Pulse, Spectrum from .utils.parameter import BareParams, BareConfig from . import utils, io, initialize, math diff --git a/src/scgenerator/physics/__init__.py b/src/scgenerator/physics/__init__.py index 392b603..6f5ed25 100644 --- a/src/scgenerator/physics/__init__.py +++ b/src/scgenerator/physics/__init__.py @@ -26,7 +26,8 @@ def material_dispersion( material: str, pressure=None, temperature=None, - ideal=False, + ideal=True, + safe=True, ): """returns the dispersion profile (beta_2) of a bulk material. @@ -41,7 +42,7 @@ def material_dispersion( 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 + whether to use the ideal gas law instead of the van der Waals equation, by default True Returns ------- @@ -50,6 +51,9 @@ def material_dispersion( """ w = units.m(wavelengths) + + order = np.argsort(w) + material_dico = load_material_dico(material) if ideal: n_gas_2 = materials.sellmeier(wavelengths, material_dico, pressure, temperature) + 1 @@ -59,12 +63,19 @@ def material_dispersion( ) 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)) + if safe: + disp = np.zeros(len(w)) + ind = w > 0 + disp[ind] = material_dispersion( + units.To.m(w[ind]), material, pressure, temperature, ideal, False + ) + return disp + else: + return fiber.beta2(w[order], np.sqrt(n_gas_2[order]))[order] def find_optimal_depth( - spectrum: T, w_c: np.ndarray, w0: float, material: str = "silica", max_z: float = 1.0 + spectrum: T, w_c: np.ndarray, w0: float, material: str, max_z: float = 1.0 ) -> tuple[T, pulse.OptimizeResult]: """finds the optimal silica depth to compress a pulse @@ -86,12 +97,26 @@ def find_optimal_depth( float distance in m """ - silica_disp = material_dispersion(units.To.m(w_c + w0), material) + w = w_c + w0 + disp = np.zeros(len(w)) + ind = w > (w0 / 10) + disp[ind] = material_dispersion(units.To.m(w[ind]), material) - propagate = lambda z: spectrum * np.exp(-0.5j * silica_disp * w_c ** 2 * z) + propagate = lambda z: spectrum * np.exp(-0.5j * disp * w_c ** 2 * z) + integrate = lambda z: math.abs2(np.fft.ifft(propagate(z))) def score(z): - return 1 / np.max(math.abs2(np.fft.ifft(propagate(z)))) + return -np.nansum(integrate(z) ** 6) - opti = minimize_scalar(score, bracket=(0, max_z)) - return opti.x + # import matplotlib.pyplot as plt + + # to_test = np.linspace(0, max_z, 200) + # scores = [score(z) for z in to_test] + # fig, ax = plt.subplots() + # ax.plot(to_test, scores / np.min(scores)) + # plt.show() + # plt.close(fig) + # ama = np.argmin(scores) + + opti = minimize_scalar(score, method="bounded", bounds=(0, max_z)) + return propagate(opti.x), opti diff --git a/src/scgenerator/physics/fiber.py b/src/scgenerator/physics/fiber.py index 8867aa2..d8b40b6 100644 --- a/src/scgenerator/physics/fiber.py +++ b/src/scgenerator/physics/fiber.py @@ -647,7 +647,7 @@ def PCF_dispersion(lambda_, pitch, ratio_d, w0=None, n2=None, A_eff=None): if A_eff is None: V_eff = pi2a / lambda_ * np.sqrt(n_co ** 2 - n_FSM2) w_eff = a_eff * (0.65 + 1.619 / V_eff ** 1.5 + 2.879 / V_eff ** 6) - A_eff = interp1d(lambda_, w_eff, kind="linear")(units.m.inv(w0)) ** pipi + A_eff = interp1d(lambda_, w_eff, kind="linear")(units.m.inv(w0)) ** 2 * pi if n2 is None: n2 = 2.6e-20 diff --git a/src/scgenerator/physics/materials.py b/src/scgenerator/physics/materials.py index f7bd06c..da091bc 100644 --- a/src/scgenerator/physics/materials.py +++ b/src/scgenerator/physics/materials.py @@ -1,8 +1,13 @@ +from typing import Any, Callable import numpy as np +import scipy.special +from scipy.integrate import cumulative_trapezoid + +from scgenerator import math from ..logger import get_logger from . import units -from .units import NA, c, kB, me, e +from .units import NA, c, kB, me, e, hbar def pressure_from_gradient(ratio, p0, p1): @@ -174,3 +179,57 @@ def non_linear_refractive_index(material_dico, pressure=None, temperature=None): def adiabadicity(w: np.ndarray, I: float, field: np.ndarray) -> np.ndarray: return w * np.sqrt(2 * me * I) / (e * np.abs(field)) + + +def free_electron_density( + t: np.ndarray, field: np.ndarray, N0: float, rate_func: Callable[[np.ndarray], np.ndarray] +) -> np.ndarray: + return N0 * (1 - np.exp(-cumulative_trapezoid(rate_func(field), t, initial=0))) + + +def ionization_rate_ADK( + ionization_energy: float, atomic_number +) -> Callable[[np.ndarray], np.ndarray]: + Z = -(atomic_number - 1) * e + + omega_p = ionization_energy / hbar + nstar = Z * np.sqrt(2.1787e-18 / ionization_energy) + omega_t = lambda field: e * np.abs(field) / np.sqrt(2 * me * ionization_energy) + Cnstar = 2 ** (2 * nstar) / (scipy.special.gamma(nstar + 1) ** 2) + omega_C = omega_p / 4 * math.abs2(Cnstar) + + def rate(field: np.ndarray) -> np.ndarray: + opt4 = 4 * omega_t(field) / om + return omega_C * (4 * omega_p / ot) ** (2 * nstar - 1) * np.exp(-4 * omega_p / (3 * ot)) + + return rate + + +class Plasma: + def __init__(self, t: np.ndarray, ionization_energy: float, atomic_number: int): + self.t = t + self.Ip = ionization_energy + self.atomic_number = atomic_number + self.rate = ionization_rate_ADK(self.Ip, self.atomic_number) + + def __call__(self, field: np.ndarray, N0: float) -> np.ndarray: + """returns the number density of free electrons as function of time + + Parameters + ---------- + field : np.ndarray + electric field in V/m + N0 : float + total number density of matter + + Returns + ------- + np.ndarray + number density of free electrons as function of time + """ + Ne = free_electron_density(self.t, field, N0, self.rate) + return cumulative_trapezoid( + np.gradient(Ne, self.t) * self.Ip / field, self.t, initial=0 + ) + e ** 2 / me * cumulative_trapezoid( + cumulative_trapezoid(Ne * field, self.t, initial=0), self.t, initial=0 + ) diff --git a/src/scgenerator/physics/pulse.py b/src/scgenerator/physics/pulse.py index 0b903a1..1538b1e 100644 --- a/src/scgenerator/physics/pulse.py +++ b/src/scgenerator/physics/pulse.py @@ -572,7 +572,7 @@ def peak_ind(values, mam=None): ) + 3 except IndexError: right_ind = len(values) - 1 - return left_ind, right_ind + return max(0, left_ind), min(len(values) - 1, right_ind) def setup_splines(x_axis, values, mam=None): @@ -685,8 +685,7 @@ def find_lobe_limits(x_axis, values, debug="", already_sorted=True): dd_roots, fwhm_pos, peak_pos, - folder_name, - file_name, + out_path, fig, ax, color, @@ -716,7 +715,7 @@ def find_lobe_limits(x_axis, values, debug="", already_sorted=True): c=color[5], ) ax.legend() - fig.savefig(os.path.join(folder_name, file_name), bbox_inches="tight") + fig.savefig(out_path, bbox_inches="tight") plt.close(fig) else: @@ -818,9 +817,7 @@ def _detailed_find_lobe_limits( # if measurement of the peak is not straightforward, we plot the situation to see # if the final measurement is good or not - folder_name, file_name, fig, ax = plot_setup( - file_name=f"it_{iterations}_{debug}", folder_name="measurements_errors_plots" - ) + out_path, fig, ax = plot_setup(out_path=f"measurement_errors_plots/it_{iterations}_{debug}") new_fwhm_pos = np.array([np.max(left_pos), np.min(right_pos)]) @@ -842,8 +839,7 @@ def _detailed_find_lobe_limits( dd_roots, fwhm_pos, peak_pos, - folder_name, - file_name, + out_path, fig, ax, color, @@ -946,7 +942,7 @@ def measure_field(t: np.ndarray, field: np.ndarray) -> Tuple[float, float, float def remove_2nd_order_dispersion( - spectrum: T, w_c: np.ndarray, beta2: float, max_z: float = -1.0 + spectrum: T, w_c: np.ndarray, beta2: float, max_z: float = -100.0 ) -> tuple[T, OptimizeResult]: """attempts to remove 2nd order dispersion from a complex spectrum @@ -964,11 +960,49 @@ def remove_2nd_order_dispersion( np.ndarray, shape (n, ) spectrum with 2nd order dispersion removed """ - # makeshift_t = np.linspace(0, 1, len(w_c)) propagate = lambda z: spectrum * np.exp(-0.5j * beta2 * w_c ** 2 * z) def score(z): - return 1 / np.max(abs2(np.fft.ifft(propagate(z)))) + return -np.max(abs2(np.fft.ifft(propagate(z)))) opti = minimize_scalar(score, bracket=(max_z, 0)) return propagate(opti.x), opti + + +def remove_2nd_order_dispersion2( + spectrum: T, w_c: np.ndarray, max_gdd: float = 1000e-30 +) -> tuple[T, OptimizeResult]: + """attempts to remove 2nd order dispersion from a complex spectrum + + 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) + + Returns + ------- + np.ndarray, shape (n, ) + spectrum with 2nd order dispersion removed + """ + propagate = lambda gdd: spectrum * np.exp(-0.5j * w_c ** 2 * 1e-30 * gdd) + integrate = lambda gdd: abs2(np.fft.ifft(propagate(gdd))) + + def score(gdd): + return -np.sum(integrate(gdd) ** 6) + + # def score(gdd): + # return -np.max(integrate(gdd)) + + # to_test = np.linspace(-max_gdd, max_gdd, 200) + # scores = [score(g) for g in to_test] + # fig, ax = plt.subplots() + # ax.plot(to_test, scores / np.min(scores)) + # plt.show() + # plt.close(fig) + # ama = np.argmin(scores) + + opti = minimize_scalar(score, bounds=(-max_gdd * 1e30, max_gdd * 1e30)) + opti["x"] *= 1e-30 + return propagate(opti.x * 1e30), opti diff --git a/src/scgenerator/physics/simulate.py b/src/scgenerator/physics/simulate.py index 07d395b..c071627 100644 --- a/src/scgenerator/physics/simulate.py +++ b/src/scgenerator/physics/simulate.py @@ -1,5 +1,6 @@ import multiprocessing import os +import random from datetime import datetime from pathlib import Path from typing import Dict, List, Tuple, Type @@ -704,7 +705,7 @@ def new_simulation( if prev_sim_dir is not None: config_dict["prev_sim_dir"] = str(prev_sim_dir) - task_id = np.random.randint(1e9, 1e12) + task_id = random.randint(1e9, 1e12) if prev_sim_dir is None: param_seq = initialize.ParamSequence(config_dict) @@ -718,7 +719,7 @@ def new_simulation( def resume_simulations(sim_dir: Path, method: Type[Simulations] = None) -> Simulations: - task_id = np.random.randint(1e9, 1e12) + task_id = random.randint(1e9, 1e12) config = io.load_toml(sim_dir / "initial_config.toml") io.set_data_folder(task_id, sim_dir) param_seq = initialize.RecoveryParamSequence(config, task_id) diff --git a/src/scgenerator/plotting.py b/src/scgenerator/plotting.py index 0f5b03c..103ac96 100644 --- a/src/scgenerator/plotting.py +++ b/src/scgenerator/plotting.py @@ -30,6 +30,7 @@ def plot_setup( - an axis """ out_path = defaults["name"] if out_path is None else out_path + out_path = Path(out_path) plot_name = out_path.name.replace(f".{file_type}", "") out_dir = out_path.resolve().parent