solver2d.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_OPTICAL_MODAL_FOURIER_REFLECTION_2D_H
#define PLASK__SOLVER_OPTICAL_MODAL_FOURIER_REFLECTION_2D_H
 
#include <plask/plask.hpp>
 
#include "../solver.hpp"
#include "expansion2d.hpp"
 
namespace plask { namespace optical { namespace modal {
 
struct PLASK_SOLVER_API FourierSolver2D: public ModalSolver<SolverWithMesh<Geometry2DCartesian,MeshAxis>> {
 
    typedef ModalSolver<SolverWithMesh<Geometry2DCartesian,MeshAxis>> BaseType;
 
    friend struct ExpansionPW2D;
 
    std::string getClassName() const override { return "optical.Fourier2D"; }
 
    enum What {
        WHAT_WAVELENGTH,    
        WHAT_K0,            
        WHAT_NEFF,          
        WHAT_KTRAN,         
        WHAT_BETA           
    };
 
 
    struct Mode {
        Expansion::Component symmetry;          
        Expansion::Component polarization;      
        double lam0;                            
        dcomplex k0;                            
        dcomplex beta;                          
        dcomplex ktran;                         
        double power;                           
        double tolx;                            
 
        Mode(const ExpansionPW2D& expansion, double tolx):
            symmetry(expansion.symmetry),
            polarization(expansion.polarization),
            lam0(expansion.lam0),
            k0(expansion.k0),
            beta(expansion.beta),
            ktran(expansion.ktran),
            power(1.),
            tolx(tolx) {}
 
        bool operator==(const Mode& other) const {
            return is_equal(k0, other.k0) && is_equal(beta, other.beta) && is_equal(ktran, other.ktran)
                && symmetry == other.symmetry && polarization == other.polarization &&
                ((isnan(lam0) && isnan(other.lam0)) || lam0 == other.lam0)
            ;
        }
 
        bool operator==(const ExpansionPW2D& other) const {
            return is_equal(k0, other.k0) && is_equal(beta, other.beta) && is_equal(ktran, other.ktran)
                && symmetry == other.symmetry && polarization == other.polarization &&
                ((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;
        }
    };
 
    enum FourierType {
        FOURIER_DISCRETE,
        FOURIER_ANALYTIC
    };
 
  protected:
 
    dcomplex beta,                     
             ktran;                     
 
    Expansion::Component symmetry;      
    Expansion::Component polarization;  
 
    size_t size;
 
    void onInitialize() override;
 
    void onInvalidate() override;
 
    void computeIntegrals() override {
        expansion.computeIntegrals();
    }
 
    int dct;
 
    FourierType ftt;
 
  public:
 
    ExpansionPW2D expansion;
 
    std::vector<Mode> modes;
 
    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.getBeta() != getBeta()) { changed = true; expansion.setBeta(getBeta()); }
        if (expansion.getKtran() != getKtran()) { changed = true; expansion.setKtran(getKtran()); }
        if (expansion.getSymmetry() != getSymmetry()) { changed = true; expansion.setSymmetry(getSymmetry()); }
        if (expansion.getPolarization() != getPolarization()) { changed = true; expansion.setPolarization(getPolarization()); }
        return changed;
    }
 
    size_t refine;
 
    PML pml;
 
    plask::optional<std::pair<double,double>> mirrors;
 
    ProviderFor<ModeEffectiveIndex>::Delegate outNeff;
 
    FourierSolver2D(const std::string& name="");
 
    void loadConfiguration(XMLReader& reader, Manager& manager) override;
 
    void setSimpleMesh() {
        writelog(LOG_DETAIL, "Creating simple mesh");
        setMesh(plask::make_shared<OrderedMesh1DSimpleGenerator>());
    }
 
    size_t findMode(FourierSolver2D::What what, dcomplex start);
 
    size_t getSize() const { return size; }
    void setSize(size_t n) {
        size = n;
        invalidate();
    }
 
    int getDCT() const { return dct; }
    void setDCT(int n) {
        if (n != 1 && n != 2)
            throw BadInput(getId(), "bad DCT type (can be only 1 or 2)");
        if (dct != n) {
            dct = n;
            if (symmetric() && ftt == FOURIER_DISCRETE) invalidate();
        }
    }
    bool dct2() const { return dct == 2; }
 
 
    FourierType getFourierType() const { return ftt; }
    void setFourierType(FourierType ft) {
        if (ft != ftt) {
            ftt = ft;
            invalidate();
        }
    }
 
    dcomplex getKtran() const { return ktran; }
 
    void setKtran(dcomplex k)  {
        if (k != 0. && (symmetric() || symmetry != Expansion::E_UNSPECIFIED)) {
            Solver::writelog(LOG_WARNING, "Resetting mode symmetry");
            symmetry = Expansion::E_UNSPECIFIED;
            invalidate();
        }
        if (k != ktran && transfer) transfer->fields_determined = Transfer::DETERMINED_NOTHING;
        ktran = k;
    }
 
    void setBeta(dcomplex k)  {
        if (k != 0. && (separated() || polarization != Expansion::E_UNSPECIFIED)) {
            Solver::writelog(LOG_WARNING, "Resetting polarizations separation");
            polarization = Expansion::E_UNSPECIFIED;
            invalidate();
        }
        if (k != beta && transfer) transfer->fields_determined = Transfer::DETERMINED_NOTHING;
        beta = k;
    }
 
    dcomplex getBeta() const { return beta; }
 
    Expansion::Component getSymmetry() const { return symmetry; }
 
    void setSymmetry(Expansion::Component sym) {
        if (sym != Expansion::E_UNSPECIFIED && geometry && !geometry->isSymmetric(Geometry2DCartesian::DIRECTION_TRAN))
            throw BadInput(getId(), "symmetry not allowed for asymmetric structure");
        if ((symmetry == Expansion::E_UNSPECIFIED) != (sym == Expansion::E_UNSPECIFIED))
            invalidate();
        if (ktran != 0. && sym != Expansion::E_UNSPECIFIED) {
            Solver::writelog(LOG_WARNING, "Resetting ktran to 0.");
            ktran = 0.;
            expansion.setKtran(0.);
        }
        symmetry = sym;
    }
 
    Expansion::Component getPolarization() const { return polarization; }
 
    void setPolarization(Expansion::Component pol) {
        if (polarization != pol)
            invalidate();
        if (beta != 0. && pol != Expansion::E_UNSPECIFIED) {
            Solver::writelog(LOG_WARNING, "Resetting beta to 0.");
            beta = 0.;
            expansion.setBeta(0.);
        }
        polarization = pol;
    }
 
    bool symmetric() const { return symmetry != Expansion::E_UNSPECIFIED; }
 
    bool separated() const { return polarization != Expansion::E_UNSPECIFIED; }
 
    Expansion& getExpansion() override { return expansion; }
 
    // RegularAxis getMesh() const { return *expansion.mesh->tran(); }
 
    cvector incidentVector(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam=NAN);
 
    using BaseType::incidentVector;
 
    cvector incidentGaussian(Transfer::IncidentDirection side, Expansion::Component polarization,
                             double sigma, double center=0., dcomplex lam=NAN);
 
  private:
 
    size_t initIncidence(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam=NAN);
 
  public:
 
    LazyData<Vec<3,dcomplex>> getScatteredFieldE(const cvector& incident,
                                                 Transfer::IncidentDirection side,
                                                 shared_ptr<const MeshD<2>> dst_mesh,
                                                 InterpolationMethod method,
                                                 PropagationDirection part = PROPAGATION_TOTAL) {
        if (!Solver::initCalculation()) setExpansionDefaults();
        // expansion.setK0(2e3*PI / wavelength);
        if (expansion.separated())
            expansion.setPolarization(polarization);
        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,
                                                 shared_ptr<const MeshD<2>> dst_mesh,
                                                 InterpolationMethod method,
                                                 PropagationDirection part = PROPAGATION_TOTAL) {
        if (!Solver::initCalculation()) setExpansionDefaults();
        // expansion.setK0(2e3*PI / wavelength);
        if (expansion.separated())
            expansion.setPolarization(polarization);
        if (!transfer) initTransfer(expansion, true);
        return transfer->getScatteredFieldH(incident, side, dst_mesh, method, part);
    }
 
    LazyData<double> getScatteredFieldMagnitude(const cvector& incident,
                                                Transfer::IncidentDirection side,
                                                shared_ptr<const MeshD<2>> dst_mesh,
                                                InterpolationMethod method) {
        if (!Solver::initCalculation()) setExpansionDefaults();
        // expansion.setK0(2e3*PI / wavelength);
        if (expansion.separated())
            expansion.setPolarization(polarization);
        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) {
            writelog(LOG_WARNING, "Mode fields are not normalized");
            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);
        outNeff.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;
    }
 
    dcomplex getEffectiveIndex(size_t n) {
        if (n >= modes.size()) throw NoValue(ModeEffectiveIndex::NAME);
        return modes[n].beta / modes[n].k0;
    }
 
    double getMirrorLosses(double 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 = (n-n1) / (n+n1); R1 *= R1;
            R2 = (n-n2) / (n+n2); R2 *= R2;
        }
        return 0.5 * std::log(R1*R2) / L;
    }
 
    void applyMode(const Mode& mode) {
        writelog(LOG_DEBUG, "Current mode <lam: {:.2f}nm, neff: {}, ktran: {}/um, polarization: {}, symmetry: {}>",
                 real(2e3*PI/mode.k0),
                 str(mode.beta/mode.k0, "{:.3f}{:+.3g}j"),
                 str(mode.ktran, "({:.3g}{:+.3g}j)", "{:.3g}"),
                 (mode.polarization == Expansion::E_LONG)? "El" : (mode.polarization == Expansion::E_TRAN)? "Et" : "none",
                 (mode.symmetry == Expansion::E_LONG)? "El" : (mode.symmetry == Expansion::E_TRAN)? "Et" : "none"
                );
        if (mode != expansion) {
            expansion.setLam0(mode.lam0);
            expansion.setK0(mode.k0);
            expansion.beta = mode.beta;
            expansion.ktran = mode.ktran;
            expansion.symmetry = mode.symmetry;
            expansion.polarization = mode.polarization;
            clearFields();
        }
    }
 
    double getWavelength(size_t n) override;
};
 
 
}}} // namespace
 
#endif