eim.hpp

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.
 */
#ifndef PLASK__MODULE_OPTICAL_EIM_HPP
#define PLASK__MODULE_OPTICAL_EIM_HPP
 
#include <limits>
 
#include <plask/plask.hpp>
 
#include "bisection.hpp"
#include "rootdigger.hpp"
 
namespace plask { namespace optical { namespace effective {
 
struct PLASK_SOLVER_API EffectiveIndex2D : public SolverWithMesh<Geometry2DCartesian, RectangularMesh<2>> {
    enum Symmetry { SYMMETRY_DEFAULT, SYMMETRY_POSITIVE, SYMMETRY_NEGATIVE, SYMMETRY_NONE };
 
    enum Polarization {
        TE,
        TM,
    };
 
    enum Emission { FRONT, BACK };
 
    struct Field {
        dcomplex F, B;
        Field() = default;
        Field(dcomplex f, dcomplex b) : F(f), B(b) {}
        Field operator*(dcomplex a) const { return Field(F * a, B * a); }
        Field operator/(dcomplex a) const { return Field(F / a, B / a); }
        Field operator*=(dcomplex a) {
            F *= a;
            B *= a;
            return *this;
        }
        Field operator/=(dcomplex a) {
            F /= a;
            B /= a;
            return *this;
        }
    };
 
    struct Matrix {
        dcomplex ff, fb, bf, bb;
        Matrix() = default;
        Matrix(dcomplex t1, dcomplex t2, dcomplex t3, dcomplex t4) : ff(t1), fb(t2), bf(t3), bb(t4) {}
        static Matrix eye() { return Matrix(1., 0., 0., 1.); }
        static Matrix diag(dcomplex f, dcomplex b) { return Matrix(f, 0., 0., b); }
        Matrix operator*(const Matrix& T) {
            return Matrix(ff * T.ff + fb * T.bf, ff * T.fb + fb * T.bb, bf * T.ff + bb * T.bf, bf * T.fb + bb * T.bb);
        }
        Field solve(const Field& v) { return Field(bb * v.F - fb * v.B, -bf * v.F + ff * v.B) / (ff * bb - fb * bf); }
    };
 
    struct Mode {
        EffectiveIndex2D* solver;                                 
        Symmetry symmetry;                                        
        dcomplex neff;                                            
        bool have_fields;                                         
        std::vector<Field, aligned_allocator<Field>> xfields;     
        std::vector<double, aligned_allocator<double>> xweights;  
        double power;                                             
 
        Mode(EffectiveIndex2D* solver, Symmetry sym)
            : solver(solver), have_fields(false), xfields(solver->xend), xweights(solver->xend), power(1.) {
            setSymmetry(sym);
        }
 
        void setSymmetry(Symmetry sym) {
            if (solver->geometry->isSymmetric(Geometry::DIRECTION_TRAN)) {
                if (sym == SYMMETRY_DEFAULT)
                    sym = SYMMETRY_POSITIVE;
                else if (sym == SYMMETRY_NONE)
                    throw BadInput(solver->getId(), "for symmetric geometry specify positive or negative symmetry");
            } else {
                if (sym == SYMMETRY_DEFAULT)
                    sym = SYMMETRY_NONE;
                else if (sym != SYMMETRY_NONE)
                    throw BadInput(solver->getId(), "for non-symmetric geometry no symmetry may be specified");
            }
            symmetry = sym;
        }
 
        bool operator==(const Mode& other) const { return symmetry == other.symmetry && is_zero(neff - other.neff); }
 
        double loss() const { return -2e7 * imag(neff * solver->k0); }
    };
 
  protected:
    friend struct RootDigger;
 
    size_t xbegin,  
        xend,       
        ybegin,     
        yend;       
 
    DataLog<dcomplex, dcomplex> log_value;
 
    std::vector<std::vector<dcomplex, aligned_allocator<dcomplex>>> nrCache;
 
    std::vector<Field, aligned_allocator<Field>> yfields;
 
    std::vector<double, aligned_allocator<double>> yweights;
 
    std::vector<dcomplex, aligned_allocator<dcomplex>> epsilons;
 
    double stripex;  
 
    Polarization polarization;  
 
    bool recompute_neffs;  
 
  public:
    Emission emission;  
 
    dcomplex vneff;  
 
    plask::optional<std::pair<double, double>> mirrors;
 
    RootDigger::Params root;
 
    RootDigger::Params stripe_root;
 
    std::vector<Mode> modes;
 
    ReceiverFor<Temperature, Geometry2DCartesian> inTemperature;
 
    ReceiverFor<Gain, Geometry2DCartesian> inGain;
 
    ReceiverFor<CarriersConcentration, Geometry2DCartesian> inCarriersConcentration;
 
    typename ProviderFor<ModeEffectiveIndex>::Delegate outNeff;
 
    typename ProviderFor<ModeLightMagnitude, Geometry2DCartesian>::Delegate outLightMagnitude;
 
    typename ProviderFor<ModeLightE, Geometry2DCartesian>::Delegate outLightE;
 
    typename ProviderFor<RefractiveIndex, Geometry2DCartesian>::Delegate outRefractiveIndex;
 
    typename ProviderFor<Heat, Geometry2DCartesian>::Delegate outHeat;
 
    EffectiveIndex2D(const std::string& name = "");
 
    virtual ~EffectiveIndex2D() {
        inTemperature.changedDisconnectMethod(this, &EffectiveIndex2D::onInputChange);
        inGain.changedDisconnectMethod(this, &EffectiveIndex2D::onInputChange);
        inCarriersConcentration.changedDisconnectMethod(this, &EffectiveIndex2D::onInputChange);
    }
 
    std::string getClassName() const override { return "optical.EffectiveIndex2D"; }
 
    std::string getClassDescription() const override {
        return "Calculate optical modes and optical field distribution using the effective index method "
               "in Cartesian two-dimensional space.";
    }
 
    void loadConfiguration(plask::XMLReader& reader, plask::Manager& manager) override;
 
    double getStripeX() const { return stripex; }
 
    void setStripeX(double x) {
        stripex = x;
        invalidate();
    }
 
    Polarization getPolarization() const { return polarization; }
 
    void setPolarization(Polarization polar) {
        polarization = polar;
        invalidate();
    }
 
    dcomplex getWavelength() const { return 2e3 * PI / k0; }
 
    void setWavelength(dcomplex wavelength) {
        k0 = 2e3 * PI / wavelength;
        invalidate();
    }
 
    void setSimpleMesh() {
        writelog(LOG_DETAIL, "Creating simple mesh");
        setMesh(plask::make_shared<RectangularMesh2DSimpleGenerator>());
    }
 
    void setHorizontalMesh(shared_ptr<MeshAxis> meshx) {  // TODO pointer to mesh is held now, is this fine?
        writelog(LOG_DETAIL, "Setting horizontal mesh");
        if (!geometry) throw NoChildException();
        auto meshxy = RectangularMesh2DSimpleGenerator().generate_t<RectangularMesh<2>>(geometry->getChild());
        meshxy->setAxis0(meshx);
        setMesh(meshxy);
    }
 
    std::vector<dcomplex> searchVNeffs(plask::dcomplex neff1 = 0.,
                                       plask::dcomplex neff2 = 0.,
                                       size_t resteps = 256,
                                       size_t imsteps = 64,
                                       dcomplex eps = dcomplex(1e-6, 1e-9));
 
    size_t findMode(dcomplex neff, Symmetry symmetry = SYMMETRY_DEFAULT);
 
    std::vector<size_t> findModes(dcomplex neff1 = 0.,
                                  dcomplex neff2 = 0.,
                                  Symmetry symmetry = SYMMETRY_DEFAULT,
                                  size_t resteps = 256,
                                  size_t imsteps = 64,
                                  dcomplex eps = dcomplex(1e-6, 1e-9));
 
    dcomplex getVertDeterminant(dcomplex neff) {
        updateCache();
        size_t stripe = mesh->tran()->findIndex(stripex);
        if (stripe < xbegin)
            stripe = xbegin;
        else if (stripe >= xend)
            stripe = xend - 1;
        return detS1(neff, nrCache[stripe]);
    }
 
    dcomplex getDeterminant(dcomplex neff, Symmetry sym = SYMMETRY_DEFAULT) {
        stageOne();
        Mode mode(this, sym);
        dcomplex det = detS(neff, mode);
        return det;
    }
 
    size_t setMode(dcomplex neff, Symmetry sym = SYMMETRY_DEFAULT);
 
    void clearModes() { modes.clear(); }
 
    double getTotalAbsorption(Mode& mode);
 
    double getTotalAbsorption(size_t num);
 
    dcomplex getDeltaNeff(double x) {
        stageOne();
        size_t i(clamp(mesh->tran()->findIndex(geometry->isSymmetric(Geometry::DIRECTION_TRAN) ? abs(x) : x), xbegin, xend - 1));
        return sqrt(epsilons[i]);
    }
 
  protected:
    void onInputChange(ReceiverBase&, ReceiverBase::ChangeReason) { invalidate(); }
 
    void onInitialize() override;
 
    void onInvalidate() override;
 
    dcomplex k0;
 
    bool need_gain;
 
    double getMirrorLosses(dcomplex n) {
        double L = geometry->getExtrusion()->getLength();
        if (isinf(L)) return 0.;
        const double lambda = real(2e3 * PI / k0);
        double R1, R2;
        if (mirrors) {
            std::tie(R1, R2) = *mirrors;
        } else {
            const double n1 = real(geometry->getFrontMaterial()->Nr(lambda, 300.)),
                         n2 = real(geometry->getBackMaterial()->Nr(lambda, 300.));
            R1 = abs((n - n1) / (n + n1));
            R2 = abs((n - n2) / (n + n2));
        }
        return lambda * std::log(R1 * R2) / (4e3 * PI * L);
    }
 
    void updateCache();
 
    void stageOne();
 
    void computeWeights(size_t stripe);
 
    void normalizeFields(Mode& mode, const std::vector<dcomplex, aligned_allocator<dcomplex>>& kx);
 
    dcomplex detS1(const dcomplex& x, const std::vector<dcomplex, aligned_allocator<dcomplex>>& NR, bool save = false);
 
    dcomplex detS(const dcomplex& x, Mode& mode, bool save = false);
 
    size_t insertMode(const Mode& mode) {
        for (size_t i = 0; i != modes.size(); ++i)
            if (modes[i] == mode) return i;
        modes.push_back(mode);
        outNeff.fireChanged();
        outLightMagnitude.fireChanged();
        outLightE.fireChanged();
        return modes.size() - 1;
    }
 
    size_t nmodes() const { return modes.size(); }
 
    dcomplex getEffectiveIndex(size_t n) {
        if (n >= modes.size()) throw NoValue(ModeEffectiveIndex::NAME);
        return modes[n].neff;
    }
 
    template <typename T> struct FieldDataBase;
    template <typename T> struct FieldDataInefficient;
    template <typename T> struct FieldDataEfficient;
    struct HeatDataImpl;
 
    const LazyData<double> getLightMagnitude(std::size_t num,
                                             shared_ptr<const plask::MeshD<2>> dst_mesh,
                                             plask::InterpolationMethod = INTERPOLATION_DEFAULT);
 
    const LazyData<Vec<3, dcomplex>> getElectricField(std::size_t num,
                                                      shared_ptr<const plask::MeshD<2>> dst_mesh,
                                                      plask::InterpolationMethod = INTERPOLATION_DEFAULT);
 
    const LazyData<dcomplex> getRefractiveIndex(RefractiveIndex::EnumType component,
                                                shared_ptr<const MeshD<2>> dst_mesh,
                                                dcomplex lam,
                                                InterpolationMethod = INTERPOLATION_DEFAULT);
 
    const LazyData<double> getHeat(shared_ptr<const MeshD<2>> dst_mesh, InterpolationMethod method = INTERPOLATION_DEFAULT);
};
 
}}}  // namespace plask::optical::effective
 
#endif  // PLASK__MODULE_OPTICAL_EIM_HPP