00001 #include <iomanip>
00002 #include <fstream>
00003 #include <ctime>
00004 #include "tools.h"
00005 #include "./Grating1D.h"
00006
00007 #include "tbci/solver/bicgstab.h"
00008 #include "tbci/solver/cgs.h"
00009 #include "tbci/solver/diagprecond.h"
00010
00011
00012 #include "tbci/solver/noprecond.h"
00013 #include "tbci/solver/lu_solver.h"
00014
00015
00016 #ifndef FIND_NEXT
00017 #define FIND_NEXT
00018 #define find_next(file,symbol) {while( file.get() != symbol && !file.eof() ){}}
00019 #endif
00020
00021 struct cplx{double re; double im;};
00022
00023
00024 extern "C" void zgbsv_(int *, int *, int *, int *, cplx *,int *, int *, cplx *, int *, int *);
00025
00026 #ifndef STD
00027 #define STD
00028 using namespace std;
00029 #endif
00030
00031
00032
00039 Grating1D::Grating1D(int truncation, int layers, int nSegments)
00040 {
00041 stacks = 1;
00042 gratingSegs=nSegments;
00043 nPadSegs = 1;
00044 nPermSegs=nSegments+2;
00045 nLayers = layers;
00046 nOrds = truncation; maxOrd=(nOrds-1)/2; minOrd=-maxOrd;
00047 material.resize(nLayers,nPermSegs);
00048 indices.resize(nLayers,nPermSegs);
00049 transpt.resize(nLayers,nPermSegs);
00050 slabThick.resize(nLayers);
00051 nUnknowns = 2*truncation*layers*stacks;
00052 layer.assign(nLayers, GratingLayer(nOrds, nPermSegs));
00053 r.resize(nOrds);
00054 t.resize(nOrds);
00055 rEFF.resize(nOrds);
00056 tEFF.resize(nOrds);
00057 incid.resize(nOrds);
00058 I = complex<double>(0.,1.);
00059 refl.resize(nOrds);
00060 tran.resize(nOrds);
00061 gamma.resize(nOrds);
00062 polMode = "TE";
00063 }
00064
00065 Grating1D::~Grating1D()
00066 {
00067 }
00068
00072 void Grating1D::setTheta(double angle)
00073 {
00074 theta = angle*M_PI/180.;
00075 for (int i=0;i<nLayers;i++) layer[i].setAngle(angle);
00076 }
00080 double Grating1D::getTheta() {return theta*180./M_PI;}
00085 void Grating1D::setLambda(double wavelength)
00086 {
00087 lambda = wavelength;
00088 for (int i=0;i<nLayers;i++) layer[i].setLambda(lambda);
00089 calcGamma();
00090 }
00092 double Grating1D::getLambda(){return lambda;}
00097 void Grating1D::setPeriod(double size)
00098 {
00099 period =size;
00100 calcGamma();
00101 for (int i=0;i<nLayers;i++) layer[i].setPeriod(period);
00102 }
00104 double Grating1D::getPeriod(void){return period;}
00106 double Grating1D::getGratingSize(void){return gratingSize;}
00109 void Grating1D::setPadding(double pad)
00110 {
00111 padding = pad;
00112 setPeriod(gratingSize+2.*padding);
00113 for (int i=0;i<nLayers;i++)
00114 {
00115 for (int j=1;j<nPermSegs-1;j++) {transpt(i,j)=(padding+(gratingSize*j)/gratingSegs)/period;}
00116 transpt(i,0)=padding/period;
00117 transpt(i,nPermSegs-1)=1.;
00118 }
00119
00120 for (int i=0; i<nLayers; i++) layer[i].setParameters(lambda, theta, n0, period, slabThick(i), transpt.get_row(i), indices.get_row(i), polMode);
00121 }
00122
00124 double Grating1D::getPadding(void){return padding;}
00128 void Grating1D::setnPermSegs(int segments)
00129 {
00130 nPermSegs=segments;
00131 for (int i=0;i<nLayers;i++) layer[i].setnPermSegs(nPermSegs);
00132 for (int i=0; i<nLayers; i++) layer[i].setParameters(lambda, theta, n0, period, slabThick(i), transpt.get_row(i), indices.get_row(i), polMode);
00133 }
00134
00136 int Grating1D::getnPermSegs(void){return nPermSegs;}
00140 void Grating1D::setPolMode(string mode)
00141 {
00142 polMode = mode;
00143 for (int i=0;i<nLayers;i++) layer[i].setPolMode(polMode);
00144 }
00145
00147 string Grating1D::getPolMode(){return polMode;}
00149 void Grating1D::setN0(complex<double> index)
00150 {
00151 n0 = index;
00152 for (int i=0;i<nLayers;i++) layer[i].setN0(n0);
00153 }
00155 complex<double>Grating1D::getN0(){return n0;}
00157 void Grating1D::setN3(complex<double> index){n3 = index;}
00159 complex<double>Grating1D::getN3(){return n3;}
00161 void Grating1D::setNPad(complex<double> index){nPad = index;}
00163 complex<double>Grating1D::getNPad(){return nPad;}
00165 void Grating1D::setMatPad(int mp){matPad=mp;}
00167 int Grating1D::getnOrds(){return nOrds;}
00169 void Grating1D::setIncid(Vector< complex<double> >incident){for (int i=0; i<nOrds; i++) {incid[i]=incident[i];}}
00174 complex<double>Grating1D::getIncid(int i){return incid[i];}
00179 complex<double>Grating1D::getRhsVec(int i){return rhsVec[i];}
00180 complex<double> Grating1D::getA(int layer, int comp){return solVec[layer*2*nOrds+2*(comp+maxOrd)];}
00181 complex<double> Grating1D::getAdash(int layer, int comp){return solVec[layer*2*nOrds+2*(comp+maxOrd)+1];}
00182 double Grating1D::getEfficiencyR(int i){return rEFF[i];}
00183 double Grating1D::getEfficiencyT(int i){return tEFF[i];}
00184 complex<double> Grating1D::getTrans(int i){return tran[i];}
00185 complex<double> Grating1D::getRef(int i){return refl[i];}
00186 int Grating1D::getnLayers(void){return nLayers;}
00187 double Grating1D::getEmbSize(){return embPeriod;}
00188 complex<double>Grating1D::getGamma(int n){return gamma[n];}
00189 void Grating1D::calcGamma(void)
00190 {
00191 complex<double> gamma0;
00192
00193 gamma0=n0*2.*M_PI/lambda*sin(theta);
00194 for(int i=0;i<nOrds;i++)
00195 {
00196 gamma[i] = gamma0 + 2*M_PI*(i+minOrd)/period;
00197 }
00198 }
00199
00203 complex<double> Grating1D::getLhsMat(int i, int j)
00204 {
00205 return lhsMatBand(i,j);
00206 }
00207
00220 void Grating1D::setParameters(double theta, double lambda, double size, complex<double> n0, complex<double>n3, complex<double>np, double pad, Vector < complex<double> >incident, string pol="TE")
00221 {
00222 setTheta(theta);
00223 setPeriod(size+2.*pad);
00224 gratingSize=size;
00225 setPadding(pad);
00226 setLambda(lambda);
00227 setN0(n0);
00228 setN3(n3);
00229 setNPad(np);
00230 setIncid(incident);
00231 setPolMode(pol);
00232 calcGamma();
00233 }
00234
00235 void Grating1D::calcIncid(Vector <complex<double> > field, int dim)
00236 {
00237 complex <double> c;
00238 for (int p=-(nOrds/2);p<nOrds/2;p++)
00239 {
00240 c=0.;
00241 for (int j=0;j<dim;j++)
00242 {
00243 if (p!=0) c+= I/(2.*M_PI*p)*field[j]*(exp((-I*2.*M_PI*(p*(j+1.)))/(1.*dim))-exp((-I*2.*(M_PI*p*j))/(1.*dim)));
00244 else c+= field[j];
00245 }
00246 c/=period;
00247 if (abs(c)<1e-15) c=0.;
00248 incid[p+nOrds/2]=c;
00249 }
00250 }
00251
00256 void Grating1D::importMaterials(string mFile)
00257 {
00258 nMat = nLines(mFile);
00259 char *pEnd;
00260 string newline;
00261 refIndMat.resize(nMat);
00262 double NCoeffs[nMat][8];
00263 ifstream fMat;
00264 fMat.open(mFile.c_str());
00265 for (int i=0;i<nMat;i++)
00266 {
00267 newline='#';
00268 while (newline[0] == '#' || string(newline).find_first_not_of(" \t\v\r\n")==string::npos) getline(fMat,newline);
00269 for (int j=0;j<8;j++) NCoeffs[i][j] = strtod(((j==0)? newline.c_str() : pEnd),&pEnd);
00270 refIndMat[i] = NCoeffs[i][0];
00271 for(int j=0;j<3;j++)
00272 {
00273 refIndMat[i] += (NCoeffs[i][2*j+2]*pow(lambda, 2.0))/(pow(lambda, 2.0) - pow(NCoeffs[i][2*j+3],2.0));
00274 }
00275 refIndMat[i]=complex<double>(real(refIndMat[i]),NCoeffs[i][1]);
00276 }
00277 fMat.close();
00278
00279 }
00280
00287 void Grating1D::initLayers(string initfile, string initStyle, string matFile)
00288 {
00289 if (initStyle == "material") importMaterials(matFile.c_str());
00290 double ind;
00291 int mat;
00292
00293 ifstream finit;
00294 finit.open(initfile.c_str());
00295
00296 for (int i=0; i<nLayers; i++)
00297 {
00298 finit>>slabThick[i]; find_next(finit, ',');
00299 for (int segment=1;segment<nPermSegs-1;segment++)
00300 {
00301 if (initStyle == "index")
00302 {
00303 finit >> ind;
00304 indices(i,segment)=complex<double>(ind,0.);
00305 }
00306 else if (initStyle == "material")
00307 {
00308 finit >> mat;
00309 material(i,segment)=mat;
00310 indices(i,segment)=refIndMat[mat];
00311 }
00312 transpt(i,segment)=(padding+(gratingSize*segment)/(gratingSegs))/period;
00313 }
00314 transpt(i,0)=padding/period;
00315 transpt(i,nPermSegs-1)=1.;
00316 indices(i,0)=nPad;
00317 indices(i,nPermSegs-1)=nPad;
00318 material(i,0)=-1;
00319 material(i,nPermSegs-1)=-1;
00320 }
00321 finit.close();
00322
00323 for (int i=0; i<nLayers; i++) layer[i].setParameters(lambda, theta, n0, period, slabThick(i), transpt.get_row(i), indices.get_row(i), polMode);
00324 }
00325
00332 void Grating1D::initLayers(F_Matrix<int> element, Vector<double> thickness, string matFile)
00333 {
00334 importMaterials(matFile.c_str());
00335 for (int i=0; i<nLayers; i++)
00336 {
00337 slabThick[i]=thickness[i];
00338 for (int segment=0;segment<nPermSegs;segment++)
00339 {
00340 if (segment<nPadSegs || segment >=nPermSegs-nPadSegs) indices(i,segment)=nPad;
00341 else if (segment>=nPadSegs && segment<nPermSegs-nPadSegs) indices(i,segment)=refIndMat[element(i,segment-nPadSegs)];
00342 transpt(i,segment)=(segment+1.)/nPermSegs;
00343 }
00344 }
00345 for (int i=0; i<nLayers; i++) layer[i].setParameters(lambda, theta, n0, period, slabThick(i), transpt.get_row(i), indices.get_row(i), polMode);
00346 }
00347
00353 void Grating1D::initLayers(F_Matrix< complex<double> > element, Vector<double> thickness)
00354 {
00355 for (int i=0; i<nLayers; i++)
00356 {
00357 slabThick[i]=thickness[i];
00358 for (int segment=1;segment<nPermSegs-1;segment++)
00359 {
00360 indices(i,segment)=element(i,segment-1);
00361 transpt(i,segment)=(padding+(gratingSize*segment)/(gratingSegs))/period;
00362 }
00363 transpt(i,0)=padding/period;
00364 transpt(i,nPermSegs-1)=1.;
00365 indices(i,0)=nPad;
00366 indices(i,nPermSegs-1)=nPad;
00367 }
00368 for (int i=0; i<nLayers; i++) layer[i].setParameters(lambda, theta, n0, period, slabThick(i), transpt.get_row(i), indices.get_row(i), polMode);
00369 }
00370
00377 void Grating1D::changeRefInd(int slab, int segment, complex<double> ind)
00378 {
00379 indices(slab,segment+nPadSegs) = ind;
00380 layer[slab].setRefInd(indices.get_row(slab));
00381 }
00388 void Grating1D::changeMaterial(int slab, int segment, int mat)
00389 {
00390 material(slab,segment+nPadSegs) = mat;
00391 indices(slab,segment+nPadSegs) = refIndMat[mat];
00392 layer[slab].setRefInd(indices.get_row(slab));
00393 }
00399 void Grating1D::addFrame(complex<double> ind, double frame)
00400 {
00401 F_Matrix< complex<double> > indicesNew(nLayers,nPermSegs+2);
00402 F_Matrix<double> transptNew(nLayers,nPermSegs+2);
00403 F_Matrix<int> materialNew(nLayers,nPermSegs+2);
00404 nPermSegs+=2;
00405 setPeriod(period+2.*frame);
00406 for (int i=0;i<nLayers;i++)
00407 {
00408 for (int j=2;j<nPermSegs-2;j++)
00409 {
00410 indicesNew(i,j)=indices(i,j-1);
00411 transptNew(i,j)=(padding+frame+(transpt(i,j-1)-transpt(i,0))*(period-2.*frame))/period;
00412 }
00413
00414 transptNew(i,0)=padding/period;
00415 transptNew(i,nPermSegs-1)=1.;
00416
00417 transptNew(i,1)=(padding+frame)/period;
00418 transptNew(i,nPermSegs-2)=(period-padding)/period;
00419
00420 indicesNew(i,0)=nPad;
00421 indicesNew(i,nPermSegs-1)=nPad;
00422 indicesNew(i,1)=ind;
00423 indicesNew(i,nPermSegs-2)=ind;
00424
00425 materialNew(i,0)=-1;
00426 materialNew(i,nPermSegs-1)=-1;
00427 materialNew(i,1)=-1;
00428 materialNew(i,nPermSegs-2)=-1;
00429 }
00430
00431 gratingSize+=2.*frame;
00432
00433
00434 material.resize(nLayers,nPermSegs);
00435 indices.resize(nLayers,nPermSegs);
00436 transpt.resize(nLayers,nPermSegs);
00437
00438 indices=indicesNew;
00439 transpt=transptNew;
00440 setnPermSegs(nPermSegs);
00441 for (int i=0; i<nLayers; i++) layer[i].setParameters(lambda, theta, n0, period, slabThick(i), transpt.get_row(i), indices.get_row(i), polMode);
00442 }
00443
00448 void Grating1D::solve(string method)
00449 {
00450 if (method=="T")
00451 {
00452 if ((signed int) rhsVec.size()!=nUnknowns) rhsVec.resize(nUnknowns);
00453 if ((signed int) lhsMatBand.size()!=(8*nLayers-4)*nOrds*nOrds) lhsMatBand.resize(nUnknowns,3*nOrds-1,3*nOrds-1);
00454 if ((signed int) solVec.size()!=nUnknowns) solVec.resize(nUnknowns);
00455 setEquations();
00456 solveEquations();
00457 findEfficiencyT();
00458 }
00459 else if (method=="S")
00460 {
00461 if ((signed int) F.size()!=nOrds*nOrds) F.resize(nOrds, nOrds);
00462 if ((signed int) G.size()!=nOrds*nOrds) G.resize(nOrds, nOrds);
00463 if ((signed int) R.size()!=nOrds*nOrds) R.resize(nOrds, nOrds);
00464 if ((signed int) T.size()!=nOrds*nOrds) T.resize(nOrds, nOrds);
00465 if ((signed int) tau.size()!=nOrds*nOrds) tau.resize(nOrds, nOrds);
00466 if ((signed int) T_tilde.size()!=nOrds*nOrds) T_tilde.resize(nOrds, nOrds);
00467 if ((signed int) R_tilde.size()!=nOrds*nOrds) R_tilde.resize(nOrds, nOrds);
00468 if ((signed int) Q1.size()!=nOrds*nOrds) Q1.resize(nOrds, nOrds);
00469 if ((signed int) Q2.size()!=nOrds*nOrds) Q2.resize(nOrds, nOrds);
00470 setSMatrix();
00471 findEfficiencyS();
00472 }
00473 }
00480 void Grating1D::setEquations()
00481 {
00482 if ("TE"==polMode) {setLhsMatTE();}
00483 else if ("TM" == polMode) setLhsMatTM();
00484 else cout<<"Polarization not properly set\n";
00485 }
00486
00487 void Grating1D::setLhsMatTE(void)
00488 {
00489 int row;
00490 int column;
00491 complex<double> em,c1, c2, c3, c4, z, alpha;
00492
00493
00494 za = pow(n0*2.*M_PI/lambda, 2);
00495 zb = pow(n3*2.*M_PI/lambda, 2);
00496
00497 for(int i=0;i<nOrds;i++)
00498 {
00499 z = pow(layer[0].getGamma(i),2);
00500 r[i]=sqrt(za - z);
00501 if(real(r[i])==0.0 && imag(r[i])==0.0)
00502 cout << "WARNING: r(" << minOrd+i << ") is zero" << endl ;
00503
00504 t[i]=sqrt(zb - z);
00505 if(real(t[i])==0.0 && imag(t[i])==0.0)
00506 cout << "WARNING: t(" << minOrd+i <<") is zero" << endl;
00507 }
00508
00509
00510
00511
00512
00513
00514
00515 for (int slab=0;slab<nLayers;slab++)
00516 {
00517 if (slab == 0)
00518 {
00519 for (int p = 0; p < nOrds; p++)
00520 {
00521 for (int m = 0; m < nOrds; m++)
00522 {
00523 lhsMatBand.setval((1.+layer[slab].getAlpha(m)/r[p])*layer[slab].getEigVecE(p,m), p, 2*m );
00524
00525 lhsMatBand.setval((1.-layer[slab].getAlpha(m)/r[p])*exp(I*layer[slab].getAlpha(m)*layer[slab].getThickness())*layer[slab].getEigVecE(p,m), p,2*m+1);
00526
00527 }
00528 }
00529 }
00530 if (slab == nLayers-1)
00531 {
00532 row = nOrds+slab*2*nOrds+(stacks-1)*nLayers*2*nOrds;
00533 column = slab*2*nOrds+(stacks-1)*nLayers*2*nOrds;
00534 for (int p = 0; p < nOrds; p++)
00535 {
00536 for (int m = 0; m < nOrds; m++)
00537 {
00538 lhsMatBand.setval((1.-layer[slab].getAlpha(m)/t[p])*exp(I*layer[slab].getAlpha(m)*layer[slab].getThickness())*layer[slab].getEigVecE(p,m), p+row,2*m+column);
00539
00540 lhsMatBand.setval((1.+layer[slab].getAlpha(m)/t[p])*layer[slab].getEigVecE(p,m) , p+row,2*m+1+column);
00541
00542 }
00543 }
00544 }
00545
00546
00547
00548 row = nOrds*(2*slab+1);
00549 column = slab*2*nOrds;
00550 for (int p = 0; p < nOrds; p++){
00551 for (int m = 0; m < nOrds; m++){
00552 alpha=layer[slab].getAlpha(m);
00553 c1 = exp(I*alpha*layer[slab].getThickness())*layer[slab].getEigVecE(p,m);
00554 c2 = layer[slab].getEigVecE(p,m);
00555 if (slab!=nLayers-1)
00556 {
00557
00558
00559
00560 lhsMatBand.setval(c1, 2*p+row, 2*m+column);
00561 lhsMatBand.setval(c2, 2*p+row, 2*m+1+column);
00562
00563 lhsMatBand.setval(alpha*c1, 2*p+1+row, 2*m+column);
00564 lhsMatBand.setval(-alpha*c2, 2*p+1+row, 2*m+1+column);
00565 }
00566
00567 if (slab!=0)
00568 {
00569
00570 lhsMatBand.setval(-c2, 2*p+row-layerLines, 2*m+column);
00571 lhsMatBand.setval(-c1, 2*p+row-layerLines, 2*m+1+column);
00572
00573 lhsMatBand.setval(-c2*alpha, 2*p+1+row-layerLines, 2*m+column);
00574 lhsMatBand.setval(alpha*c1, 2*p+1+row-layerLines, 2*m+1+column);
00575
00576 }
00577 }
00578 }
00579 }
00580 }
00581
00582
00583
00584 void Grating1D::setLhsMatTM(void)
00585 {
00586 int row;
00587 int column;
00588 complex<double> em,c1, c2, c3, c4, z,sum;
00589
00590
00591 za = pow(n0*2.*M_PI/lambda, 2);
00592 zb = pow(n3*2.*M_PI/lambda, 2);
00593
00594 for(int i=0;i<nOrds;i++)
00595 {
00596 z = pow(layer[0].getGamma(i),2);
00597 r[i]=sqrt(za - z);
00598 if(real(r[i])==0.0 && imag(r[i])==0.0)
00599 cout << "WARNING: r(" << minOrd+i << ") is zero" << endl ;
00600
00601 t[i]=sqrt(zb - z);
00602 if(real(t[i])==0.0 && imag(t[i])==0.0)
00603 cout << "WARNING: t(" << minOrd+i <<") is zero" << endl;
00604 }
00605
00606
00607
00608
00609 for (int slab=0;slab<nLayers;slab++)
00610 {
00611
00612
00613 if (slab == 0)
00614 {
00615 for (int p = 0; p < nOrds; p++)
00616 {
00617 for (int m = 0; m < nOrds; m++)
00618 {
00619 sum=0.;
00620 for (int q=0;q<nOrds;q++)
00621 {
00622 sum += layer[slab].getPermCoeff(p-q)*r[q]*layer[slab].getEigVecE(q,m);
00623 }
00624 lhsMatBand.setval(n0*n0*layer[slab].getAlpha(m)*layer[slab].getEigVecE(p,m)+sum, p, 2*m);
00625 lhsMatBand.setval((-n0*n0*layer[slab].getAlpha(m)*layer[slab].getEigVecE(p,m)+sum)*exp(I*layer[slab].getAlpha(m)*layer[slab].getThickness()), p, 2*m+1);
00626
00627 }
00628 }
00629 }
00630 if (slab == nLayers-1)
00631 {
00632 row = nOrds+slab*2*nOrds+(stacks-1)*nLayers*2*nOrds;
00633 column = slab*2*nOrds+(stacks-1)*nLayers*2*nOrds;
00634 for (int p = 0; p < nOrds; p++)
00635 {
00636 for (int m = 0; m < nOrds; m++)
00637 {
00638 sum=0.;
00639 for (int q=0;q<nOrds;q++)
00640 {
00641 sum += layer[slab].getPermCoeff(p-q)*t[q]*layer[slab].getEigVecE(q,m);
00642 }
00643 lhsMatBand.setval((-n3*n3*layer[slab].getAlpha(m)*layer[slab].getEigVecE(p,m)+sum)*exp(I*layer[slab].getAlpha(m)*layer[slab].getThickness()), p+row, 2*m+column);
00644 lhsMatBand.setval(n3*n3*layer[slab].getAlpha(m)*layer[slab].getEigVecE(p,m)+sum, p+row, 2*m+1+column);
00645 }
00646 }
00647 }
00648
00649
00650 row = nOrds*(2*slab+1);
00651 column = slab*2*nOrds;
00652 for (int p = 0; p < nOrds; p++){
00653 for (int m = 0; m < nOrds; m++)
00654 {
00655 c1 = exp(I*layer[slab].getAlpha(m)*layer[slab].getThickness())*layer[slab].getEigVecE(p,m);
00656 c2 = layer[slab].getEigVecE(p,m);
00657
00658 c3=0.;
00659
00660 for (int i=0;i<nOrds;i++)
00661 {
00662 if (slab!=nLayers-1) c3 += layer[slab+1].getPermCoeff(p-i)*layer[slab].getEigVecE(i,m);
00663 else if (slab==nLayers-1) c3 += layer[0].getPermCoeff(p-i)*layer[slab].getEigVecE(i,m);
00664 }
00665 c4=0.;
00666
00667 for (int i=0;i<nOrds;i++)
00668 {
00669 if (slab!=nLayers-1) c4 += layer[slab].getPermCoeff(p-i)*layer[slab+1].getEigVecE(i,m);
00670 else if (slab==nLayers-1) c4 += layer[slab].getPermCoeff(p-i)*layer[0].getEigVecE(i,m);
00671 }
00672
00673 for (int i=0;i<(slab == nLayers-1? stacks-1: stacks); i++)
00674 {
00675
00676
00677
00678 lhsMatBand.setval(c1, 2*p+row+i*stackLines, 2*m+column+i*stackColumns);
00679 lhsMatBand.setval(c2, 2*p+row+i*stackLines, 2*m+1+column+i*stackColumns);
00680
00681
00682 lhsMatBand.setval(c3*layer[slab].getAlpha(m)*exp(I*layer[slab].getAlpha(m)*layer[slab].getThickness()), 2*p+1+row+i*stackLines, 2*m+column+i*stackColumns);
00683 lhsMatBand.setval(-c3*layer[slab].getAlpha(m), 2*p+1+row+i*stackLines, 2*m+1+column+i*stackColumns);
00684
00685 lhsMatBand.setval(-c4*layer[slab+1].getAlpha(m), 2*p+1+row+i*stackLines, 2*m+column+i*stackColumns+2*nOrds);
00686 lhsMatBand.setval(c4*layer[slab+1].getAlpha(m)*exp(I*layer[slab+1].getAlpha(m)*layer[slab+1].getThickness()), 2*p+1+row+i*stackLines, 2*m+1+column+i*stackColumns+2*nOrds);
00687
00688 }
00689
00690 for (int i=(slab==0?1:0); i<stacks; i++)
00691 {
00692
00693 lhsMatBand.setval(-c2, 2*p+row-layerLines+i*stackLines, 2*m+column+i*stackColumns);
00694 lhsMatBand.setval(-c1, 2*p+row-layerLines+i*stackLines, 2*m+1+column+i*stackColumns);
00695 }
00696 }
00697 }
00698 }
00699 }
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00733 void Grating1D::solveEquations(void)
00734 {
00735 int i, nrhs=1, ifail=0;
00736 cplx *lhsMatNag, *rhsVecNag;
00737
00738 rhsVecNag = new cplx[nUnknowns];
00739
00740
00741 if ("TE"== polMode)
00742 {
00743 for (int i=0; i<nUnknowns; i++){rhsVec[i]=(i<nOrds ? 2.*incid[i] : 0.);}
00744 }
00745 if ("TM" == polMode)
00746 {
00747 for (int i=0;i<nUnknowns;i++)
00748 {
00749 rhsVec[i]=0.;
00750 if (i<nOrds)
00751 {
00752 for (int j=0;j<nOrds;j++) rhsVec[i]+= 2.*layer[0].getPermCoeff(i-j)*r[j]*incid[j];
00753 }
00754 }
00755 }
00756 for (int i=0; i < nUnknowns; i++){
00757 rhsVecNag[i].re = real(rhsVec[i]);
00758 rhsVecNag[i].im = imag(rhsVec[i]);
00759 }
00760
00761 int *IPIV;
00762 int KL = 3*nOrds-1;
00763 int KU = 3*nOrds-1;
00764
00765 int LDAB = 2*KL+KU+1;
00766 lhsMatNag = new cplx[nUnknowns*LDAB];
00767 IPIV = new int[nUnknowns];
00768
00769 for (i=0;i<nUnknowns;i++)
00770 {
00771 for (int j=((i-KL)<0? 0 : i-KL); j<((i+KU)>=nUnknowns? nUnknowns : i+KU); j++)
00772 {
00773 lhsMatNag[2*(3*nOrds-1)+i-j+j*(3*(3*nOrds-1)+1)].re=real(lhsMatBand(i,j));
00774 lhsMatNag[2*(3*nOrds-1)+i-j+j*(3*(3*nOrds-1)+1)].im=imag(lhsMatBand(i,j));
00775 }
00776 }
00777
00778
00779 zgbsv_(&nUnknowns, &KL, &KU, &nrhs, lhsMatNag, &LDAB, IPIV, rhsVecNag, &nUnknowns, &ifail);
00780
00781 for (int i=0;i<nUnknowns;i++){solVec[i]=complex<double> (rhsVecNag[i].re,rhsVecNag[i].im);}
00782
00783
00784
00785 if(ifail!=0) {
00786 cout << "Failed to solve linear system. Error ifail = "<< ifail << endl;
00787 }
00788
00789 delete [] rhsVecNag;
00790 delete [] lhsMatNag;
00791 delete [] IPIV;
00792 }
00793
00798 void Grating1D::setSMatrix(void)
00799 {
00800 complex<double> c;
00801 if ((signed int)expTerm.size()!=nOrds) expTerm.resize(nOrds);
00802 if ((signed int)T.size()!=nOrds*nOrds) T.resize(nOrds);
00803 if ((signed int)R.size()!=nOrds*nOrds) R.resize(nOrds);
00804
00805 GratingLayer lastBoundary(nOrds, 1);
00806 Vector<double> trans(1.,1);
00807 Vector< complex<double> >ind(n3,1);
00808 lastBoundary.setParameters(lambda, getTheta(), n0, period, 0., trans, ind, polMode);
00809 lastBoundary.solve("S");
00810
00811
00812 T.setunit();
00813 R.fill(0.);
00814
00815 for (int i=-1;i<nLayers;i++)
00816 {
00817 if (i==nLayers-1)
00818 {
00819 ind[0]=n0;
00820 lastBoundary.setParameters(lambda, getTheta(), n0, period, 0., trans, ind, polMode);
00821 lastBoundary.solve("S");
00822
00823 }
00824
00825 Q1.fill(0.); Q2.fill(0.);
00826 if (i==nLayers-1)
00827 {
00828 Q1=layer[i].eigVecE;
00829 Q2=lastBoundary.invEigVecH*layer[i].eigVecH;
00830 }
00831 else if(i==-1)
00832 {
00833 Q1=layer[i+1].invEigVecE*lastBoundary.eigVecE;
00834 Q2=layer[i+1].invEigVecH*lastBoundary.eigVecH;
00835 }
00836 else
00837 {
00838 Q1=layer[i+1].invEigVecE*layer[i].eigVecE;
00839 Q2=layer[i+1].invEigVecH*layer[i].eigVecH;
00840 }
00841
00842
00843 for (int j=0;j<nOrds;j++)
00844 {
00845 expTerm[j]=(i>-1 ? exp(I*layer[i].getAlpha(j)*layer[i].getThickness()) : 1.);
00846 }
00847 for (int j=0;j<nOrds;j++)
00848 {
00849 R_tilde.set_col(expTerm[j]*R.get_col(j),j);
00850 R_tilde.set_row(expTerm[j]*R.get_row(j),j);
00851 }
00852
00853 F=Q1*(unit(F.rows())+R_tilde);
00854 G=Q1*(unit(G.rows())-R_tilde);
00855 tau=inv(F+G);
00856
00857
00858
00859 for (int j=0;j<nOrds;j++)
00860 {
00861 T_tilde.set_col(T.get_col(j)*expTerm[j],j);
00862 }
00863 T=complex<double>(2.)*T_tilde*tau;
00864 R=unit(R.columns())-complex<double>(2.)*G*tau;
00865 }
00866 }
00867
00874 void Grating1D::findEfficiency(string method)
00875 {
00876 if (method=="S") findEfficiencyS();
00877 else if (method=="T") findEfficiencyT();
00878 }
00879
00880 void Grating1D::findEfficiencyS(void)
00881 {
00882 double incidPowFactor = 0.0;
00883 complex<double> z;
00884
00885 za = pow(n0*2.*M_PI/lambda, 2);
00886 zb = pow(n3*2.*M_PI/lambda, 2);
00887
00888 for(int i=0;i<nOrds;i++)
00889 {
00890 z = pow(layer[0].getGamma(i),2);
00891 r[i]=sqrt(za - z);
00892 if(real(r[i])==0.0 && imag(r[i])==0.0)
00893 cout << "WARNING: r(" << minOrd+i << ") is zero" << endl ;
00894
00895 t[i]=sqrt(zb - z);
00896 if(real(t[i])==0.0 && imag(t[i])==0.0)
00897 cout << "WARNING: t(" << minOrd+i <<") is zero" << endl;
00898 }
00899
00900
00901 for(int i=0;i<nOrds;i++){
00902 if(real(incid[i]) != 0. || imag(incid[i]) != 0.){
00903 if(imag(r[i]) != 0.)
00904 cout << "ERROR: evanescent incident wave "<<endl;
00905 incidPowFactor += real(r[i])*norm(incid[i]);
00906 }
00907 }
00908 for (int i=0;i<nOrds;i++)
00909 {
00910 tran[i]=0.; refl[i]=0.;
00911 tran=T*incid;
00912 refl=R*incid;
00913 rEFF[i] = norm(refl[i])*real(r[i])/incidPowFactor;
00914 tEFF[i] = norm(tran[i])*real(t[i])/incidPowFactor;
00915 }
00916 }
00917
00918
00919
00920 void Grating1D::findEfficiencyT(void)
00921 {
00922 incidPowFactor = 0.0;
00923 complex <double> em;
00924
00925
00926 for(int i=0;i<nOrds;i++){
00927 if(real(incid[i]) != 0. || imag(incid[i]) != 0.){
00928 if(imag(r[i]) != 0.)
00929 cout << "ERROR: evanescent incident wave "<<endl;
00930 incidPowFactor += real(r[i])*norm(incid[i]);
00931 }
00932 }
00933
00934 for(int i=0;i<nOrds;i++){
00935 refl[i] = -incid[i];
00936 tran[i] = complex<double>(0., 0.);
00937
00938
00939 for(int m=0;m<nOrds;m++){
00940 em = exp(I*layer[0].getAlpha(m) * layer[0].getThickness());
00941 refl[i] += layer[0].getEigVecE(i,m) * (solVec[2*m] + solVec[2*m+1]* em);
00942
00943 em = exp(I*layer[nLayers-1].getAlpha(m) *layer[nLayers-1].getThickness());
00944 tran[i] +=layer[nLayers-1].getEigVecE(i,m) * (solVec[nUnknowns-2*nOrds+2*m+1] + solVec[nUnknowns-2*nOrds+2*m] * em);
00945 }
00946
00947
00948 rEFF[i] = norm(refl[i])*real(r[i])/incidPowFactor;
00949 tEFF[i] = norm(tran[i])*real(t[i])/incidPowFactor;
00950 if ("TM"==polMode) tEFF[i]*=(norm(n0)/norm(n3));
00951 }
00952 }
00953
00958 void Grating1D::calcIntField(int X, int Y)
00959 {
00960 dimX = X;
00961 dimY = Y;
00962 double thickTot=0.;
00963 fieldAmplitude.resize(dimX,dimY);
00964 string filewithending;
00965
00966 calcGamma();
00967
00968 for (int i=0;i<nLayers;i++) thickTot+=layer[i].getThickness();
00969
00970 for (int j=0;j<dimX;j++)
00971 {
00972 for (int k=0;k<dimY;k++)
00973 {
00974 fieldAmplitude(j,k)=abs(getIntField((j+1.)/(1.*dimX), (k+1)/(1.*dimY)));
00975 }
00976 }
00977 }
00978
00979
00987 void Grating1D::calcExtField(int X, int Y, double propBehind, double start=0., double stop=1.)
00988 {
00989 dimX = X;
00990 dimY = Y;
00991 fieldAmplitude.resize(dimX,dimY);
00992 propAfter=propBehind;
00993 plotSize=(stop-start)*period;
00994
00995
00996 for (int j=0;j<dimX;j++)
00997 {
00998 for (int k=0;k<dimY;k++)
00999 {
01000 fieldAmplitude(j,k)=norm(getExtField((j/(1.*dimX)*(stop-start)+start), ((1.*k+1)/dimY)*propBehind));
01001 }
01002 }
01003 }
01004
01012 void Grating1D::calcEmbField(int X, int Y, double propBehind, double start=0., double stop=1.)
01013 {
01014 dimX = X;
01015 dimY = Y;
01016 fieldAmplitude.resize(dimX,dimY);
01017 string filewithending;
01018
01019 for (int j=0;j<dimX;j++)
01020 {
01021 for (int k=0;k<=dimY;k++)
01022 {
01023 fieldAmplitude(j,k)=norm(getEmbExtField((j/(1.*dimX)*(stop-start)+start), ((1.*k)/dimY)*propBehind));
01024 }
01025 }
01026 }
01030 void Grating1D::writeField(string file)
01031 {
01032 ofstream out;
01033 out.open(file.c_str());
01034 for (int k=0;k<dimX-1;k++)
01035 {
01036 for (int j=0;j<dimY-1;j++)
01037 {
01038 out << fieldAmplitude(k,j) <<" ";
01039 }
01040 out <<endl;
01041 }
01042 out.close();
01043 }
01044
01050 complex<double> Grating1D::getIntField(double x, double y)
01051 {
01052 complex<double> field=0.;
01053 complex<double> XTerm;
01054 double totThickness=0., thick=0.;
01055 int n=0;
01056
01057 for (int i=0;i<nLayers*stacks;i++) totThickness += layer[i%nLayers].getThickness();
01058 totThickness*=stacks;
01059 thick=layer[0].getThickness();
01060 while (y*totThickness>thick)
01061 {
01062 thick += layer[(n+1)%nLayers].getThickness();
01063
01064 n++;
01065 }
01066 thick-= layer[n].getThickness();
01067 y = (y*totThickness-thick)/layer[n%nLayers].getThickness();
01068 for (int m=minOrd;m<=maxOrd;m++)
01069 {
01070 XTerm=0.;
01071
01072 for (int p=0;p<nOrds;p++)
01073 {
01074 XTerm += layer[n%nLayers].getEigVecE(p,m+maxOrd)*exp(I*layer[n%nLayers].getGamma(p)*x*period);
01075 }
01076
01077 field += (getA(n,m)*exp(I*layer[n%nLayers].getAlpha(m+maxOrd)*y*(layer[n%nLayers].getThickness()))
01078 + getAdash(n,m)*exp(-I*layer[n%nLayers].getAlpha(m+maxOrd)*(y-1.)*(layer[n%nLayers].getThickness())))
01079 * XTerm;
01080 }
01081 return field;
01082 }
01083
01089 complex<double> Grating1D::getExtField(double x, double y)
01090 {
01091 complex<double> field = 0.;
01092 for (int m=0;m<nOrds;m++)
01093 {
01094 field+=tran[m]*exp(I*(gamma[m]*x*period+t[m]*y));
01095 }
01096 return field;
01097 }
01098
01103 void Grating1D::calcTransEmb(double embSize, int embOrd)
01104 {
01105 complex<double> c=0.,d=0., z;
01106 nOrdsEmb=embOrd+1-embOrd%2;
01107 embedding=embSize*gratingSize;
01108 embPeriod=2.*embedding+period;
01109
01110 if ((signed int) TransEmb.size()!=nOrdsEmb) TransEmb.resize(nOrdsEmb);
01111 if ((signed int) tEmb.size()!=nOrdsEmb) tEmb.resize(nOrdsEmb);
01112 if ((signed int) gammaEmb.size()!=nOrdsEmb)gammaEmb.resize(nOrdsEmb);
01113
01114
01115 complex<double> gamma0=n0*2.*M_PI/lambda*sin(theta);
01116 for(int i=0;i<nOrdsEmb;i++)
01117 {
01118 gammaEmb[i] = gamma0 + 2*M_PI*(i-(nOrdsEmb-1)/2)/embPeriod;
01119 }
01120
01121 complex<double> zb = pow(n3*2.*M_PI/lambda, 2);
01122 for(int i=0;i<nOrdsEmb;i++)
01123 {
01124 z = pow(gammaEmb[i],2);
01125 tEmb[i]=sqrt(zb - z);
01126 if(real(tEmb[i])==0.0 && imag(tEmb[i])==0.0)
01127 cout << "WARNING: t(" << minOrd+i <<") is zero" << endl;
01128 }
01129
01130 for(int m=0;m<nOrdsEmb;m++)
01131 {
01132 TransEmb[m]=0.;
01133 for (int p=0;p<nOrds;p++)
01134 {
01135 d=tran[p]*exp(-I*gamma[p]*embedding);
01136 if (gammaEmb[m]-gamma[p]!=0.) d*=(I/(gammaEmb[m]-gamma[p]))
01137 *(exp(I*(gamma[p]-gammaEmb[m])*(period+embedding))-exp(I*(gamma[p]-gammaEmb[m])*embedding));
01138 else d*=period;
01139 TransEmb[m]+=d;
01140 }
01141 TransEmb[m]/=embPeriod;
01142 }
01143
01144 }
01145
01152 complex<double> Grating1D::getEmbExtField(double x, double y)
01153 {
01154 complex<double> field = 0.;
01155 for (int m=0;m<nOrdsEmb;m++)
01156 {
01157 field+=TransEmb[m]*exp(I*(gammaEmb[m]*x*embPeriod+tEmb[m]*y));
01158 }
01159 return field;
01160 }
01161
01165 void Grating1D::writeParameters(string filename)
01166 {
01167 complex<double> c;
01168 ofstream out;
01169 out.open(filename.c_str());
01170 out << "Polarization = "<<polMode<<endl;
01171 out << "Lambda = "<<lambda<<endl;
01172 out << "Incident Angle = "<<getTheta()<<endl;
01173 out <<"Number of layers = "<<getnLayers()<<endl;
01174 out <<"Number of segments per layer = "<<layer[0].getnPermSegs()<<endl;
01175 out <<"Size of grating = "<<getGratingSize()<<endl;
01176 out <<"Thickness of layers ";
01177 for (int i=0;i<nLayers;i++) out <<layer[i].getThickness()<<" , ";
01178 out <<endl;
01179 out <<"Number of Fourier orders = "<<nOrds<<endl;
01180 out <<"Refractive index in region 1 = "<<n0<<endl;
01181 out <<"Refractive index in region 3 = "<<n3<<endl;
01182 out <<"Padding size = "<<getPadding()<<endl;
01183 out <<"Index of padding = "<<nPad<<endl;
01184 out <<"Size of embedding behind grating = "<<embPeriod/period<<endl;
01185 out <<"Number of Fourier coefficients for field calculations behind grating = "<<nOrdsEmb<<endl;
01186 double TransTot=0., RefTot=0.;
01187 for (int i=0;i<nOrds;i++)
01188 {
01189 TransTot+=getEfficiencyT(i);
01190 RefTot+=getEfficiencyR(i);
01191 }
01192 out <<"Total reflection = "<<RefTot<<endl;
01193 out <<"Total transmission = "<<TransTot<<endl;
01194
01195 out <<"\n\n Order , transmission efficiency , reflection efficiency\n";
01196 for (int j=0;j<nOrds;j++)
01197 {
01198 out << j-(nOrds-1)/2 <<" , "<<getEfficiencyT(j)<<" , "<<getEfficiencyR(j)<<endl;
01199 }
01200 out <<endl;
01201
01202
01203
01204
01205
01206
01207
01208 out << endl;
01209 out<<"Grating structure (real part)\n";
01210 out << endl;
01211 for (int i=0;i<nLayers;i++)
01212 {
01213 out << layer[i].getThickness()<< " , ";
01214 for (int j=nPadSegs;j<nPermSegs-nPadSegs;j++)
01215 {
01216 out << real(layer[i].getRefInd(j))<<" ";
01217 }
01218 out << endl;
01219 }
01220 out << endl;
01221 out<<"Grating structure (complex)\n";
01222 out << endl;
01223 for (int i=0;i<nLayers;i++)
01224 {
01225 out << layer[i].getThickness()<< " , ";
01226 for (int j=0;j<nPermSegs;j++)
01227 {
01228 out << layer[i].getRefInd(j)<<" ";
01229 }
01230 out << endl;
01231 }
01232 out << endl;
01233 out<<"Transition points\n";
01234 out << endl;
01235 for (int i=0;i<nLayers;i++)
01236 {
01237 for (int j=0;j<nPermSegs;j++)
01238 {
01239 out << layer[i].getTransPt(j)<<" ";
01240 }
01241 out << endl;
01242 }
01243 out << endl;
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294 out.close();
01295
01296 }
01297
01298 int Grating1D::nLines(string filename)
01299 {
01300 string line;
01301 int n=0;
01302 ifstream f;
01303 f.open(filename.c_str());
01304 while (!f.eof())
01305 {
01306 line="#";
01307 while (((line[0] == '#') || line.find_first_not_of(" \t\v\r\n")==string::npos)&& !f.eof()) {getline(f,line);}
01308 n++;
01309 }
01310 f.close();
01311 return --n;
01312 }
01313
01314 double Grating1D::delta(int i, int j) {return (i==j?1.:0.);}
