Source code for pynlo.medium

# -*- coding: utf-8 -*-
"""
Optical modes in waveguides or other media.

"""

__all__ = ["Mode"]


# %% Imports

import bisect
import collections
import copy

import numpy as np
from scipy.constants import c, pi


# %% Collections

_LinearOperator = collections.namedtuple("LinearOperator", ["u", "gain", "phase", "phase_raw"])
_LinearZ = collections.namedtuple("LinearZ", ["any", "alpha", "beta"])
_NonlinearZ = collections.namedtuple("NonlinearZ", ["any", "g2", "pol", "g3", "r3"])


# %% Single Mode

[docs] class Mode(): """ An optical mode. An optical mode is defined in the frequency domain. All properties must be input as effective values, i.e. those found by integrating out the transverse spatial dependence of the medium and mode. If a parameter is z dependent, it can be input as a function who's first argument is the propagation distance. Parameters ---------- v_grid : array_like of float The frequency grid. beta : array_like of float or callable The phase coefficient, the real part of the complex wavenumber. alpha : array_like of float or callable, optional The gain coefficient, twice the imaginary part of the complex wavenumber. g2 : array_like of complex or callable, optional The effective 2nd-order nonlinearity. g2_inv : array_like of float, optional The location of all poled domain inversion boundaries. g3 : array_like of complex or callable, optional The effective 3rd-order nonlinearity. rv_grid : array_like of float, optional An origin-contiguous frequency grid associated with the 3rd-order nonlinear response function. r3 : array_like of complex or callable, optional The effective 3rd-order nonlinear response function containing both the Raman and instantaneous nonlinearities. z : float, optional The initial position within the mode. The default is 0. Notes ----- Forward traveling waves of a mode are defined using the following conventions: .. math:: E, H \\sim a \\, e^{i(\\omega t - \\kappa z)} + \\text{c.c.} \\\\ \\kappa = \\beta + i \\frac{\\alpha}{2}, \\quad \\beta = n \\frac{\\omega}{c} """ def __init__(self, v_grid, beta, alpha=None, g2=None, g2_inv=None, g3=None, rv_grid=None, r3=None, z=0.0): #---- Position self._z = z #---- Frequency Grid self._v_grid = np.asarray(v_grid, dtype=float) self._w_grid = 2*pi*self._v_grid #---- Refractive Index if callable(beta): assert (len(beta(z)) == len(v_grid)), "The length of beta must match v_grid." self._beta = beta else: assert (len(beta) == len(v_grid)), "The length of beta must match v_grid." self._beta = np.asarray(beta, dtype=float) #---- Gain if (alpha is None) or callable(alpha): self._alpha = alpha else: self._alpha = np.asarray(alpha, dtype=float) #---- 2nd-Order Nonlinearity if (g2 is None) or callable(g2): self._g2 = g2 else: self._g2 = np.asarray(g2, dtype=complex) if g2_inv is None: self._g2_inv = None self._g2_inv_sorted = [] else: assert (g2 is not None) and (g2_inv is not None), ( "Poling can only be defined when g2 is defined") self._g2_inv_sorted = sorted(g2_inv) self._g2_inv = {z:(idx + 1) % 2 for idx, z in enumerate(self._g2_inv_sorted)} #---- 3rd-Order Nonlinearity if (g3 is None) or callable(g3): self._g3 = g3 else: self._g3 = np.asarray(g3, dtype=complex) # Nonlinear Response Function if (rv_grid is not None) and (r3 is not None): assert (g3 is not None) and (r3 is not None), ( "Raman nonlinearity can only be defined when g3 is defined") self._rv_grid = np.asarray(rv_grid, dtype=float) if callable(r3): self._r3 = r3 else: assert (len(r3) == len(rv_grid)), "The length of r3 must match rv_grid." self._r3 = np.asarray(r3, dtype=complex) else: assert (rv_grid is None) and (r3 is None), ( "rv_grid and r3 must both be defined at the same time or not at all.") self._rv_grid = None self._r3 = None #---- Z Dependence self._z_linear = _LinearZ( any=callable(alpha) or callable(beta), alpha=callable(alpha), beta=callable(beta)) self._z_nonlinear = _NonlinearZ( any=callable(g2) or callable(g3) or callable(r3), g2=callable(g2), pol=g2_inv is not None, g3=callable(g3), r3=callable(r3)) self._z_mode = self.z_linear.any or self.z_nonlinear.any or self.z_nonlinear.pol #---- General Properties @property def z(self): """ The position within the mode, with units of ``m``. Returns ------- float """ return self._z @z.setter def z(self, z): self._z = z @property def v_grid(self): """ The frequency grid, with units of ``Hz``. Returns ------- ndarray of float """ return self._v_grid @property def rv_grid(self): """ The origin-contiguous frequency grid associated with the Raman response. Units are in ``Hz``. Returns ------- None or ndarray of float """ return self._rv_grid @property def z_mode(self): """ The z dependence of the mode. Returns ------- bool """ return self._z_mode @property def z_linear(self): """ The z dependence of the linear terms. Returns ------- any : bool Whether there is any z dependence of the linearity. alpha : bool Z-dependent gain coefficient. beta : bool Z-dependent phase coefficient. """ return self._z_linear @property def z_nonlinear(self): """ The z dependence of the nonlinear terms. Returns ------- any : bool Whether there is any z dependence of the nonlinearity (excluding poling). g2 : bool Z-dependent 2nd-order nonlinear parameter. pol : bool Poled 2nd-order nonlinearity. g3 : bool Z-dependent 3rd-order nonlinear parameter. r3 : bool Z-dependent Raman response. """ return self._z_nonlinear #---- 1st-Order Properties @property def alpha(self): """ The gain coefficient, with units of ``1/m``. Positive values correspond to gain and negative values to loss. Returns ------- None or ndarray of float """ return self._alpha(self.z) if callable(self._alpha) else self._alpha @property def beta(self): """ The phase coefficient, or angular wavenumber, with units of ``1/m``. Returns ------- ndarray of float """ return self._beta(self.z) if callable(self._beta) else self._beta @property def n(self): """ The refractive index. Returns ------- ndarray of float """ return self.beta*c/self._w_grid @property def beta1(self): """ The group walk-off parameter, with units of ``s/m``. Returns ------- ndarray of float """ return np.gradient(self.beta, self._w_grid, edge_order=2)
[docs] def d_12(self, v0=None): """ The group velocity mismatch, with units of ``s/m``. Parameters ---------- v0 : float, optional The target reference frequency. The central frequency is selected by default. Returns ------- ndarray of float """ if v0 is None: v0 = self.v_grid[self.v_grid.size//2] v0_idx = np.argmin(np.abs(v0 - self.v_grid)) beta1 = self.beta1 return beta1[v0_idx] - beta1
@property def n_g(self): """ The group index. Returns ------- ndarray of float """ return c*self.beta1 @property def v_g(self): """ The group velocity, with units of ``m/s``. Returns ------- ndarray of float """ return 1/self.beta1 @property def beta2(self): """ The group velocity dispersion (GVD), with units of ``s**2/m``. Returns ------- ndarray of float """ return np.gradient(self.beta1, self._w_grid, edge_order=2) @property def D(self): """ The dispersion parameter D, with units of ``s/m**2``. Returns ------- ndarray of float """ return -2*pi/c * self.v_grid**2 * self.beta2 #TODO: test against chi1 helper functions
[docs] def linear_operator(self, dz, v0=None): """ The linear operator which advances a pulse over a distance `dz`. The linear operator acts on the analytic spectrum through multiplication. Parameters ---------- dz : float The step size. v0 : float, optional The target reference frequency of the comoving frame. The default selects the central frequency. Returns ------- u : ndarray of complex The forward evolution operator. gain : float or ndarray The accumulated gain or loss (multiplicative). phase : ndarray of float The accumulated phase in the comoving frame (additive). phase_raw : ndarray of float The raw accumulated phase. """ #---- Gain alpha = self.alpha if alpha is None: alpha = 0.0 gain = np.exp(alpha*dz) #---- Phase beta_raw = self.beta # Comoving frame if v0 is None: v0 = self.v_grid[self.v_grid.size//2] v0_idx = np.argmin(np.abs(v0 - self.v_grid)) beta_cm = beta_raw - self.beta1[v0_idx]*self._w_grid #---- Propagation Constant kappa = beta_cm + 0.5j*alpha #---- Linear Operator operator = np.exp(-1j*kappa*dz) lin_operator = _LinearOperator( u=operator, gain=gain, phase=dz*beta_cm, phase_raw=dz*beta_raw) return lin_operator
#---- 2nd-Order Properties @property def g2(self): """ The magnitude of the effective 2nd-order nonlinear parameter, with units of ``1/(W**0.5*m*Hz)``. Returns ------- None or ndarray of complex """ return self._g2(self.z) if callable(self._g2) else self._g2 @property def g2_inv(self): """ The location of all 2nd-order domain inversion boundaries within the mode. A value of 1 indicates the start of an inverted domain. A value of 0 indicates the start of an unpoled region. Returns ------- None or dict of int """ return self._g2_inv @property def g2_pol(self): """ The poling status at the current z position. A value of 1 indicates that the current position is in a region with an inverted domain. A value of 0 indicates an unpoled region. Returns ------- int """ poled = bisect.bisect_right(self._g2_inv_sorted, self.z) % 2 return poled #---- 3rd-Order Properties @property def g3(self): """ The effective 3rd-order nonlinear parameter, with units of ``1/(W*m*Hz)``. Returns ------- None or ndarray of complex """ return self._g3(self.z) if callable(self._g3) else self._g3 @property def gamma(self): """ The nonlinear parameter :math:`\\gamma`, with units of ``1/(W*m)``. Returns ------- None or ndarray of complex """ g3 = self.g3 if g3 is not None and len(g3.shape)>=2: g3 = g3[0] * np.sum(g3[1:]**3, axis=0) return 3/2*self._w_grid*g3 if g3 is not None else None #TODO: test against chi3 helper functions @property def r3(self): """ The effective 3rd-order nonlinear response function containing both the Raman and instantaneous nonlinearities. Returns ------- None or ndarray of complex """ return self._r3(self.z) if callable(self._r3) else self._r3 #---- Misc
[docs] def copy(self): """Copy the mode.""" return copy.deepcopy(self)
# class GaussianMode(): # """ # Fundamental Gaussian modes for simulating single-mode free-space # propagation. # - effective area based on distance to nominal waist location # - could also include convenience functions for setting up the beam # through focusing, propagation, etc. # - check Boyd 2.10 "Nonlinear Optical Interactions with Focused Gaussian # Beams" for complicating factors. # # """ # %% Multimode # class Waveguide(): # """ # Collection of modes and the nonlinear interactions between them # """ # def __init__(self, modes, coupling): # pass # class FreeSpace(Waveguide): # """ # Collection of Hermite–Gaussian or Laguerre–Gaussian modes for simulating # free space propagation of arbitrary distribution. # - effective area based on distance to nominal waist location # - could also include convenience functions for setting up the beam # through focusing, propagation, etc. # # """