threshold.py

Go to the documentation of this file.

# This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
# Copyright (c) 2022 Lodz University of Technology
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 3.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
 
# coding: utf8
# Copyright (C) 2014 Photonics Group, Lodz University of Technology
 
import numpy as np
 
import electrical.diffusion
import electrical.shockley
import gain.freecarrier
import plask
import thermal.static
 
try:
    import scipy
except ImportError:
    plask.print_log('warning', "scipy could not be imported."
                               " You will not be able to run some algorithms."
                               " Install scipy to resolve this issue.")
else:
    import scipy.optimize
 
from .thermoelectric import attribute, h5open, ThermoElectric
from . import _doc
 
 
class ThresholdSearch(ThermoElectric):
 
    _OPTICAL_ROOTS = {}
 
    _Diffusion = None
    _Gain = None
    _Optical = None
    _optarg = 'lam'
    _lam0 = 'lam0'
 
    tfreq = 6.0
    vmin = None
    vmax = None
    optical_resolution = (800, 600)
    vtol = 1e-5
    quick = False
    maxiter = 50
    skip_thermal = False
 
    def __init__(self, name):
        super().__init__(name)
        self.diffusion = self._Diffusion(name)
        self.gain = self._Gain(name)
        self.optical = self._Optical(name)
        self.__reconnect__reconnect()
        self.threshold_voltage = None
        self.threshold_current = None
        self._invalidate = None
        self.modeno = None
        self._sn = 0
 
    def __reconnect(self):
        self.diffusion.inTemperature = self.thermal.outTemperature
        self.diffusion.inCurrentDensity = self.electrical.outCurrentDensity
        self.gain.inTemperature = self.thermal.outTemperature
        self.gain.inCarriersConcentration = self.diffusion.outCarriersConcentration
        self.optical.inTemperature = self.thermal.outTemperature
        self.optical.inGain = self.gain.outGain
 
    def reconnect(self):
        """
        Reconnect all internal solvers.
 
        This method should be called if some of the internal solvers were changed manually.
        """
        super().reconnect()
        self.__reconnect__reconnect()
 
    def _parse_xpl(self, tag, manager):
        if tag == 'root':
            self.vmin = tag.get('vmin', self.vmin)
            self.vmax = tag.get('vmax', self.vmax)
            self.ivb = tag['bcond']
            self.vtolvtol = tag.get('vtol', self.vtolvtol)
            self.maxitermaxiter = tag.get('maxiter', self.maxitermaxiter)
            self.quickquick = tag.get('quick', self.quickquick)
        elif tag == 'diffusion':
            self._read_attr(tag, 'fem-method', self.diffusion, str, 'fem_method')
            self._read_attr(tag, 'accuracy', self.diffusion, float)
            self._read_attr(tag, 'abs-accuracy', self.diffusion, float, 'abs_accuracy')
            self._read_attr(tag, 'maxiters', self.diffusion, int)
            self._read_attr(tag, 'maxrefines', self.diffusion, int)
            self._read_attr(tag, 'interpolation', self.diffusion, str)
        elif tag == 'gain':
            self._read_attr(tag, 'lifetime', self.gain, float)
            self._read_attr(tag, 'matrix-elem', self.gain, float, 'matrix_element')
            self._read_attr(tag, 'strained', self.gain, bool)
        elif tag.name in self._OPTICAL_ROOTS:
            root = getattr(self.optical, self._OPTICAL_ROOTS[tag.name])
            self._read_attr(tag, 'method', root, str)
            self._read_attr(tag, 'tolx', root, float)
            self._read_attr(tag, 'tolf-min', root, float, 'tolf_min')
            self._read_attr(tag, 'tolf-max', root, float, 'tolf_max')
            self._read_attr(tag, 'maxstep', root, float)
            self._read_attr(tag, 'maxiter', root, int)
            self._read_attr(tag, 'alpha', root, float)
            self._read_attr(tag, 'lambda', root, float)
            self._read_attr(tag, 'initial-range', root, tuple, 'initial_range')
        elif tag == 'output':
            self.optical_resolutionoptical_resolution = tag.get('optical-res-x', self.optical_resolutionoptical_resolution[0]), \
                                      tag.get('optical-res-y', self.optical_resolutionoptical_resolution[1])
        else:
            if tag == 'geometry':
                self.optical.geometry = self.diffusion.geometry = self.gain.geometry = \
                    tag.getitem(manager.geo, 'optical')
            elif tag == 'mesh':
                if 'diffusion' in tag: self.diffusion.mesh = tag.getitem(manager.msh, 'diffusion')
                if 'optical' in tag: self.optical.mesh = tag.getitem(manager.msh, 'optical')
                if 'gain' in tag: self.gain.mesh = tag.getitem(manager.msh, 'gain')
            elif tag == 'loop':
                self.skip_thermalskip_thermal = tag.get('skip-thermal', self.skip_thermalskip_thermal)
                self._read_attr(tag, 'inittemp', self.gain, float, 'T0')
            super()._parse_xpl(tag, manager)
 
    def on_initialize(self):
        if not self.skip_thermalskip_thermal:
            self.thermal.initialize()
        self.electrical.initialize()
        self.diffusion.initialize()
        self.gain.initialize()
        self.optical.initialize()
 
    def on_invalidate(self):
        if not self.skip_thermalskip_thermal:
            self.thermal.invalidate()
        self.electrical.invalidate()
        self.diffusion.invalidate()
        self.optical.invalidate()
 
    def _optargs(self):
        return {}
 
    def _get_lam(self):
        if callable(self.get_lamget_lam):
            return self.get_lamget_lam()
        else:
            return self.get_lamget_lam
 
    def get_lam(self):
        raise NotImplemented('get_lam')
 
    def compute_thermoelectric(self):
        """
        Perform thermo-electric calculations.
 
        This method may be called manually to perform thermo-electric calculations.
        Afterwards, one may investigate gain spectrum or verify settings of the optical
        solver.
        """
        self.initialize()
        if self.skip_thermalskip_thermal:
            self.electrical.compute()
        else:
            if self._invalidate:
                self.thermal.invalidate()
                self.electrical.invalidate()
                self.diffusion.invalidate()
                self.gain.invalidate()
            verr = 2. * self.electrical.maxerr
            terr = 2. * self.thermal.maxerr
            while terr > self.thermal.maxerr or verr > self.electrical.maxerr:
                verr = self.electrical.compute(self.tfreqtfreq)
                terr = self.thermal.compute(1)
        self.diffusion.compute()
 
    def step(self, volt, save=False):
        """
        Function performing one step of the threshold search.
 
        Args:
            volt (float): Voltage on a specified boundary condition (V).
 
            save (bool): If `True` the computed fields are saved to the
                HDF5 file after each computations step.
 
        Returns:
            float: Loss of a specified mode
        """
        plask.print_log('detail', "ThresholdSearch: V = {0:.4f} V".format(volt))
        self.electrical.voltage_boundary[self.ivb].value = volt
        self.compute_thermoelectric()
        self.optical.invalidate()
        optstart = getattr(self, 'get_'+self._optarg)()
        if self._optarg is None:
            self.modeno = self.optical.find_mode(optstart, **self._optargs())
        else:
            _optargs = self._optargs().copy()
            _optargs[self._optarg] = optstart
            self.modeno = self.optical.find_mode(**_optargs)
        val = self.optical.modes[self.modeno].loss
        plask.print_log('result', "ThresholdSearch: V = {:.4f} V, loss = {:g} / cm".format(volt, val))
        if save:
            self.savesave(None if save is True else save, 'ThresholdSearchStep/{:03}[{:.4f}V]'.format(self._sn, volt))
            self._sn += 1
        return val
 
    def _quickstep(self, arg):
        """
        Function performing one step of the quick threshold search.
 
        Args:
            arg (array): Array containing voltage on a specified boundary condition (V) and wavelength.
 
        Returns:
            array: Imaginary and real part of a specified mode
        """
        volt, lam = arg * self._quickscale
        plask.print_log('detail', "ThresholdSearch: V = {0:.4f} V, lam = {1:g} nm".format(volt, lam))
        self.electrical.voltage_boundary[self.ivb].value = volt
        self.compute_thermoelectric()
        if self._optarg is None:
            det = self.optical.get_determinant(lam, **self._optargs())
        else:
            _optargs = self._optargs().copy()
            _optargs[self._optarg] = lam
            det = self.optical.get_determinant(**_optargs)
        plask.print_log('result', "ThresholdSearch: V = {0:.4f} V, lam = {1:g} nm, det = {2.real}{2.imag:+g}j"
                        .format(volt, lam, det))
        return np.array([det.real, det.imag])
 
    def get_optical_determinant(self, lam):
        """
        Function computing determinant of the optical solver.
 
        Args:
             lam (float or array): Wavelength to compute the determinant for (nm).
 
        Returns:
            float or array: Optical determinant.
        """
        self.compute_thermoelectric()
        if self._optarg is None:
            return self.optical.get_determinant(lam, **self._optargs())
        else:
            _optargs = self._optargs().copy()
            _optargs[self._optarg] = lam
            return self.optical.get_determinant(**_optargs)
 
    def _get_in_junction(self, func, axis=None):
        if axis is None: axis = self.optical.mesh
        results = {}
        for i, (lb, msh) in enumerate(self._iter_levels(self.diffusion.geometry, axis, 'QW', 'QD', 'gain')):
            if not lb:
                lb = 0
            else:
                try: lb = int(lb)
                except ValueError: pass
            results[lb] = func(msh)
        return results
 
    def get_junction_concentrations(self, interpolation='linear'):
        """
        Get carriers concentration at the active regions.
 
        Args:
            interpolation (str): Interpolation used when retrieving current density.
 
        Return:
            dict: Dictionary of junction current density data.
                  Keys are the junction number.
        """
        return self._get_in_junction(lambda msh: self.diffusion.outCarriersConcentration(msh, interpolation),
                                     self.diffusion.mesh)
 
    def get_junction_gains(self, axis=None, interpolation='linear'):
        """
        Get gain at the active regions.
 
        Args:
            axis (mesh or sequence): Points along horizontal axis to plot gain at.
                                     Defaults to the optical mesh.
 
           interpolation (str): Interpolation used when retrieving current density.
 
        Return:
            dict: Dictionary of junction current density data.
                  Keys are the junction number.
        """
        lam = getattr(self.optical, self._lam0).real
        if lam is None: lam = self._get_lam().real
        return self._get_in_junction(lambda msh: self.gain.outGain(msh, lam, interpolation), axis)
 
    def get_gain_spectrum(self, lams, pos=0., junction=0):
        """
        Get gain spectrum for specified junction.
 
        Args:
            lams (array of floats): Wavelengths for which the spectrum should be plotted.
 
            pos (float): Lateral position fo the point in which the spectrum is plotted.
 
            junction (int): Junction number to take gain from.
 
        Return:
            Data: Gain spectrum.
        """
        level = list(self._iter_levels(self.diffusion.geometry, [pos]))[junction][1]
        return self.gain.spectrum(level[0])(lams)
 
    def _plot_in_junction(self, func, axis, bounds, kwargs, label):
        if axis is None: axis = self.optical.mesh
        i = 0
        for i, (lb, msh) in enumerate(self._iter_levels(self.diffusion.geometry, axis, 'QW', 'QD', 'gain')):
            values = np.array(func(msh))
            if label is None:
                lab = "Junction {:s}".format(lb)
            elif isinstance(label, tuple) or isinstance(label, tuple):
                lab = label[i]
            else:
                lab = label
            plask.plot(msh.axis0, values, label=lab, **kwargs)
        if i > 1:
            plask.legend(loc='best')
        plask.xlabel(u"${}$ (µm)".format(plask.config.axes[-2]))
        if bounds:
            self._plot_hbounds(self.optical)
 
    def plot_junction_concentration(self, bounds=True, interpolation='linear', label=None, **kwargs):
        """
        Plot carriers concentration at the active region.
 
        Args:
            bounds (bool): If *True* then the geometry objects boundaries are
                           plotted.
 
            interpolation (str): Interpolation used when retrieving current density.
 
            label (str or sequence): Label for each junction. It can be a sequence of
                                     consecutive labels for each junction, or a string
                                     in which case the same label is used for each
                                     junction. If omitted automatic label is generated.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        self._plot_in_junction(lambda msh: self.diffusion.outCarriersConcentration(msh, interpolation),
                               self.diffusion.mesh, bounds, kwargs, label)
        plask.ylabel(u"Carriers Concentration (1/cm\xb3)")
        plask.window_title("Carriers Concentration")
 
    def plot_junction_gain(self, axis=None, bounds=True, interpolation='linear', label=None, **kwargs):
        """
        Plot gain at the active region.
 
        Args:
            axis (mesh or sequence): Points along horizontal axis to plot gain at.
                                     Defaults to the optical mesh.
 
            bounds (bool): If *True* then the geometry objects boundaries are
                           plotted.
 
            interpolation (str): Interpolation used when retrieving current density.
 
            label (str or sequence): Label for each junction. It can be a sequence of
                                     consecutive labels for each junction, or a string
                                     in which case the same label is used for each
                                     junction. If omitted automatic label is generated.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        lam = getattr(self.optical, self._lam0).real
        if lam is None: lam = self._get_lam().real
        self._plot_in_junction(lambda msh: self.gain.outGain(msh, lam, interpolation),
                               axis, bounds, kwargs, label)
        plask.ylabel(u"Gain (1/cm)")
        plask.window_title("Gain Profile")
 
    def plot_gain_spectrum(self, lams, pos=0., junction=0, comp=None, **kwargs):
        """
        Plot gain spectrum for specified junction.
 
        Args:
            lams (array of floats): Wavelengths for which the spectrum should be plotted.
 
            pos (float): Lateral position fo the point in which the spectrum is plotted.
 
            junction (int): Junction number to take gain from.
 
            comp (int or str): Spectrum component to plot
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        level = list(self._iter_levels(self.diffusion.geometry, [pos], 'QW', 'QD', 'gain'))[junction][1]
        spectrum = self.gain.spectrum(level[0])
        if comp is None:
            plask.plot(lams, spectrum(lams), **kwargs)
        else:
            _comp = comp
            try:
                if isinstance(comp, str):
                    comp = plask.config.axes.index(comp)
                    if comp != 0: comp -= 1
                if comp < 0 or comp > 1:
                    raise ValueError(comp)
            except (ValueError, TypeError) as err:
                    raise ValueError("Bad spectrum component '{}'".format(_comp)) from err
            plask.plot(lams, np.array(spectrum(lams))[:,comp], **kwargs)
        plask.xlabel(u"Wavelength (nm)")
        plask.ylabel(u"Gain (1/cm)")
        plask.window_title("Gain Spectrum")
 
    def plot_optical_determinant(self, lams, **kwargs):
        """
        Function plotting determinant of the optical solver.
 
        Args:
            lams (array): Wavelengths to plot the determinant for (nm).
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        vals = self.get_optical_determinant(lams)
        plask.plot(lams, abs(vals))
        plask.yscale('log')
        plask.xlabel("Wavelength (nm)")
        plask.ylabel("Determinant [ar.u.]")
 
    def compute(self, save=True, invalidate=False, group='ThresholdSearch', stepsave=False):
        """
        Execute the algorithm.
 
        In the beginning the solvers are invalidated and next, the self-
        consistent loop of thermal, electrical, gain, and optical calculations
        are run within the root-finding algorithm until the mode is found
        with zero optical losses.
 
        Args:
            save (bool or str): If `True` the computed fields are saved to the
                HDF5 file named after the script name with the suffix denoting
                either the batch job id or the current time if no batch system
                is used. The filename can be overridden by setting this parameter
                as a string.
 
            invalidate (bool): If this flag is set, solvers are invalidated
                               in the beginning of the computations.
 
            group (str): HDF5 group to save the data under.
 
            stepsave (bool): If `True` the computed fields are saved to the
                HDF5 file after each computations step.
 
        Returns:
            The voltage set to ``ivolt`` boundary condition for the threshold.
            The threshold current can be then obtained by calling:
 
            >>> solver.get_total_current()
            123.0
        """
 
        if invalidate:
            if not self.skip_thermalskip_thermal:
                self.thermal.invalidate()
            self.electrical.invalidate()
            self.diffusion.invalidate()
            self.optical.invalidate()
        self._invalidate = invalidate
        self.initialize()
 
        self._sn = 0
        self._max_concentration = None
        self._max_gain = None
 
        if (self.vmin is None) != (self.vmax is None):
            raise ValueError("Both 'vmin' and 'vmax' must be either None or a float")
        if self.vmin is None:
            volt = self.electrical.voltage_boundary[self.ivb].value
            if self.quickquick:
                if stepsave:
                    plask.print_log('warning', "Fields cannot be saved in each step if the quick method is used")
                self.compute_thermoelectric()
                lam = self._get_lam().real
                self._quickscale = np.array([volt, lam])
                result = scipy.optimize.root(self._quickstep, np.array([1., 1.]),
                                             tol=1e-3 * min(self.vtolvtol/volt, 1e-6),
                                             options=dict(maxfev=2*self.maxitermaxiter, eps=1e-6))
                if not result.status:
                    raise plask.ComputationError(result.message)
                self.threshold_voltage, lam = result.x * self._quickscale
                if self._optarg is None:
                    self.modeno = self.optical.set_mode(lam, **self._optargs())
                else:
                    _optargs = self._optargs().copy()
                    _optargs[self._optarg] = lam
                    self.modeno = self.optical.set_mode(**_optargs)
            else:
                self.threshold_voltage = scipy.optimize.newton(self.step, volt, args=(stepsave,),
                                                               tol=self.vtolvtol/volt, maxiter=self.maxitermaxiter)
        else:
            if self.quickquick:
                raise RuntimeError("Quick computation only allowed with a single starting point")
            self.threshold_voltage = scipy.optimize.brentq(self.step, self.vmin, self.vmax, args=(stepsave,),
                                                           xtol=self.vtolvtol * 2. / (self.vmin+self.vmax),
                                                           maxiter=self.maxitermaxiter, disp=True)
 
        self.threshold_current = abs(self.electrical.get_total_current())
 
        infolines = self._get_info_get_info()
 
        if save:
            filename = self.savesave(None if save is True else save)
            if filename.endswith('.h5'): filename = filename[:-3]
            plask.print_log('info', "Results saved to file '{}.txt'".format(filename))
            with open(filename+'.txt', 'a', encoding='utf-8') as out:
                out.writelines(line + '\n' for line in infolines)
                out.write("\n")
 
        plask.print_log('important', "Threshold Search Finished")
        for line in infolines:
            plask.print_log('important', "  " + line)
 
        return self.threshold_voltage
 
    def _get_info(self):
        result = self._get_defines_info() + [
            "Threshold voltage (V):         {:8.3f}".format(self.threshold_voltage),
            "Threshold current (mA):        {:8.3f}".format(self.threshold_current),
            "Maximum temperature (K):       {:8.3f}".format(max(self.thermal.outTemperature(self.thermal.mesh)))
        ]
        levels = list(self._iter_levels(self.diffusion.geometry, self.diffusion.mesh, 'QW', 'QD', 'gain'))
        max_concentration = []
        for no, mesh in levels:
            value = self.diffusion.outCarriersConcentration(mesh)
            max_concentration.append(max(value))
        result.append("Maximum concentration (1/cm3):    {}"
                      .format(', '.join('{:.3e}'.format(c) for c in max_concentration)))
        lam = getattr(self.optical, self._lam0).real
        if lam is None: lam = self._get_lam().real
        max_gain = []
        for no, mesh in levels:
            value = self.gain.outGain(mesh, lam).array[:,:,0]
            max_gain.append(max(value))
        result.append("Maximum gain (1/cm):          {}"
                      .format(', '.join('{:9.3f}'.format(g[0]) for g in max_gain)))
        return result
 
    def save(self, filename=None, group='ThresholdSearch', optical_resolution=None):
        """
        Save the computation results to the HDF5 file.
 
        Args:
            filename (str): The file name to save to.
                If omitted, the file name is generated automatically based on
                the script name with suffix denoting either the batch job id or
                the current time if no batch system is used.
 
            group (str): HDF5 group to save the data under.
 
            optical_resolution (tuple of ints): Number of points in horizontal and vertical directions
                for optical field.
        """
        if optical_resolution is None: optical_resolution = self.optical_resolutionoptical_resolution
        h5file, group, filename, close = h5open(filename, group)
        self._save_thermoelectric(h5file, group)
        levels = list(self._iter_levels(self.diffusion.geometry, self.diffusion.mesh, 'QW', 'QD', 'gain'))
        for no, mesh in levels:
            value = self.diffusion.outCarriersConcentration(mesh)
            plask.save_field(value, h5file, group + '/Junction'+no+'CarriersConcentration')
        lam = getattr(self.optical, self._lam0).real
        if lam is None: lam = self._get_lam().real
        for no, mesh in levels:
            value = self.gain.outGain(mesh, lam)
            plask.save_field(value, h5file, group + '/Junction'+no+'Gain')
        if self.modeno is not None:
            obox = self.optical.geometry.bbox
            oaxis = plask.mesh.Regular(obox.left, obox.right, optical_resolution[0])
            omesh = plask.mesh.Rectangular2D(oaxis,
                                             plask.mesh.Regular(obox.bottom, obox.top, optical_resolution[1]))
            ofield = self.optical.outLightMagnitude(self.modeno, omesh)
            plask.save_field(ofield/max(ofield), h5file, group + '/LightMagnitude')
            try:
                nrfield = self.optical.outRefractiveIndex(omesh)
                plask.save_field(nrfield, h5file, group + '/RefractiveIndex')
            except AttributeError:
                epsfield = self.optical.outEpsilon(omesh)
                plask.save_field(epsfield, h5file, group + '/Epsilon')
            rmesh = next(self._iter_levels(self.optical.geometry, oaxis, 'QW', 'QD', 'gain'))[1]
            orfield = self.optical.outLightMagnitude(self.modeno, rmesh)
            plask.save_field(orfield/max(orfield), h5file, group + '/HorizontalLightMagnitude')
 
        if close:
            h5file.close()
        plask.print_log('info', "Fields saved to file '{}' in group '{}'".format(filename, group))
        return filename
 
    def get_optical_field(self, resolution=None):
        """
        Get computed optical mode field at threshold.
 
        Args:
            resolution (tuple of ints): Number of points in horizontal and vertical directions.
       """
        if resolution is None: resolution = self.optical_resolutionoptical_resolution
        box = self.optical.geometry.bbox
        intensity_mesh = plask.mesh.Rectangular2D(plask.mesh.Regular(box.left, box.right, resolution[0]),
                                                  plask.mesh.Regular(box.bottom, box.top, resolution[1]))
        field = self.optical.outLightMagnitude(self.modeno, intensity_mesh)
        return field
 
    def get_optical_field_horizontal(self, resolution=None, interpolation='linear'):
        """
        Get horizontal distribution of the computed optical mode field at threshold.
 
        Args:
            resolution (int): Number of points in horizontal direction.
 
            interpolation (str): Interpolation used when retrieving current density.
        """
        if resolution is None:
            resolution = self.optical_resolutionoptical_resolution[0]
        box = self.optical.geometry.bbox
        raxis = plask.mesh.Regular(box.left, box.right, resolution)
        rmesh = next(self._iter_levels(self.optical.geometry, raxis, 'QW', 'QD', 'gain'))[1]
        field = self.optical.outLightMagnitude(self.modeno, rmesh, interpolation)
        return field
 
    def get_optical_field_vertical(self, pos=0.01, offset=0.5, resolution=None, interpolation='linear'):
        """
        Plot vertical distribution of the computed optical mode field at threshold and
        refractive index profile.
 
        Args:
            resolution (int): Number of points in horizontal direction.
 
            pos (float): Horizontal position to get the field at.
 
            offset (float): Distance above and below geometry boundary to include into
                            the plot.
 
            interpolation (str): Interpolation used when retrieving current density.
        """
        if resolution is None:
            resolution = self.optical_resolutionoptical_resolution[1]
        box = self.optical.geometry.bbox
        zaxis = plask.mesh.Regular(box.bottom - offset, box.top + offset, 10 * resolution)
        zmesh = plask.mesh.Rectangular2D([pos], zaxis)
        field = self.optical.outLightMagnitude(self.modeno, zmesh, interpolation)
        return field
 
    def plot_optical_field(self, resolution=None, geometry_color='0.75', geometry_alpha=0.35, geometry_lw=1.0, **kwargs):
        """
        Plot computed optical mode field at threshold.
 
        Args:
            resolution (tuple of ints): Number of points in horizontal and vertical directions.
 
            geometry_color (str or ``None``): Matplotlib color specification
                for the geometry. If ``None``, structure is not plotted.
 
            geometry_alpha (float): Geometry opacity (1 — fully opaque, 0 – invisible).
 
            geometry_lw (float): Line width for geometry.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        field = self.get_optical_field(resolution)
        plask.plot_field(field, **kwargs)
        plask.plot_geometry(self.optical.geometry, color=geometry_color, lw=geometry_lw)
        plask.window_title("Light Intensity")
 
    def plot_optical_field_horizontal(self, resolution=None, bounds=True, interpolation='linear', **kwargs):
        """
        Plot horizontal distribution of the computed optical mode field at threshold.
 
        Args:
            resolution (int): Number of points in horizontal direction.
 
            bounds (bool): If *True* then the geometry objects boundaries are
                           plotted.
 
            interpolation (str): Interpolation used when retrieving current density.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        field = self.get_optical_field_horizontal(resolution, interpolation)
        color = kwargs.pop('color', None)
        if color is None:
            try:
                color = plask.rcParams['axes.prop_cycle']
            except KeyError:
                color = plask.rcParams['axes.color_cycle'][1]
            else:
                color = color.by_key()['color'][1]
        plask.plot_profile(field/max(field), color=color, **kwargs)
        if bounds:
            self._plot_hbounds(self.optical)
        plask.ylabel("Light Intensity [arb.u.]")
        plask.window_title("Radial Light Intensity")
 
    def plot_optical_field_vertical(self, pos=0.01, offset=0.5, resolution=None, interpolation='linear', **kwargs):
        """
        Plot vertical distribution of the computed optical mode field at threshold and
        refractive index profile.
 
        Args:
            resolution (int): Number of points in horizontal direction.
 
            pos (float): Horizontal position to get the field at.
 
            offset (float): Distance above and below geometry boundary to include into
                            the plot.
 
            interpolation (str): Interpolation used when retrieving current density.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        field = self.get_optical_field_vertical(pos, offset, resolution, interpolation)
        try:
            cc = plask.rcParams['axes.prop_cycle']
        except KeyError:
            cc = plask.rcParams['azmeshxes.color_cycle']
        else:
            cc = cc.by_key()['color']
        color2 = cc[0]
        color = kwargs.pop('color', None)
        if color is None:
            color = cc[1]
        ax1 = plask.gca()
        plask.plot_profile(field/max(field), color=color, **kwargs)
        plask.ylabel("Light Intensity [arb.u.]")
        ax2 = plask.twinx()
        try:
            plask.plot_profile(self.optical.outRefractiveIndex(field.mesh).real, comp=0, color=color2)
        except AttributeError:
            eps = self.optical.outEpsilon(field.mesh)
            nr = plask.Data((np.array(eps)[:,0,0]**0.5).real, eps.mesh)
            plask.plot_profile(nr, color=color2)
        plask.ylabel("Refractive Index")
        ax1.set_zorder(ax2.get_zorder()+1)
        ax1.patch.set_visible(False)
        ax2.patch.set_visible(True)
        plask.xlim(field.mesh.axis1[0], field.mesh.axis1[-1])
        plask.window_title("Vertical Light Intensity")
 
 
class ThresholdSearchCyl(ThresholdSearch):
    """
    Solver for threshold search of semiconductor laser.
 
    This solver performs thermo-electrical computations followed by
    determination ot threshold current and optical analysis in order to
    determine the threshold of a semiconductor laser. The search is
    performed by ``scipy`` root finding algorithm in order to determine
    the voltage and electric current ensuring no optical loss in the
    laser cavity.
 
    The computations can be executed using `compute` method, after which
    the results may be save to the HDF5 file with `save` or presented visually
    using ``plot_...`` methods. If ``save`` parameter of the :meth:`compute` method
    is *True* the fields are saved automatically after the computations.
    The file name is based on the name of the executed script with suffix denoting
    either the launch time or the identifier of a batch job if a batch system
    (like SLURM, OpenPBS, or SGE) is used.
    """
 
    _Thermal = thermal.static.StaticCyl
    _Electrical = electrical.shockley.ShockleyCyl
    _Diffusion = electrical.diffusion.DiffusionCyl
    _Gain = gain.freecarrier.FreeCarrierCyl
 
    _OPTICAL_ROOTS = {'optical-root': 'root', 'optical-stripe-root': 'stripe_root'}
 
    outTemperature = property(lambda self: self.thermalthermal.outTemperature, doc=_Thermal.outTemperature.__doc__)
    outHeatFlux = property(lambda self: self.thermalthermal.outHeatFlux, doc=_Thermal.outHeatFlux.__doc__)
 
    outThermalConductivity = property(lambda self: self.thermalthermal.outThermalConductivity,
                                      doc=_Thermal.outThermalConductivity.__doc__)
    outVoltage = property(lambda self: self.electricalelectrical.outVoltage, doc=_Electrical.outVoltage.__doc__)
    outCurrentDensity = property(lambda self: self.electricalelectrical.outCurrentDensity,
                                 doc=_Electrical.outCurrentDensity.__doc__)
    outHeat = property(lambda self: self.electricalelectrical.outHeat, doc=_Electrical.outHeat.__doc__)
    outConductivity = property(lambda self: self.electricalelectrical.outConductivity, doc=_Electrical.outConductivity.__doc__)
    outCarriersConcentration = property(lambda self: self.diffusiondiffusion.outCarriersConcentration,
                                        doc=_Diffusion.outCarriersConcentration.__doc__)
    outGain = property(lambda self: self.gaingain.outGain, doc=_Gain.outGain.__doc__)
    outLightMagnitude = property(lambda self: self.opticaloptical.outLightMagnitude, doc=_doc.outLightMagnitude)
    outLoss = property(lambda self: self.opticaloptical.outLoss, doc=_doc.outLoss)
    outWavelength = property(lambda self: self.opticaloptical.outWavelength, doc=_doc.outWavelength)
    outRefractiveIndex = property(lambda self: self.opticaloptical.outRefractiveIndex, doc=_doc.outRefractiveIndex)
    outLightE = property(lambda self: self.opticaloptical.outLightE, doc=_doc.outLightE)
 
    thermal = attribute(_Thermal.__name__+"()")
    ":class:`thermal.static.StaticCyl` solver used for thermal calculations."
 
    electrical = attribute(_Electrical.__name__+"()")
    ":class:`electrical.shockley.ShockleyCyl` solver used for electrical calculations."
 
    diffusion = attribute(_Diffusion.__name__+"()")
    ":class:`electrical.diffusion.DiffusionCyl` solver used for electrical calculations."
 
    gain = attribute(_Gain.__name__+"()")
    ":class:`gain.freecarrier.FreeCarrierCyl` solver used for gain calculations."
 
    optical = attribute("EffectiveFrequencyCyl()")
    ":class:`optical.effective.EffectiveFrequencyCyl` solver used for optical calculations."
 
    tfreq = 6.0
    """
    Number of electrical iterations per single thermal step.
 
    As temperature tends to converge faster, it is reasonable to repeat thermal
    solution less frequently.
    """
 
    vmin = None
    """
    Minimum voltage to search threshold for.
 
    It should be below the threshold.
    """
 
    vmax = None
    """
    Maximum voltage to search threshold for.
 
    It should be above the threshold.
    """
 
    vtol = 1e-5
    "Tolerance on voltage in the root search."
 
    maxiter = 50
    "Maximum number of root finding iterations."
 
    maxlam = attribute("optical.lam0")
    "Maximum wavelength considered for the optical mode search."
 
    dlam = 0.02
    """
    Wavelength step.
 
    Step, by which the wavelength is swept while searching for the approximate mode.
    """
 
    lpm = 0
    """
    Angular mode number $m$.
 
    0 for LP0x, 1 for LP1x, etc.
    """
 
    lpn = 1
    """
    Radial mode number $n$.
 
    1 for LPx1, 2 for LPx2, etc.
    """
 
    optical_resolution = (800, 600)
    """
    Number of points along the horizontal and vertical axes for the saved
    and plotted optical field.
    """
 
    skip_thermal = False
    """
    Skip thermal computations.
 
    The structure is assumed to have a constant temperature.
    This can be used to look for the threshold under pulse laser operation.
    """
 
    def __init__(self, name=''):
        from optical.effective import EffectiveFrequencyCyl
        self._Optical_Optical = EffectiveFrequencyCyl
        super().__init__(name)
        self.maxlam = None
 
    # def on_initialize(self):
    #     super().on_initialize()
 
    def get_lam(self):
        """
        Get approximate wavelength for optical computations.
 
        This method returns approximate wavelength for optical computations.
        By default if browses the wavelength range starting from :attr:`maxlam`,
        decreasing it by :attr:`dlam` until radial mode :attr:`lpn` is found.
 
        You can override this method or set it to a a fixed value to use custom
        mode approximation.
 
        Example:
             >>> solver = ThresholdSearchCyl()
             >>> solver.get_lam = 980.
             >>> solver.compute()
        """
 
        lam = self.maxlam if self.maxlam is not None else self.opticaloptical.lam0
        n = 0
        prev = 0.
        decr = False
        while n < self.lpnlpn and lam.real > 0.:
            curr = abs(self.opticaloptical.get_determinant(lam=lam, m=self.lpmlpm))
            if decr and curr > prev:
                n += 1
            decr = curr < prev
            prev = curr
            lam -= self.dlamdlam
        if n == self.lpnlpn:
            return lam + 2. * self.dlamdlam
        raise ValueError("Approximation of mode LP{0.lpm}{0.lpn} not found".format(self))
 
    def _optargs(self):
        return dict(m=self.lpmlpm)
 
    def _parse_xpl(self, tag, manager):
        if tag == 'optical':
            self._read_attr(tag, 'lam0', self.opticaloptical, float)
            self._read_attr(tag, 'vlam', self.opticaloptical, float)
            self._read_attr(tag, 'vat', self.opticaloptical, float)
            self._read_attr(tag, 'emission', self.opticaloptical)
            maxlam = tag.get('maxlam')
            if maxlam is not None:
                try:
                    self.maxlam = complex(maxlam)
                except ValueError:
                    raise plask.XMLError("{}: attribute maxlam has illegal value '{}'".format(tag, maxlam))
            self.dlamdlam = float(tag.get('dlam', self.dlamdlam))
            self.lpmlpm = int(tag.get('m', self.lpmlpm))
            self.lpnlpn = int(tag.get('n', self.lpnlpn))
        else:
            if tag == 'optical-root':
                self._read_attr(tag, 'determinant', self.opticaloptical, str, 'determinant_mode')
            super()._parse_xpl(tag, manager)
 
    def get_vert_optical_determinant(self, vlam):
        """
        Function computing ‘vertical determinant’ of the optical solver.
 
        Args:
             vlam (float or array): ‘Vertical wavelength’ to compute the vertical
                                     determinant for (nm).
 
        Returns:
            float or array: Optical vertical determinant.
        """
        self.compute_thermoelectric()
        if self._optarg is None:
            return self.opticaloptical.get_determinant(vlam, **self._optargs_optargs())
        else:
            _optargs = self._optargs_optargs().copy()
            _optargs[self._optarg] = vlam
            return self.opticaloptical.get_determinant(**_optargs)
 
    def plot_vert_optical_determinant(self, vlams, **kwargs):
        """
        Function plotting ‘vertical determinant’ of the optical solver.
 
        Args:
            vlams (array): ‘Vertical wavelengths’ to plot the determinant for (nm).
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        vals = self.get_vert_optical_determinant(vlams)
        plask.plot(vlams, abs(vals))
        plask.yscale('log')
        plask.xlabel("Vertical Wavelength (nm)")
        plask.ylabel("Determinant [ar.u.]")
 
    def _get_info(self):
        return super()._get_info() + [
            "LP{}{} mode wavelength (nm):     {:8.3f}".format(self.lpmlpm, self.lpnlpn, self.opticaloptical.modes[self.modeno].lam.real)
        ]
 
 
class ThresholdSearchBesselCyl(ThresholdSearch):
    """
    Solver for threshold search of semiconductor laser with vector optical solver.
 
    This solver performs thermo-electrical computations followed by
    determination ot threshold current and optical analysis in order to
    determine the threshold of a semiconductor laser. The search is
    performed by ``scipy`` root finding algorithm in order to determine
    the voltage and electric current ensuring no optical loss in the
    laser cavity.
 
    This solver uses vector optical solver :class:`~plask.optical.modal.BesselCyl`.
 
    The computations can be executed using `compute` method, after which
    the results may be save to the HDF5 file with `save` or presented visually
    using ``plot_...`` methods. If ``save`` parameter of the :meth:`compute` method
    is *True* the fields are saved automatically after the computations.
    The file name is based on the name of the executed script with suffix denoting
    either the launch time or the identifier of a batch job if a batch system
    (like SLURM, OpenPBS, or SGE) is used.
    """
 
    _Thermal = thermal.static.StaticCyl
    _Electrical = electrical.shockley.ShockleyCyl
    _Diffusion = electrical.diffusion.DiffusionCyl
    _Gain = gain.freecarrier.FreeCarrierCyl
 
    _OPTICAL_ROOTS = {'optical-root': 'root'}
 
    outTemperature = property(lambda self: self.thermalthermal.outTemperature, doc=_Thermal.outTemperature.__doc__)
    outHeatFlux = property(lambda self: self.thermalthermal.outHeatFlux, doc=_Thermal.outHeatFlux.__doc__)
 
    outThermalConductivity = property(lambda self: self.thermalthermal.outThermalConductivity,
                                      doc=_Thermal.outThermalConductivity.__doc__)
    outVoltage = property(lambda self: self.electricalelectrical.outVoltage, doc=_Electrical.outVoltage.__doc__)
    outCurrentDensity = property(lambda self: self.electricalelectrical.outCurrentDensity,
                                 doc=_Electrical.outCurrentDensity.__doc__)
    outHeat = property(lambda self: self.electricalelectrical.outHeat, doc=_Electrical.outHeat.__doc__)
    outConductivity = property(lambda self: self.electricalelectrical.outConductivity, doc=_Electrical.outConductivity.__doc__)
    outCarriersConcentration = property(lambda self: self.diffusiondiffusion.outCarriersConcentration,
                                        doc=_Diffusion.outCarriersConcentration.__doc__)
    outGain = property(lambda self: self.gaingain.outGain, doc=_Gain.outGain.__doc__)
    outLightMagnitude = property(lambda self: self.opticaloptical.outLightMagnitude, doc=_doc.outLightMagnitude)
    outLoss = property(lambda self: self.opticaloptical.outLoss, doc=_doc.outLoss)
    outWavelength = property(lambda self: self.opticaloptical.outWavelength, doc=_doc.outWavelength)
    outEpsilon = property(lambda self: self.opticaloptical.outEpsilon, doc=_doc.outEpsilon)
    outLightE = property(lambda self: self.opticaloptical.outLightE, doc=_doc.outLightE)
 
    thermal = attribute(_Thermal.__name__ + "()")
    ":class:`thermal.static.StaticCyl` solver used for thermal calculations."
 
    electrical = attribute(_Electrical.__name__ + "()")
    ":class:`electrical.shockley.ShockleyCyl` solver used for electrical calculations."
 
    diffusion = attribute(_Diffusion.__name__ + "()")
    ":class:`electrical.diffusion.DiffusionCyl` solver used for electrical calculations."
 
    gain = attribute(_Gain.__name__ + "()")
    ":class:`gain.freecarrier.FreeCarrierCyl` solver used for gain calculations."
 
    optical = attribute("BesselCyl()")
    ":class:`optical.modal.BesselFrequencyCyl` solver used for optical calculations."
 
    tfreq = 6.0
    """
    Number of electrical iterations per single thermal step.
 
    As temperature tends to converge faster, it is reasonable to repeat thermal
    solution less frequently.
    """
 
    vmin = None
    """
    Minimum voltage to search threshold for.
 
    It should be below the threshold.
    """
 
    vmax = None
    """
    Maximum voltage to search threshold for.
 
    It should be above the threshold.
    """
 
    vtol = 1e-5
    "Tolerance on voltage in the root search."
 
    maxiter = 50
    "Maximum number of root finding iterations."
 
    maxlam = attribute("optical.lam0")
    "Maximum wavelength considered for the optical mode search."
 
    dlam = 0.05
    """
    Wavelength step.
 
    Step, by which the wavelength is swept while searching for the approximate mode.
    """
 
    hem = 1
    """
      Angular mode number $m$.
 
      1 for HE1x, 2 for HE2x, etc.
    """
 
    hen = 1
    """
      Radial mode number $n$.
 
      1 for HEx1, 2 for HEx2, etc.
    """
 
    lam = None
    """
    Initial wavelength for optical search.
 
    If this value is set, the computations are started from this value. If this
    value is set, the radial mode number :attr:`hen` is ignored.
 
    Note that it is safer to leave this empty and allow the solver to look for it
    automatically, however, it may increase the time of optical computations.
    """
 
    optical_resolution = (800, 600)
    """
    Number of points along the horizontal and vertical axes for the saved
    and plotted optical field.
    """
 
    skip_thermal = False
    """
    Skip thermal computations.
 
    The structure is assumed to have a constant temperature.
    This can be used to look for the threshold under pulse laser operation.
    """
 
    def __init__(self, name=''):
        from optical.modal import BesselCyl
        self._Optical_Optical = BesselCyl
        super().__init__(name)
        self.maxlam = None
 
    # def on_initialize(self):
    #     super().on_initialize()
 
    # def on_invalidate(self):
    #     super().on_invalidate()
 
    def get_lam(self):
        """
        Get approximate wavelength for optical computations.
 
        This method returns approximate wavelength for optical computations.
        By default if browses the wavelength range starting from :attr:`maxlam`,
        decreasing it by :attr:`dlam` until radial mode :attr:`hen` is found.
 
        You can override this method or set it to a a fixed value to use custom
        mode approximation.
 
        Example:
             >>> solver = ThresholdSearchBesselCyl()
             >>> solver.get_lam = 980.
             >>> solver.compute()
        """
 
        if self.lam is not None:
            return self.lam
 
        lam = self.maxlam if self.maxlam is not None else self.opticaloptical.lam0
        n = 0
        prev = 0.
        decr = False
        while n < self.henhen and lam.real > 0.:
            curr = abs(self.opticaloptical.get_determinant(lam=lam, m=self.hemhem))
            if decr and curr > prev:
                n += 1
            decr = curr < prev
            prev = curr
            lam -= self.dlamdlam
        if n == self.henhen:
            return lam + 2. * self.dlamdlam
        raise ValueError("Approximation of mode HE{0.hem}{0.hen} not found".format(self))
 
    def _optargs(self):
        return dict(m=self.hemhem)
 
    def _parse_xpl(self, tag, manager):
        if tag == 'optical':
            self._read_attr(tag, 'lam0', self.opticaloptical, float)
            self._read_attr(tag, 'update-gain', self.opticaloptical, bool, 'update_gain')
            self._read_attr(tag, 'domain', self.opticaloptical)
            self._read_attr(tag, 'size', self.opticaloptical, int, 'size')
            self._read_attr(tag, 'group-layers', self.opticaloptical, bool, 'group_layers')
            self._read_attr(tag, 'k-method', self.opticaloptical, pyattr='kmethod')
            self._read_attr(tag, 'k-scale', self.opticaloptical, float, 'kscale')
            self._read_attr(tag, 'transfer', self.opticaloptical)
            self._read_attr(tag, 'emission', self.opticaloptical)
            maxlam = tag.get('maxlam')
            if maxlam is not None:
                try:
                    self.maxlam = complex(maxlam)
                except ValueError:
                    raise plask.XMLError("{}: attribute maxlam has illegal value '{}'".format(tag, maxlam))
            self.dlamdlam = float(tag.get('dlam', self.dlamdlam))
            self.lam = tag.get('lam', self.lam)
            if self.lam is not None: self.lam = complex(self.lam)
            self.hemhem = int(tag.get('m', self.hemhem))
            self.henhen = int(tag.get('n', self.henhen))
        elif tag == 'optical-interface':
            attrs = {key: val for (key, val) in ((key, tag.get(key)) for key in ('position', 'object', 'path'))
                     if val is not None}
            if len(attrs) > 1 and (len(attrs) > 2 or 'position' in attrs):
                raise plask.XMLError("{}: conflicting attributes '{}'".format(tag, "' and '".join(attrs.keys())))
            elif 'position' in attrs:
                self.opticaloptical.set_interface(attrs['position'])
            elif 'object' in attrs:
                path = attrs.get('path')
                if path is not None:
                    self.opticaloptical.set_interface(manager.geo[attrs['object']], manager.pth[path])
                else:
                    self.opticaloptical.set_interface(manager.geo[attrs['object']])
        elif tag == 'optical-vpml':
            self._read_attr(tag, 'factor', self.opticaloptical.vpml, complex, 'factor')
            self._read_attr(tag, 'dist', self.opticaloptical.vpml, float, 'dist')
            self._read_attr(tag, 'size', self.opticaloptical.vpml, float, 'size')
        elif tag == 'optical-pml':
            self._read_attr(tag, 'factor', self.opticaloptical.pml, complex, 'factor')
            self._read_attr(tag, 'shape', self.opticaloptical.pml, float, 'shape')
            self._read_attr(tag, 'dist', self.opticaloptical.pml, float, 'dist')
            self._read_attr(tag, 'size', self.opticaloptical.pml, float, 'size')
        else:
            super()._parse_xpl(tag, manager)
 
    def _get_info(self):
        return super()._get_info() + [
            "HE{}{} mode wavelength (nm):     {:8.3f}".format(self.hemhem, self.henhen, self.opticaloptical.modes[self.modeno].lam.real)
        ]
 
 
class ThresholdSearch2D(ThresholdSearch):
    """
    Solver for threshold search of semiconductor laser.
 
    This solver performs thermo-electrical computations followed by
    determination ot threshold current and optical analysis in order to
    determine the threshold of a semiconductor laser. The search is
    performed by ``scipy`` root finding algorithm in order to determine
    the voltage and electric current ensuring no optical loss in the
    laser cavity.
 
    The computations can be executed using `compute` method, after which
    the results may be save to the HDF5 file with `save` or presented visually
    using ``plot_...`` methods. If ``save`` parameter of the :meth:`compute` method
    is *True* the fields are saved automatically after the computations.
    The file name is based on the name of the executed script with suffix denoting
    either the launch time or the identifier of a batch job if a batch system
    (like SLURM, OpenPBS, or SGE) is used.
    """
 
    _optarg = 'neff'
    _lam0 = 'wavelength'
 
    _Thermal = thermal.static.Static2D
    _Electrical = electrical.shockley.Shockley2D
    _Diffusion = electrical.diffusion.Diffusion2D
    _Gain = gain.freecarrier.FreeCarrier2D
 
    _OPTICAL_ROOTS = {'optical-root': 'root', 'optical-stripe-root': 'stripe_root'}
 
    outTemperature = property(lambda self: self.thermalthermal.outTemperature, doc=_Thermal.outTemperature.__doc__)
    outHeatFlux = property(lambda self: self.thermalthermal.outHeatFlux, doc=_Thermal.outHeatFlux.__doc__)
 
    outThermalConductivity = property(lambda self: self.thermalthermal.outThermalConductivity,
                                      doc=_Thermal.outThermalConductivity.__doc__)
    outVoltage = property(lambda self: self.electricalelectrical.outVoltage, doc=_Electrical.outVoltage.__doc__)
    outCurrentDensity = property(lambda self: self.electricalelectrical.outCurrentDensity,
                                 doc=_Electrical.outCurrentDensity.__doc__)
    outHeat = property(lambda self: self.electricalelectrical.outHeat, doc=_Electrical.outHeat.__doc__)
    outConductivity = property(lambda self: self.electricalelectrical.outConductivity, doc=_Electrical.outConductivity.__doc__)
    outCarriersConcentration = property(lambda self: self.diffusiondiffusion.outCarriersConcentration,
                                        doc=_Diffusion.outCarriersConcentration.__doc__)
    outGain = property(lambda self: self.gaingain.outGain, doc=_Gain.outGain.__doc__)
    outLightMagnitude = property(lambda self: self.opticaloptical.outLightMagnitude, doc=_doc.outLightMagnitude)
    outLoss = property(lambda self: self.opticaloptical.outLoss, doc=_doc.outLoss)
    outNeff = property(lambda self: self.opticaloptical.outNeff, doc=_doc.outNeff)
    outRefractiveIndex = property(lambda self: self.opticaloptical.outRefractiveIndex, doc=_doc.outRefractiveIndex)
    outLightE = property(lambda self: self.opticaloptical.outLightE, doc=_doc.outLightE)
 
    thermal = attribute(_Thermal.__name__+"()")
    ":class:`thermal.static.Static2D` solver used for thermal calculations."
 
    electrical = attribute(_Electrical.__name__+"()")
    ":class:`electrical.shockley.Shockley2D` solver used for electrical calculations."
 
    diffusion = attribute(_Diffusion.__name__+"()")
    ":class:`electrical.diffusion.Diffusion2D` solver used for electrical calculations."
 
    gain = attribute(_Gain.__name__+"()")
    ":class:`gain.freecarrier.FreeCarrier2D` solver used for gain calculations."
 
    optical = attribute("EffectiveIndex2D()")
    ":class:`optical.effective.EffectiveIndex2D` solver used for optical calculations."
 
    tfreq = 6.0
    """
    Number of electrical iterations per single thermal step.
 
    As temperature tends to converge faster, it is reasonable to repeat thermal
    solution less frequently.
    """
 
    vmin = None
    """
    Minimum voltage to search threshold for.
 
    It should be below the threshold.
    """
 
    vmax = None
    """
    Maximum voltage to search threshold for.
 
    It should be above the threshold.
    """
 
    vtol = 1e-5
    "Tolerance on voltage in the root search."
 
    maxiter = 50
    "Maximum number of root finding iterations."
 
    wavelength = None
    "Emission wavelength (nm)."
 
    dneff = 0.02
    """
    Effective index step.
 
    Step, by which the effective index is swept while searching for the approximate mode.
    """
 
    mn = 1
    """
    Lateral mode number $n$.
    """
 
    optical_resolution = (800, 600)
    """
    Number of points along the horizontal and vertical axes for the saved
    and plotted optical field.
    """
 
    skip_thermal = False
    """
    Skip thermal computations.
 
    The structure is assumed to have a constant temperature.
    This can be used to look for the threshold under pulse laser operation.
    """
 
    def __init__(self, name=''):
        from optical.effective import EffectiveIndex2D
        self._Optical_Optical = EffectiveIndex2D
        super().__init__(name)
 
    def on_initialize(self):
        super().on_initialize()
        points = plask.mesh.Rectangular2D.SimpleGenerator()(self.opticaloptical.geometry).elements.mesh
        self._maxneff = max(self.opticaloptical.geometry.get_material(point).Nr(self.opticaloptical.wavelength.real).real
                            for point in points)
 
    def get_nng(self):
        """
        Get approximate effective index for optical computations.
 
        This method returns approximate wavelength for optical computations.
        By default if browses the wavelength range starting from :attr:`maxneff`,
        decreasing it by :attr:`dneff` until lateral mode :attr:`mn` is found.
 
        You can override this method to use custom mode approximation.
 
        Example:
             >>> solver = ThresholdSearch2D()
             >>> solver.get_nng = lambda: 3.5
             >>> solver.compute()
        """
 
        neff = self._maxneff
        n = 0
        prev = 0.
        decr = False
        while n < self.mnmn and neff.real > 0.:
            curr = abs(self.opticaloptical.get_determinant(neff))
            if decr and curr > prev:
                n += 1
            decr = curr < prev
            prev = curr
            neff -= self.dneffdneff
        if n == self.mnmn:
            return neff + 2. * self.dneffdneff
        raise ValueError("Mode approximation not found")
 
    def _parse_xpl(self, tag, manager):
        if tag == 'optical':
            self._read_attr(tag, 'lam', self.opticaloptical, float)
            self._read_attr(tag, 'polarization', self.opticaloptical, float)
            self._read_attr(tag, 'vneff', self.opticaloptical, float)
            self._read_attr(tag, 'vat', self.opticaloptical, float)
            self._read_attr(tag, 'emission', self.opticaloptical)
            self.dneffdneff = float(tag.get('dneff', self.dneffdneff))
            self.mnmn = int(tag.get('mn', self.mnmn))
        else:
            if tag == 'optical-root':
                self._read_attr(tag, 'determinant', self.opticaloptical, str)
            super()._parse_xpl(tag, manager)
 
    def get_optical_determinant(self, neff):
        """
        Function computing determinant of the optical solver.
 
        Args:
             neff (float or array): Effective index to compute the determinant for.
 
        Returns:
            float or array: Optical determinant.
        """
        self.compute_thermoelectric()
        if self._optarg_optarg is None:
            return self.opticaloptical.get_determinant(neff, **self._optargs())
        else:
            _optargs = self._optargs().copy()
            _optargs[self._optarg_optarg] = neff
            return self.opticaloptical.get_determinant(**_optargs)
 
    def plot_optical_determinant(self, neffs, **kwargs):
        """
        Function plotting determinant of the optical solver.
 
        Args:
            neffs (array): Array of effective indices to plot the determinant for.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        vals = self.get_optical_determinantget_optical_determinant(neffs)
        plask.plot(neffs, abs(vals))
        plask.yscale('log')
        plask.xlabel("Effective Index")
        plask.ylabel("Determinant [ar.u.]")
 
    def get_vert_optical_determinant(self, vneff):
        """
        Function computing ‘vertical determinant’ of the optical solver.
 
        Args:
             vneff (float or array): Effective index to compute the vertical
                                     determinant for.
 
        Returns:
            float or array: Optical determinant.
        """
        self.compute_thermoelectric()
        if self._optarg_optarg is None:
            return self.opticaloptical.get_determinant(vneff, **self._optargs())
        else:
            _optargs = self._optargs().copy()
            _optargs[self._optarg_optarg] = vneff
            return self.opticaloptical.get_determinant(**_optargs)
 
    def plot_vert_optical_determinant(self, vneffs, **kwargs):
        """
        Function plotting ‘vertical determinant’ of the optical solver.
 
        Args:
            vneffs (array): Array of effective indices to plot the vertical
                            determinant for.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        vals = self.get_vert_optical_determinant(vneffs)
        plask.plot(vneffs, abs(vals))
        plask.yscale('log')
        plask.xlabel("Vertical Effective Index")
        plask.ylabel("Determinant [ar.u.]")
 
    def _get_info(self):
        return super()._get_info() + [
            "Effective index:     {:8.3f}".format(self.opticaloptical.modes[self.modeno].neff.real)
        ]
 
 
class ThresholdSearchFourier2D(ThresholdSearch):
    """
    Solver for threshold search of semiconductor laser.
 
    This solver performs thermo-electrical computations followed by
    determination ot threshold current and optical analysis in order to
    determine the threshold of a semiconductor laser. The search is
    performed by ``scipy`` root finding algorithm in order to determine
    the voltage and electric current ensuring no optical loss in the
    laser cavity.
 
    This solver uses vector optical solver :class:`~plask.optical.modal.Fourier2D`.
 
    The computations can be executed using `compute` method, after which
    the results may be save to the HDF5 file with `save` or presented visually
    using ``plot_...`` methods. If ``save`` parameter of the :meth:`compute` method
    is *True* the fields are saved automatically after the computations.
    The file name is based on the name of the executed script with suffix denoting
    either the launch time or the identifier of a batch job if a batch system
    (like SLURM, OpenPBS, or SGE) is used.
    """
 
    _optarg = 'neff'
    _lam0 = 'lam'
 
    _Thermal = thermal.static.Static2D
    _Electrical = electrical.shockley.Shockley2D
    _Diffusion = electrical.diffusion.Diffusion2D
    _Gain = gain.freecarrier.FreeCarrier2D
 
    _OPTICAL_ROOTS = {'optical-root': 'root'}
 
    outTemperature = property(lambda self: self.thermalthermal.outTemperature, doc=_Thermal.outTemperature.__doc__)
    outHeatFlux = property(lambda self: self.thermalthermal.outHeatFlux, doc=_Thermal.outHeatFlux.__doc__)
 
    outThermalConductivity = property(lambda self: self.thermalthermal.outThermalConductivity,
                                      doc=_Thermal.outThermalConductivity.__doc__)
    outVoltage = property(lambda self: self.electricalelectrical.outVoltage, doc=_Electrical.outVoltage.__doc__)
    outCurrentDensity = property(lambda self: self.electricalelectrical.outCurrentDensity,
                                 doc=_Electrical.outCurrentDensity.__doc__)
    outHeat = property(lambda self: self.electricalelectrical.outHeat, doc=_Electrical.outHeat.__doc__)
    outConductivity = property(lambda self: self.electricalelectrical.outConductivity, doc=_Electrical.outConductivity.__doc__)
    outCarriersConcentration = property(lambda self: self.diffusiondiffusion.outCarriersConcentration,
                                        doc=_Diffusion.outCarriersConcentration.__doc__)
    outGain = property(lambda self: self.gaingain.outGain, doc=_Gain.outGain.__doc__)
    outLightMagnitude = property(lambda self: self.opticaloptical.outLightMagnitude, doc=_doc.outLightMagnitude)
    outNeff = property(lambda self: self.opticaloptical.outNeff, doc=_doc.outNeff)
    outEpsilon = property(lambda self: self.opticaloptical.outEpsilon, doc=_doc.outEpsilon)
    outLightE = property(lambda self: self.opticaloptical.outLightE, doc=_doc.outLightE)
 
    thermal = attribute(_Thermal.__name__+"()")
    ":class:`thermal.static.Static2D` solver used for thermal calculations."
 
    electrical = attribute(_Electrical.__name__+"()")
    ":class:`electrical.shockley.Shockley2D` solver used for electrical calculations."
 
    diffusion = attribute(_Diffusion.__name__+"()")
    ":class:`electrical.diffusion.Diffusion2D` solver used for electrical calculations."
 
    gain = attribute(_Gain.__name__+"()")
    ":class:`gain.freecarrier.FreeCarrier2D` solver used for gain calculations."
 
    optical = attribute("Fourier2D()")
    ":class:`optical.modal.Fourier2D` solver used for optical calculations."
 
    tfreq = 6.0
    """
    Number of electrical iterations per single thermal step.
 
    As temperature tends to converge faster, it is reasonable to repeat thermal
    solution less frequently.
    """
 
    vmin = None
    """
    Minimum voltage to search threshold for.
 
    It should be below the threshold.
    """
 
    vmax = None
    """
    Maximum voltage to search threshold for.
 
    It should be above the threshold.
    """
 
    vtol = 1e-5
    "Tolerance on voltage in the root search."
 
    maxiter = 50
    "Maximum number of root finding iterations."
 
    wavelength = None
    "Emission wavelength (nm)."
 
    dneff = 0.02
    """
    Effective index step.
 
    Step, by which the effective index is swept while searching for the approximate mode.
    """
 
    mn = 1
    """
    Lateral mode number $n$.
    """
 
    optical_resolution = (800, 600)
    """
    Number of points along the horizontal and vertical axes for the saved
    and plotted optical field.
    """
 
    skip_thermal = False
    """
    Skip thermal computations.
 
    The structure is assumed to have a constant temperature.
    This can be used to look for the threshold under pulse laser operation.
    """
 
    def __init__(self, name=''):
        from optical.modal import Fourier2D
        self._Optical_Optical = Fourier2D
        super().__init__(name)
 
    def on_initialize(self):
        super().on_initialize()
        points = plask.mesh.Rectangular2D.SimpleGenerator()(self.opticaloptical.geometry).elements.mesh
        self._maxneff = max(self.opticaloptical.geometry.get_material(point).Nr(self.opticaloptical.wavelength.real).real
                            for point in points)
 
    def get_nng(self):
        """
        Get approximate effective index for optical computations.
 
        This method returns approximate wavelength for optical computations.
        By default if browses the wavelength range starting from :attr:`maxneff`,
        decreasing it by :attr:`dneff` until lateral mode :attr:`mn` is found.
 
        You can override this method to use custom mode approximation.
 
        Example:
             >>> solver = ThresholdSearchFourier2D()
             >>> solver.get_nng = lambda: 3.5
             >>> solver.compute()
        """
 
        neff = self._maxneff
        n = 0
        prev = 0.
        decr = False
        while n < self.mnmn and neff.real > 0.:
            curr = abs(self.opticaloptical.get_determinant(neff=neff))
            if decr and curr > prev:
                n += 1
            decr = curr < prev
            prev = curr
            neff -= self.dneffdneff
        if n == self.mnmn:
            return neff + 2. * self.dneffdneff
        raise ValueError("Mode approximation not found")
 
    def _parse_xpl(self, tag, manager):
        if tag == 'optical':
            self._read_attr(tag, 'size', self.opticaloptical, int)
            self._read_attr(tag, 'refine', self.opticaloptical, int)
            self._read_attr(tag, 'smooth', self.opticaloptical, float)
            self._read_attr(tag, 'symmetry', self.opticaloptical)
            self._read_attr(tag, 'group-layers', self.opticaloptical, bool, 'group_layers')
            self._read_attr(tag, 'transfer', self.opticaloptical)
            self._read_attr(tag, 'lam', self.opticaloptical, float)
            self.dneffdneff = float(tag.get('dneff', self.dneffdneff))
            self.mnmn = int(tag.get('mn', self.mnmn))
        elif tag == 'optical-interface':
            attrs = {key: val for (key, val) in ((key, tag.get(key)) for key in ('position', 'object', 'path'))
                     if val is not None}
            if len(attrs) > 1 and (len(attrs) > 2 or 'position' in attrs):
                raise plask.XMLError("{}: conflicting attributes '{}'".format(tag, "' and '".join(attrs.keys())))
            elif 'position' in attrs:
                self.opticaloptical.set_interface(attrs['position'])
            elif 'object' in attrs:
                path = attrs.get('path')
                if path is not None:
                    self.opticaloptical.set_interface(manager.geo[attrs['object']], manager.pth[path])
                else:
                    self.opticaloptical.set_interface(manager.geo[attrs['object']])
        elif tag == 'optical-vpml':
            self._read_attr(tag, 'factor', self.opticaloptical.vpml, complex, 'factor')
            self._read_attr(tag, 'dist', self.opticaloptical.vpml, float, 'dist')
            self._read_attr(tag, 'size', self.opticaloptical.vpml, float, 'size')
        elif tag == 'optical-pml':
            self._read_attr(tag, 'factor', self.opticaloptical.pml, complex, 'factor')
            self._read_attr(tag, 'shape', self.opticaloptical.pml, float, 'shape')
            self._read_attr(tag, 'dist', self.opticaloptical.pml, float, 'dist')
            self._read_attr(tag, 'size', self.opticaloptical.pml, float, 'size')
        else:
            super()._parse_xpl(tag, manager)
 
    def get_optical_determinant(self, neff):
        """
        Function computing determinant of the optical solver.
 
        Args:
             neff (float or array): Effective index to compute the determinant for.
 
        Returns:
            float or array: Optical determinant.
        """
        self.compute_thermoelectric()
        if self._optarg_optarg is None:
            return self.opticaloptical.get_determinant(neff, **self._optargs())
        else:
            _optargs = self._optargs().copy()
            _optargs[self._optarg_optarg] = neff
            return self.opticaloptical.get_determinant(**_optargs)
 
    def plot_optical_determinant(self, neffs, **kwargs):
        """
        Function plotting determinant of the optical solver.
 
        Args:
            neffs (array): Array of effective indices to plot the determinant for.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        vals = self.get_optical_determinantget_optical_determinant(neff=neffs)
        plask.plot(neffs, abs(vals))
        plask.yscale('log')
        plask.xlabel("Effective Index")
        plask.ylabel("Determinant [ar.u.]")
 
    def plot_vert_optical_determinant(self, vneffs, **kwargs):
        """
        Function plotting ‘vertical determinant’ of the optical solver.
 
        Args:
            vneffs (array): Array of effective indices to plot the vertical
                            determinant for.
 
            **kwargs: Keyword arguments passed to the plot function.
        """
        vals = self.get_vert_optical_determinant(vneffs)
        plask.plot(vneffs, abs(vals))
        plask.yscale('log')
        plask.xlabel("Vertical Effective Index")
        plask.ylabel("Determinant [ar.u.]")
 
    def _get_info(self):
        return super()._get_info() + [
            "Effective index:     {:8.3f}".format(self.opticaloptical.modes[self.modeno].neff.real)
        ]
 
 
__all__ = 'ThresholdSearchCyl', 'ThresholdSearchBesselCyl', 'ThresholdSearch2D', 'ThresholdSearchFourier2D'