/dbfs1/fh22/linux/src/FMMLIB/GratingLayer2D.h

Go to the documentation of this file.

#pragma once

//   Copyright (C) 2007 by Florian Hudelist
//  fh22@.hw.ac.uk
#ifndef GRATINGLAYER2D_H
#define GRATINGLAYER2D_H

#include <math.h>
#include <complex>
#include <fstream>
#include <string>

#include "tbci/std_cplx.h"
#include "tbci/vector.h"
#include "tbci/matrix.h"
#include "tbci/f_matrix.h"
#include "tbci/solver/lu_solver.h"
#include "tbci/lapack_stdcplx.h"


using std::complex;
using std::ofstream;
using std::string;

using TBCI::F_Matrix;
using TBCI::Matrix;
using TBCI::Vector;
using TBCI::lu_invert;


class GratingLayer2D
{
private:
        double mu;
        // unneeded arrays are deleted if memsave = true
        // for iterative algorithm this might decrease the performance as memory
        // has to be allocated and delocated in each iteration step
        bool memorysave;
        // grating parameters
        Matrix< complex<double> > refInd;
        double periodX;
        double periodY;
        double thickness;
        int nPixelsX, nPixelsY;
        complex<double> n1;
        Vector<double> transPtX;
        Vector<double> transPtY;
        // incident wave parameters
        double lambda;
        double theta;
        double phi;
        double k0;
        
        // for calculation of characteristic field
        // nOrds=number of spatial Fourier orders used to calculate the field
        // nPermCoeffs= number of Fourier modes for permitivity (must be twice as much as nOrds)
        int nOrds, nOrdsSq, nPermCoeffs, maxPermCoeff;  
        Matrix< complex<double> > epsilon, invEpsilon, invEpsilonMatrix, invEpsilonMatrix_inv;
        Vector< complex<double> > alpha, beta;
        // matrices of [3]eq5-9: t-> top closed bracket (eq6), b-> bottom closed bracket (eq7)
        // bt -> bottom-top brackets (eq8) tb -> top-bottom brackets (eq9)
        // inv before (t/b) indicates the inverse of the Fourier matrix, inv after (t/b/tb/bt) 
        // indicates that the Fourier expansion of the reciprocal of a function (1/f(x)) is calculated
        // invEps is the inverse of the epsilon matrix eps (eq5)
        Matrix< complex<double> > eps, invEps, tInvEps, bInvEps, btEps, tbEps;
        // products of matrices which occur several times: A_B -> A*B
        Matrix < complex<double> > invEps_TbEps, invEps_BtEps;
        Vector < complex<double> > eigVal, gamma;
        Matrix< complex<double> > eigMatE, F, G, FG;

public:
        void calcEig(void);


        Matrix < complex<double> > eigVecE, eigVecH, invEigVecE, invEigVecH;
        // initiate the layer
        GratingLayer2D();
        
        virtual ~GratingLayer2D();

        const GratingLayer2D &operator=(GratingLayer2D &);

        void minimiseMemory();
        // set required parameters needed for the calculations
        // 1: the number of Fourier orders in x and y (total number is then nOrds^2)
        // 2: pixel size in x [my]
        // 3: pixel size in y [my]
        // 4 :thickness of the layer [my]
        // 5: incident wavelength [my]
        // 6: incident angle in respect to the z axis: theta [degree] ([3], Fig1)
        // 7: incident angle in respect to the x axis: phi   [degree] ([3], Fig1)
        // 8: refractive index in region 0
        void setParameters(int, double, double, double, double, double, double, complex<double>);

        void setTransPt(Vector<double>, Vector<double>);
        Vector<double> getTransPtX();
        Vector<double> getTransPtY();
        double getTransPtX(int);
        double getTransPtY(int);
        void reset(Vector<double>, Vector<double>, Matrix< complex<double> >);
        void memsave(bool b){memorysave=b;}
        bool memsave(void){return memorysave;}

        // grating size in x-direction [my m]
        double getSizeX(void);          void setSizeX(double);
        // grating size in y-direction [my m]
        double getSizeY(void);          void setSizeY(double);
        // the thickness of the layer in my m
        double getThickness(void);      void setThickness(double);
        // incident wavelength in vacuum [my m]
        double getLambda(void);         void setLambda(double);
        // angle to z-axis (incident angle) [degree]
        double getTheta(void);          void setTheta(double);
        // angle to x-axis [degree]
        double getPhi(void);            void setPhi(double);
        // refractive index of region 1 (incident region)
        complex<double>getN1(void);     void setN1(complex<double>);
        complex<double>getRefInd(int, int);     
        Matrix< complex<double> >getRefInd(void);
        void setRefInd(Matrix< complex<double> >);
        // number of pixels in x-direction
        int getnPixelsX(void);
        // number of pixels in y-direction
        int getnPixelsY(void);

        int getnOrds(void);
        void setnOrds(int);
        
        // get an element of the eigenvector matrices:
        // the eigenvectors are the columns of the matrix
        // 1: row of the matrix (column)
        // 2: column of the matrix (vector)
        complex<double> getEigMatE(int, int);
        complex<double> getEigVecE(int, int);
        complex<double> getEigVecH(int, int);
        complex<double> getInvEigVecE(int, int);
        complex<double> getInvEigVecH(int, int);
        complex<double> getG(int, int);
        // returns the n_th eigenvalue
        complex<double> getEigVal(int);
        complex<double> getAlpha(int);
        complex<double> getBeta(int);
        complex<double> getGamma(int);
        
        // calculate alpha and beta ([3] equation 11: gamma_x=alpha, gamma_y=beta)
        void calcAlphaBeta(void);
        // calculate the epsilon matrices needed for the eigenvalue equation
        // [3]: equations 5, 8, 9
        // also calculated: [[epsilon]]^(-1) which is needed for the matrix F 
        // [3]: equation 33
        virtual void calcPermCoeffs(void);
        void calcEpsilon(void);
        void calcTBEpsilon(void);
        void calcBTEpsilon(void);
        // set the eigenvalue equation to calculate the eigenmodes of the electric field
        // according to [3]eq.35; by calculating the matrix product by hand and setting zeta=0
        // several factors will cancel out
        void setEigMat(void);
        // solve the eigenmatrix obtained by setEigMat. The eigenvectors are stored
        // in eigVecE and the eigenvalues in eigVal
        // the magnetic part of the eigenvectors (eigVecH) are calculated with [3]eq.36
        void solveEigMat(void);
        // calculates the solution analytically for a homogeneous layer
        void solveHomogeneous(void);
        // Performs : calcPermCoeffs(); setEigMat(); solveEigMat() if the layer is not homogeneous
        // calculates the eigenvectors and eigenvalues analytically if the layer is homogeneous
        // this is necessary to make the algorithm numerically stable (in the homogeneous case, rounding errors can
        // have a huge impact on the result)
        virtual void solve(void);
        
};


#endif

---