/dbfs1/fh22/linux/src/FMMLIB/GratingLayer1D.cpp

Go to the documentation of this file.

#include <fstream>
#include <iomanip>
#include "tools.h"
#include "GratingLayer1D.h"


#ifndef STD
#define STD
using namespace std;
#endif
using TBCI::F_Matrix;
using TBCI::F_BandMatrix;
using TBCI::Vector;


GratingLayer1D::GratingLayer(int truncate, int segments){
        nOrds = truncate; maxOrd=(nOrds-1)/2; minOrd=-maxOrd;
        nPermSegs = segments;
        minPermCoeff = -(nOrds-1);
        maxPermCoeff = nOrds-1;
        nPermCoeffs = maxPermCoeff-minPermCoeff+1;
        permCoeff.resize(2*nOrds-1);
        permCoeffInv.resize(2*nOrds-1);
        eigVal.resize(nOrds);
        alpha.resize(nOrds);
        transPt.resize(segments);
        refInd.resize(segments);
        gamma.resize(nOrds);
        I = complex<double>(0.,1.);
        B.resize(nOrds, nOrds);
        epsMat.resize(nOrds, nOrds);
        invEpsMat.resize(nOrds, nOrds);
        epsMat_inv.resize(nOrds, nOrds);
        invEpsMat_inv.resize(nOrds, nOrds);
        eigMat.resize(nOrds, nOrds);
        polMode="TM";
        eigVecE.resize(nOrds, nOrds);
        eigVecH.resize(nOrds, nOrds);
        invEigVecE.resize(nOrds, nOrds);
        invEigVecH.resize(nOrds, nOrds);
}

GratingLayer1D::~GratingLayer()
{}


// set and get data members
void GratingLayer1D::setPolMode(string polarisation){polMode=polarisation;}
string GratingLayer1D::getPolMode(void){return polMode;}

void GratingLayer1D::setThickness(double layerThickness) {thickness=layerThickness;}
double GratingLayer1D::getThickness() {return thickness;}
void GratingLayer1D::setLambda(double wavelength){lambda=wavelength;}
double GratingLayer1D::getLambda(){return lambda;}

void GratingLayer1D::setN0(complex<double> n){n0=n;}
complex<double> GratingLayer1D::getN0(void){return n0;}
void GratingLayer1D::setPeriod(double per){period = per;}
double GratingLayer1D::getPeriod(void){return period;}
int GratingLayer1D::getnOrds(void){return nOrds;}

void GratingLayer1D::setnPermSegs(int segments) 
{
        nPermSegs = segments;
        refInd.resize(segments);
        transPt.resize(segments);
}
int GratingLayer1D::getnPermSegs() {return nPermSegs;}
void GratingLayer1D::setTransPt(Vector<double> point) {for (int i=0; i<nPermSegs; i++){transPt[i]=point[i];}}
double GratingLayer1D::getTransPt(int comp) {return transPt[comp];}

void GratingLayer1D::setAngle(double angle){theta = angle*M_PI/180.;}

void GratingLayer1D::setRefInd(Vector < complex<double> > index){for (int i = 0; i<nPermSegs; i++){refInd[i]=index[i];}}
complex<double> GratingLayer1D::getRefInd(int comp){return refInd[comp];}

complex<double> GratingLayer1D::getPermCoeff(int order){return permCoeff[order+maxPermCoeff];}
complex<double> GratingLayer1D::getAlpha(int order){return alpha[order];}
complex<double> GratingLayer1D::getEigVal(int order){return eigVal[order];}
complex<double> GratingLayer1D::getEigMat(int row, int column){return eigMat(row,column);}
complex<double> GratingLayer1D::getGamma(int order){return gamma[order];}
complex<double> GratingLayer1D::getEigVecE(int vec, int comp){return eigVecE(vec,comp);}
complex<double> GratingLayer1D::getEigVecH(int vec, int comp){return eigVecH(vec,comp);}
complex<double> GratingLayer1D::getInvEigVecE(int vec, int comp){return invEigVecE(vec,comp);}
complex<double> GratingLayer1D::getInvEigVecH(int vec, int comp){return invEigVecH(vec,comp);}

void GratingLayer1D::setParameters(double wavelength, double angle, complex<double>ind0, double period, double height, Vector <double> points, Vector < complex<double> > indices, string pol){
        setLambda(wavelength);
        setPeriod(period);
        setThickness(height);
        setTransPt(points);
        setN0(ind0);
        setAngle(angle);
        setRefInd(indices);
        setPolMode(pol);
}

void GratingLayer1D::solve(string method="T")
{
        CalculatePermitivity(); 
        createEigMat();
        solveEigMat();
        if (method=="S")
        {
                calcSMatrix();
        }
}

void GratingLayer1D::CalculatePermitivity(void)
{
        // p is the actual index of the coefficient: -M..M is saved in array 0 .. 2*M+1
        int p = minPermCoeff;
        complex<double> prev, curr;
        for(int i=0; i<nPermCoeffs; i++)
        {
                permCoeff[i] = complex<double>(0., 0.);
                permCoeffInv[i] = complex<double>(0., 0.);
                prev = complex<double>((p!=0 ? 1. : 0.), 0.);
                // numerical integration to calculate epsilon[i]: sum of integrals with constant epsilon
                for(int j=0; j<nPermSegs; j++){
                        if (p!=0) {curr = 1.*exp(-I*2.*M_PI*(1.*p)*transPt[j]) ;}
                        else {curr = complex<double>(transPt[j], 0.);}
                        
                        permCoeff[i] += pow(refInd[j],2) * (curr -  prev);
                        if (polMode == "TM")  permCoeffInv[i] += (1./pow(refInd[j],2.)) * (curr -  prev);
                        prev = curr;
                }
                if(p!=0){
                        permCoeff[i] *= complex<double>(0., 1./(2.*M_PI*p));
                        if (polMode == "TM") permCoeffInv[i] *= complex<double>(0., 1./(2.*M_PI*p));
                }
                p++;
        }
}



void GratingLayer1D::createEigMat(void)
{
        int  l, m;
        complex<double> z;
        double ksq = pow((2.*M_PI/lambda), 2);
        
        // calculate gamma[i] (Ben Layet, p.14)
        gamma0=n0*2.*M_PI/lambda*sin(theta);            
        for(int i=0;i<nOrds;i++) {gamma[i] = gamma0 + 2*M_PI*(i+minOrd)/period;}

        // calculate the eigenmatrix elements for TE case (Ben Layet, eq 2.43)
        if (polMode == "TE")
        {       
                l= minOrd;
                for(int i=0;i<nOrds;i++){
                        m=minOrd;
                        for(int j=0;j<nOrds;j++)
                        { 
                                z=permCoeff[nOrds-1+l-m] * ksq;
                                if(l==m) {z -= pow(gamma[m+maxOrd],2);}
                                eigMat(i,j) = z;
                                if (norm(eigMat(i,j))<1e-12) eigMat(i,j)=0.;
                                m++;
                        }
                        l++;
                }
        }
        
        // calculate the eigenmatrix elements for TM case (Granet, Guizal, 1995, eq 13)
        else if (polMode == "TM")
        {
                for (int i=0;i<nOrds;i++)
                {
                        for (int j=0;j<nOrds;j++)
                        {
                                epsMat(i,j)=permCoeff[nOrds-1+i-j];
                                invEpsMat(i,j)=permCoeffInv[nOrds-1+i-j];
                        }
                }
                epsMat_inv=inv(epsMat);
                invEpsMat_inv=inv(invEpsMat);
                
                for (int i=0;i<nOrds;i++)
                {
                        for (int j=0;j<nOrds;j++)
                        {
                                B(i,j)=complex<double>(((i==j) ? ksq : 0.),0.) - epsMat_inv(i,j)*gamma[i]*gamma[j];
                        }
                }
                eigMat=invEpsMat_inv*B;
        }
}


void GratingLayer1D::solveEigMat(void){
        tiny=1e-12;
        for (int i=0;i<nOrds;i++)
        {
                for (int j=0;j<nOrds;j++)
                {
                        if (fabs(real(eigMat(i,j)))<tiny) {eigMat(i,j)=complex<double>(0.,imag(eigMat(i,j)));}
                        if (fabs(imag(eigMat(i,j)))<tiny) {eigMat(i,j)=complex<double>(real(eigMat(i,j)), 0.); }
                }
        }
        if (eig(eigMat,eigVal,eigVecE)!=0) cout<<"Problems with Eigenvector calculation\n";;
        /*Use eigenvalues to set alpha*/
        for(int i=0;i<nOrds;i++)
        {
                alpha[i] = sqrt(eigVal[i]);
                alpha[i] = complex<double>(fabs(real(alpha[i])),fabs(imag(alpha[i])));
        }
}

void GratingLayer1D::calcSMatrix(void)
{
        complex<double> z;
        // eigVecE contains the eigenVectors as columns eigVec[component][vector]
        if (polMode=="TE")
        {
                for (int i=0;i<nOrds; i++)
                {
                        eigVecH.set_col(alpha[i]*eigVecE.get_col(i),i);
                }
        }
        if (polMode=="TM")
        {
                eigVecH=invEpsMat*eigVecE;
                for (int i=0;i<nOrds;i++)
                {
                                eigVecH.set_col(eigVecH.get_col(i)*alpha[i],i);
                }
        }
        // calculate the inverse of the eigenvector matrices.
        // Matrices consist of orthogonal columns, so 
        // inv(M)=conjugate(transpose(M))
        double normal;
        for (int i=0;i<nOrds;i++)
        {
                normal=0.;
                for (int j=0;j<nOrds;j++) normal+=norm(eigVecH(j,i));
                for (int j=0;j<nOrds;j++)
                {
                        invEigVecE(i,j)=conj(eigVecE(j,i));
                        invEigVecH(i,j)=conj(eigVecH(j,i))/normal;
                }
        }       
}


void GratingLayer1D::writeIndex(int orders=-1, string filename="index.dat")
{
        if (orders==-1) orders=nOrds;
        int size = 1000;
        Vector < complex<double> > index(size);
        if (orders%2==0) orders++;
        ofstream out;
        out.open("filename");
        for (int i=0; i<size;i++)
        {
                index[i]=0.;
                for (int p=-(orders-1)/2;p<=(orders-1)/2; p++)
                {
                        index[i]+=getPermCoeff(p)*exp(I*(2*M_PI*p*i/size));
                }
                out<<index[i]<<"\n";
        }
        out.close();
}

---