diagonalizer.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_DIAGONALIZER_H
#define PLASK__SOLVER_OPTICAL_MODAL_DIAGONALIZER_H
 
#include <utility>
 
#ifdef OPENMP_FOUND
#   include <omp.h>
#endif
 
#include <plask/plask.hpp>
 
#include "matrices.hpp"
#include "expansion.hpp"
 
namespace plask { namespace optical { namespace modal {
 
class Diagonalizer
{
  protected:
    Expansion* src;                     
    std::vector<bool> diagonalized;     
 
  public:
    const std::size_t lcount;           
 
    Diagonalizer(Expansion* src);
 
    virtual ~Diagonalizer();
 
    virtual std::size_t matrixSize() const = 0;
 
    inline const Expansion* source() const { return src; }
 
    inline Expansion* source() { return src; }
 
    virtual void initDiagonalization() = 0;
 
    virtual bool diagonalizeLayer(size_t layer) = 0;
 
    bool isDiagonalized(size_t layer) {
        return diagonalized[layer];
    }
 
    virtual const cdiagonal& Gamma(size_t layer) const = 0;
 
    virtual const cmatrix& TE(size_t layer) const = 0;
 
    virtual const cmatrix& TH(size_t layer) const = 0;
 
    virtual const cmatrix& invTE(size_t layer) const = 0;
 
    virtual const cmatrix& invTH(size_t layer) const = 0;
 
    // Diagonalization function to compute the smallest eigenvalue of the provided matrix
//     virtual dcomplex smallest_eigevalue(const cmatrix& M) = 0;
};
 
 
class SimpleDiagonalizer : public Diagonalizer
{
  protected:
    std::vector<cdiagonal> gamma;       
    std::vector<cmatrix> Te, Th;        
    std::vector<cmatrix> Te1, Th1;      
 
    void sqrtGamma(cdiagonal& gam) {
        const size_t N = src->matrixSize();
        for (std::size_t j = 0; j < N; j++) {
            dcomplex g = sqrt(gam[j]);
            if (g == 0.) g = SMALL; // Ugly hack to avoid singularity!
            if (real(g) < -SMALL) g = -g;
            if (imag(g) > SMALL) g = -g;
            gam[j] = g;
        }
    }
 
  public:
    SimpleDiagonalizer(Expansion* g);
    ~SimpleDiagonalizer() {}
 
    std::size_t matrixSize() const override;
 
    void initDiagonalization() override;
 
    bool diagonalizeLayer(size_t layer) override;
 
    // Functions returning references to calculated matrices
    const cdiagonal& Gamma(size_t layer) const override { return gamma[layer]; }
    const cmatrix& TE(size_t layer) const override { return Te[layer]; }
    const cmatrix& TH(size_t layer) const override { return Th[layer]; }
    const cmatrix& invTE(size_t layer) const override { return Te1[layer]; }
    const cmatrix& invTH(size_t layer) const override { return Th1[layer]; }
};
 
}}} // namespace plask::optical::modal
#endif // PLASK__SOLVER_OPTICAL_MODAL_DIAGONALIZER_H