efm.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_EFM_HPP
#define PLASK__MODULE_OPTICAL_EFM_HPP
 
#include <limits>
 
#include <camos/camos.h>
#include <plask/plask.hpp>
 
#include "bisection.hpp"
#include "rootdigger.hpp"
 
namespace plask { namespace optical { namespace effective {
 
static constexpr int MH = 2;  // Hankel function type (1 or 2)
 
struct PLASK_SOLVER_API EffectiveFrequencyCyl : public SolverWithMesh<Geometry2DCylindrical, RectangularMesh<2>> {
    struct FieldZ {
        dcomplex F, B;
        FieldZ() = default;
        FieldZ(dcomplex f, dcomplex b) : F(f), B(b) {}
        FieldZ operator*(dcomplex a) const { return FieldZ(F * a, B * a); }
        FieldZ operator/(dcomplex a) const { return FieldZ(F / a, B / a); }
        FieldZ operator*=(dcomplex a) {
            F *= a;
            B *= a;
            return *this;
        }
        FieldZ operator/=(dcomplex a) {
            F /= a;
            B /= a;
            return *this;
        }
    };
 
    struct MatrixZ {
        dcomplex ff, fb, bf, bb;
        MatrixZ() = default;
        MatrixZ(dcomplex t1, dcomplex t2, dcomplex t3, dcomplex t4) : ff(t1), fb(t2), bf(t3), bb(t4) {}
        static MatrixZ eye() { return MatrixZ(1., 0., 0., 1.); }
        static MatrixZ diag(dcomplex f, dcomplex b) { return MatrixZ(f, 0., 0., b); }
        MatrixZ operator*(const MatrixZ& T) {
            return MatrixZ(ff * T.ff + fb * T.bf, ff * T.fb + fb * T.bb, bf * T.ff + bb * T.bf, bf * T.fb + bb * T.bb);
        }
        FieldZ operator*(const FieldZ& v) { return FieldZ(ff * v.F + fb * v.B, bf * v.F + bb * v.B); }
        FieldZ solve(const FieldZ& v) { return FieldZ(bb * v.F - fb * v.B, -bf * v.F + ff * v.B) / (ff * bb - fb * bf); }
    };
 
    struct FieldR {
        dcomplex J, H;
        FieldR() = default;
        FieldR(dcomplex j, dcomplex h) : J(j), H(h) {}
        FieldR operator*(dcomplex a) const { return FieldR(a * J, a * H); }
        FieldR operator/(dcomplex a) const { return FieldR(J / a, H / a); }
        FieldR& operator*=(dcomplex a) {
            J *= a;
            H *= a;
            return *this;
        }
        FieldR& operator/=(dcomplex a) {
            J /= a;
            H /= a;
            return *this;
        }
    };
 
    struct MatrixR {
        dcomplex JJ, JH, HJ, HH;
        MatrixR(dcomplex jj, dcomplex jh, dcomplex hj, dcomplex hh) : JJ(jj), JH(jh), HJ(hj), HH(hh) {}
        static MatrixR eye() { return MatrixR(1., 0., 0., 1.); }
        static MatrixR diag(dcomplex j, dcomplex h) { return MatrixR(j, 0., 0., h); }
        MatrixR operator*(dcomplex c) { return MatrixR(c * JJ, c * JH, c * HJ, c * HH); }
        friend MatrixR operator*(dcomplex c, const MatrixR& M) { return MatrixR(c * M.JJ, c * M.JH, c * M.HJ, c * M.HH); }
        MatrixR operator/(dcomplex d) {
            dcomplex c = 1. / d;
            return MatrixR(c * JJ, c * JH, c * HJ, c * HH);
        }
        MatrixR& operator*=(dcomplex c) {
            JJ *= c;
            JH *= c;
            HJ *= c;
            HH *= c;
            return *this;
        }
        MatrixR& operator/=(dcomplex d) {
            dcomplex c = 1. / d;
            JJ *= c;
            JH *= c;
            HJ *= c;
            HH *= c;
            return *this;
        }
        FieldR operator*(const FieldR& v) { return FieldR(JJ * v.J + JH * v.H, HJ * v.J + HH * v.H); }
        MatrixR operator*(const MatrixR& o) {
            return MatrixR(JJ * o.JJ + JH * o.HJ, JJ * o.JH + JH * o.HH, HJ * o.JJ + HH * o.HJ, HJ * o.JH + HH * o.HH);
        }
        FieldR solve(const FieldR& v) { return FieldR(HH * v.J - JH * v.H, -HJ * v.J + JJ * v.H) / (JJ * HH - JH * HJ); }
        MatrixR solve(const MatrixR& o) {
            MatrixR result(HH * o.JJ - JH * o.HJ, HH * o.JH - JH * o.HH, -HJ * o.JJ + JJ * o.HJ, -HJ * o.JH + JJ * o.HH);
            result /= (JJ * HH - JH * HJ);
            return result;
        }
    };
 
    enum Emission {
        TOP,    
        BOTTOM  
    };
 
    enum Determinant {
        DETERMINANT_INWARDS,   
        DETERMINANT_OUTWARDS,  
        DETERMINANT_FULL       
    };
 
    struct Mode {
        EffectiveFrequencyCyl* solver;                            
        int m;                                                    
        bool have_fields;                                         
        std::vector<FieldR, aligned_allocator<FieldR>> rfields;   
        std::vector<double, aligned_allocator<double>> rweights;  
        dcomplex lam;                                             
        double power;                                             
 
        Mode(EffectiveFrequencyCyl* solver, int m = 0)
            : solver(solver), m(m), have_fields(false), rfields(solver->rsize), rweights(solver->rsize), power(1.) {}
 
        bool operator==(const Mode& other) const { return m == other.m && is_zero(lam - other.lam); }
 
        dcomplex rField(double r) const {
            double Jr, Ji, Hr, Hi;
            long nz, ierr;
            size_t ir = solver->mesh->axis[0]->findIndex(r);
            if (ir > 0) --ir;
            if (ir >= solver->veffs.size()) ir = solver->veffs.size() - 1;
            dcomplex x = r * solver->k0 * sqrt(solver->nng[ir] * (solver->veffs[ir] - solver->freqv(lam)));
            if (real(x) < 0.) x = -x;
            if (imag(x) > SMALL) x = -x;
            if (ir == solver->rsize - 1) {
                Jr = Ji = 0.;
            } else {
                zbesj(x.real(), x.imag(), m, 1, 1, &Jr, &Ji, nz, ierr);
                if (ierr != 0) throw ComputationError(solver->getId(), "could not compute J({0}, {1}) @ r = {2}um", m, str(x), r);
            }
            if (ir == 0) {
                Hr = Hi = 0.;
            } else {
                zbesh(x.real(), x.imag(), m, 1, MH, 1, &Hr, &Hi, nz, ierr);
                if (ierr != 0) throw ComputationError(solver->getId(), "could not compute H({0}, {1}) @ r = {2}um", m, str(x), r);
            }
            return rfields[ir].J * dcomplex(Jr, Ji) + rfields[ir].H * dcomplex(Hr, Hi);
        }
 
        double loss() const { return imag(4e7 * PI / lam); }
    };
 
    dcomplex freqv(dcomplex lam) { return 2. - 4e3 * PI / lam / k0; }
 
    dcomplex lambda(dcomplex freq) { return 2e3 * PI / (k0 * (1. - freq / 2.)); }
 
  protected:
    friend struct RootDigger;
 
    DataLog<dcomplex, dcomplex> log_value;
 
    size_t rsize,  
        zbegin,    
        zsize;     
 
    std::vector<std::vector<dcomplex, aligned_allocator<dcomplex>>> nrCache;
 
    std::vector<std::vector<dcomplex, aligned_allocator<dcomplex>>> ngCache;
 
    std::vector<FieldZ> zfields;
 
    std::vector<double, aligned_allocator<double>> zintegrals;
 
    std::vector<dcomplex, aligned_allocator<dcomplex>> veffs;
 
    std::vector<dcomplex, aligned_allocator<dcomplex>> nng;
 
    dcomplex old_k0;
 
    Emission emission;
 
    void onInputChange(ReceiverBase&, ReceiverBase::ChangeReason) { cache_outdated = true; }
 
    int rstripe;
 
  public:
    Determinant determinant;
 
    int getStripe() const { return rstripe; }
 
    void setStripe(int stripe) {
        if (!mesh) setSimpleMesh();
        if (stripe < 0 || std::size_t(stripe) >= mesh->axis[0]->size()) throw BadInput(getId(), "wrong stripe number specified");
        rstripe = stripe;
        invalidate();
    }
 
    double getStripeR() const {
        if (rstripe == -1 || !mesh) return NAN;
        return mesh->axis[0]->at(rstripe);
    }
 
    void setStripeR(double r = 0.) {
        if (!mesh) setSimpleMesh();
        if (r < 0) throw BadInput(getId(), "radial position cannot be negative");
        rstripe = int(std::lower_bound(mesh->axis[0]->begin() + 1, mesh->axis[0]->end(), r) - mesh->axis[0]->begin() - 1);
        invalidate();
    }
 
    void useAllStripes() {
        rstripe = -1;
        invalidate();
    }
 
    dcomplex getDeltaNeff(double r) {
        stageOne();
        if (r < 0) throw BadInput(getId(), "radial position cannot be negative");
        size_t ir = mesh->axis[0]->findIndex(r);
        if (ir > 0) --ir;
        if (ir >= veffs.size()) ir = veffs.size() - 1;
        return sqrt(nng[ir] * veffs[ir]);
    }
 
    dcomplex getNNg(double r) {
        stageOne();
        if (r < 0) throw BadInput(getId(), "radial position cannot be negative");
        size_t ir = mesh->axis[0]->findIndex(r);
        if (ir > 0) --ir;
        if (ir >= veffs.size()) ir = veffs.size() - 1;
        return sqrt(nng[ir]);
    }
 
    // Parameters for rootdigger
    RootDigger::Params root;         
    RootDigger::Params stripe_root;  
 
    double perr;
 
    dcomplex k0;
 
    dcomplex vlam;
 
    std::vector<Mode> modes;
 
    ReceiverFor<Temperature, Geometry2DCylindrical> inTemperature;
 
    ReceiverFor<Gain, Geometry2DCylindrical> inGain;
 
    ReceiverFor<CarriersConcentration, Geometry2DCylindrical> inCarriersConcentration;
 
    typename ProviderFor<ModeWavelength>::Delegate outWavelength;
 
    typename ProviderFor<ModeLoss>::Delegate outLoss;
 
    typename ProviderFor<ModeLightMagnitude, Geometry2DCylindrical>::Delegate outLightMagnitude;
 
    typename ProviderFor<ModeLightE, Geometry2DCylindrical>::Delegate outLightE;
 
    typename ProviderFor<RefractiveIndex, Geometry2DCylindrical>::Delegate outRefractiveIndex;
 
    typename ProviderFor<Heat, Geometry2DCylindrical>::Delegate outHeat;
 
    EffectiveFrequencyCyl(const std::string& name = "");
 
    virtual ~EffectiveFrequencyCyl() {
        inTemperature.changedDisconnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
        inGain.changedDisconnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
        inCarriersConcentration.changedDisconnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
    }
 
    std::string getClassName() const override { return "optical.EffectiveFrequencyCyl"; }
 
    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;
 
    Emission getEmission() const { return emission; }
 
    void setEmission(Emission emis) {
        emission = emis;
        for (auto& mode : modes) mode.have_fields = false;
    }
 
    void setSimpleMesh() {
        writelog(LOG_DETAIL, "Creating simple mesh");
        setMesh(plask::make_shared<RectangularMesh2DSimpleGenerator>());
    }
 
    void setHorizontalMesh(shared_ptr<MeshAxis> meshx) {
        writelog(LOG_DETAIL, "Setting horizontal mesh");
        if (!geometry) throw NoChildException();
        auto meshxy = makeGeometryGrid(geometry->getChild());
        meshxy->setTran(meshx);
        setMesh(meshxy);
    }
 
    bool getAsymptotic() const { return asymptotic; }
 
    void setAsymptotic(bool value) {
        asymptotic = value;
        invalidate();
    }
 
    size_t findMode(dcomplex lambda, int m = 0);
 
    std::vector<size_t> findModes(plask::dcomplex lambda1 = 0.,
                                  plask::dcomplex lambda2 = 0.,
                                  int m = 0,
                                  size_t resteps = 256,
                                  size_t imsteps = 64,
                                  dcomplex eps = dcomplex(1e-6, 1e-9));
 
    dcomplex getVertDeterminant(dcomplex vlambda) {
        updateCache();
        if (rstripe < 0) throw BadInput(getId(), "this works only for the weighted approach");
        if (vlam == 0. && isnan(k0.real())) throw BadInput(getId(), "no reference wavelength `lam0` specified");
        dcomplex v = freqv(vlambda);
        return this->detS1(v, nrCache[rstripe], ngCache[rstripe]);
    }
 
    dcomplex getDeterminant(dcomplex lambda, int m = 0) {
        if (isnan(k0.real())) throw BadInput(getId(), "no reference wavelength `lam0` specified");
        stageOne();
        Mode mode(this, m);
        dcomplex det = detS(lambda, mode);
        // log_value(v, det);
        return det;
    }
 
    size_t setMode(dcomplex clambda, int m = 0);
 
    inline size_t setMode(double lambda, double loss, int m = 0) {
        return setMode(dcomplex(lambda, -lambda * lambda / (4e7 * PI) * loss), m);
    }
 
    void clearModes() { modes.clear(); }
 
    double getTotalAbsorption(Mode& mode);
 
    double getTotalAbsorption(size_t num);
 
    double getGainIntegral(Mode& mode);
 
    double getGainIntegral(size_t num);
 
  protected:
    bool need_gain;
 
    bool cache_outdated;
 
    bool have_veffs;
 
    bool asymptotic;
 
    void onInitialize() override;
 
    void onInvalidate() override;
 
    void updateCache();
 
    void stageOne();
 
    dcomplex detS1(const dcomplex& v,
                   const std::vector<dcomplex, aligned_allocator<dcomplex>>& NR,
                   const std::vector<dcomplex, aligned_allocator<dcomplex>>& NG,
                   std::vector<FieldZ>* saveto = nullptr);
 
    void computeStripeNNg(size_t stripe, bool save_integrals = false);
 
    double integrateBessel(Mode& mode);
 
    void computeBessel(size_t i, dcomplex v, const Mode& mode, dcomplex* J1, dcomplex* H1, dcomplex* J2, dcomplex* H2);
 
    dcomplex detS(const plask::dcomplex& lam, Mode& mode, bool save = false);
 
    size_t getMainStripe() {
        if (rstripe < 0) {
            size_t stripe = 0;
            // Look for the innermost stripe with not constant refractive index
            bool all_the_same = true;
            while (all_the_same) {
                dcomplex same_nr = nrCache[stripe].front();
                dcomplex same_ng = ngCache[stripe].front();
                for (auto nr = nrCache[stripe].begin(), ng = ngCache[stripe].begin(); nr != nrCache[stripe].end(); ++nr, ++ng)
                    if (*nr != same_nr || *ng != same_ng) {
                        all_the_same = false;
                        break;
                    }
                if (all_the_same) ++stripe;
            }
            writelog(LOG_DETAIL, "Vertical field distribution taken from stripe {0}", stripe);
            return stripe;
        } else {
            return rstripe;
        }
    }
 
    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);
        outWavelength.fireChanged();
        outLoss.fireChanged();
        outLightMagnitude.fireChanged();
        outLightE.fireChanged();
        return modes.size() - 1;
    }
 
    size_t nmodes() const { return modes.size(); }
 
    double getWavelength(size_t n) {
        if (n >= modes.size()) throw NoValue(ModeWavelength::NAME);
        return real(modes[n].lam);
    }
 
    double getModalLoss(size_t n) {
        if (n >= modes.size()) throw NoValue(ModeLoss::NAME);
        return imag(4e7 * PI / modes[n].lam);  // 2e4  2/µm -> 2/cm
    }
 
    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,
                                             const shared_ptr<const MeshD<2>>& dst_mesh,
                                             InterpolationMethod = INTERPOLATION_DEFAULT);
 
    const LazyData<Vec<3, dcomplex>> getElectricField(std::size_t num,
                                                      const shared_ptr<const plask::MeshD<2>>& dst_mesh,
                                                      plask::InterpolationMethod = INTERPOLATION_DEFAULT);
 
    const LazyData<dcomplex> getRefractiveIndex(RefractiveIndex::EnumType component,
                                                const shared_ptr<const MeshD<2>>& dst_mesh,
                                                dcomplex lam,
                                                InterpolationMethod = INTERPOLATION_DEFAULT);
 
    const LazyData<double> getHeat(const shared_ptr<const MeshD<2>>& dst_mesh, InterpolationMethod method = INTERPOLATION_DEFAULT);
};
 
}}}  // namespace plask::optical::effective
 
#endif  // PLASK__MODULE_OPTICAL_EFM_HPP