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

Go to the documentation of this file.

//   Copyright (C) 2007 by Florian Hudelist
//  fh22@.hw.ac.uk
//
//References: 
//===========
// S-MATRIX APPOACH
// [1] Li, JOSA A, Vol.13, No.5, 1996 - Formulation and comparison of two recursive matrix algorithms for modeling diffraction gratings
// [2] Li, JOSA A, Vol.20, No.4, 2003 - Note on the S-matrix algorithm
//
// FMM: CROSSED GRATINGS
// [3] Li, JOSA A, Vol.14, No.10, 1997 - New formulation of the Fourier modal method for crossed surface-relief gratings 
//

#include "GratingLayer2D.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include "tools.h"
using std::endl;
using std::ofstream;
using std::fstream;
using std::cout;

extern "C" void zgeev_(char *, char *, int *, cmplx*, int *, cmplx*, cmplx*, int*, cmplx*, int*, cmplx*, int*, double*, int*);


GratingLayer2D::GratingLayer2D(void)
{
        mu=1.;
        memorysave=true;
}

GratingLayer2D::~GratingLayer2D(){}

void GratingLayer2D::minimiseMemory()
{
        epsilon.resize(0);
        invEpsilon.resize(0);
        invEpsilonMatrix.resize(0);
        invEpsilonMatrix_inv.resize(0);
        eps.resize(0);
        invEps.resize(0);
        tInvEps.resize(0);
        bInvEps.resize(0);
        btEps.resize(0);
        tbEps.resize(0);
        invEps_TbEps.resize(0);
        invEps_BtEps.resize(0);
        eigVal.resize(0);
        gamma.resize(0);
        eigMatE.resize(0);
        F.resize(0);
        G.resize(0);
        FG.resize(0);
        eigVecE.resize(0);
        eigVecH.resize(0);
        invEigVecE.resize(0);
        invEigVecH.resize(0);
}


void GratingLayer2D::setParameters(int ords, double pX, double pY, double thick, double lam, double angleZ, double angleX, complex<double> ref0)
{
        setnOrds(ords);
        periodX=pX;
        periodY=pY;
        lambda=lam;
        theta=M_PI/180.*angleZ;
        phi=M_PI/180.*angleX;
        n1=ref0;
        k0=2.*M_PI*real(n1)/lambda;
        thickness = thick;
        calcAlphaBeta();
}

double GratingLayer2D::getSizeX(void){return periodX;}          
void GratingLayer2D::setSizeX(double d){periodX=d;}
double GratingLayer2D::getSizeY(void){return periodY;}          
void GratingLayer2D::setSizeY(double d){periodY=d;}
double GratingLayer2D::getThickness(void){return thickness;}    
void GratingLayer2D::setThickness(double d){thickness=d;}
double GratingLayer2D::getLambda(void){return lambda;}          

void GratingLayer2D::setLambda(double wavelength)
{
        lambda=wavelength;
        k0=2.*M_PI*real(n1)/lambda;
}
double GratingLayer2D::getTheta(void){return theta*180./M_PI;}  
void GratingLayer2D::setTheta(double angle)
{
        theta=M_PI/180.*angle;
}
double GratingLayer2D::getPhi(void){return phi*180./M_PI;}      
void GratingLayer2D::setPhi(double angle)
{
        phi=M_PI/180.*angle;
}
complex<double>GratingLayer2D::getN1(void){return n1;}          
void GratingLayer2D::setN1(complex<double>index)
{
        n1=index;
        k0=2.*M_PI*real(n1)/lambda;
}

complex<double>GratingLayer2D::getRefInd(int i, int j){return refInd(i,j);}
Matrix< complex<double> > GratingLayer2D::getRefInd(void){return refInd;}
void GratingLayer2D::setRefInd(Matrix < complex<double> >ref)
{
        nPixelsX=ref.columns();
        nPixelsY=ref.rows();
        refInd.resize(ref.rows(),ref.columns());
        refInd=ref;
}


void GratingLayer2D::setTransPt(Vector<double> tX, Vector<double> tY)
{
        //need to add 1 element to set transition point 0 to 0 (is needed later in the calculations)
        nPixelsX=tX.size();
        nPixelsY=tY.size();
        transPtX.resize(0.,tX.size()+1);
        transPtY.resize(0.,tY.size()+1);
        for (unsigned int i=0;i<transPtX.size()-1;i++)  transPtX[i+1]=tX[i];
        for (unsigned int i=0;i<transPtY.size()-1;i++)  transPtY[i+1]=tY[i];
}

Vector<double> GratingLayer2D::getTransPtX(void){return transPtX.slice(1,transPtX.size());}
Vector<double> GratingLayer2D::getTransPtY(void){return transPtY.slice(1,transPtY.size());}
double GratingLayer2D::getTransPtX(int n){return transPtX[n+1];}
double GratingLayer2D::getTransPtY(int n){return transPtY[n+1];}

void GratingLayer2D::reset(Vector<double> pX, Vector<double> pY, Matrix< complex<double> > N)
{
        nPixelsX=pX.size();
        nPixelsY=pY.size();
        setTransPt(pX,pY);
        setRefInd(N);
}



int GratingLayer2D::getnPixelsX(void) {return nPixelsX;}
int GratingLayer2D::getnPixelsY(void) {return nPixelsY;}
int GratingLayer2D::getnOrds(void) {return nOrds;}
void GratingLayer2D::setnOrds(int ords)
{
        nOrds=ords;
        nOrdsSq=nOrds*nOrds;
        nPermCoeffs=2*nOrds-1;
        maxPermCoeff=(nPermCoeffs-1)/2;
}

complex<double> GratingLayer2D::getEigMatE(int a, int b){return eigMatE(a,b);}
complex<double> GratingLayer2D::getEigVecE(int a, int b){return eigVecE(a,b);}
complex<double> GratingLayer2D::getEigVecH(int a, int b){return eigVecH(a,b);}
complex<double> GratingLayer2D::getInvEigVecE(int a, int b){return invEigVecE(a,b);}
complex<double> GratingLayer2D::getInvEigVecH(int a, int b){return invEigVecH(a,b);}    
complex<double> GratingLayer2D::getG(int a, int b){return G(a,b);}      
complex<double> GratingLayer2D::getEigVal(int a){return eigVal[a];}
complex<double> GratingLayer2D::getAlpha(int a){return alpha[a];}
complex<double> GratingLayer2D::getBeta(int a){return beta[a];}
complex<double> GratingLayer2D::getGamma(int a){return gamma[a];}


const GratingLayer2D &GratingLayer2D::operator=(GratingLayer2D &in)
{
        //assign all parameters
        setnOrds(in.getnOrds());
        setN1(in.getN1());
        setLambda(in.getLambda());
        setPhi(in.getPhi());
        setTheta(in.getTheta());
        setSizeX(in.getSizeX());
        setSizeY(in.getSizeY());
        setThickness(in.getThickness());
        nPixelsX=in.getnPixelsX();
        nPixelsY=in.getnPixelsY();
        k0=in.k0;
        mu=in.mu;

        // assign all matrices
        refInd.resize(in.refInd);
        transPtX.resize(in.transPtX);
        transPtY.resize(in.transPtY);

        alpha.resize(in.alpha);
        beta.resize(in.beta);
        epsilon.resize(in.epsilon);
        invEpsilon.resize(in.invEpsilon);
        invEpsilonMatrix.resize(in.invEpsilonMatrix);
        invEpsilonMatrix_inv.resize(in.invEpsilonMatrix_inv);
        eps.resize(in.eps);
        invEps.resize(in.invEps);
        tInvEps.resize(in.invEps);
        bInvEps.resize(in.bInvEps);
        btEps.resize(in.btEps);
        tbEps.resize(in.tbEps);
        invEps_TbEps.resize(in.invEps_TbEps);
        invEps_BtEps.resize(in.invEps_BtEps);
        eigVal.resize(in.eigVal);
        gamma.resize(in.gamma);
        eigMatE.resize(in.eigMatE);
        F.resize(in.F);
        G.resize(in.G);
        FG.resize(in.FG);
        eigVecE.resize(in.eigVecE);
        eigVecH.resize(in.eigVecH);
        invEigVecE.resize(in.invEigVecE);
        invEigVecH.resize(in.invEigVecH);

        return *this;
}


// calculate the exponents of the base functions ([3], eq11)
void GratingLayer2D::calcAlphaBeta(void)
{
        if ((signed int)alpha.size()!=nOrds) alpha.resize(nOrds);
        if ((signed int)beta.size()!=nOrds) beta.resize(nOrds);
        for (int i=0;i<nOrds;i++)
        {
                alpha[i]=n1*2.*M_PI/lambda*sin(theta)*cos(phi)+2*M_PI*(i-(nOrds-1)/2)/periodX;
                beta[i]=n1*2.*M_PI/lambda*sin(theta)*sin(phi)+2*M_PI*(i-(nOrds-1)/2)/periodY;
        }
}



void GratingLayer2D::calcPermCoeffs(void)
{
        calcEpsilon();
        calcBTEpsilon();
        calcTBEpsilon();
}

void GratingLayer2D::calcEpsilon(void)
{
        if ((signed int)epsilon.size()!=nPermCoeffs*nPermCoeffs) epsilon.resize(nPermCoeffs, nPermCoeffs);
        if ((signed int)eps.size()!=pow(nOrds,4)) eps.resize(nOrds*nOrds, nOrds*nOrds);         
        if ((signed int)invEps.size()!=pow(nOrds,4)) invEps.resize(nOrds*nOrds, nOrds*nOrds);
        double k2, k1;
        complex<double> cx,cy, c;
        k1=2.*M_PI/periodX;
        k2=2.*M_PI/periodY;

        //calc epsilon_n,m: [3] equation 4
        for (int i=0;i<nPermCoeffs;i++)
        {
                for (int j=0;j<nPermCoeffs;j++)
                {
                        epsilon(i,j)=0.;
                        for (int p=0;p<nPixelsX;p++)
                        {
                                for (int q=0;q<nPixelsY;q++)
                                {
                                        // if i-maxPermCoeff=0 (i.e. zero order coefficient) then the exponential function vanishes
                                        // and in the integral is only a constant
                                        cx = ((i-maxPermCoeff==0) ? (transPtX[p+1]-transPtX[p]) : (exp(-I*k1*(1.*i-maxPermCoeff)*transPtX[p+1])-exp(-I*k1*(1.*i-maxPermCoeff)*(transPtX[p])))*I*(periodX/((i-maxPermCoeff)*k1)));
                                        // if j-maxPermCoeff=0 (i.e. zero order coefficient) then the exponential function vanishes
                                        // and in the integral is only a constant
                                        cy = ((j-maxPermCoeff==0) ? (transPtY[q+1]-transPtY[q]) : (exp(-I*k2*(1.*j-maxPermCoeff)*transPtY[q+1])-exp(-I*k2*(1.*j-maxPermCoeff)*(transPtY[q])))*I*(periodY/((j-maxPermCoeff)*k2)));

                                        epsilon(i,j)+=cx*cy*pow(refInd(q,p),2.);
                                }
                        }
                        epsilon(i,j)/=(periodX*periodY);
                }
        }

         //[3] eq5: [[epsilon]]_(mn,jl)=epsilon[m-j][n-l]
         //[[epsilon]]_(mn,jl)=[[epsilon]][m*nOrds+n][j*nOrds+l]
        for (int m=0;m<nOrds;m++)
        {
                for (int n=0;n<nOrds;n++)
                {
                        for (int j=0;j<nOrds;j++)
                        {
                                for (int l=0;l<nOrds;l++)
                                {
                                        eps(m*nOrds+n,j*nOrds+l)=epsilon(nOrds-1+m-j,nOrds-1+n-l);
                                }
                        }
                }
        }
        invEps=invert(eps);


        if (memorysave)
        {
                epsilon.resize(0);
                eps.resize(0);
        }
}


void GratingLayer2D::calcTBEpsilon(void)
{
        if ((signed int)bInvEps.size()!=nPixelsX*nPermCoeffs) bInvEps.resize(nPixelsX, nPermCoeffs);
        if ((signed int)tbEps.size()!=pow(nOrds,4)) tbEps.resize(0.,nOrds*nOrds, nOrds*nOrds);
        Matrix< complex<double> > M1(nOrds);
        double k2, k1;
        k1=2.*M_PI/periodX;
        k2=2.*M_PI/periodY;
        // bInvEps: [3] equation 7
        for (int i=0;i<nPixelsX;i++)
        {
                for(int j=0;j<nPermCoeffs;j++)
                {
                        bInvEps(i,j)=0.;
                        for (int k=0;k<nPixelsY;k++)
                        {
                                if (j-maxPermCoeff!=0)
                                {
                                        bInvEps(i,j)+=pow(1./refInd(k,i),2.)/(2.*M_PI*(j-maxPermCoeff))*I*(exp(-I*k2*(1.*j-maxPermCoeff)*transPtY[k+1])-exp(-I*k2*(1.*j-maxPermCoeff)*transPtY[k]));
                                }
                                else bInvEps(i,j)+=pow(1./refInd(k,i),2.)/periodY*(transPtY[k+1]-transPtY[k]);
                        }
                }
        }


        // calc tbEps according to Li, [3],eq9
        // M1 is the Toeplitz natrix formed by the coefficients of invTinvEps 
        // [3]eq7 for 1/epsilon to be used in eq9, invM1 is the inverse of M1
        for (int x=0;x<nPixelsX;x++)
        {       
                for (int i=0;i<nOrds;i++)
                {
                        for (int j=0;j<nOrds;j++)
                        {
                                M1(i,j) = bInvEps(x,nOrds-1+i-j);
                        }
                }
                // invM1 is the inverse of bInvEps needed in eq 9
                M1=invert(M1);

                // tbEps[i][k] = tbEps_(m-j , n-l)=matrix(m*nOrds+n , j*nOrds+l)
                for (int m=0;m<nOrds;m++)
                {
                        for (int n=0;n<nOrds;n++)
                        {
                                for (int j=0;j<nOrds;j++)
                                {
                                        for (int l=0;l<nOrds;l++)
                                        {
                                                if (m-j!=0)
                                                {
                                                        tbEps(m*nOrds+n,j*nOrds+l)+= M1(n,l)*I/((1.*m-j)*2.*M_PI)*(exp(-I*k1*(1.*m-j)*transPtX[x+1])-exp(-I*k1*(1.*m-j)*transPtX[x]));
                                                }
                                                else tbEps(m*nOrds+n,j*nOrds+l)+=M1(n,l)/periodX*(transPtX[x+1]-transPtX[x]);
                                        }
                                }
                        }
                }
        }
        if (memorysave)
        {
                bInvEps.resize(0);
        }
}

void GratingLayer2D::calcBTEpsilon(void)
{
        if ((signed int)tInvEps.size()!=nPixelsY*nPermCoeffs) tInvEps.resize(nPixelsY, nPermCoeffs);
        if ((signed int)btEps.size()!=pow(nOrds,4)) btEps.resize(0.,nOrds*nOrds, nOrds*nOrds);
        Matrix< complex<double> > M1(nOrds);
        double k2, k1;
        k1=2.*M_PI/periodX;
        k2=2.*M_PI/periodY;

        // tInvEps: [3] equation 6
        for (int i=0;i<nPixelsY;i++)
        {
                for(int j=0;j<nPermCoeffs;j++)
                {
                        tInvEps(i,j)=0.;
                        for (int k=0;k<nPixelsX;k++)
                        {
                                if (j-maxPermCoeff!=0)
                                {
                                        tInvEps(i,j)+=pow(1./refInd(i,k),2.)/(2.*M_PI*(j-maxPermCoeff))*I*(exp(-I*k1*(1.*j-maxPermCoeff)*transPtX[k+1])-exp(-I*k1*(1.*j-maxPermCoeff)*(transPtX[k])));
                                }
                                else tInvEps(i,j)+=pow(1./refInd(i,k),2.)*(transPtX[k+1]-transPtX[k])/periodX;
                        }
                }
        }
        
        // calc btEps [3], eq8
        // M1 is the Toeplitz natrix formed by the coefficients of invTinvEps 
        // [3]eq6 for 1/epsilon to be used in eq8, invM1 is the inverse of M1
        for (int y=0;y<nPixelsY;y++)
        {       
                for (int i=0;i<nOrds;i++)
                {
                        for (int j=0;j<nOrds;j++)
                        {
                                M1(i,j) = tInvEps(y,nOrds-1+i-j);
                        }
                }

                // invM1 is the inverse of tInvEps needed in [3], eq8
                M1=invert(M1);
                // btEps[i][k] = btEps_(m-j , n-l)=matrix(m*nOrds+n , j*nOrds+l)
                for (int m=0;m<nOrds;m++)
                {
                        for (int n=0;n<nOrds;n++)
                        {
                                for (int j=0;j<nOrds;j++)
                                {
                                        for (int l=0;l<nOrds;l++)
                                        {
                                                if (n-l!=0)
                                                {
                                                        btEps(m*nOrds+n,j*nOrds+l) +=M1(m,j)*I/((1.*n-l)*2.*M_PI)*(exp(-I*k2*(1.*n-l)*transPtY[y+1])-exp(-I*k2*(1.*n-l)*(transPtY[y])));
                                                }
                                                else btEps(m*nOrds+n,j*nOrds+l) +=M1(m,j)*(transPtY[y+1]-transPtY[y])/periodY;
                                        }
                                }
                        }
                }
        }
        if (memorysave)
        {
                tInvEps.resize(0);
        }
}



void GratingLayer2D::solve(void)
{
        int hom=0;
        for (int i=0;i<nPixelsX;i++) for (int j=0;j<nPixelsY;j++) 
        {
                if (refInd(j,i)!=refInd(0,0)) hom++;
        }
        if (hom==0) solveHomogeneous();
        else
        {
                calcPermCoeffs();
                setEigMat();
                solveEigMat();
        }
}

void GratingLayer2D::solveHomogeneous()
{
        if ((signed) eigVal.size()!=2*nOrdsSq) {eigVal.resize(2*nOrdsSq); }
        if ((signed) invEigVecE.size()!=pow(2*nOrdsSq,2)) {invEigVecE.resize(2*nOrdsSq,2*nOrdsSq); }
        if ((signed) invEigVecH.size()!=pow(2*nOrdsSq,2)) invEigVecH.resize(2*nOrdsSq,2*nOrdsSq);
        if ((signed int) eigVecE.size()!=pow(2*nOrdsSq,2)) eigVecE.resize(0.,2*nOrdsSq, 2*nOrdsSq);
        if ((signed int) eigVecH.size()!=pow(2*nOrdsSq,2)) eigVecH.resize(0.,2*nOrdsSq, 2*nOrdsSq);
        if ((signed)gamma.size()!=2*nOrdsSq) gamma.resize(2*nOrdsSq);

        for (int i=0;i<2*nOrdsSq;i++)
        {
                eigVal[i]=sqrt(-alpha[(i%nOrdsSq)/nOrds]*alpha[(i%nOrdsSq)/nOrds]-beta[i%nOrds]*beta[i%nOrds]+(refInd(0,0)*refInd(0,0)*mu*k0*k0));
        }
        /*Use eigenvalues to set gamma*/
        for(int i=0;i<2*nOrdsSq;i++){
                gamma(i) = sqrt(eigVal(i));
                if(real(gamma(i))+imag(gamma(i))<0) gamma(i)*=-1.;
        }
        for (int i=0;i<nOrdsSq;i++)
        {
                eigVecE(i,i)=1.;
                eigVecE(i+nOrdsSq,i+nOrdsSq)=1.;
                invEigVecE(i,i)=1.;
                invEigVecE(i+nOrdsSq,i+nOrdsSq)=1.;
                eigVecH(i,i)=-alpha[i/nOrds]*beta[i%nOrds]/(mu*k0*eigVal[i]);
                eigVecH(i+nOrdsSq,i+nOrdsSq)=(alpha[i/nOrds]*beta[i%nOrds])/(mu*k0*eigVal[i]);
                eigVecH(i,i+nOrdsSq)=(alpha[i/nOrds]*alpha[i/nOrds]-mu*k0*k0*refInd(0,0)*refInd(0,0))/(mu*k0*eigVal[i]);
                eigVecH(i+nOrdsSq,i)=(mu*k0*k0*refInd(0,0)*refInd(0,0)-beta[i%nOrds]*beta[i%nOrds])/(mu*k0*eigVal[i]);
        }
        invEigVecH=invert(eigVecH);
}

void GratingLayer2D::setEigMat(void)
{
        if (G.size()!=pow(2*nOrdsSq,2)) G.resize(2*nOrdsSq,2*nOrdsSq);
        if (eigMatE.size()!=pow(2*nOrdsSq,2)) eigMatE.resize(2*nOrds*nOrds, 2*nOrds*nOrds);
        complex<double> c;
        // set F*G matrix of [3]eq35. For ceta=0 (i.e. rectangular lattice vectors) a few terms cancel out
        // after executing the matrix product
        //cout <<"\nEIGMATE\n\n";
        for (int i=0;i<2*nOrdsSq;i++)
        {
                for (int j=0;j<2*nOrdsSq;j++)
                {
                        eigMatE(i,j)=0.;
                        // set (FG)_(1,1) submatrix
                        if ((i<nOrdsSq) && (j<nOrdsSq))
                        {       
                                eigMatE(i,j) += mu*k0*k0*btEps(i,j);
                                if (i==j) eigMatE(i,j) -= beta[i%nOrds]*beta[i%nOrds];
                                // calculate the (i,j) element of ([alpha]*[eps])*([alpha]*[btEps])
                                c=0.;
                                for (int k=0;k<nOrdsSq;k++)
                                {
                                        c+=alpha[i/nOrds]*invEps(i,k)*alpha[k/nOrds]*btEps(k,j);
                                }
                                eigMatE(i,j) -= c;
                        }
                        //set (FG)_(1,2) submatrix
                        if ((i<nOrdsSq) && (j>=nOrdsSq))
                        {
                                // calculate the (i,j) element of ([alpha]*[eps])*([beta]*[tbEps])
                                c=0.;
                                for (int k=0;k<nOrdsSq;k++)
                                {
                                        c+=alpha[i/nOrds]*invEps(i,k)*beta[k%nOrds]*tbEps(k,j-nOrdsSq);
                                }
                                eigMatE(i,j) = ((i==j-nOrdsSq?alpha[i/nOrds]*beta[i%nOrds]:0.)-c);
                        }
                        
                        //set (FG)_(2,1) submatrix
                        if ((i>=nOrdsSq) && (j<nOrdsSq))
                        {
                                // calculate the (i,j) element of ([alpha]*[eps])*([beta]*[tbEps])
                                c=0.;
                                for (int k=0;k<nOrdsSq;k++)
                                {
                                        c+=beta[i%nOrds]*invEps(i-nOrdsSq,k)*alpha[k/nOrds]*btEps(k,j);
                                }
                                eigMatE(i,j) += -c+(i==j+nOrdsSq?alpha[(i%nOrdsSq)/nOrds]*beta[i%nOrds]:0.);
                        }
                        // set (FG)_(2,2) submatrix
                        if ((i>=nOrdsSq) && (j>=nOrdsSq))
                        {
                                eigMatE(i,j) += mu*k0*k0*tbEps(i-nOrdsSq,j-nOrdsSq);
                                if (i==j) eigMatE(i,j) -= alpha[(i%nOrdsSq)/nOrds]*alpha[(i%nOrdsSq)/nOrds];
                                // calculate the (i,j) element of ([alpha]*[eps])*([alpha]*[btEps])
                                c=0.;
                                for (int k=0;k<nOrdsSq;k++)
                                {
                                        c+=beta[i%nOrds]*invEps(i-nOrdsSq,k)*beta[k%nOrds]*tbEps(k,j-nOrdsSq);
                                }
                                eigMatE(i,j) -= c;
                        }
                        
                        
                        // set G matrix ([3]eq34) to calculate the magnetic eigenvectors ([3]eq36 )
                        // the angle zeta=0, so all terms with sin(ceta) cancel out, cos(zeta)=1
                        // set (G)_(1,1) submatrix
                        if ((i<nOrdsSq) && (j<nOrdsSq))
                        {       
                                G(i,j)=(i==j ? -alpha[i/nOrds]*beta[j%nOrds] : 0.);
                        }
                        //set (G)_(1,2) submatrix
                        if ((i<nOrdsSq) && (j>=nOrdsSq))
                        {
                                G(i,j)= (i==j-nOrdsSq ? alpha[i/nOrds]*alpha[i/nOrds] : 0.) - mu*k0*k0*tbEps(i,j-nOrdsSq);
                        }
                        //set (G)_(2,1) submatrix
                        if ((i>=nOrdsSq) && (j<nOrdsSq))
                        {
                                G(i,j)=mu*k0*k0*btEps(i-nOrdsSq,j) - (i-nOrdsSq==j ? beta[i%nOrds]*beta[i%nOrds] : 0.);
                        }
                        // set (G)_(2,2) submatrix
                        if ((i>=nOrdsSq) && (j>=nOrdsSq))
                        {
                                G(i,j) = (i==j ? alpha[(i%nOrdsSq)/nOrds]*beta[j%nOrds] : 0.);
                        }
                }
        }

        if (memorysave)
        {
        invEps.resize(0);
                btEps.resize(0);
                tbEps.resize(0);
        }
}



void GratingLayer2D::solveEigMat(void)
{

        if ((signed int) eigVecE.size()!=pow(2*nOrdsSq,2)) eigVecE.resize(0.,2*nOrdsSq, 2*nOrdsSq);
        if ((signed int) eigVecH.size()!=pow(2*nOrdsSq,2)) eigVecH.resize(0.,2*nOrdsSq, 2*nOrdsSq);
        if ((signed int) eigVal.size()!=2*nOrdsSq) eigVal.resize(0.,2*nOrdsSq);
        if ((signed int) gamma.size()!=2*nOrdsSq) gamma.resize(2*nOrdsSq);
        F_Matrix< complex<double> >eigVecE_F(eigVecE.rows(), eigVecE.columns());


        TBCI::eig(F_Matrix< complex<double> >(eigMatE),eigVal,eigVecE_F);
        eigVecE=Matrix< complex<double> >(eigVecE_F);


        /*Use eigenvalues to set gamma*/
        for(int i=0;i<2*nOrdsSq;i++){
                gamma(i) = sqrt(eigVal(i));
                if(real(gamma(i))+imag(gamma(i))<0) gamma(i)*=-1.;
        }

        eigVecH=G*eigVecE;
        for (int i=0;i<2*nOrdsSq;i++)
        {
                for (int j=0;j<2*nOrdsSq;j++)
                        eigVecH(i,j)/=(mu*k0*gamma[j]);
        }

        if ((signed int) invEigVecE.size()!=pow(2*nOrdsSq,2)) {invEigVecE.resize(2*nOrdsSq,2*nOrdsSq); }
        if ((signed int) invEigVecH.size()!=pow(2*nOrdsSq,2)) invEigVecH.resize(2*nOrdsSq,2*nOrdsSq);
        invEigVecH=invert(eigVecH);
        invEigVecE=invert(eigVecE);

        if (memorysave)
        {
                eigMatE.resize(0);
                G.resize(0);
        }
}       

void GratingLayer2D::calcEig(void)
{

}

---