solver3d.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_3D_H
#define PLASK__SOLVER_OPTICAL_MODAL_FOURIER_REFLECTION_3D_H
 
#include <plask/plask.hpp>
 
#include "../solver.hpp"
#include "../reflection.hpp"
#include "expansion3d.hpp"
 
#ifdef minor
#   undef minor
#endif
 
namespace plask { namespace optical { namespace modal {
 
struct PLASK_SOLVER_API FourierSolver3D: public ModalSolver<SolverOver<Geometry3D>> {
 
    typedef ModalSolver<SolverOver<Geometry3D>> BaseType;
 
    friend struct ExpansionPW3D;
 
    std::string getClassName() const override { return "optical.Fourier3D"; }
 
    enum What {
        WHAT_WAVELENGTH,    
        WHAT_K0,            
        WHAT_KLONG,         
        WHAT_KTRAN          
    };
 
    enum ExpansionRule {
        RULE_OLD,
        RULE_DIRECT,
        RULE_INVERSE,
        RULE_COMBINED
    };
 
    struct Mode {
        Expansion::Component symmetry_long;     
        Expansion::Component symmetry_tran;     
        double lam0;                            
        dcomplex k0;                            
        dcomplex klong;                         
        dcomplex ktran;                         
        double power;                           
        double tolx;                            
 
        Mode(const ExpansionPW3D& expansion, double tolx):
            symmetry_long(expansion.symmetry_long),
            symmetry_tran(expansion.symmetry_tran),
            lam0(expansion.lam0),
            k0(expansion.k0),
            klong(expansion.klong),
            ktran(expansion.ktran),
            power(1.),
            tolx(tolx) {}
 
        bool operator==(const Mode& other) const {
            return is_equal(k0, other.k0) && is_equal(klong, other.klong) && is_equal(ktran, other.ktran)
                && symmetry_long == other.symmetry_long && symmetry_tran == other.symmetry_tran &&
                ((isnan(lam0) && isnan(other.lam0)) || lam0 == other.lam0)
            ;
        }
 
        bool operator==(const ExpansionPW3D& other) const {
            return is_equal(k0, other.k0) && is_equal(klong, other.klong) && is_equal(ktran, other.ktran)
                && symmetry_long == other.symmetry_long && symmetry_tran == other.symmetry_tran &&
                ((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;
        }
    };
 
    size_t size_long;
    size_t size_tran;
 
  protected:
 
    dcomplex klong,                             
             ktran;                             
 
    Expansion::Component symmetry_long,         
                         symmetry_tran;         
 
    void onInitialize() override;
 
    void onInvalidate() override;
 
    void computeIntegrals() override {
        expansion.computeIntegrals();
    }
 
    int dct;
 
    ExpansionRule expansion_rule;
 
    size_t initIncidence(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam);
 
  public:
 
    ExpansionPW3D 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.getKlong() != getKlong()) { changed = true; expansion.setKlong(getKlong()); }
        if (expansion.getKtran() != getKtran()) { changed = true; expansion.setKtran(getKtran()); }
        if (expansion.getSymmetryLong() != getSymmetryLong()) { changed = true; expansion.setSymmetryLong(getSymmetryLong()); }
        if (expansion.getSymmetryTran() != getSymmetryTran()) { changed = true; expansion.setSymmetryTran(getSymmetryTran()); }
        return changed;
    }
 
    size_t refine_long;
    size_t refine_tran;
 
    double grad_smooth;
 
    PML pml_long;
    PML pml_tran;
 
    ProviderFor<GradientFunctions, Geometry3D>::Delegate outGradients;
 
    FourierSolver3D(const std::string& name="");
 
    void loadConfiguration(XMLReader& reader, Manager& manager) override;
 
    size_t findMode(What what, dcomplex start);
 
    size_t getLongSize() const { return size_long; }
 
    size_t getTranSize() const { return size_tran; }
 
    void setLongSize(size_t n) {
        size_long = n;
        invalidate();
    }
    void setTranSize(size_t n) {
        size_tran = n;
        invalidate();
    }
 
    void setSizes(size_t nl, size_t nt) {
        size_long = nl;
        size_tran = nt;
        invalidate();
    }
 
    Expansion::Component getSymmetryLong() const { return symmetry_long; }
 
    void setSymmetryLong(Expansion::Component symmetry) {
        if (symmetry != Expansion::E_UNSPECIFIED && geometry && !geometry->isSymmetric(Geometry3D::DIRECTION_LONG))
            throw BadInput(getId(), "longitudinal symmetry not allowed for asymmetric structure");
        if ((symmetry_long == Expansion::E_UNSPECIFIED) != (symmetry == Expansion::E_UNSPECIFIED))
            invalidate();
        if (klong != 0. && symmetry != Expansion::E_UNSPECIFIED) {
            Solver::writelog(LOG_WARNING, "Resetting klong to 0.");
            klong = 0.;
            expansion.setKlong(0.);
        }
        symmetry_long = symmetry;
    }
 
    Expansion::Component getSymmetryTran() const { return symmetry_tran; }
 
    void setSymmetryTran(Expansion::Component symmetry) {
        if (symmetry != Expansion::E_UNSPECIFIED && geometry && !geometry->isSymmetric(Geometry3D::DIRECTION_TRAN))
            throw BadInput(getId(), "transverse symmetry not allowed for asymmetric structure");
        if ((symmetry_tran == Expansion::E_UNSPECIFIED) != (symmetry == Expansion::E_UNSPECIFIED))
            invalidate();
        if (ktran != 0. && symmetry != Expansion::E_UNSPECIFIED) {
            Solver::writelog(LOG_WARNING, "Resetting ktran to 0.");
            ktran = 0.;
            expansion.setKtran(0.);
        }
        symmetry_tran = symmetry;
    }
 
    bool symmetricLong() const { return expansion.symmetric_long(); }
 
    bool symmetricTran() const { return expansion.symmetric_tran(); }
 
    dcomplex getKlong() const { return klong; }
 
    void setKlong(dcomplex k)  {
        if (k != 0. && (expansion.symmetric_long() || symmetry_long != Expansion::E_UNSPECIFIED)) {
            Solver::writelog(LOG_WARNING, "Resetting longitudinal mode symmetry");
            symmetry_long = Expansion::E_UNSPECIFIED;
            invalidate();
        }
        klong = k;
    }
 
    dcomplex getKtran() const { return ktran; }
 
    void setKtran(dcomplex k)  {
        if (k != 0. && (expansion.symmetric_tran() || symmetry_tran != Expansion::E_UNSPECIFIED)) {
            Solver::writelog(LOG_WARNING, "Resetting transverse mode symmetry");
            symmetry_tran = Expansion::E_UNSPECIFIED;
            invalidate();
        }
        ktran = k;
    }
 
    double getGradSmooth() const { return grad_smooth; }
 
    void setGradSmooth(double value) {
        bool changed = grad_smooth != value;
        grad_smooth = value;
        if (changed) this->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 (expansion.symmetric_long() || expansion.symmetric_tran()) invalidate();
        }
    }
    bool dct2() const { return dct == 2; }
 
    ExpansionRule getRule() const { return expansion_rule; }
    void setRule(ExpansionRule rule) {
        if (rule != expansion_rule) {
            expansion_rule = rule;
            invalidate();
        }
    }
 
    // /// Get mesh at which material parameters are sampled along longitudinal axis
    // RegularAxis getLongMesh() const { return expansion.mesh->lon(); }
    //
    // /// Get mesh at which material parameters are sampled along transverse axis
    // RegularAxis getTranMesh() const { return expansion.mesh->tran(); }
 
    Expansion& getExpansion() override { return expansion; }
 
    size_t minor() const { return expansion.Nl; }
 
    cvector incidentVector(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam=NAN);
 
    using BaseType::incidentVector;
 
    cvector incidentGaussian(Transfer::IncidentDirection side, Expansion::Component polarization,
                             double sigma_long, double sigma_tran, double center_long=0., double center_tran=0., dcomplex lam=NAN);
 
    LazyData<Vec<3,dcomplex>> getScatteredFieldE(const cvector& incident,
                                                 Transfer::IncidentDirection side,
                                                 const shared_ptr<const MeshD<3>>& 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<3>>& 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<3>>& 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) {
            writelog(LOG_WARNING, "Mode fields are not normalized unless emission is set to 'top' or 'bottom'");
            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);
        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].klong / modes[n].k0;
    }
 
    void applyMode(const Mode& mode) {
        writelog(LOG_DEBUG, "Current mode <lam: {}nm, klong: {}/um, ktran: {}/um, symmetry: ({},{})>",
                 str(2e3*PI/mode.k0, "({:.3f}{:+.3g}j)", "{:.3f}"),
                 str(mode.klong, "({:.3f}{:+.3g}j)", "{:.3f}"),
                 str(mode.ktran, "({:.3f}{:+.3g}j)", "{:.3f}"),
                 (mode.symmetry_long == Expansion::E_LONG)? "El" : (mode.symmetry_long == Expansion::E_TRAN)? "Et" : "none",
                 (mode.symmetry_tran == Expansion::E_LONG)? "El" : (mode.symmetry_tran == Expansion::E_TRAN)? "Et" : "none"
                );
        if (mode != expansion) {
            expansion.setLam0(mode.lam0);
            expansion.setK0(mode.k0);
            expansion.klong = mode.klong;
            expansion.ktran = mode.ktran;
            expansion.symmetry_long = mode.symmetry_long;
            expansion.symmetry_tran = mode.symmetry_tran;
            clearFields();
        }
    }
 
    double getWavelength(size_t n) override;
 
    LazyData<double> getGradients(GradientFunctions::EnumType what,
                                  const shared_ptr<const MeshD<3>>& dst_mesh,
                                  InterpolationMethod interp);
 
};
 
 
}}} // namespace
 
#endif