dedebenui/scgenerator
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

replaced material_dico with Gas object

Browse Source
This commit is contained in:
Benoît Sierro
2023-01-24 11:22:53 +01:00
parent 239eaddd7c
commit 1baef68656
24 changed files with 212 additions and 613 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/scgenerator/__init__.py
View File

@@ -1,10 +1,11 @@
# flake8: noqa
from scgenerator import math, operators
from scgenerator.evaluator import Evaluator
from scgenerator.helpers import *
from scgenerator.legacy import convert_sim_folder
from scgenerator.math import abs2, argclosest, normalized, span, tspace, wspace
from scgenerator.parameter import FileConfiguration, Parameters
from scgenerator.physics import fiber, materials, pulse, simulate, units, plasma
from scgenerator.physics import fiber, materials, plasma, pulse, simulate, units
from scgenerator.physics.simulate import RK4IP, parallel_RK4IP, run_simulation
from scgenerator.physics.units import PlotRange
from scgenerator.plotting import (
@@ -25,4 +26,3 @@ from scgenerator.variationer import (
Variationer,
VariationSpecsError,
)
from scgenerator.helpers import *

1
src/scgenerator/cache.py
View File

@@ -1,6 +1,7 @@
from collections import namedtuple
from copy import copy
from functools import wraps
import numpy as np
CacheInfo = namedtuple("CacheInfo", "hits misses size")

4
src/scgenerator/const.py
View File

@@ -1,7 +1,7 @@
__version__ = "0.2.6dev"
from typing import Any
__version__ = "0.2.7dev"
ONE_2 = 1 / 2
ONE_3 = 1 / 3
ONE_FOURTH = 1 / 4

2
src/scgenerator/defaults.py
View File

@@ -1,6 +1,6 @@
import matplotlib.pyplot as plt
from pathlib import Path
import matplotlib.pyplot as plt
default_plotting = dict(
figsize=(10, 7),

1
src/scgenerator/env.py
View File

@@ -2,7 +2,6 @@ import os
from pathlib import Path
from typing import Any, Dict, Literal, Optional, Set
ENVIRON_KEY_BASE = "SCGENERATOR_"
TMP_FOLDER_KEY_BASE = ENVIRON_KEY_BASE + "SC_TMP_"
PREFIX_KEY_BASE = ENVIRON_KEY_BASE + "PREFIX_"

13
src/scgenerator/evaluator.py
View File

@@ -5,12 +5,11 @@ from typing import Any, Callable, Optional, Union
import numpy as np
from . import math, operators, utils
from .const import MANDATORY_PARAMETERS
from .errors import EvaluatorError, NoDefaultError
from .physics import fiber, materials, pulse, units
from .utils import _mock_function, func_rewrite, get_arg_names, get_logger
from scgenerator import math, operators, utils
from scgenerator.const import MANDATORY_PARAMETERS
from scgenerator.errors import EvaluatorError, NoDefaultError
from scgenerator.physics import fiber, materials, pulse, units
from scgenerator.utils import _mock_function, func_rewrite, get_arg_names, get_logger
class Rule:
@@ -336,7 +335,7 @@ default_rules: list[Rule] = [
# Fiber Dispersion
Rule("w_for_disp", units.m, ["wl_for_disp"]),
Rule("hr_w", fiber.delayed_raman_w),
Rule("gas_info", materials.GasInfo),
Rule("gas_info", materials.Gas),
Rule("chi_gas", lambda gas_info, wl_for_disp: gas_info.sellmeier.chi(wl_for_disp)),
Rule("n_gas_2", materials.n_gas_2),
Rule("n_eff", fiber.n_eff_hasan, conditions=dict(model="hasan")),

6
src/scgenerator/helpers.py
View File

@@ -5,7 +5,7 @@ series of helper functions
import numpy as np
from scgenerator.math import all_zeros
from scgenerator.physics.fiber import beta2, n_eff_marcatili, n_eff_hasan
from scgenerator.physics.fiber import beta2, n_eff_hasan, n_eff_marcatili
from scgenerator.physics.materials import n_gas_2
from scgenerator.physics.units import c
@@ -40,7 +40,7 @@ def capillary_dispersion(
pressure = 101325
if temperature is None:
temperature = 293.15
n = n_eff_marcatili(wl, n_gas_2(wl, gas_name.lower(), pressure, temperature, False), radius)
n = n_eff_marcatili(wl, n_gas_2(wl, gas_name.lower(), pressure, temperature), radius)
w = 2 * np.pi * c / wl
return beta2(w, n)[2:-2]
@@ -117,7 +117,7 @@ def revolver_dispersion(
temperature = 293.15
n = n_eff_hasan(
wl,
n_gas_2(wl, gas_name.lower(), pressure, temperature, False),
n_gas_2(wl, gas_name.lower(), pressure, temperature),
core_radius,
capillary_num,
capillary_nested,

10
src/scgenerator/legacy.py
View File

@@ -8,11 +8,11 @@ import numpy as np
import tomli
import tomli_w
from .const import SPEC1_FN, SPEC1_FN_N, SPECN_FN1
from .parameter import FileConfiguration, Parameters
from .pbar import PBars
from .utils import save_parameters
from .variationer import VariationDescriptor
from scgenerator.const import SPEC1_FN, SPEC1_FN_N, SPECN_FN1
from scgenerator.parameter import FileConfiguration, Parameters
from scgenerator.pbar import PBars
from scgenerator.utils import save_parameters
from scgenerator.variationer import VariationDescriptor
def load_config(path: os.PathLike) -> dict[str, Any]:

2
src/scgenerator/logger.py
View File

@@ -1,6 +1,6 @@
import logging
from .env import log_file_level, log_print_level
from scgenerator.env import log_file_level, log_print_level
lvl_map: dict[str, int] = dict(

4
src/scgenerator/math.py
View File

@@ -7,7 +7,7 @@ from scipy.interpolate import interp1d, lagrange
from scipy.special import jn_zeros
from functools import cache
from .cache import np_cache
from scgenerator.cache import np_cache
pi = np.pi
c = 299792458.0
@@ -83,7 +83,7 @@ def _power_fact_array(x, n):
def abs2(z: np.ndarray) -> np.ndarray:
return z.real ** 2 + z.imag ** 2
return z.real**2 + z.imag**2
def normalized(z: np.ndarray) -> np.ndarray:

52
src/scgenerator/operators.py
View File

@@ -6,15 +6,14 @@ from __future__ import annotations
from abc import ABC, abstractmethod
from dataclasses import dataclass
from typing import Any, Callable
from typing import Callable
import numpy as np
from scipy.interpolate import interp1d
from . import math
from .logger import get_logger
from .physics import fiber, materials, pulse, units, plasma
from .utils import load_material_dico
from scgenerator import math
from scgenerator.logger import get_logger
from scgenerator.physics import fiber, materials, plasma, pulse, units
class CurrentState:
@@ -212,8 +211,7 @@ class NoOpFreq(Operator):
class AbstractGas(Operator):
gas_name: str
material_dico: dict[str, Any]
gas: materials.Gas
def __call__(self, state: CurrentState) -> np.ndarray:
return self.square_index(state)
@@ -274,18 +272,10 @@ class ConstantGas(AbstractGas):
ideal_gas: bool,
wl_for_disp: np.ndarray,
):
self.material_dico = load_material_dico(gas_name)
self.gas_name = gas_name
self.gas = materials.Gas(gas_name)
self.pressure_const = pressure
if ideal_gas:
self.number_density_const = materials.number_density_van_der_waals(
pressure=pressure, temperature=temperature, material_dico=self.material_dico
)
else:
self.number_density_const = self.pressure_const / (units.kB * temperature)
self.n_gas_2_const = materials.n_gas_2(
wl_for_disp, gas_name, self.pressure_const, temperature, ideal_gas
)
self.number_density_const = self.gas.number_density(temperature, pressure, ideal_gas)
self.n_gas_2_const = self.gas.sellmeier.n_gas_2(wl_for_disp, temperature, pressure)
def pressure(self, state: CurrentState = None) -> float:
return self.pressure_const
@@ -311,10 +301,9 @@ class PressureGradientGas(AbstractGas):
ideal_gas: bool,
wl_for_disp: np.ndarray,
):
self.gas_name = gas_name
self.p_in = pressure_in
self.p_out = pressure_out
self.material_dico = load_material_dico(gas_name)
self.gas = materials.Gas(gas_name)
self.temperature = temperature
self.ideal_gas = ideal_gas
self.wl_for_disp = wl_for_disp
@@ -323,23 +312,10 @@ class PressureGradientGas(AbstractGas):
return materials.pressure_from_gradient(state.z_ratio, self.p_in, self.p_out)
def number_density(self, state: CurrentState) -> float:
if self.ideal_gas:
return self.pressure(state) / (units.kB * self.temperature)
else:
return materials.number_density_van_der_waals(
pressure=self.pressure(state),
temperature=self.temperature,
material_dico=self.material_dico,
)
return self.gas.number_density(self.temperature, self.pressure(state), self.ideal_gas)
def square_index(self, state: CurrentState) -> np.ndarray:
return materials.fast_n_gas_2(
self.wl_for_disp,
self.pressure(state),
self.temperature,
self.ideal_gas,
self.material_dico,
)
return self.gas.sellmeier.n_gas_2(self.wl_for_disp, self.temperature, self.pressure(state))
##################################################
@@ -829,9 +805,7 @@ class VariableScalarGamma(AbstractGamma):
self.arr = np.ones(t_num)
def __call__(self, state: CurrentState) -> np.ndarray:
n2 = materials.non_linear_refractive_index(
self.gas_op.material_dico, self.gas_op.pressure(state), self.temperature
)
n2 = self.gas_op.gas.n2(self.temperature, self.gas_op.pressure(state))
return self.arr * fiber.gamma_parameter(n2, self.w0, self.A_eff)
@@ -847,7 +821,7 @@ class Plasma(Operator):
def __init__(self, dt: float, w: np.ndarray, gas_op: AbstractGas):
self.gas_op = gas_op
self.mat_plasma = plasma.Plasma(dt, self.gas_op.material_dico["ionization_energy"])
self.mat_plasma = plasma.Plasma(dt, self.gas_op.gas["ionization_energy"])
self.factor_out = -w / (2.0 * units.c**2 * units.epsilon0)
def __call__(self, state: CurrentState) -> np.ndarray:

20
src/scgenerator/parameter.py
View File

@@ -9,19 +9,19 @@ from dataclasses import dataclass, field, fields
from functools import lru_cache, wraps
from math import isnan
from pathlib import Path
from typing import Any, Callable, ClassVar, Iterable, Iterator, Set, Type, TypeVar, Union
from typing import Any, Callable, ClassVar, Iterable, Iterator, Set, Type, TypeVar
import numpy as np
from . import env, legacy, utils
from .const import MANDATORY_PARAMETERS, PARAM_FN, VALID_VARIABLE, __version__
from .errors import EvaluatorError
from .evaluator import Evaluator
from .logger import get_logger
from .operators import AbstractConservedQuantity, LinearOperator, NonLinearOperator
from .solver import Integrator
from .utils import DebugDict, fiber_folder, update_path_name
from .variationer import VariationDescriptor, Variationer
from scgenerator import env, legacy, utils
from scgenerator.const import MANDATORY_PARAMETERS, PARAM_FN, VALID_VARIABLE, __version__
from scgenerator.errors import EvaluatorError
from scgenerator.evaluator import Evaluator
from scgenerator.logger import get_logger
from scgenerator.operators import AbstractConservedQuantity, LinearOperator, NonLinearOperator
from scgenerator.solver import Integrator
from scgenerator.utils import DebugDict, fiber_folder, update_path_name
from scgenerator.variationer import VariationDescriptor, Variationer
T = TypeVar("T")

2
src/scgenerator/pbar.py
View File

@@ -10,7 +10,7 @@ from typing import Iterable, Union
from tqdm import tqdm
from .env import pbar_policy
from scgenerator.env import pbar_policy
T_ = typing.TypeVar("T_")

44
src/scgenerator/physics/__init__.py
View File

@@ -8,10 +8,9 @@ from typing import TypeVar
import numpy as np
from scipy.optimize import minimize_scalar
from .. import math
from . import fiber, materials, units, pulse
from ..cache import np_cache
from ..utils import load_material_dico
from scgenerator import math
from scgenerator.cache import np_cache
from scgenerator.physics import fiber, materials, pulse, units
T = TypeVar("T")
@@ -28,8 +27,6 @@ def material_dispersion(
material: str,
pressure=None,
temperature=None,
ideal=True,
safe=True,
):
"""returns the dispersion profile (beta_2) of a bulk material.
@@ -54,26 +51,11 @@ def material_dispersion(
w = units.m(wavelengths)
if safe:
disp = np.zeros(len(w))
ind = w > 0
disp[ind] = material_dispersion(
units.m.inv(w[ind]), material, pressure, temperature, ideal, False
)
return disp
else:
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
order = np.argsort(w)
unorder = np.argsort(order)
return fiber.beta2(w[order], np.sqrt(n_gas_2[order]))[unorder]
sellmeier = materials.Sellmeier.load(material)
n_gas_2 = sellmeier.n_gas_2(wavelengths, temperature, pressure)
order = np.argsort(w)
unorder = np.argsort(order)
return fiber.beta2(w[order], np.sqrt(n_gas_2[order]))[unorder]
def find_optimal_depth(
@@ -105,7 +87,7 @@ def find_optimal_depth(
disp[ind] = material_dispersion(units.m.inv(w[ind]), material)
def propagate(z):
return spectrum * np.exp(-0.5j * disp * w_c ** 2 * z)
return spectrum * np.exp(-0.5j * disp * w_c**2 * z)
def integrate(z):
return math.abs2(np.fft.ifft(propagate(z)))
@@ -118,7 +100,7 @@ def find_optimal_depth(
def propagate_field(
t: np.ndarray, field: np.ndarray, z: float, material: str, center_wl_nm: float = 1540.0
t: np.ndarray, field: np.ndarray, z: float, material: str, center_wl_nm: float
) -> np.ndarray:
"""propagates a field through bulk material
@@ -132,8 +114,8 @@ def propagate_field(
distance to propagate in m
material : str
material name
center_wl_nm : float, optional
center wavelength of the grid in nm, by default 1540
center_wl_nm : float
center wavelength of the grid in nm
Returns
-------
@@ -143,4 +125,4 @@ def propagate_field(
w_c = math.wspace(t)
l = units.m(w_c + units.nm(center_wl_nm))
disp = material_dispersion(l, material)
return np.fft.ifft(np.fft.fft(field) * np.exp(0.5j * disp * w_c ** 2 * z))
return np.fft.ifft(np.fft.fft(field) * np.exp(0.5j * disp * w_c**2 * z))

249
src/scgenerator/physics/fiber.py
View File

@@ -1,16 +1,17 @@
from typing import Any, Iterable, Literal, TypeVar
from typing import Iterable, TypeVar
import numpy as np
from numpy import e
from numpy.fft import fft
from numpy.polynomial.chebyshev import Chebyshev, cheb2poly
from scipy.interpolate import interp1d
from scgenerator import utils
from scgenerator.cache import np_cache
from scgenerator.math import argclosest, u_nm
from scgenerator.physics import materials as mat
from scgenerator.physics import units
from scgenerator.physics.units import c, pi
from scipy.interpolate import interp1d
pipi = 2 * pi
T = TypeVar("T")
@@ -255,7 +256,7 @@ def resonance_thickness(
float
thickness in m
"""
n_si_2 = mat.n_gas_2(wl_for_disp, "silica", None, None, True)
n_si_2 = mat.n_gas_2(wl_for_disp, "silica", None, None)
return (
wl_for_disp
/ (4 * np.sqrt(n_si_2))
@@ -268,7 +269,7 @@ def resonance_strength(
wl_for_disp: np.ndarray, core_radius: float, capillary_thickness: float, order: int
) -> float:
a = 1.83 + (2.3 * capillary_thickness / core_radius)
n_si = np.sqrt(mat.n_gas_2(wl_for_disp, "silica", None, None, True))
n_si = np.sqrt(mat.n_gas_2(wl_for_disp, "silica", None, None))
return (
capillary_thickness
/ (n_si * core_radius) ** (2.303 * a / n_si)
@@ -345,7 +346,7 @@ def n_eff_hasan(
n_eff_2 = n_gas_2 - (u * wl_for_disp / (pipi * R_eff)) ** 2
chi_sil = mat.sellmeier(wl_for_disp, utils.load_material_dico("silica"))
chi_sil = mat.Sellmeier.load("silica").chi(wl_for_disp)
with np.errstate(divide="ignore", invalid="ignore"):
for m, strength in enumerate(capillary_resonance_strengths):
@@ -466,143 +467,6 @@ def A_eff_from_V(core_radius: float, V_eff: T) -> T:
return out
def HCPCF_find_with_given_ZDW(
variable: Literal["pressure", "temperature"],
target: float,
search_range: tuple[float, float],
material_dico: dict[str, Any],
model="marcatili",
model_params={},
pressure=None,
temperature=None,
ideal=False,
):
"""finds the parameters (pressure or temperature) to yield the target ZDW. assign the string value 'vary' to the parameter
Parameters
----------
variable : str {"pressure", "temperature"}
which parameter to vary
target : float
the ZDW target, in m
search_range : array, shape (2,)
(min, max) of the search range
other parameters : see HCPCF_dispersion. Pressure or temperature is used as initial value if it is variable
Returns
-------
the parameter that satisfies the ZDW
"""
from scipy import optimize
l_search = [120e-9, 6000e-9]
#
fixed = [material_dico, model, model_params, ideal]
if variable == "pressure":
fixed.append(temperature)
x0 = 1e5 if pressure is None else pressure
def zdw(x, *args):
current_ZDW = HCPF_ZDW(
l_search,
args[0],
model=args[1],
model_params=args[2],
pressure=x,
temperature=args[4],
ideal=args[3],
)
out = current_ZDW - target
return out
elif variable == "temperature":
fixed.append(pressure)
x0 = 273.15 if temperature is None else temperature
def zdw(x, *args):
current_ZDW = HCPF_ZDW(
l_search,
args[0],
model=args[1],
model_params=args[2],
pressure=args[4],
temperature=x,
ideal=args[3],
)
out = current_ZDW - target
return out
else:
raise AttributeError(f"'variable' arg must be 'pressure' or 'temperature', not {variable}")
optimized = optimize.root_scalar(
zdw, x0=x0, args=tuple(fixed), method="brentq", bracket=search_range
)
return optimized.root
def HCPF_ZDW(
search_range,
material_dico,
model="marcatili",
model_params={},
pressure=None,
temperature=None,
ideal=False,
max_iter=10,
threshold=1e-36,
):
"""finds one Zero Dispersion Wavelength (ZDW) of a given HC-PCF fiber
Parameters
----------
see HCPCF_dispersion for description of most arguments
max_iter : float
How many iterations are allowed at most to reach the threashold
threshold : float
upper bound of what counts as beta2 == 0 (in si units)
Returns
-------
float:
the ZDW in m
"""
prev_find = np.inf
l = np.linspace(*search_range, 50)
core_radius = model_params["core_radius"]
zdw_ind = 0
for i in range(max_iter):
beta2 = HCPCF_dispersion(
l,
material_dico,
model=model,
model_params=model_params,
pressure=pressure,
temperature=temperature,
ideal=ideal,
)
zdw_ind = argclosest(beta2, 0)
if beta2[zdw_ind] < threshold:
break
elif beta2[zdw_ind] < prev_find:
l = np.linspace(
l[zdw_ind] - (100 / (i + 1)) * 1e-9, l[zdw_ind] + (100 / (i + 1)) * 1e-9, 50
)
prev_find = beta2[zdw_ind]
else:
raise RuntimeError(
f"Could not find a ZDW with parameters {1e6*core_radius} um, {1e-5 * pressure} bar, {temperature} K."
)
else:
print(f"Could not get to threshold in {max_iter} iterations")
return l[zdw_ind]
def beta(w_for_disp: np.ndarray, n_eff: np.ndarray) -> np.ndarray:
return n_eff * w_for_disp / c
@@ -634,12 +498,11 @@ def frame_velocity(beta1_arr: np.ndarray, w0_ind: int) -> float:
def HCPCF_dispersion(
wl_for_disp,
material_dico=None,
gas_name: str,
model="marcatili",
model_params={},
pressure=None,
temperature=None,
ideal=False,
**model_params,
):
"""returns the dispersion profile (beta_2) of a hollow-core photonic crystal fiber.
@@ -647,8 +510,8 @@ def HCPCF_dispersion(
----------
wl_for_disp : ndarray, shape (n, )
wavelengths over which to calculate the dispersion
material_dico : dict
material dictionary respecting standard format explained in FIXME
gas_name : str
name of the filling gas in lower case
model : string {"marcatili", "marcatili_adjusted", "hasan"}
which model of effective refractive index to use
model_params : tuple
@@ -666,17 +529,7 @@ def HCPCF_dispersion(
"""
w = units.m(wl_for_disp)
if material_dico is None:
n_gas_2 = np.ones_like(wl_for_disp)
else:
if ideal:
n_gas_2 = mat.sellmeier(wl_for_disp, material_dico, pressure, temperature) + 1
else:
N_1 = mat.number_density_van_der_waals(
pressure=pressure, temperature=temperature, material_dico=material_dico
)
N_0 = mat.number_density_van_der_waals(material_dico=material_dico)
n_gas_2 = mat.sellmeier(wl_for_disp, material_dico) * N_1 / N_0 + 1
n_gas_2 = mat.Sellmeier.load(gas_name).n_gas_2(wl_for_disp, temperature, pressure)
n_eff_func = dict(
marcatili=n_eff_marcatili, marcatili_adjusted=n_eff_marcatili_adjusted, hasan=n_eff_hasan
@@ -686,81 +539,6 @@ def HCPCF_dispersion(
return beta2(w, n_eff)
def dynamic_HCPCF_dispersion(
wl_for_disp: np.ndarray,
pressure_values: list[float],
core_radius: float,
fiber_model: str,
model_params: dict[str, Any],
temperature: float,
ideal_gas: bool,
w0: float,
interp_range: tuple[float, float],
material_dico: dict[str, Any],
deg: int,
):
"""returns functions for beta2 coefficients and gamma instead of static values
Parameters
----------
wl_for_disp : wavelength array
params : dict
flattened parameter dictionary
material_dico : dict
material dictionary (see README for details)
Returns
-------
beta2_coef : func(r), r is the relative position in the fiber
a function that returns an array of coefficients as function of the relative position in the fiber
to be used in disp_op
gamma : func(r), r is the relative position in the fiber
a function that returns a float corresponding to the nonlinear parameter at the relative position
in the fiber
"""
A_eff = 1.5 * core_radius**2
# defining function instead of storing every possilble value
def pressure(r):
return mat.pressure_from_gradient(r, *pressure_values)
def beta2(r):
return HCPCF_dispersion(
wl_for_disp,
core_radius,
material_dico,
fiber_model,
model_params,
pressure(r),
temperature,
ideal_gas,
)
def n2(r):
return mat.non_linear_refractive_index(material_dico, pressure(r), temperature)
ratio_range = np.linspace(0, 1, 256)
gamma_grid = np.array([gamma_parameter(n2(r), w0, A_eff) for r in ratio_range])
gamma_interp = interp1d(ratio_range, gamma_grid)
beta2_grid = np.array(
[dispersion_coefficients(wl_for_disp, beta2(r), w0, interp_range, deg) for r in ratio_range]
)
beta2_interp = [
interp1d(ratio_range, beta2_grid[:, i], assume_sorted=True) for i in range(deg + 1)
]
def beta2_func(r):
return [beta2_interp[i](r)[()] for i in range(deg + 1)]
def gamma_func(r):
return gamma_interp(r)[()]
return beta2_func, gamma_func
def gamma_parameter(n2: float, w0: float, A_eff: T) -> T:
if isinstance(A_eff, np.ndarray):
A_eff_term = np.sqrt(A_eff * A_eff[0])
@@ -815,8 +593,7 @@ def n_eff_pcf(wl_for_disp: np.ndarray, pitch: float, pitch_ratio: float) -> np.n
n_eff2 = (wl_for_disp * W / (pi2a)) ** 2 + n_FSM2
n_eff = np.sqrt(n_eff2)
material_dico = utils.load_material_dico("silica")
chi_mat = mat.sellmeier(wl_for_disp, material_dico)
chi_mat = mat.Sellmeier.load("silica").chi(wl_for_disp)
return n_eff + np.sqrt(chi_mat + 1)
@@ -1121,7 +898,7 @@ def capillary_loss(wl: np.ndarray, he_mode: tuple[int, int], core_radius: float)
np.ndarray
loss in 1/m
"""
chi_silica = abs(mat.sellmeier(wl, utils.load_material_dico("silica")))
chi_silica = abs(mat.Sellmeier.load("silica").chi(wl))
nu_n = 0.5 * (chi_silica + 2) / np.sqrt(chi_silica)
return nu_n * (u_nm(*he_mode) * wl / pipi) ** 2 * core_radius**-3

318
src/scgenerator/physics/materials.py
View File

@@ -1,10 +1,11 @@
from __future__ import annotations
import functools
from dataclasses import dataclass, field
from typing import Any, TypeVar
from functools import cache
from typing import TypeVar
import numpy as np
from scgenerator import utils
from scgenerator.cache import np_cache
from scgenerator.logger import get_logger
@@ -24,6 +25,7 @@ class Sellmeier:
constant: float = 0
@classmethod
@cache
def load(cls, name: str) -> Sellmeier:
mat_dico = utils.load_material_dico(name)
s = mat_dico.get("sellmeier", {})
@@ -38,7 +40,7 @@ class Sellmeier:
}
)
def chi(self, wl: T) -> T:
def chi(self, wl: T, temperature: float | None = None, pressure: float | None = None) -> T:
"""n^2 - 1"""
if isinstance(wl, np.ndarray):
chi = np.zeros_like(wl) # = n^2 - 1
@@ -59,32 +61,32 @@ class Sellmeier:
else:
raise ValueError(f"kind {self.kind} is not recognized.")
# if temperature is not None:
# chi *= self.temperature_ref / temperature
if temperature is not None:
chi *= self.temperature_ref / temperature
# if pressure is not None:
# chi *= pressure / self.pressure_ref
if pressure is not None:
chi *= pressure / self.pressure_ref
return chi
def n_gas_2(self, wl: T) -> T:
return self.chi(wl) + 1
def n_gas_2(self, wl: T, temperature: float | None = None, pressure: float | None = None) -> T:
return self.chi(wl, temperature, pressure) + 1
def n(self, wl: T) -> T:
return np.sqrt(self.n_gas_2(wl))
def n(self, wl: T, temperature: float | None = None, pressure: float | None = None) -> T:
return np.sqrt(self.n_gas_2(wl, temperature, pressure))
class GasInfo:
class Gas:
name: str
sellmeier: Sellmeier
n2: float
atomic_number: int
atomic_mass: float
ionization_energy: float
_n2: float
ionization_energy: float | None
def __init__(self, gas_name: str):
self.name = gas_name
self.mat_dico = utils.load_material_dico(gas_name)
self.n2 = self.mat_dico["kerr"]["n2"]
self._n2 = self.mat_dico["kerr"]["n2"]
self.atomic_mass = self.mat_dico["atomic_mass"]
self.atomic_number = self.mat_dico["atomic_number"]
self.ionization_energy = self.mat_dico.get("ionization_energy")
@@ -130,11 +132,9 @@ class GasInfo:
/ temperature
)
else:
return number_density_van_der_waals(
self.get("a"), self.get("b"), pressure, temperature
) / number_density_van_der_waals(
self.get("a"),
self.get("b"),
return self.number_density_van_der_waals(
pressure, temperature
) / self.number_density_van_der_waals(
self.sellmeier.pressure_ref,
self.sellmeier.temperature_ref,
)
@@ -161,7 +161,68 @@ class GasInfo:
if ideal_gas:
return pressure / temperature / kB
else:
return number_density_van_der_waals(self.get("a"), self.get("b"), pressure, temperature)
return self.number_density_van_der_waals(pressure, temperature)
def number_density_van_der_waals(self, pressure: float = None, temperature: float = None):
"""returns the number density of a gas
Parameters
----------
pressure : float, optional
pressure in Pa, by default the reference pressure
temperature : float, optional
temperature in K, by default the reference temperature
Returns
----------
the numbers density (/m^3)
Raises
----------
ValueError : Since the Van der Waals equation is a cubic one, there could be more than one real, positive solution
"""
logger = get_logger(__name__)
if pressure == 0:
return 0
a = self.mat_dico.get("a", 0)
b = self.mat_dico.get("b", 0)
pressure = self.mat_dico["sellmeier"].get("P0", 101325) if pressure is None else pressure
temperature = (
self.mat_dico["sellmeier"].get("T0", 273.15) if temperature is None else temperature
)
ap = a / NA**2
bp = b / NA
# setup van der Waals equation for the number density
p3 = -ap * bp
p2 = ap
p1 = -(pressure * bp + kB * temperature)
p0 = pressure
# filter out unwanted matches
roots = np.roots([p3, p2, p1, p0])
roots = roots[np.isreal(roots)].real
roots = roots[roots > 0]
if len(roots) != 1:
s = f"Van der Waals eq with parameters P={pressure}, T={temperature}, a={a}, b={b}"
s += f", There is more than one possible number density : {roots}."
s += f", {np.min(roots)} was returned"
logger.warning(s)
return np.min(roots)
def n2(self, temperature: float | None = None, pressure: float | None = None) -> float:
"""nonlinear refractive index"""
# if pressure and/or temperature are specified, adjustment is made according to number density ratio
if pressure is not None or temperature is not None:
N0 = self.number_density_van_der_waals()
N = self.number_density_van_der_waals(pressure, temperature)
ratio = N / N0
else:
ratio = 1
return ratio * self._n2
@property
def ionic_charge(self):
@@ -178,32 +239,9 @@ class GasInfo:
@np_cache
def n_gas_2(
wl_for_disp: np.ndarray, gas_name: str, pressure: float, temperature: float, ideal_gas: bool
):
def n_gas_2(wl_for_disp: np.ndarray, gas_name: str, pressure: float, temperature: float):
"""Returns the sqare of the index of refraction of the specified gas"""
material_dico = utils.load_material_dico(gas_name)
n_gas_2 = fast_n_gas_2(wl_for_disp, pressure, temperature, ideal_gas, material_dico)
return n_gas_2
def fast_n_gas_2(
wl_for_disp: np.ndarray,
pressure: float,
temperature: float,
ideal_gas: bool,
material_dico: dict[str, Any],
):
if ideal_gas:
n_gas_2 = sellmeier(wl_for_disp, material_dico, pressure, temperature) + 1
else:
N_1 = number_density_van_der_waals(
pressure=pressure, temperature=temperature, material_dico=material_dico
)
N_0 = number_density_van_der_waals(material_dico=material_dico)
n_gas_2 = sellmeier(wl_for_disp, material_dico) * N_1 / N_0 + 1
return n_gas_2
return Sellmeier.load(gas_name).n_gas_2(wl_for_disp, temperature, pressure)
def pressure_from_gradient(ratio, p0, p1):
@@ -220,190 +258,25 @@ def pressure_from_gradient(ratio, p0, p1):
return np.sqrt(p0**2 - ratio * (p0**2 - p1**2))
def number_density_van_der_waals(
a=None, b=None, pressure=None, temperature=None, material_dico=None
):
"""returns the number density of a gas
Parameters
----------
P : pressure
T : temperature
for pressure and temperature, the default
a : Van der Waals a coefficient
b : Van der Waals b coefficient
material_dico : optional. If passed, will compute the number density at given reference values found in material_dico
Returns
----------
the numbers density (/m^3)
Raises
----------
ValueError : Since the Van der Waals equation is a cubic one, there could be more than one real, positive solution
"""
logger = get_logger(__name__)
if pressure == 0:
return 0
if material_dico is not None:
a = material_dico.get("a", 0) if a is None else a
b = material_dico.get("b", 0) if b is None else b
pressure = material_dico["sellmeier"].get("P0", 101325) if pressure is None else pressure
temperature = (
material_dico["sellmeier"].get("T0", 273.15) if temperature is None else temperature
)
else:
a = 0 if a is None else a
b = 0 if b is None else b
pressure = 101325 if pressure is None else pressure
temperature = 273.15 if temperature is None else temperature
ap = a / NA**2
bp = b / NA
# setup van der Waals equation for the number density
p3 = -ap * bp
p2 = ap
p1 = -(pressure * bp + kB * temperature)
p0 = pressure
# filter out unwanted matches
roots = np.roots([p3, p2, p1, p0])
roots = roots[np.isreal(roots)].real
roots = roots[roots > 0]
if len(roots) != 1:
s = f"Van der Waals eq with parameters P={pressure}, T={temperature}, a={a}, b={b}"
s += f", There is more than one possible number density : {roots}."
s += f", {np.min(roots)} was returned"
logger.warning(s)
return np.min(roots)
@functools.singledispatch
def sellmeier(
wl_for_disp: np.ndarray,
material_dico: dict[str, Any],
pressure: float = None,
temperature: float = None,
) -> np.ndarray:
"""reads a file containing the Sellmeier values corresponding to the choses material and
returns the real susceptibility pressure and temperature adjustments are made according to
ideal gas law.
Parameters
----------
wl_for_disp : array, shape (n,)
wavelength vector over which to compute the refractive index
material_dico : dict[str, Any]
material dictionary as explained in scgenerator.utils.load_material_dico
pressure : float, optional
pressure in Pa, by default None
temperature : float, optional
temperature of the gas in Kelvin
Returns
-------
array, shape (n,)
chi = n^2 - 1
"""
logger = get_logger(__name__)
WL_THRESHOLD = 8.285e-6
ind = wl_for_disp < WL_THRESHOLD
temp_l = wl_for_disp[ind]
B = material_dico["sellmeier"]["B"]
C = material_dico["sellmeier"]["C"]
const = material_dico["sellmeier"].get("const", 0)
P0 = material_dico["sellmeier"].get("P0", 1e5)
t0 = material_dico["sellmeier"].get("t0", 273.15)
kind = material_dico["sellmeier"].get("kind", 1)
chi = np.zeros_like(wl_for_disp) # = n^2 - 1
if kind == 1:
logger.debug("materials : using Sellmeier 1st kind equation")
for b, c_ in zip(B, C):
chi[ind] += temp_l**2 * b / (temp_l**2 - c_)
elif kind == 2: # gives n-1
logger.debug("materials : using Sellmeier 2nd kind equation")
for b, c_ in zip(B, C):
chi[ind] += b / (c_ - 1 / temp_l**2)
chi += const
chi = (chi + 1) ** 2 - 1
elif kind == 3: # Schott formula
logger.debug("materials : using Schott equation")
for i, b in reversed(list(enumerate(B))):
chi[ind] += b * temp_l ** (-2 * (i - 1))
chi[ind] = chi[ind] - 1
else:
raise ValueError(f"kind {kind} is not recognized.")
if temperature is not None:
chi *= t0 / temperature
if pressure is not None:
chi *= pressure / P0
logger.debug(f"computed chi between {np.min(chi):.2e} and {np.max(chi):.2e}")
return chi
@sellmeier.register
def sellmeier_scalar(
wavelength: float,
material_dico: dict[str, Any],
pressure: float = None,
temperature: float = None,
) -> float:
"""n^2 - 1"""
return float(sellmeier(np.array([wavelength]), material_dico, pressure, temperature)[0])
def delta_gas(w, material_dico):
def delta_gas(w: np.ndarray, gas: Gas) -> np.ndarray:
"""returns the value delta_t (eq. 24 in Markos(2017))
Parameters
----------
w : angular frequency array
material_dico : material dictionary as explained in scgenerator.utils.load_material_dico
w : np.ndarray
angular frequency array
gas : Gas
Returns
----------
delta_t
since 2 gradients are computed, it is recommended to exclude the 2 extremum values
"""
chi = sellmeier(units.m.inv(w), material_dico)
N0 = number_density_van_der_waals(material_dico=material_dico)
chi = gas.sellmeier.chi(units.m.inv(w))
N0 = gas.number_density_van_der_waals()
dchi_dw = np.gradient(chi, w)
return 1 / (N0 * c) * (dchi_dw + w / 2 * np.gradient(dchi_dw, w))
def non_linear_refractive_index(material_dico, pressure=None, temperature=None):
"""returns the non linear refractive index n2 adjusted for pressure and temperature
NOTE : so far, there is no adjustment made for wavelength
Parameters
----------
lambda_ : wavelength array
material_dico :
pressure : pressure in Pa
temperature : temperature in Kelvin
Returns
----------
n2
"""
n2_ref = material_dico["kerr"]["n2"]
# if pressure and/or temperature are specified, adjustment is made according to number density ratio
if pressure is not None or temperature is not None:
N0 = number_density_van_der_waals(material_dico=material_dico)
N = number_density_van_der_waals(
pressure=pressure, temperature=temperature, material_dico=material_dico
)
ratio = N / N0
else:
ratio = 1
return ratio * n2_ref
def gas_n2(gas_name: str, pressure: float, temperature: float) -> float:
"""returns the nonlinear refractive index
@@ -421,7 +294,7 @@ def gas_n2(gas_name: str, pressure: float, temperature: float) -> float:
float
n2 in m2/W
"""
return non_linear_refractive_index(utils.load_material_dico(gas_name), pressure, temperature)
return Gas(gas_name).n2(temperature, pressure)
def gas_chi3(gas_name: str, wavelength: float, pressure: float, temperature: float) -> float:
@@ -441,11 +314,10 @@ def gas_chi3(gas_name: str, wavelength: float, pressure: float, temperature: flo
float
[description]
"""
mat_dico = utils.load_material_dico(gas_name)
chi = sellmeier_scalar(wavelength, mat_dico, pressure=pressure, temperature=temperature)
return n2_to_chi3(
non_linear_refractive_index(mat_dico, pressure, temperature), np.sqrt(chi + 1)
)
gas = Gas(gas_name)
n = gas.sellmeier.n(wavelength, temperature, pressure)
n2 = gas.n2(temperature, pressure)
return n2_to_chi3(n2, n)
def n2_to_chi3(n2: float, n0: float) -> float:

12
src/scgenerator/physics/plasma.py
View File

@@ -1,12 +1,12 @@
from typing import Any, Callable, NamedTuple, TypeVar
from typing import Any, Callable, NamedTuple
import numpy as np
import scipy.special
from scipy.interpolate import interp1d
from numpy.core.numeric import zeros_like
from scipy.interpolate import interp1d
from ..math import cumulative_simpson, expm1_int
from .units import e, hbar, me
from scgenerator.math import cumulative_simpson, expm1_int
from scgenerator.physics.units import e, hbar, me
class PlasmaInfo(NamedTuple):
@@ -99,7 +99,7 @@ class Plasma:
loss_term = np.zeros_like(field)
loss_term[nzi] = dn_dt[nzi] * self.ionization_energy / field[nzi]
phase_term = self.dt * e ** 2 / me * cumulative_simpson(electron_density * field)
phase_term = self.dt * e**2 / me * cumulative_simpson(electron_density * field)
dp_dt = loss_term + phase_term
polarization = self.dt * cumulative_simpson(dp_dt)
@@ -117,4 +117,4 @@ def free_electron_density(rate: np.ndarray, dt: float, N0: float) -> np.ndarray:
def barrier_suppression(ionpot, Z):
Ip_au = ionpot / 4.359744650021498e-18
ns = Z / np.sqrt(2 * Ip_au)
return Z ** 3 / (16 * ns ** 4) * 5.14220670712125e11
return Z**3 / (16 * ns**4) * 5.14220670712125e11

9
src/scgenerator/physics/pulse.py
View File

@@ -15,18 +15,17 @@ from dataclasses import astuple, dataclass
from pathlib import Path
from typing import Literal, Tuple, TypeVar
import numba
import matplotlib.pyplot as plt
import numpy as np
from numpy import pi
from numpy.fft import fft, fftshift, ifft
from scipy.interpolate import UnivariateSpline, interp1d
from scipy.optimize import minimize_scalar
from scipy.optimize.optimize import OptimizeResult
from scipy.optimize._optimize import OptimizeResult
from ..defaults import default_plotting
from .. import math
from . import units
from scgenerator import math
from scgenerator.defaults import default_plotting
from scgenerator.physics import units
c = 299792458.0
hbar = 1.05457148e-34

12
src/scgenerator/physics/simulate.py
View File

@@ -1,20 +1,20 @@
import multiprocessing
import multiprocessing.connection
import os
import warnings
from collections import defaultdict
from datetime import datetime
from logging import Logger
from pathlib import Path
from typing import Any, Generator, Iterator, Optional, Type, Union
import warnings
import numpy as np
from .. import solver, utils
from ..logger import get_logger
from ..operators import CurrentState
from ..parameter import FileConfiguration, Parameters
from ..pbar import PBars, ProgressBarActor, progress_worker
from scgenerator import solver, utils
from scgenerator.logger import get_logger
from scgenerator.operators import CurrentState
from scgenerator.parameter import FileConfiguration, Parameters
from scgenerator.pbar import PBars, ProgressBarActor, progress_worker
try:
import ray

23
src/scgenerator/plotting.py
View File

@@ -6,18 +6,18 @@ import matplotlib.gridspec as gs
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes
from matplotlib.figure import Figure
from matplotlib.colors import ListedColormap
from matplotlib.figure import Figure
from matplotlib.transforms import offset_copy
from scipy.interpolate import UnivariateSpline, interp1d
from scipy.interpolate import UnivariateSpline
from . import math
from .const import PARAM_SEPARATOR
from .defaults import default_plotting as defaults
from .math import abs2, linear_interp_2d, span
from .parameter import Parameters
from .physics import pulse, units
from .physics.units import PlotRange, sort_axis
from scgenerator import math
from scgenerator.const import PARAM_SEPARATOR
from scgenerator.defaults import default_plotting as defaults
from scgenerator.math import abs2, linear_interp_2d, span
from scgenerator.parameter import Parameters
from scgenerator.physics import pulse, units
from scgenerator.physics.units import PlotRange, sort_axis
RangeType = tuple[float, float, Union[str, Callable]]
NO_LIM = object()
@@ -563,7 +563,6 @@ def transform_mean_values(
np.ndarray, shape (m, n)
all the values
"""
AAA
if values.ndim != 2:
print(f"Shape was {values.shape}. Can only plot 2D arrays")
return
@@ -1133,7 +1132,5 @@ def annotate_fwhm(
else:
offset = -arrow_length_pts
x, y = (left, v_max / 2)
trans = offset_copy(
ax.transData, ax.get_figure(), offset, 0, "points"
)
trans = offset_copy(ax.transData, ax.get_figure(), offset, 0, "points")
ax.text(x, y, arrow_label, transform=trans, **text_kwargs)

13
src/scgenerator/solver.py
View File

@@ -1,23 +1,22 @@
from __future__ import annotations
from collections import defaultdict
import logging
from abc import abstractmethod
from collections import defaultdict
from typing import Iterator, Type
import numba
import numpy as np
from .logger import get_logger
from .operators import (
from scgenerator.logger import get_logger
from scgenerator.operators import (
AbstractConservedQuantity,
CurrentState,
LinearOperator,
NonLinearOperator,
ValueTracker,
)
from .utils import get_arg_names
import warnings
from scgenerator.utils import get_arg_names
class IntegratorError(Exception):
@@ -455,8 +454,8 @@ def rk4ip_step(
@numba.jit(nopython=True)
def compute_diff(coarse_spec: np.ndarray, fine_spec: np.ndarray) -> float:
diff = coarse_spec - fine_spec
diff2 = diff.imag ** 2 + diff.real ** 2
return np.sqrt(diff2.sum() / (fine_spec.real ** 2 + fine_spec.imag ** 2).sum())
diff2 = diff.imag**2 + diff.real**2
return np.sqrt(diff2.sum() / (fine_spec.real**2 + fine_spec.imag**2).sum())
def get_integrator(integration_scheme: str):

18
src/scgenerator/spectra.py
View File

@@ -7,21 +7,21 @@ from typing import Callable, Iterator, Optional, Union
import matplotlib.pyplot as plt
import numpy as np
from . import math
from .const import PARAM_FN, SPEC1_FN, SPEC1_FN_N
from .logger import get_logger
from .parameter import Parameters
from .physics import pulse, units
from .physics.units import PlotRange
from .plotting import (
from scgenerator import math
from scgenerator.const import PARAM_FN, SPEC1_FN, SPEC1_FN_N
from scgenerator.legacy import translate_parameters
from scgenerator.logger import get_logger
from scgenerator.parameter import Parameters
from scgenerator.physics import pulse, units
from scgenerator.physics.units import PlotRange
from scgenerator.plotting import (
mean_values_plot,
propagation_plot,
single_position_plot,
transform_1D_values,
transform_2D_propagation,
)
from .utils import load_spectrum, simulations_list, load_toml
from .legacy import translate_parameters
from scgenerator.utils import load_spectrum, load_toml, simulations_list
class Spectrum(np.ndarray):

4
src/scgenerator/variationer.py
View File

@@ -1,6 +1,6 @@
from math import prod
import itertools
from collections.abc import MutableMapping, Sequence
from math import prod
from pathlib import Path
from typing import Any, Callable, Generator, Generic, Iterable, Iterator, Optional, TypeVar, Union
@@ -8,7 +8,7 @@ import numpy as np
from pydantic import validator
from pydantic.main import BaseModel
from .const import PARAM_SEPARATOR
from scgenerator.const import PARAM_SEPARATOR
T = TypeVar("T")

2
testing/test_resonance.py
View File

@@ -6,7 +6,7 @@ import matplotlib.pyplot as plt
def main():
capillary_thickness = 1.4e-6
wl = np.linspace(200e-9, 2000e-9, 500)
n_gas_2 = sc.materials.n_gas_2(wl, "air", 3e5, 300, True)
n_gas_2 = sc.materials.n_gas_2(wl, "air", 3e5, 300)
resonances = []
for i in range(5):
t = sc.fiber.resonance_thickness(wl, i, n_gas_2, 40e-6)