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

Go to the documentation of this file.

//   Copyright (C) 2007 by Florian Hudelist
//  fh22@.hw.ac.uk
//
//References: 
//===========
// [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
// [3] Li, JOSA A, Vol.14, No.10, 1997 - New formulation of the Fourier modal method for crossed surface-relief gratings 
//
#include "Grating2D.h"
#include <string>
#include <iostream>
#include <ctime>
#include "tools.h"

#include "tbci/f_bandmatrix.h"
#include "tbci/matrix.h"

using std::cout;
using std::complex;
using std::cin;
using std::endl;
using TBCI::Vector;
using TBCI::F_Matrix;
using TBCI::CSCMatrix;
using TBCI::lu_invert;

bool memorysave;

Grating2D::Grating2D(int layers)
{
        nLayers=layers;
        //GratingLayer2D dummy;
        while ((signed) layer.size()<nLayers) layer.push_back(GratingLayer2D());
        mu=1.;
        setMemsave(true);
}

Grating2D::~Grating2D()
{
}

const Grating2D &Grating2D::operator=(Grating2D &in)
{
        setnLayers(in.getnLayers());
        setnOrds(in.getnOrds());
        setN1(in.getN1());
        setN3(in.getN3());
        setLambda(in.getLambda());
        setPhi(in.getPhi());
        setTheta(in.getTheta());
        setSizeY(in.getSizeY());
        setSizeX(in.getSizeX());
        incid.resize(in.incid);

        setMemsave(in.getMemsave());

        R.resize(in.R);
        T.resize(in.T);
        Q12.resize(in.Q12);
        Q2.resize(in.Q2);
        tau.resize(in.tau);
        T_tilde.resize(in.T_tilde);
        R_tilde.resize(in.R_tilde);
        G.resize(in.G);
        W21.resize(in.G);
        W23.resize(in.W23);
        alpha_1.resize(in.alpha_1);
        alpha_3.resize(in.alpha_3);
        beta_1.resize(in.beta_1);
        beta_3.resize(in.beta_3);
        gamma_1.resize(in.gamma_1);
        gamma_3.resize(in.gamma_3);
        k1=in.k1;
        k3=in.k3;

        trans.resize(in.trans);
        ref.resize(in.ref);
        efficiencyR.resize(in.efficiencyR);
        efficiencyT.resize(in.efficiencyT);

        for (int l=0;l<getnLayers();l++) layer[l]=in.layer[l];

        return *this;
}


void Grating2D::minimiseMemory(void)
{
        Q12.resize(0);
        Q2.resize(0);
        tau.resize(0);
        T_tilde.resize(0);
        R_tilde.resize(0);
        G.resize(0);
        W21.resize(0);
        W23.resize(0);
        for (int l=0;l<nLayers;l++) layer[l].minimiseMemory();
}
void Grating2D::setParameters(int ords, double pX, double pY, Vector<double> thick, double lam, double angleZ, double angleX, complex<double> ref1, complex<double> ref3)
{
        setnOrds(ords);
        setSizeX(pX);
        setSizeY(pY);
        setLambda(lam);
        setTheta(angleZ);
        setPhi(angleX);
        setN1(ref1);
        setN3(ref3);
        for (int i=0;i<nLayers;i++)
        {
                layer[i].setParameters(nOrds, periodX, periodY, thick[i], lambda, theta, phi, n1);
        }
        calcAlphaBetaGamma();
}

void Grating2D::setTransPt(Vector<double> tX, Vector<double> tY)
{
        for (int i=0;i<nLayers;i++) layer[i].setTransPt(tX,tY);
}

void Grating2D::setTransPt()
{
        double pX,pY;
        Vector<double> tX, tY;
        for (int l=0;l<nLayers;l++) 
        {
                tX.resize(layer[l].getnPixelsX());
                tY.resize(layer[l].getnPixelsY());
                // calculate pixelsizes
                pX=periodX/tX.size();
                pY=periodY/tY.size();
                for (unsigned int i=0;i<tX.size();i++) tX[i]=pX*(i+1);
                for (unsigned int i=0;i<tY.size();i++) tY[i]=pY*(i+1);
                layer[l].setTransPt(tX,tY);
        }
}


void Grating2D::setRefInd(vector< Matrix< complex<double> > > N)
{
        for (int i=0;i<nLayers;i++) layer[i].setRefInd(N[i]);
}



void Grating2D::setMemsave(bool b)
{
        memorysave=b;
        for (int i=0;i<nLayers;i++)
        {
                layer[i].memsave(b);
        }
}
bool Grating2D::getMemsave(void) {return memorysave;}


double Grating2D::getSizeX(void){return periodX;}
void Grating2D::setSizeX(double d)
{
        periodX=d;
        for (int l=0;l<nLayers;l++) layer[l].setSizeX(periodX);
}
double Grating2D::getSizeY(void){return periodY;}               
void Grating2D::setSizeY(double d)
{
        periodY=d;
        for (int l=0;l<nLayers;l++) layer[l].setSizeY(periodY);
}

int Grating2D::getnLayers(void){return nLayers;}

void Grating2D::setnLayers(int l)
{
        nLayers=layer.size();
        // delete the last elements from the layer vector
        if (l<nLayers) while ((signed) layer.size()>l) layer.pop_back();
        // append elements to the layer vector
        if (l>nLayers) while ((signed) layer.size()<l) layer.push_back(GratingLayer2D());
        nLayers=l;
}

double Grating2D::getLambda(void){return lambda;}               
void Grating2D::setLambda(double d)
{
        lambda=d;
        k1=2.*real(n1)*M_PI/lambda;
        k3=2.*real(n3)*M_PI/lambda;
        for (unsigned int l=0;l<layer.size();l++) layer[l].setLambda(lambda);
}
double Grating2D::getTheta(void){return theta*180./M_PI;}       
void Grating2D::setTheta(double d)
{
        theta=M_PI/180.*d;
        for (unsigned int l=0;l<layer.size();l++) layer[l].setTheta(d);
}

double Grating2D::getPhi(void){return phi*180./M_PI;}           
void Grating2D::setPhi(double d)
{
        phi=M_PI/180.*d;
        for (unsigned int l=0;l<layer.size();l++) layer[l].setPhi(d);
}
complex<double>Grating2D::getN1(void){return n1;}               
void Grating2D::setN1(complex<double>c)
{
        n1=c;
        k1=2.*real(n1)*M_PI/lambda;
        for (unsigned int l=0;l<layer.size();l++) layer[l].setN1(c);
}
complex<double>Grating2D::getN3(void){return n3;}               
void Grating2D::setN3(complex<double>c)
{
        n3=c;
        k3=2.*real(n3)*M_PI/lambda;
}

int Grating2D::getnOrds(void) {return nOrds;}
void Grating2D::setnOrds(int ords)
{
        if (ords%2==0)
        {
                cout<<"Attempt to set nOrds to an even number ("<<ords<<"), an odd number is required, set nOrds to "<<ords+1<<endl;
                nOrds=ords+1;
        }
        nOrds=ords;
        nOrdsSq=nOrds*nOrds;
        for (unsigned int i=0;i<layer.size();i++) layer[i].setnOrds(nOrds);
}

complex<double> Grating2D::getReflection(string pol, int row, int column)
{
        if ((pol=="X")||(pol=="x")) return ref[row*nOrds+column];
        if ((pol=="Y")||(pol=="y")) return ref[row*nOrds+column+nOrdsSq];
        else 
        {
                //cout<<"Grating2D::getReflection(string pol, int row, int column): polarisation not recognised\n";
                return 0.;
        }
}

complex<double> Grating2D::getTransmission(string pol, int row, int column)
{
        if ((pol=="X")||(pol=="x")) return trans[row*nOrds+column];
        if ((pol=="Y")||(pol=="y")) return trans[row*nOrds+column+nOrdsSq];
        else 
        {
                //cout<<"Grating2D::getTransmission(string pol, int row, int column): polarisation not recognised\n";
                return 0.;
        }
}

complex<double> Grating2D::getIncid(string pol, int row, int column)
{
        if ((pol=="X")||(pol=="x")) return incid[row*nOrds+column];
        if ((pol=="Y")||(pol=="y")) return incid[row*nOrds+column+nOrdsSq];
        else return 0.;
}

complex<double> Grating2D::getR(int a, int b){return R(a,b);}
complex<double> Grating2D::getT(int a, int b){return T(a,b);}

complex<double> Grating2D::getAlphaT(int n){return alpha_3[n];}
complex<double> Grating2D::getBetaT(int n){return beta_3[n];}
complex<double> Grating2D::getAlphaR(int n){return alpha_1[n];}
complex<double> Grating2D::getBetaR(int n){return beta_1[n];}


void Grating2D::setIncid(Matrix< complex<double> >x, Matrix< complex<double> >y)
{
        if (signed(incid.size())!=2*nOrdsSq) incid.resize(2*nOrdsSq);
        complex<double> normalise;
        for (int i=0;i<nOrds;i++)
        {
                for (int j=0;j<nOrds;j++)
                {
                        incid[i*nOrds+j]=x(i,j);
                        incid[i*nOrds+j+nOrdsSq]=y(i,j);
                }
        }
        // normalise the incident coefficients sao that sum(incid^2)=1
        // according to [3] equation 20
        normalise = 1./gamma_1[(nOrds-1)/2*nOrds+(nOrds-1)/2]*((k1-beta_1[(nOrds-1)/2])*norm(incid[(nOrds-1)/2*nOrds+(nOrds-1)/2])+(k1-alpha_1[(nOrds-1)/2])*norm(incid[(nOrds-1)/2*nOrds+(nOrds-1)/2+nOrdsSq])+alpha_1[(nOrds-1)/2]*beta_1[(nOrds-1)/2]*(incid[(nOrds-1)/2*nOrds+(nOrds-1)/2]*conj(incid[(nOrds-1)/2*nOrds+(nOrds-1)/2+nOrdsSq])+incid[(nOrds-1)/2*nOrds+(nOrds-1)/2+nOrdsSq]*conj(incid[(nOrds-1)/2*nOrds+(nOrds-1)/2])));
        
        incid[(nOrds-1)/2*nOrds+(nOrds-1)/2]/=sqrt(normalise);
        incid[(nOrds-1)/2*nOrds+(nOrds-1)/2+nOrdsSq]/=sqrt(normalise);
}

Matrix< complex<double> > Grating2D::getIncid(string pol)
{
        Matrix< complex<double> > in(nOrds,nOrds);
        for (int i=0;i<nOrds;i++)
        {
                for (int j=0;j<nOrds;j++)
                {
                        in(i,j)=incid[i*nOrds+j+(pol=="y"?nOrdsSq:0)];
                }
        }
        return in;
}


void Grating2D::solve()
{
        for (int l=0;l<nLayers;l++) layer[l].solve();
        setSMatrix();
        calcEfficiencyR();
        calcEfficiencyT();
}


void Grating2D::calcAlphaBetaGamma(void)
{
        if ((signed int)alpha_1.size()!=nOrds) alpha_1.resize(0.,nOrds);
        if ((signed int)alpha_3.size()!=nOrds) alpha_3.resize(0.,nOrds);
        if ((signed int)beta_1.size()!=nOrds) beta_1.resize(0.,nOrds);
        if ((signed int)beta_3.size()!=nOrds) beta_3.resize(0.,nOrds);
        if ((signed int)gamma_1.size()!=nOrdsSq) gamma_1.resize(0.,nOrdsSq);
        if ((signed int)gamma_3.size()!=nOrdsSq) gamma_3.resize(0.,nOrdsSq);
        
        for (int i=0;i<nOrds;i++)
        {
                alpha_3[i]=2.*M_PI/lambda*n3*sin(theta)*cos(phi)+(i-(nOrds-1)/2)*2.*M_PI/periodX;
                alpha_1[i]=2.*M_PI/lambda*n3*sin(theta)*cos(phi)+(i-(nOrds-1)/2)*2.*M_PI/periodX;
                beta_3[i]=2.*M_PI/lambda*n3*sin(theta)*sin(phi)+(i-(nOrds-1)/2)*2.*M_PI/periodY;
                beta_1[i]=2.*M_PI/lambda*n3*sin(theta)*sin(phi)+(i-(nOrds-1)/2)*2.*M_PI/periodY;
        }

        for (int i=0;i<nOrdsSq;i++)
        {
                        gamma_3[i]=sqrt(pow(2.*M_PI/lambda*n3,2)-alpha_3[i/nOrds]*alpha_3[i/nOrds]-beta_3[i%nOrds]*beta_3[i%nOrds]);
                        gamma_1[i]=sqrt(pow(2.*M_PI/lambda*n1,2)-alpha_1[i/nOrds]*alpha_1[i/nOrds]-beta_1[i%nOrds]*beta_1[i%nOrds]);
                        if ( (real(gamma_1[i])+imag(gamma_1[i]))<0 ) gamma_1[i]*=-1.;
                        if ( (real(gamma_3[i])+imag(gamma_3[i]))<0 ) gamma_3[i]*=-1.;
        }
        
}

void Grating2D::calcBoundaryW(void)
{
        complex<double> eps=n3*n3;
        //calculate the magnetic part of the eigenvectors with [3]eq.36
        //each column represents one eigenvector
        W23.resize(0.,2*nOrdsSq, 2*nOrdsSq);
        for (int i=0;i<nOrdsSq;i++)
        {
                W23(i,i)=-alpha_3[(i%nOrdsSq)/nOrds]*beta_3[i%nOrds]/(sqrt(mu)*k1*gamma_3[i]);
                W23(i+nOrdsSq,i+nOrdsSq)=(alpha_3[(i%nOrdsSq)/nOrds]*beta_3[i%nOrds])/(sqrt(mu)*k1*gamma_3[i]);
                W23(i,i+nOrdsSq)=(alpha_3[i/nOrds]*alpha_3[i/nOrds]-mu*k1*k1*n3*n3)/(sqrt(mu)*k1*gamma_3[i]);
                W23(i+nOrdsSq,i)=(mu*k1*k1*n3*n3-beta_3[i%nOrds]*beta_3[i%nOrds])/(sqrt(mu)*k1*gamma_3[i]);
        }
        W23=invert(W23);
}

void Grating2D::setSMatrix(void)
{
        Matrix< complex<double> > unit(2*nOrdsSq,2*nOrdsSq),T_tilde(2*nOrdsSq,2*nOrdsSq);
        unit.setunit();
        Vector < complex<double> > expTerm(2*nOrdsSq);
        // set R=0 and T=1 ([1]: equation 16a)
        calcBoundaryW();        
        
        // calculate the initial values of T and R. This method is not very elegant but it works
        GratingLayer2D dummy;
        Matrix< complex<double> > dummyn;
        dummyn.resize(1,1);
        dummy.setParameters(nOrds, periodX, periodY, 0., lambda, theta, phi, n1);
        dummyn(0,0)=n1;
        dummy.setRefInd(dummyn);
        dummy.solveHomogeneous();

        if ((signed int)R.size()!=pow(2*nOrdsSq,2)) R.resize(0.,2*nOrdsSq,2*nOrdsSq);
        if ((signed int)R_tilde.size()!=pow(2*nOrdsSq,2)) R_tilde.resize(0.,2*nOrdsSq, 2*nOrdsSq);
        if ((signed int)T.size()!=pow(2*nOrdsSq,2)) T.resize(0,2*nOrdsSq, 2*nOrdsSq);
        if ((signed int)G.size()!=pow(2*nOrdsSq,2)) G.resize(0.,2*nOrdsSq, 2*nOrdsSq);
        if ((signed int)tau.size()!=pow(2*nOrdsSq,2)) tau.resize(0.,2*nOrdsSq, 2*nOrdsSq);
        if ((signed int)Q12.size()!=pow(2*nOrdsSq,2)) Q12.resize(0.,2*nOrdsSq,2*nOrdsSq);
        if ((signed int)Q2.size()!=pow(2*nOrdsSq,2)) Q2.resize(0.,2*nOrdsSq,2*nOrdsSq);
        if ((signed int)trans.size()!=2*nOrdsSq) trans.resize(0.,2*nOrdsSq);
        if ((signed int)ref.size()!=2*nOrdsSq) ref.resize(0.,2*nOrdsSq);
        T.setunit();
        R.fill(0.);
        R_tilde.fill(0.);
        Q2=layer[0].invEigVecH*dummy.eigVecH;
        //Q1=layer[0].invEigVecE*dummy.eigVecE;
        // Q12 = Q1+Q2
        Q12=layer[0].invEigVecE*dummy.eigVecE+Q2;
        tau=invert(Q12);

        // calculate G ([2]: equation 11)
        for (int j=0;j<2*nOrdsSq;j++) 
        {
                expTerm(j)=exp(I*dummy.getGamma(j)*dummy.getThickness());
        }
        
        // calculate new R_ud and T_dd matrices ([2]: equations 13a,13b)
        R=unit-complex<double>(2.,0.)*Q2*tau;
        T=complex<double>(2.)*T;
        T.mult_rows(expTerm);
        T=T*tau;
        // calculate the Transition and Reflection matrix
        for (int i=0;i<nLayers;i++)
        {
                // calculate G ([2]: equation 11)
                for (int j=0;j<2*nOrdsSq;j++) 
                {
                        expTerm(j)=exp(I*layer[i].getGamma(j)*layer[i].getThickness());
                }
                G=0.; tau=0.;
                R_tilde=multFrontBack(R,expTerm);
                // calculate Q_1/2, G, tau ([2]: equation 11)
                if (i<nLayers-1)
                {
                        G=layer[i+1].invEigVecH*layer[i].eigVecH*(unit-R_tilde);
                        tau=layer[i+1].invEigVecE*layer[i].eigVecE*(unit+R_tilde)+G;
                        tau=invert(tau);
                }
                else
                {
                        G=W23*layer[i].eigVecH*(unit-R_tilde);
                        tau=layer[i].eigVecE*(unit+R_tilde)+G;
                        tau=invert(tau);
                }
                R=unit-complex<double>(2.)*G*tau;
                T=complex<double>(2.)*multFront(T,expTerm)*tau;

        }

        trans=T*incid;
        ref=R*incid;

        
        if (memorysave)
        {
                R_tilde.resize(0);
                T_tilde.resize(0);
                Q12.resize(0);
                Q2.resize(0);
                G.resize(0);
                tau.resize(0);
        }       
}

void Grating2D::calcEfficiencyR()
{
        double e=0.;
        efficiencyR.resize(nOrds, nOrds);
        for (int m=0;m<nOrds;m++)
        {
                for (int n=0;n<nOrds;n++)
                {
                        efficiencyR(m,n)=1./gamma_1[m*nOrds+n]*(
                        (k1*k1-beta_1[n]*beta_1[n])*norm(getReflection("x",m,n))+(k1*k1-alpha_1[m]*alpha_1[m])*norm(getReflection("y",m,n))+alpha_1[m]*beta_1[n]*(getReflection("x",m,n)*conj(getReflection("y",m,n))+getReflection("y",m,n)*conj(getReflection("x",m,n))));
                        e+=abs(efficiencyR(m,n));
                }
        }
}

void Grating2D::calcEfficiencyT()
{
        double e=0.;
        efficiencyT.resize(nOrds, nOrds);
        for (int m=0;m<nOrds;m++)
        {
                for (int n=0;n<nOrds;n++)
                {
                        efficiencyT(m,n)=1./gamma_3[m*nOrds+n]*((k3*k3-beta_3[n]*beta_3[n])*norm(getTransmission("x",m,n))+(k3*k3-alpha_3[m]*alpha_3[m])*norm(getTransmission("y",m,n))+alpha_3[m]*beta_3[n]*(getTransmission("x",m,n)*conj(getTransmission("y",m,n))+getTransmission("y",m,n)*conj(getTransmission("x",m,n))));
                        e+=abs(efficiencyT(m,n));
                }
        }
}

complex<double> Grating2D::calcTransField(double x, double y, double z, string pol)
{
        complex<double> f=0, f2=0;
        if (pol=="x") 
        {
                for (int i=0;i<nOrdsSq;i++)             
                {
                        f+= trans[i]*exp(I*(alpha_3[i/nOrds]*x+beta_3[i%nOrds]*y-gamma_3[i]*z)); 
                }
        }
        else if (pol=="y") 
        {       
                for (int i=0;i<nOrdsSq;i++)     
                {
                        f+= trans[i+nOrdsSq]*exp(I*(alpha_3[i/nOrds]*x+beta_3[i%nOrds]*y-gamma_3[i]*z));
                }
        }
        return f;
}
void Grating2D::writeField(string filename, string refTrans, string pol, double z, int pointsX, int pointsY)
{
        ofstream out;
        complex<double> c;
        out.open((filename+".dat").c_str());
        for (int i=0;i<pointsX;i++)
        {
                for (int j=0;j<pointsY;j++)
                {       
                        c=calcTransField((1.*i)/pointsX*periodX, (1.*j)/pointsY*periodY, z, pol);
                        if (refTrans=="t") out << norm(c) <<"\t";
                }
                out<<endl;
        }
        out.close();
}

void Grating2D::writeParameters(string filename)
{
        ofstream out;
        out.open((filename).c_str());
        out<<"Lambda=\t"<<lambda<<endl;
        out<<"Phi=\t"<<phi<<endl;
        out<<"Theta=\t"<<theta<<endl;
        out<<"n1=\t"<<real(n1)<<"\t"<<imag(n1)<<endl;
        out<<"n3=\t"<<real(n3)<<"\t"<<imag(n3)<<endl;
        out<<"PeriodX=\t"<<periodX<<endl;
        out<<"PeriodY=\t"<<periodY<<endl;
        out<<"nOrds=\t"<<nOrds<<endl;
        out<<"nPixelsX=\t"<<layer[0].getnPixelsX()<<endl;
        out<<"nPixelsY=\t"<<layer[0].getnPixelsY()<<endl;
    out <<"\n\nREFRACTIVE INDEX\n";
        for (unsigned int i=0;i<layer[0].getRefInd().rows();i++)
        {
                for (unsigned int j=0;j<layer[0].getRefInd().columns();j++)
                {
                        out<<real(layer[0].getRefInd(i,j))<<"\t"<<imag(layer[0].getRefInd(i,j))<<endl;
                }
        }
        out<<"\n\nTRANSMISSION ORDERS\n";
        for (int i=0;i<2*nOrdsSq;i++) out<<real(trans[i])<<"\t"<<imag(trans[i])<<endl;
        out<<"\n\nREFLECTION ORDERS\n";
        for (int i=0;i<2*nOrdsSq;i++) out<<real(ref[i])<<"\t"<<imag(ref[i])<<endl;
        out.close();
}

//void Grating2D::calcEmbOrds(int nEmbOrds, double dPadd)
//{
//    nOrdsEmb=nEmb;
//    periodXEmb=size[0];
//    periodYEmb=size[1];
//    alphaEmb.resize(nOrdsEmb);
//    betaEmb.resize(nOrdsEmb);
//    gammaEmb.resize(nOrdsEmb);
//    f=Vector< complex<double> >(2*nOrdsEmb*nOrdsEmb);
//    complex<double> x[2], y[2];
//      complex<double> xIntegral, yIntegral;
//
//    x[0]=(size[0]-periodX)/2.;
//    x[1]=x[0]+periodX;
//    y[0]=(size[1]-periodY)/2.;
//    y[1]=y[0]+periodY;
//    for (int s=0;s<nOrdsEmb;s++)
//      {
//              for (int t=0;t<nOrdsEmb;t++)
//              {
//            for (int p=0;p<nOrds;p++)
//                      {
//                if (alphaEmb[p]==alpha[s]):
//                    xIntegral=x[1]-x[0];
//                else
//                              {
//                    xIntegral=1./(I*(alpha[p]-alphaEmb[s]))*(exp(1j*(alphaEmb[s]-alpha[p])*x[1])-exp(1j*(alphaEmbb[s]-alpha[p])*x[0]));
//                              }
//                for (int q=0;q<nOrds;q++)
//                              {
//                    if (betaEmb[q]==beta[t])
//                        yIntegral=y[1]-y[0];
//                    else:
//                        yIntegral=1./(1j*(self.b[q]-self.bEmb[t]))*(exp(1j*(self.bEmb[t]-self.b[q])*y[1])-exp(1j*(self.bEmb[t]-self.b[q])*y[0]))
//                    for pol in range (0,2):
//                        f[pol][s][t]+=-1./(size[0]*size[1])*self.T[pol][p][q]*exp(1j*(self.a[p]*x[0]+self.b[q]*y[0]))*xIntegral*yIntegral
//                              }
//                      }
//              }
//      }
//    return f
//
//}


double Grating2D::delta(int i, int j){return ((i==j)?1.:0.);}

---