fourier2d-python.cpp

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.
 */
#define PY_ARRAY_UNIQUE_SYMBOL PLASK_OPTICAL_SLAB_ARRAY_API
#define NO_IMPORT_ARRAY
 
#include "fourier2d-python.hpp"
#include "modal-python.hpp"
 
namespace plask { namespace optical { namespace modal { namespace python {
 
template <> inline const char* solver_compute_reflectivity_name<FourierSolver2D>() {
    return "Fourier2D.compute_reflectivity";
}
 
template <> inline const char* solver_compute_transmittivity_name<FourierSolver2D>() {
    return "Fourier2D.compute_transmittivity";
}
 
 
template <> py::object Eigenmodes<FourierSolver2D>::array(const dcomplex* data, size_t N) const {
    int dim = 2, strid = 2;
    if (solver.separated()) strid = dim = 1;
    npy_intp dims[] = {npy_intp(N / strid), npy_intp(strid)};
    npy_intp strides[] = {npy_intp(strid * sizeof(dcomplex)), npy_intp(sizeof(dcomplex))};
    PyObject* arr = PyArray_New(&PyArray_Type, dim, dims, NPY_CDOUBLE, strides, (void*)data, 0, 0, NULL);
    if (arr == nullptr) throw plask::CriticalException("cannot create array");
    return py::object(py::handle<>(arr));
}
 
inline static std::string polarization_str(Expansion::Component val) {
    AxisNames* axes = getCurrentAxes();
    switch (val) {
        case Expansion::E_TRAN: return "E" + axes->getNameForTran();
        case Expansion::E_LONG: return "E" + axes->getNameForLong();
        default: return "none";
    }
}
 
static py::object FourierSolver2D_getMirrors(const FourierSolver2D& self) {
    if (!self.mirrors) return py::object();
    return py::make_tuple(self.mirrors->first, self.mirrors->second);
}
 
static void FourierSolver2D_setMirrors(FourierSolver2D& self, py::object value) {
    if (value.is_none())
        self.mirrors.reset();
    else {
        try {
            double v = py::extract<double>(value);
            self.mirrors.reset(std::make_pair(v, v));
        } catch (py::error_already_set&) {
            PyErr_Clear();
            try {
                if (py::len(value) != 2) throw py::error_already_set();
                self.mirrors.reset(
                    std::make_pair<double, double>(double(py::extract<double>(value[0])), double(py::extract<double>(value[1]))));
            } catch (py::error_already_set&) {
                throw ValueError("none, float, or tuple of two floats required");
            }
        }
    }
}
 
static py::object FourierSolver2D_getDeterminant(py::tuple args, py::dict kwargs) {
    if (py::len(args) != 1)
        throw TypeError(u8"get_determinant() takes exactly one non-keyword argument ({0} given)", py::len(args));
    FourierSolver2D* self = py::extract<FourierSolver2D*>(args[0]);
    auto* expansion = &self->expansion;
 
    enum What { WHAT_NOTHING = 0, WHAT_WAVELENGTH, WHAT_K0, WHAT_NEFF, WHAT_KTRAN, WHAT_BETA };
    What what = WHAT_NOTHING;
    py::object array;
 
    plask::optional<dcomplex> k0, neff, ktran, beta;
    std::string _ktran, _beta;
 
    AxisNames* axes = getCurrentAxes();
    py::stl_input_iterator<std::string> begin(kwargs), end;
    for (auto i = begin; i != end; ++i) {
        if (*i == "lam") {
            if (what == WHAT_K0 || k0) throw BadInput(self->getId(), u8"'lam' and 'k0' are mutually exclusive");
            if (PyArray_Check(py::object(kwargs[*i]).ptr())) {
                if (what) throw TypeError(u8"only one key may be an array");
                what = WHAT_WAVELENGTH;
                array = kwargs[*i];
            } else
                k0.reset(2e3 * PI / py::extract<dcomplex>(kwargs[*i])());
        } else if (*i == "k0") {
            if (what == WHAT_WAVELENGTH || k0) throw BadInput(self->getId(), u8"'lam' and 'k0' are mutually exclusive");
            if (PyArray_Check(py::object(kwargs[*i]).ptr())) {
                if (what) throw TypeError(u8"only one key may be an array");
                what = WHAT_K0;
                array = kwargs[*i];
            } else
                k0.reset(py::extract<dcomplex>(kwargs[*i]));
        } else if (*i == "neff") {
            if (what == WHAT_BETA || beta) throw BadInput(self->getId(), u8"'neff' and '{}' are mutually exclusive", _beta);
            if (PyArray_Check(py::object(kwargs[*i]).ptr())) {
                if (expansion->separated())
                    throw Exception("{0}: Cannot get determinant for effective index array with polarization separation",
                                    self->getId());
                if (what) throw TypeError(u8"only one key may be an array");
                what = WHAT_NEFF;
                array = kwargs[*i];
            } else {
                neff.reset(py::extract<dcomplex>(kwargs[*i]));
                if (expansion->separated() && *neff != 0.)
                    throw Exception("{0}: Effective index must be 0 with polarization separation",
                                    self->getId());
            }
        } else if (*i == "ktran" || *i == "kt" || *i == "k" + axes->getNameForTran()) {
            if (what == WHAT_KTRAN || ktran) throw BadInput(self->getId(), u8"'{}' was already specified as '{}'", *i, _ktran);
            _ktran = *i;
            if (PyArray_Check(py::object(kwargs[*i]).ptr())) {
                if (expansion->symmetric())
                    throw Exception("{0}: Cannot get determinant for transverse wavevector array with symmetry", self->getId());
                if (what) throw TypeError(u8"only one key may be an array");
                what = WHAT_KTRAN;
                array = kwargs[*i];
            } else {
                ktran.reset(py::extract<dcomplex>(kwargs[*i]));
                if (expansion->symmetric() && *ktran != 0.)
                    throw Exception("{0}: Transverse wavevector must be 0 with symmetry",
                                    self->getId());
            }
        } else if (*i == "beta" || *i == "klong" || *i == "kl" || *i == "k"+axes->getNameForLong()) {
            if (what == WHAT_BETA || beta) throw BadInput(self->getId(), u8"'{}' was already specified as '{}'", *i, _beta);
            _beta = *i;
            if (what == WHAT_NEFF || neff) throw BadInput(self->getId(), u8"'{}' and 'neff' are mutually exclusive", *i);
            if (PyArray_Check(py::object(kwargs[*i]).ptr())) {
                if (expansion->separated())
                    throw Exception("{0}: Cannot get determinant for longitudinal wavevector array with polarization separation",
                                    self->getId());
                if (what) throw TypeError(u8"only one key may be an array");
                what = WHAT_BETA;
                array = kwargs[*i];
            } else {
                beta.reset(py::extract<dcomplex>(kwargs[*i]));
                if (expansion->separated() && *beta != 0.)
                    throw Exception("{0}: Longitudinal wavevector must be 0 with polarization separation",
                                    self->getId());
            }
        } else
            throw TypeError(u8"get_determinant() got unexpected keyword argument '{0}'", *i);
    }
 
    self->Solver::initCalculation();
 
    if (k0) expansion->setK0(*k0);
    else expansion->setK0(self->getK0());
    if (neff) { if (what != WHAT_WAVELENGTH && what != WHAT_K0) expansion->setBeta(*neff * expansion->k0); }
    else if (beta) expansion->setBeta(*beta);
    else expansion->setBeta(self->getBeta());
    if (ktran) expansion->setKtran(*ktran);
    else expansion->setKtran(self->getKtran());
    expansion->setLam0(self->getLam0());
    expansion->setSymmetry(self->getSymmetry());
    expansion->setPolarization(self->getPolarization());
 
    switch (what) {
        case WHAT_NOTHING: return py::object(self->getDeterminant());
        case WHAT_WAVELENGTH:
            return UFUNC<dcomplex>(
                [self, neff](dcomplex x) -> dcomplex {
                    self->expansion.setK0(2e3 * PI / x);
                    if (neff) self->expansion.setBeta(*neff * self->expansion.k0);
                    return self->getDeterminant();
                },
                array,
                "Fourier2D.get_determinant",
                "lam");
        case WHAT_K0:
            return UFUNC<dcomplex>(
                [self, neff](dcomplex x) -> dcomplex {
                    self->expansion.setK0(x);
                    if (neff) self->expansion.setBeta(*neff * x);
                    return self->getDeterminant();
                },
                array,
                "Fourier2D.get_determinant",
                "k0");
        case WHAT_NEFF:
            return UFUNC<dcomplex>(
                [self](dcomplex x) -> dcomplex {
                    self->expansion.setBeta(x * self->getK0());
                    return self->getDeterminant();
                },
                array,
                "Fourier2D.get_determinant",
                "neff");
        case WHAT_KTRAN:
            return UFUNC<dcomplex>(
                [self](dcomplex x) -> dcomplex {
                    self->expansion.setKtran(x);
                    return self->getDeterminant();
                },
                array,
                "Fourier2D.get_determinant",
                "ktran");
        case WHAT_BETA:
            return UFUNC<dcomplex>(
                [self](dcomplex x) -> dcomplex {
                    self->expansion.setBeta(x);
                    return self->getDeterminant();
                },
                array,
                "Fourier2D.get_determinant",
                "beta");
    }
    return py::object();
}
 
static size_t FourierSolver2D_setMode(py::tuple args, py::dict kwargs) {
    if (py::len(args) != 1) throw TypeError(u8"set_mode() takes exactly one non-keyword argument ({0} given)", py::len(args));
    FourierSolver2D* self = py::extract<FourierSolver2D*>(args[0]);
    auto* expansion = &self->expansion;
 
    AxisNames* axes = getCurrentAxes();
    plask::optional<dcomplex> k0, neff, ktran;
    py::stl_input_iterator<std::string> begin(kwargs), end;
    for (auto i = begin; i != end; ++i) {
        if (*i == "lam" || *i == "wavelength") {
            if (k0) throw BadInput(self->getId(), u8"'lam' and 'k0' are mutually exclusive");
            k0.reset(2e3 * PI / py::extract<dcomplex>(kwargs[*i])());
        } else if (*i == "k0") {
            if (k0) throw BadInput(self->getId(), u8"'lam' and 'k0' are mutually exclusive");
            k0.reset(py::extract<dcomplex>(kwargs[*i]));
        } else if (*i == "neff") {
            neff.reset(py::extract<dcomplex>(kwargs[*i]));
        } else if (*i == "ktran" || *i == "kt" || *i == "k" + axes->getNameForTran()) {
            ktran.reset(py::extract<dcomplex>(kwargs[*i]));
        } else
            throw TypeError(u8"set_mode() got unexpected keyword argument '{0}'", *i);
    }
 
    self->Solver::initCalculation();
 
    if (k0)
        expansion->setK0(*k0);
    else
        expansion->setK0(self->getK0());
    if (neff)
        expansion->setBeta(*neff * expansion->k0);
    else
        expansion->setBeta(self->getBeta());
    if (ktran)
        expansion->setKtran(*ktran);
    else
        expansion->setKtran(self->getKtran());
    expansion->setLam0(self->getLam0());
    expansion->setSymmetry(self->getSymmetry());
    expansion->setPolarization(self->getPolarization());
 
    return self->setMode();
}
 
static size_t FourierSolver2D_findMode(py::tuple args, py::dict kwargs) {
    if (py::len(args) != 1) throw TypeError(u8"find_mode() takes exactly one non-keyword argument ({0} given)", py::len(args));
    FourierSolver2D* self = py::extract<FourierSolver2D*>(args[0]);
 
    if (py::len(kwargs) != 1) throw TypeError(u8"find_mode() takes exactly one keyword argument ({0} given)", py::len(kwargs));
    std::string key = py::extract<std::string>(kwargs.keys()[0]);
    dcomplex value = py::extract<dcomplex>(kwargs[key]);
    AxisNames* axes = getCurrentAxes();
    FourierSolver2D::What what;
 
    if (key == "lam" || key == "wavelength")
        what = FourierSolver2D::WHAT_WAVELENGTH;
    else if (key == "k0")
        what = FourierSolver2D::WHAT_K0;
    else if (key == "neff")
        what = FourierSolver2D::WHAT_NEFF;
    else if (key == "ktran" || key == "kt" || key == "k" + axes->getNameForTran())
        what = FourierSolver2D::WHAT_KTRAN;
    else if (key == "beta" || key == "klong" || key == "kl" || key == "k" + axes->getNameForTran())
        what = FourierSolver2D::WHAT_BETA;
    else
        throw TypeError(u8"find_mode() got unexpected keyword argument '{0}'", key);
 
    return self->findMode(what, value);
}
 
static dcomplex FourierSolver2D_Mode_Neff(const FourierSolver2D::Mode& mode) { return mode.beta / mode.k0; }
 
static py::object FourierSolver2D_Mode__getattr__(const FourierSolver2D::Mode& mode, const std::string name) {
    auto axes = getCurrentAxes();
    if (name == "k" + axes->getNameForLong()) return py::object(mode.beta);
    if (name == "k" + axes->getNameForTran()) return py::object(mode.ktran);
    throw AttributeError(u8"'Mode' object has no attribute '{0}'", name);
    return py::object();
}
 
static std::string FourierSolver2D_Mode_str(const FourierSolver2D::Mode& self) {
    AxisNames* axes = getCurrentAxes();
    std::string pol;
    switch (self.polarization) {
        case Expansion::E_TRAN: pol = "E" + axes->getNameForTran(); break;
        case Expansion::E_LONG: pol = "E" + axes->getNameForLong(); break;
        default: pol = "none";
    }
    std::string sym;
    switch (self.symmetry) {
        case Expansion::E_TRAN: sym = "E" + axes->getNameForTran(); break;
        case Expansion::E_LONG: sym = "E" + axes->getNameForLong(); break;
        default: sym = "none";
    }
    return format(u8"<lam: {:.2f}nm, neff: {}, ktran: {}/um, polarization: {}, symmetry: {}, power: {:.2g} mW>",
                  real(2e3 * PI / self.k0), str(self.beta / self.k0, "{:.3f}{:+.3g}j"),
                  str(self.ktran, "({:.3g}{:+.3g}j)", "{:.3g}"), pol, sym, self.power);
}
static std::string FourierSolver2D_Mode_repr(const FourierSolver2D::Mode& self) {
    AxisNames* axes = getCurrentAxes();
    std::string pol;
    switch (self.polarization) {
        case Expansion::E_TRAN: pol = "'E" + axes->getNameForTran() + "'"; break;
        case Expansion::E_LONG: pol = "'E" + axes->getNameForLong() + "'"; break;
        default: pol = "None";
    }
    std::string sym;
    switch (self.symmetry) {
        case Expansion::E_TRAN: sym = "'E" + axes->getNameForTran() + "'"; break;
        case Expansion::E_LONG: sym = "'E" + axes->getNameForLong() + "'"; break;
        default: sym = "None";
    }
    return format(u8"Fourier2D.Mode(lam={0}, neff={1}, ktran={2}, polarization={3}, symmetry={4}, power={5})",
                  str(2e3 * PI / self.k0), str(self.beta / self.k0), str(self.ktran), pol, sym, self.power);
}
 
static py::object FourierSolver2D_getFieldVectorE(FourierSolver2D& self, int num, double z) {
    if (num < 0) num += int(self.modes.size());
    if (std::size_t(num) >= self.modes.size()) throw IndexError(u8"bad mode number {:d}", num);
    return arrayFromVec2D<NPY_CDOUBLE>(self.getFieldVectorE(num, z), self.separated(), 2);
}
 
static py::object FourierSolver2D_getFieldVectorH(FourierSolver2D& self, int num, double z) {
    if (num < 0) num += int(self.modes.size());
    if (std::size_t(num) >= self.modes.size()) throw IndexError(u8"bad mode number {:d}", num);
    return arrayFromVec2D<NPY_CDOUBLE>(self.getFieldVectorH(num, z), self.separated(), 2);
}
 
template <>
py::object Scattering<FourierSolver2D>::getFieldVectorE(double z) {
    return arrayFromVec2D<NPY_CDOUBLE>(solver->getScatteredFieldVectorE(incident, side, z), solver->separated(), 2);
}
 
template <>
py::object Scattering<FourierSolver2D>::getFieldVectorH(double z) {
    return arrayFromVec2D<NPY_CDOUBLE>(solver->getScatteredFieldVectorH(incident, side, z), solver->separated(), 2);
}
 
 
static py::object FourierSolver2D_incidentGaussian(FourierSolver2D& self,
                                                   Transfer::IncidentDirection side,
                                                   Expansion::Component polarization,
                                                   double sigma,
                                                   double center) {
    return arrayFromVec<NPY_CDOUBLE>(self.incidentGaussian(side, polarization, sigma, center));
}
 
static shared_ptr<Scattering<FourierSolver2D>> FourierSolver2D_scatteringGaussian(FourierSolver2D& self,
                                                                                  Transfer::IncidentDirection side,
                                                                                  Expansion::Component polarization,
                                                                                  double sigma,
                                                                                  double center) {
    return shared_ptr<Scattering<FourierSolver2D>>(
        new Scattering<FourierSolver2D>(&self, side, self.incidentGaussian(side, polarization, sigma, center)));
}
 
namespace detail {
    static Expansion::Component getPolarization2D(PyObject* obj) {
        Expansion::Component val;
        if (obj == Py_None) {
            val = Expansion::E_UNSPECIFIED;
        } else {
            AxisNames* axes = getCurrentAxes();
            try {
                std::string repr = py::extract<std::string>(obj);
                if (repr == "none" || repr == "NONE" || repr == "None")
                    val = Expansion::E_UNSPECIFIED;
                else if (repr == "Etran" || repr == "Et" || repr == "E"+axes->getNameForTran() ||
                        repr == "Hlong" || repr == "Hl" || repr == "H"+axes->getNameForLong() || repr == "TM")
                    val = Expansion::E_TRAN;
                else if (repr == "Elong" || repr == "El" || repr == "E"+axes->getNameForLong() ||
                        repr == "Htran" || repr == "Ht" || repr == "H"+axes->getNameForTran() || repr == "TE")
                    val = Expansion::E_LONG;
                else
                    throw py::error_already_set();
            } catch (py::error_already_set&) {
                throw ValueError("wrong component specification.");
            }
        }
        return val;
    }
}
 
void FourierSolver2D_setPolarization(FourierSolver2D* self, PyObject* polarization) {
    self->setPolarization(detail::getPolarization2D(polarization));
}
 
py::object FourierSolver2D_computeReflectivity_polarization(FourierSolver2D* self,
                                                            py::object wavelength,
                                                            Transfer::IncidentDirection side,
                                                            PyObject* polarization
                                                           ) {
    return Solver_computeReflectivity_polarization(self, wavelength, side, detail::getPolarization2D(polarization));
}
 
py::object FourierSolver2D_computeTransmittivity_polarization(FourierSolver2D* self,
                                                              py::object wavelength,
                                                              Transfer::IncidentDirection side,
                                                              PyObject* polarization
                                                             ) {
    return Solver_computeTransmittivity_polarization(self, wavelength, side, detail::getPolarization2D(polarization));
}
 
shared_ptr<Scattering<FourierSolver2D>> FourierSolver2D_Scattering_from_polarization(FourierSolver2D* parent,
                                                                                     Transfer::IncidentDirection side,
                                                                                     PyObject* polarization
                                                                                    ) {
    return Scattering<FourierSolver2D>::from_polarization(parent, side, detail::getPolarization2D(polarization));
}
 
 
void export_FourierSolver2D() {
    py_enum<typename FourierSolver2D::FourierType>()
        .value("DISCRETE", FourierSolver2D::FOURIER_DISCRETE)
        .value("ANALYTIC", FourierSolver2D::FOURIER_ANALYTIC);
 
    CLASS(FourierSolver2D, "Fourier2D",
          u8"Optical Solver using Fourier expansion in 2D.\n\n"
          u8"It calculates optical modes and optical field distribution using Fourier modal method\n"
          u8"and reflection transfer in two-dimensional Cartesian space.")
    export_base(solver);
    PROVIDER(outNeff, "Effective index of the last computed mode.");
    RW_PROPERTY(size, getSize, setSize, "Orthogonal expansion size.");
    RW_PROPERTY(symmetry, getSymmetry, setSymmetry, "Mode symmetry.");
    solver.add_property("polarization", &__Class__::getPolarization, &FourierSolver2D_setPolarization, "Mode polarization.");
    solver.add_property("lam", &__Class__::getLam, &Solver_setLam<__Class__>,
                        u8"Wavelength of the light (nm).\n\n"
                        u8"Use this property only if you are looking for anything else than\n"
                        u8"the wavelength, e.g. the effective index of lateral wavevector.\n");
    solver.add_property("wavelength", &__Class__::getLam, &Solver_setLam<__Class__>, u8"Alias for :attr:`lam`");
    solver.add_property("k0", &__Class__::getK0, &Solver_setK0<__Class__>,
                        u8"Normalized frequency of the light (1/µm).\n\n"
                        u8"Use this property only if you are looking for anything else than\n"
                        u8"the wavelength,e.g. the effective index of lateral wavevector.\n");
    RW_PROPERTY(klong, getBeta, setBeta,
                u8"Longitudinal propagation constant of the light (1/µm).\n\n"
                u8"Use this property only if you are looking for anything else than\n"
                u8"the longitudinal component of the propagation vector and the effective index.\n");
    RW_PROPERTY(beta, getBeta, setBeta,
                u8"Alias for :attr:`klong`\n");
    RW_PROPERTY(ktran, getKtran, setKtran,
                u8"Transverse propagation constant of the light (1/µm).\n\n"
                u8"Use this property only  if you are looking for anything else than\n"
                u8"the transverse component of the propagation vector.\n");
    RW_FIELD(refine, "Number of refinement points for refractive index averaging.");
    solver.add_property("ft", &__Class__::getFourierType, &__Class__::setFourierType,
                        u8"Type of the Fourier transform. Analytic transform is faster and more precise,\n"
                        u8"however it ignores temperature and gain distributions.\n");
    solver.add_property("dct", &__Class__::getDCT, &__Class__::setDCT,
                        "Type of discrete cosine transform for symmetric expansion.");
    RW_FIELD(emission,
             "Direction of the useful light emission.\n\n"
             u8"Necessary for the over-threshold model to correctly compute the output power.\n"
             u8"Currently the fields are normalized only if this parameter is set to\n"
             u8"``top`` or ``bottom``. Otherwise, it is ``undefined`` (default) and the fields\n"
             u8"are not normalized.");
    solver.def("get_determinant", py::raw_function(FourierSolver2D_getDeterminant),
               u8"Compute discontinuity matrix determinant.\n\n"
               u8"Arguments can be given through keywords only.\n\n"
               u8"Args:\n"
               u8"    lam (complex): Wavelength (nm).\n"
               u8"    k0 (complex): Normalized frequency (1/µm).\n"
               u8"    neff (complex): Longitudinal effective index.\n"
               u8"    ktran (complex): Transverse wavevector (1/µm).\n");
    solver.def("find_mode", py::raw_function(FourierSolver2D_findMode),
               u8"Compute the mode near the specified effective index.\n\n"
               u8"Only one of the following arguments can be given through a keyword.\n"
               u8"It is the starting point for search of the specified parameter.\n\n"
               u8"Args:\n"
               u8"    lam (complex): Wavelength (nm).\n"
               u8"    k0 (complex): Normalized frequency (1/µm).\n"
               u8"    neff (complex): Longitudinal effective index.\n"
               u8"    ktran (complex): Transverse wavevector (1/µm).\n");
    solver.def("set_mode", py::raw_function(FourierSolver2D_setMode),
               u8"Set the mode for specified parameters.\n\n"
               u8"This method should be used if you have found a mode manually and want to insert\n"
               u8"it into the solver in order to determine the fields. Calling this will raise an\n"
               u8"exception if the determinant for the specified parameters is too large.\n\n"
               u8"Arguments can be given through keywords only.\n\n"
               u8"Args:\n"
               u8"    lam (complex): Wavelength (nm).\n"
               u8"    k0 (complex): Normalized frequency (1/µm).\n"
               u8"    neff (complex): Longitudinal effective index.\n"
               u8"    ktran (complex): Transverse wavevector (1/µm).\n");
    solver.def("compute_reflectivity", &FourierSolver2D_computeReflectivity_polarization,
               (py::arg("lam"), "side", "polarization"));
    solver.def("compute_reflectivity", &Solver_computeReflectivity_index<FourierSolver2D>, (py::arg("lam"), "side", "index"));
    solver.def("compute_reflectivity", &Solver_computeReflectivity_array<FourierSolver2D>, (py::arg("lam"), "side", "coffs"),
               u8"Compute reflection coefficient on planar incidence [%].\n\n"
               u8"Args:\n"
               u8"    lam (float or array of floats): Incident light wavelength (nm).\n"
               u8"    side (`top` or `bottom`): Side of the structure where the incident light is\n"
               u8"        present.\n"
               u8"    polarization: Specification of the incident light polarization.\n"
               u8"        It should be a string of the form 'E\\ *#*\\ ', where *#* is the axis\n"
               u8"        name of the non-vanishing electric field component.\n"
               u8"    idx: Eigenmode number.\n"
               u8"    coeffs: expansion coefficients of the incident vector.\n");
    solver.def("compute_transmittivity", &FourierSolver2D_computeTransmittivity_polarization,
               (py::arg("lam"), "side", "polarization"));
    solver.def("compute_transmittivity", &Solver_computeTransmittivity_index<FourierSolver2D>, (py::arg("lam"), "side", "index"));
    solver.def("compute_transmittivity", &Solver_computeTransmittivity_array<FourierSolver2D>, (py::arg("lam"), "side", "coffs"),
               u8"Compute transmission coefficient on planar incidence [%].\n\n"
               u8"Args:\n"
               u8"    lam (float or array of floats): Incident light wavelength (nm).\n"
               u8"    side (`top` or `bottom`): Side of the structure where the incident light is\n"
               u8"        present.\n"
               u8"    polarization: Specification of the incident light polarization.\n"
               u8"        It should be a string of the form 'E\\ *#*\\ ', where *#* is the axis name\n"
               u8"        of the non-vanishing electric field component.\n"
               u8"    idx: Eigenmode number.\n"
               u8"    coeffs: expansion coefficients of the incident vector.\n");
    solver.add_property("mirrors", FourierSolver2D_getMirrors, FourierSolver2D_setMirrors,
                        u8"Mirror reflectivities. If None then they are automatically estimated from the\n"
                        u8"Fresnel equations.");
    solver.add_property("pml", py::make_function(&Solver_getPML<FourierSolver2D>, py::with_custodian_and_ward_postcall<0, 1>()),
                        &Solver_setPML<FourierSolver2D>, u8"Side Perfectly Matched Layers boundary conditions.\n\n" PML_ATTRS_DOC);
    RO_FIELD(modes, "Computed modes.");
    solver.def("scattering", FourierSolver2D_Scattering_from_polarization, py::with_custodian_and_ward_postcall<0, 1>(),
               (py::arg("side"), "polarization"));
    solver.def("scattering", Scattering<FourierSolver2D>::from_index, py::with_custodian_and_ward_postcall<0, 1>(),
               (py::arg("side"), "idx"));
    solver.def("scattering", Scattering<FourierSolver2D>::from_array, py::with_custodian_and_ward_postcall<0, 1>(),
               (py::arg("side"), "coeffs"),
               u8"Access to the reflected field.\n\n"
               u8"Args:\n"
               u8"    side (`top` or `bottom`): Side of the structure where the incident light is\n"
               u8"        present.\n"
               u8"    polarization: Specification of the incident light polarization.\n"
               u8"        It should be a string of the form 'E\\ *#*\\ ', where *#* is the axis name\n"
               u8"        of the non-vanishing electric field component.\n"
               u8"    idx: Eigenmode number.\n"
               u8"    coeffs: expansion coefficients of the incident vector.\n\n"
               u8":rtype: Fourier2D.Scattering\n");
    solver.def("get_raw_E", FourierSolver2D_getFieldVectorE, (py::arg("num"), "level"),
               u8"Get Fourier expansion coefficients for the electric field.\n\n"
               u8"This is a low-level function returning $E_l$ and/or $E_t$ Fourier\n"
               u8"expansion coefficients. Please refer to the detailed solver description for their\n"
               u8"interpretation.\n\n"
               u8"Args:\n"
               u8"    num (int): Computed mode number.\n"
               u8"    level (float): Vertical level at which the coefficients are computed.\n\n"
               u8":rtype: numpy.ndarray\n");
    solver.def("get_raw_H", FourierSolver2D_getFieldVectorH, (py::arg("num"), "level"),
               u8"Get Fourier expansion coefficients for the magnetic field.\n\n"
               u8"This is a low-level function returning $H_l$ and/or $H_t$ Fourier\n"
               u8"expansion coefficients. Please refer to the detailed solver description for their\n"
               u8"interpretation.\n\n"
               u8"Args:\n"
               u8"    num (int): Computed mode number.\n"
               u8"    level (float): Vertical level at which the coefficients are computed.\n\n"
               u8":rtype: numpy.ndarray\n");
    solver.def("layer_eigenmodes", &Eigenmodes<FourierSolver2D>::fromZ, py::arg("level"),
               u8"Get eignemodes for a layer at specified level.\n\n"
               u8"This is a low-level function to access diagonalized eigenmodes for a specific\n"
               u8"layer. Please refer to the detailed solver description for the interpretation\n"
               u8"of the returned values.\n\n"
               u8"Args:\n"
               u8"    level (float): Vertical level at which the coefficients are computed.\n\n"
               u8":rtype: :class:`~optical.modal.Fourier2D.Eigenmodes`\n",
               py::with_custodian_and_ward_postcall<0, 1>());
    solver.def("gaussian", &FourierSolver2D_incidentGaussian, (py::arg("side"), "polarization", "sigma", py::arg("center") = 0.),
               u8"Create coefficients vector with Gaussian profile.\n\n"
               u8"This method is intended to use for :py:meth:`scattering` method.\n\n"
               u8"Args:\n"
               u8"    side (`top` or `bottom`): Side of the structure where the incident light is\n"
               u8"        present.\n"
               u8"    polarization: Specification of the incident light polarization.\n"
               u8"        It should be a string of the form 'E\\ *#*\\ ', where *#* is the axis name\n"
               u8"        of the non-vanishing electric field component.\n"
               u8"    sigma (float): Gaussian standard deviation (µm).\n"
               u8"    center (float): Position of the beam center (µm).\n\n"
               u8"Example:\n"
               u8"   >>> scattered = fourier.scattering('top', \n"
               u8"   ...     fourier.gaussian('top', 'Ex', 0.2))\n");
    solver.def("scattering_gaussian", &FourierSolver2D_scatteringGaussian,
               (py::arg("side"), "polarization", "sigma", py::arg("center") = 0.),
               u8"Helper function to Access reflected fields for access incidence.\n\n"
               u8"This method is equivalent to calling:\n\n"
               u8"   >>> fourier.scattering(side,\n"
               u8"   ...     fourier.gaussian(side, polarization, sigma, center))\n\n"
               u8"Args:\n"
               u8"    side (`top` or `bottom`): Side of the structure where the incident light is\n"
               u8"        present.\n"
               u8"    polarization: Specification of the incident light polarization.\n"
               u8"        It should be a string of the form 'E\\ *#*\\ ', where *#* is the axis name\n"
               u8"        of the non-vanishing electric field component.\n"
               u8"    sigma (float): Gaussian standard deviation (µm).\n"
               u8"    center (float): Position of the beam center (µm).\n\n");
 
    py::scope scope = solver;
    (void)scope;  // don't warn about unused variable scope
 
    register_vector_of<FourierSolver2D::Mode>("Modes");
    py::class_<FourierSolver2D::Mode>("Mode", u8"Detailed information about the mode.", py::no_init)
        .def_readonly("symmetry", &FourierSolver2D::Mode::symmetry, u8"Mode horizontal symmetry.")
        .def_readonly("polarization", &FourierSolver2D::Mode::polarization, u8"Mode polarization.")
        .add_property("lam", &getModeWavelength<FourierSolver2D::Mode>, u8"Mode wavelength (nm).")
        .add_property("wavelength", &getModeWavelength<FourierSolver2D::Mode>, u8"Mode wavelength (nm).")
        .def_readonly("k0", &FourierSolver2D::Mode::k0, u8"Mode normalized frequency (1/µm).")
        .def_readonly("beta", &FourierSolver2D::Mode::beta, u8"Mode longitudinal wavevector (1/µm).")
        .add_property("neff", &FourierSolver2D_Mode_Neff, u8"Mode longitudinal effective index (-).")
        .def_readonly("ktran", &FourierSolver2D::Mode::ktran, u8"Mode transverse wavevector (1/µm).")
        .def_readwrite("power", &FourierSolver2D::Mode::power, u8"Total power emitted into the mode (mW).")
        .def("__str__", &FourierSolver2D_Mode_str)
        .def("__repr__", &FourierSolver2D_Mode_repr)
        .def("__getattr__", &FourierSolver2D_Mode__getattr__);
 
    Scattering<FourierSolver2D>::registerClass("2D", "Fourier");
    Eigenmodes<FourierSolver2D>::registerClass("Fourier2D", "2D");
}
 
}}}}  // namespace plask::optical::modal::python