solvercyl.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__SOLVER__SLAB_SOLVERCYL_H
#define PLASK__SOLVER__SLAB_SOLVERCYL_H
 
#include <plask/plask.hpp>
 
#include "../solver.hpp"
#include "../reflection.hpp"
#include "expansioncyl-fini.hpp"
#include "expansioncyl-infini.hpp"
 
 
namespace plask { namespace optical { namespace modal {
 
struct PLASK_SOLVER_API BesselSolverCyl: public ModalSolver<SolverWithMesh<Geometry2DCylindrical,MeshAxis>> {
 
    typedef ModalSolver<SolverWithMesh<Geometry2DCylindrical,MeshAxis>> BaseType;
 
    friend struct ExpansionBessel;
    friend struct ExpansionBesselFini;
    friend struct ExpansionBesselInfini;
 
    enum BesselDomain {
        DOMAIN_FINITE,
        DOMAIN_INFINITE
    };
 
    enum InfiniteWavevectors {
        WAVEVECTORS_UNIFORM,
        WAVEVECTORS_NONUNIFORM,
        // WAVEVECTORS_LEGENDRE,
        WAVEVECTORS_LAGUERRE,
        WAVEVECTORS_MANUAL
    };
 
    enum Rule {
        RULE_DIRECT,
        RULE_COMBINED_1,
        RULE_COMBINED_2,
        RULE_OLD
    };
 
    const char* ruleName() {
        switch (rule) {
            case RULE_DIRECT: return "inverse";
            case RULE_COMBINED_1: return "inverse (mod 1)";
            case RULE_COMBINED_2: return "inverse (mod 2)";
            case RULE_OLD: return "direct";
        }
        return "unknown";
    }
 
    std::string getClassName() const override { return "optical.BesselCyl"; }
 
    struct Mode {
        double lam0;                    
        dcomplex k0;                    
        int m;                          
        double power;                   
        double tolx;                    
 
        Mode(const std::unique_ptr<ExpansionBessel>& expansion, double tolx):
            lam0(expansion->lam0), k0(expansion->k0), m(expansion->m), power(1.), tolx(tolx) {}
 
        bool operator==(const Mode& other) const {
            return m == other.m && is_equal(k0, other.k0) && is_equal(lam0, other.lam0) &&
                   ((isnan(lam0) && isnan(other.lam0)) || lam0 == other.lam0);
        }
 
        bool operator==(const std::unique_ptr<ExpansionBessel>& other) const {
            return m == other->m && is_equal(k0, other->k0) && is_equal(lam0, other->lam0) &&
                   ((isnan(lam0) && isnan(other->lam0)) || lam0 == other->lam0);
        }
 
        template <typename T>
        bool operator!=(const T& other) const {
            return !(*this == other);
        }
 
      private:
 
        template <typename T>
        bool is_equal(T a, T b) const {
            return abs(a-b) <= tolx;
        }
    };
 
  protected:
 
    BesselDomain domain;
 
    int m;
 
    size_t size;
 
    void onInitialize() override;
 
    void onInvalidate() override;
 
    void computeIntegrals() override {
        expansion->computeIntegrals();
    }
 
    Rule rule;
 
    double kscale;
 
    double kmax;
 
    InfiniteWavevectors kmethod;
 
  public:
 
    std::vector<double> klist;
 
    boost::optional<std::vector<double>> kweights;
 
    std::unique_ptr<ExpansionBessel> expansion;
 
    std::vector<Mode> modes;
 
    std::vector<double> getKlist() {
        initCalculation();
        computeIntegrals();
        return expansion->getKpts();
    }
 
    void setKlist(const std::vector<double>& values) {
        if (kmethod != WAVEVECTORS_MANUAL) {
            invalidate();
            this->writelog(LOG_WARNING, "Setting Hankel transform method to Manual");
            kmethod = WAVEVECTORS_MANUAL;
        }
        klist = values;
    }
 
    std::vector<double> getKweights() {
        if (domain == DOMAIN_INFINITE) {
            initCalculation();
            computeIntegrals();
            ExpansionBesselInfini* ex = dynamic_cast<ExpansionBesselInfini*>(expansion.get());
            if (ex) return std::vector<double>(ex->kdelts.begin(), ex->kdelts.end());
        }
        return std::vector<double>();
    }
 
    void setKweights(const std::vector<double>& values) {
        if (kmethod != WAVEVECTORS_MANUAL) {
            invalidate();
            this->writelog(LOG_WARNING, "Setting Hankel transform method to Manual");
            kmethod = WAVEVECTORS_MANUAL;
        }
        kweights.reset(values);
    }
 
    void clearKweights() {
        kweights.reset();
    }
 
    void clearModes() override {
        modes.clear();
    }
 
    bool setExpansionDefaults(bool with_k0=true) override {
        bool changed = false;
        if (expansion->getLam0() != getLam0()) { changed = true; expansion->setLam0(getLam0()); }
        if (with_k0) {
            if (expansion->getK0() != getK0()) { changed = true; expansion->setK0(getK0()); }
        }
        if (expansion->getM() != getM()) { changed = true; expansion->setM(getM()); }
        return changed;
    }
 
    double integral_error;
 
    size_t max_integration_points;
 
    PML pml;
 
    typename ProviderFor<ModeLoss>::Delegate outLoss;
 
    BesselSolverCyl(const std::string& name="");
 
    void loadConfiguration(XMLReader& reader, Manager& manager) override;
 
    size_t findMode(dcomplex start, int m=1);
 
    size_t getSize() const { return size; }
    void setSize(size_t n) {
        size = n;
        invalidate();
    }
 
    double getKscale() const { return kscale; }
    void setKscale(double s) { kscale = s; }
 
    double getKmax() const { return kmax; }
    void setKmax(double s) { kmax = s; }
 
    InfiniteWavevectors getKmethod() const { return kmethod; }
    void setKmethod(InfiniteWavevectors k) { kmethod = k; }
 
    BesselDomain getDomain() const { return domain; }
    void setDomain(BesselDomain dom) {
        domain = dom;
        invalidate();
    }
 
    Rule getRule() const { return rule; }
    void setRule(Rule r) {
        rule = r;
        invalidate();
    }
 
    unsigned getM() const { return m; }
    void setM(unsigned n) { m = n; }
 
    Expansion& getExpansion() override { return *expansion; }
 
    LazyData<Vec<3,dcomplex>> getScatteredFieldE(const cvector& incident,
                                                 Transfer::IncidentDirection side,
                                                 const shared_ptr<const MeshD<2>>& dst_mesh,
                                                 InterpolationMethod method,
                                                 PropagationDirection part = PROPAGATION_TOTAL) {
        if (!Solver::initCalculation()) setExpansionDefaults(false);
        if (!transfer) initTransfer(*expansion, true);
        return transfer->getScatteredFieldE(incident, side, dst_mesh, method, part);
    }
 
    LazyData<Vec<3,dcomplex>> getScatteredFieldH(const cvector& incident,
                                                 Transfer::IncidentDirection side,
                                                 const shared_ptr<const MeshD<2>>& dst_mesh,
                                                 InterpolationMethod method,
                                                 PropagationDirection part = PROPAGATION_TOTAL) {
        if (!Solver::initCalculation()) setExpansionDefaults(false);
        if (!transfer) initTransfer(*expansion, true);
        return transfer->getScatteredFieldH(incident, side, dst_mesh, method, part);
    }
 
    LazyData<double> getScatteredFieldMagnitude(const cvector& incident,
                                                Transfer::IncidentDirection side,
                                                const shared_ptr<const MeshD<2>>& dst_mesh,
                                                InterpolationMethod method) {
        if (!Solver::initCalculation()) setExpansionDefaults(false);
        if (!transfer) initTransfer(*expansion, true);
        return transfer->getScatteredFieldMagnitude(incident, side, dst_mesh, method);
    }
 
    cvector getFieldVectorE(size_t num, double z) {
        applyMode(modes[num]);
        return transfer->getFieldVectorE(z);
    }
 
    cvector getFieldVectorH(size_t num, double z) {
        applyMode(modes[num]);
        return transfer->getFieldVectorH(z);
    }
 
    double getIntegralEE(size_t num, double z1, double z2) {
        applyMode(modes[num]);
        return transfer->getFieldIntegral(FIELD_E, z1, z2, modes[num].power);
    }
 
    double getIntegralHH(size_t num, double z1, double z2) {
        applyMode(modes[num]);
        return transfer->getFieldIntegral(FIELD_H, z1, z2, modes[num].power);
    }
 
    cvector getScatteredFieldVectorE(const cvector& incident, Transfer::IncidentDirection side, double z) {
        if (!Solver::initCalculation()) setExpansionDefaults(false);
        if (!transfer) initTransfer(*expansion, true);
        return transfer->getScatteredFieldVectorE(incident, side, z, PROPAGATION_TOTAL);
    }
 
    cvector getScatteredFieldVectorH(const cvector& incident, Transfer::IncidentDirection side, double z) {
        if (!Solver::initCalculation()) setExpansionDefaults(false);
        if (!transfer) initTransfer(*expansion, true);
        return transfer->getScatteredFieldVectorH(incident, side, z, PROPAGATION_TOTAL);
    }
 
    double getScatteredIntegralEE(const cvector& incident, Transfer::IncidentDirection side, double z1, double z2) {
        if (!Solver::initCalculation()) setExpansionDefaults(false);
        if (!transfer) initTransfer(*expansion, true);
        return transfer->getScatteredFieldIntegral(FIELD_E, incident, side, z1, z2);
    }
 
    double getScatteredIntegralHH(const cvector& incident, Transfer::IncidentDirection side, double z1, double z2) {
        if (!Solver::initCalculation()) setExpansionDefaults(false);
        if (!transfer) initTransfer(*expansion, true);
        return transfer->getScatteredFieldIntegral(FIELD_H, incident, side, z1, z2);
    }
 
    size_t setMode() {
        if (abs2(this->getDeterminant()) > root.tolf_max*root.tolf_max)
            throw BadInput(this->getId(), "cannot set the mode, determinant too large");
        return insertMode();
    }
 
  protected:
    size_t insertMode() {
        static bool warn = true;
        if (warn && ((emission != EMISSION_TOP && emission != EMISSION_BOTTOM) || domain == DOMAIN_INFINITE)) {
            if (domain == DOMAIN_INFINITE)
                writelog(LOG_WARNING, "Mode fields are not normalized (infinite domain)");
            else
                writelog(LOG_WARNING, "Mode fields are not normalized (emission direction not specified)");
            warn = false;
        }
        Mode mode(expansion, root.tolx);
        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();
        outLightH.fireChanged();
        return modes.size()-1;
    }
 
    size_t nummodes() const override { return modes.size(); }
 
    double applyMode(size_t n) override {
        if (n >= modes.size()) throw BadInput(this->getId(), "mode {0} has not been computed", n);
        applyMode(modes[n]);
        return modes[n].power;
    }
 
    void applyMode(const Mode& mode) {
        writelog(LOG_DEBUG, "Current mode <m: {:d}, lam: {}nm>", mode.m, str(2e3*PI/mode.k0, "({:.3f}{:+.3g}j)"));
        expansion->setLam0(mode.lam0);
        expansion->setK0(mode.k0);
        expansion->setM(mode.m);
    }
 
    double getModalLoss(size_t n) {
        if (n >= modes.size()) throw NoValue(ModeLoss::NAME);
        return 2e4 * modes[n].k0.imag();  // 2e4  2/µm -> 2/cm
    }
 
    double getWavelength(size_t n) override;
 
#ifndef NDEBUG
  public:
    cmatrix epsV_k(size_t layer);
    cmatrix epsTss(size_t layer);
    cmatrix epsTsp(size_t layer);
    cmatrix epsTps(size_t layer);
    cmatrix epsTpp(size_t layer);
    cmatrix muV_k();
    cmatrix muTss();
    cmatrix muTsp();
    cmatrix muTps();
    cmatrix muTpp();
#endif
 
};
 
 
 
}}} // # namespace plask::optical::modal
 
#endif // PLASK__SOLVER__SLAB_SOLVERCYL_H