### added source files from https://metric.computational-photonics.eu/metric_unix.zip

parent 726a2cd7
 /* * METRIC-Application: * A three layer waveguide, attenuating core, * dispersion curves, perturbational computation of mode attenuation */ #include #include #include #include"metric.h" #define Pol TE // light polarization #define Wgpns 1.45 // substrate refractive index #define Wgpnf 1.98 // film refractive index #define Wgpnc 1.0 // cover: air double Wgpt = 0.5; // slab thickness /mum #define Wavel 1.55 // vacuum wavelength /mum #define WgpAtt 1.0e-4 // attenuation constant of the core, // intensity attenuation (in 1/mum = 10^4/cm) // of a plane wave in the bulk material #define MaxOrd 100 // highest order of a recorded mode #define Pbeg 0.01 // Parameter variations: minimum value, #define Pend 2.0 // maximum value, #define NumP 200 // number of points /* ------------------------------------------------------------------------ */ /* waveguide definition */ Waveguide wgdef() { Waveguide g(1); g.hx(0) = 0.0; g.hx(1) = Wgpt; g.n(0) = Wgpns; g.n(1) = Wgpnf; g.n(2) = Wgpnc; g.lambda = Wavel; return g; } /* ------------------------------------------------------------------------ */ /* calculate effective indices and relative attenuation levels (intensity attenuation of the guided modes relative to the intensity attenuation of plane waves in the bulk core material) versus the layer thickness parameter */ int main() { ModeArray ma; char name = "t___neff"; char attname = "t___att"; name = polchr(Pol); attname = polchr(Pol); Dmatrix neff(MaxOrd+1, NumP+1); Dmatrix att(MaxOrd+1, NumP+1); Dvector pval(NumP+1); int omax = -1; neff.init(0.0); fprintf(stderr, "\nDispersion curves:\n"); for(int pi=0; pi<=NumP; ++pi) { pval(pi) = Pbeg+((double)pi)/((double)NumP)*(Pend-Pbeg); Wgpt = pval(pi); // define the waveguide Waveguide wg = wgdef(); // the permittivity perturbation, loss in the core layer Perturbation pert = attenuation(WgpAtt, wg.n(1), wg.lambda, wg.layer(1)); // find the guided modes modeanalysis(wg, Pol, ma, 1); // store effective mode indices, evaluate the perturbation for(int j=0; j<=ma.num-1; ++j) { int o = ma(j).nodes(); if(o <= MaxOrd) { neff(o, pi) = ma(j).neff; // mode intensity attenuation, // relative to the bulk attenuation att(o, pi) = -ma(j).phaseshift(pert, FORW).im *2.0/WgpAtt; if(o > omax) omax = o; } } int s = NumP/60; if(pi % s == 0) fprintf(stderr, "."); } fprintf(stderr, "\n"); // write the data to files, individually per mode order, for ... Dvector oarg(NumP+1); Dvector oval(NumP+1); Dvector aval(NumP+1); for(int o=0; o<=omax; ++o) { int num = 0; name = dig10(o); name = dig1(o); attname = dig10(o); attname = dig1(o); for(int pi=0; pi<=NumP; ++pi) { if(neff(o, pi) > 0.0) { oarg(num) = pval(pi); oval(num) = neff(o, pi); aval(num) = att(o, pi); ++num; } } if(num >= 1) { // ... effective mode indices, toxyf(name, oarg.subvector(num, 0), oval.subvector(num, 0)); // ... and attenuation levels toxyf(attname, oarg.subvector(num, 0), aval.subvector(num, 0)); } } fprintf(stderr, "Ok.\n"); return 0; }
 /* * METRIC-Application: * Propagation of Gaussian beams in a homogeneous dielectric medium * - crossing & interference of wide beams originating from the * left and lower border */ #include #include #include #include"metric.h" #define Pol TE // light polarization #define FPn 1.45 // refractive index of the medium #define Wavel 1.55 // vacuum wavelength /mum #define CW 20.0 // computational window, x- and z-extension #define NumM 50 // number of spectral terms in the field discretizations #define DW 20.0 // display window, x- and z-extension /* ------------------------------------------------------------------------ */ /* specification of the incoming Gaussian beams */ #define BPwidth 3.0 // half with #define BPx0 -1.0 // center position #define BPtheta 10.0 // tilt angle, in degrees Gaussianbeam Beam1(Pol, BPwidth, BPx0, BPtheta, Wavel, FPn); Gaussianbeam Beam2(Pol, BPwidth, -BPx0, -BPtheta, Wavel, FPn); /* the computational domain, homogeneous space */ Circuit circuitdef() { Circuit cr; cr = Circuit(1, 1); cr.hx(0) = -CW/6.0; cr.hx(1) = CW/6.0; cr.hz(0) = -CW/6.0; cr.hz(1) = CW/6.0; cr.n(0, 0) = FPn; cr.n(0, 1) = FPn; cr.n(0, 2) = FPn; cr.n(1, 0) = FPn; cr.n(1, 1) = FPn; cr.n(1, 2) = FPn; cr.n(2, 0) = FPn; cr.n(2, 1) = FPn; cr.n(2, 2) = FPn; cr.lambda = Wavel; return cr; } /* ------------------------------------------------------------------------ */ /* QUEP simulation */ int main() { fprintf(stderr, "\nPropagation of Gaussian beams:\n"); fprintf(stderr, "------------------------------\n"); fprintf(stderr, "lambda = %g mum\n", Wavel); fprintf(stderr, "n = %g\n", FPn); fprintf(stderr, "w = %g mum\n", BPwidth); fprintf(stderr, "x0 = %g mum\n", BPx0); fprintf(stderr, "theta = %g^o\n", BPtheta); // the structure under investigation Circuit fsp = circuitdef(); // vertical and horizontal computational window Interval cw(-CW/2.0, CW/2.0); // QUEP setup QuepField qfld(fsp, Pol, cw, NumM, cw, NumM); // input: the beams as defined above, at the left and lower borders qfld.input(LEFT, Beam1, CC1); qfld.input(BOTTOM, Beam2, CC1); // QUEP solution of the Helmholtz problem qfld.quepsim(); double pl, pr, pb, pt; fprintf(stderr, "\nTotal incoming power:\n"); pl = qfld.Pin(LEFT); pr = qfld.Pin(RIGHT); pb = qfld.Pin(BOTTOM); pt = qfld.Pin(TOP); fprintf(stderr, "P_left = %g\n", pl); fprintf(stderr, "P_right = %g\n", pr); fprintf(stderr, "P_bottom = %g\n", pb); fprintf(stderr, "P_top = %g\n", pt); fprintf(stderr, "Sum P_in = %g\n", pl+pr+pb+pt); fprintf(stderr, "\nTotal outgoing power:\n"); pl = qfld.Pout(LEFT); pr = qfld.Pout(RIGHT); pb = qfld.Pout(BOTTOM); pt = qfld.Pout(TOP); fprintf(stderr, "P_left = %g\n", pl); fprintf(stderr, "P_right = %g\n", pr); fprintf(stderr, "P_bottom = %g\n", pb); fprintf(stderr, "P_top = %g\n", pt); fprintf(stderr, "Sum P_out = %g\n", pl+pr+pb+pt); char pc = polchr(Pol); Fcomp fc = principalcomp(Pol); Mlop_Print = NO; Mlop_Colour = YES; qfld.viewer(-0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 't', pc); qfld.plot(fc, MOD, -0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 't', pc); qfld.movie(-0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 30, 't', pc); fprintf(stderr, "\nOk.\n"); return 0; }
 /* * METRIC-Application: * Grating in a symmetric three layer waveguide, * (infinite sequence of equally spaced identical dielectric squares), * parameters: J.D. Joannopoulos et al., "Photonic Crystals: * Molding the flow of light", Princeton University Press, 2008; * band structure analysis, BEP. */ #include #include #include #include"metric.h" #define Pol TE // light polarization #define GPnc 1.0 // cover refractive index #define GPnf (sqrt(12.0)) // film refractive index #define GPp 1.0 // grating period /mum #define GPt 0.4 // film thickness /mum #define GPs 0.6 // width of the holes /mum double Wavel = 0.09; // vacuum wavelength /mum double CWbot = -5.123; // vertical computational window double CWtop = 5.123; // int NumMx = 70; // number of spectral terms in the field discretization, #define NumMxMin 20 // adapted to a minimum to ensure convergence of the eigensolver #define NSDTPWL 20 // number of initial (spectral) discretiz. terms per wavelength #define Pmin 0.01 // Parameter scan: range and step size #define Pmax 0.50 // #define Pstep 0.005 // /* periodic dielectric waveguide: one Period */ SegWgStruct grdef() { Waveguide slab(1); slab.hx(1) = GPt/2.0; slab.hx(0) = -GPt/2.0; slab.n(2) = GPnc; slab.n(1) = GPnf; slab.n(0) = GPnc; slab.lambda = Wavel; Waveguide etched(1); etched.hx(1) = GPt/2.0; etched.hx(0) = -GPt/2.0; etched.n(2) = GPnc; etched.n(1) = GPnc; etched.n(0) = GPnc; etched.lambda = Wavel; SegWgStruct gr(3); // the period extends from gr.hz(0) to gr.hz(3) double z = 0.0; int l = 0; gr(l) = slab; gr.hz(l) = z; ++l; gr(l) = slab; z += 0.5*(GPp-GPs); gr.hz(l) = z; ++l; gr(l) = etched; z += GPs; gr.hz(l) = z; ++l; gr(l) = slab; z += 0.5*(GPp-GPs); gr.hz(l) = z; ++l; gr(l) = slab; return gr; } int main() { for(double pv = Pmin; pv <= Pmax; pv += Pstep) { Wavel = GPp/pv; fprintf(stderr, "\n [lambda = %g, Lambda/lambda = %g] \n", Wavel, GPp/Wavel); // the grating, one period SegWgStruct gr = grdef(); // vertical computational window, adapt to vacuum wavelength double cwt = CWtop; double cwb = CWbot; if((cwt-cwb) < 10.0*Wavel) { double cw0 = 0.5*(cwt+cwb); cwt = cw0+5.0123*Wavel; cwb = cw0-5.0123*Wavel; } NumMx = (int) ceil((cwt-cwb)/Wavel*((double)NSDTPWL)); Interval vcw(cwb, cwt); //invoke the BEP-Floquet-Bloch solver FBMAXATTLVL = 0.01; FBVOPLVL = 1.0e-4; FBModeArray ma; fbmsolve(gr, Pol, vcw, NumMx, 1, ma, 0); // data output: FB wavenumbers and directional power levels for(int j=0; j<=ma.num-1; ++j) { apptoxyf("bands", ma(j).beta*GPp/2.0/PI, GPp/Wavel); if(ma(j).beta > 0) { apptoxyf("pf", GPp/Wavel, ma(j).fld.Pout(SBRIGHT)); apptoxyf("pb", GPp/Wavel, ma(j).fld.Pout(SBLEFT)); } } } fprintf(stderr, "\nOk.\n"); return 0; }
 /* * METRIC-Application: * Propagation of Gaussian beams in a homogeneous dielectric medium * - two narrow, rapidly sperading sources at the top border * "double-slit" configuration */ #include #include #include #include"metric.h" #define Pol TE // light polarization #define FPn 1.45 // refractive index of the medium #define Wavel 1.55 // vacuum wavelength /mum #define CW 20.0 // computational window, x- and z-extension #define NumM 50 // number of spectral terms in the field discretizations #define DW 20.0 // display window, x- and z-extension /* ------------------------------------------------------------------------ */ /* specification of the incoming Gaussian beams */ #define BPwidth 0.5 // half with #define BPx0 3.0 // center position #define BPtheta 0.0 // tilt angle, in degrees Gaussianbeam Beam1(Pol, BPwidth, BPx0, BPtheta, Wavel, FPn); Gaussianbeam Beam2(Pol, BPwidth, -BPx0, BPtheta, Wavel, FPn); /* the computational domain, homogeneous space */ Circuit circuitdef() { Circuit cr; cr = Circuit(1, 1); cr.hx(0) = -CW/6.0; cr.hx(1) = CW/6.0; cr.hz(0) = -CW/6.0; cr.hz(1) = CW/6.0; cr.n(0, 0) = FPn; cr.n(0, 1) = FPn; cr.n(0, 2) = FPn; cr.n(1, 0) = FPn; cr.n(1, 1) = FPn; cr.n(1, 2) = FPn; cr.n(2, 0) = FPn; cr.n(2, 1) = FPn; cr.n(2, 2) = FPn; cr.lambda = Wavel; return cr; } /* ------------------------------------------------------------------------ */ /* QUEP simulation */ int main() { fprintf(stderr, "\nPropagation of Gaussian beams:\n"); fprintf(stderr, "------------------------------\n"); fprintf(stderr, "lambda = %g mum\n", Wavel); fprintf(stderr, "n = %g\n", FPn); fprintf(stderr, "w = %g mum\n", BPwidth); fprintf(stderr, "x0 = %g mum\n", BPx0); fprintf(stderr, "theta = %g^o\n", BPtheta); // the structure under investigation Circuit fsp = circuitdef(); // vertical and horizontal computational window Interval cw(-CW/2.0, CW/2.0); // QUEP setup QuepField qfld(fsp, Pol, cw, NumM, cw, NumM); // input: the beams as defined above, from the top border qfld.input(TOP, Beam1, CC1); qfld.input(TOP, Beam2, CC1); // QUEP solution of the Helmholtz problem qfld.quepsim(); double pl, pr, pb, pt; fprintf(stderr, "\nTotal incoming power:\n"); pl = qfld.Pin(LEFT); pr = qfld.Pin(RIGHT); pb = qfld.Pin(BOTTOM); pt = qfld.Pin(TOP); fprintf(stderr, "P_left = %g\n", pl); fprintf(stderr, "P_right = %g\n", pr); fprintf(stderr, "P_bottom = %g\n", pb); fprintf(stderr, "P_top = %g\n", pt); fprintf(stderr, "Sum P_in = %g\n", pl+pr+pb+pt); fprintf(stderr, "\nTotal outgoing power:\n"); pl = qfld.Pout(LEFT); pr = qfld.Pout(RIGHT); pb = qfld.Pout(BOTTOM); pt = qfld.Pout(TOP); fprintf(stderr, "P_left = %g\n", pl); fprintf(stderr, "P_right = %g\n", pr); fprintf(stderr, "P_bottom = %g\n", pb); fprintf(stderr, "P_top = %g\n", pt); fprintf(stderr, "Sum P_out = %g\n", pl+pr+pb+pt); char pc = polchr(Pol); Fcomp fc = principalcomp(Pol); Mlop_Print = NO; Mlop_Colour = YES; qfld.viewer(-0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 't', pc); qfld.plot(fc, MOD, -0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 't', pc); qfld.movie(-0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 30, 't', pc); fprintf(stderr, "\nOk.\n"); return 0; }
 /* * METRIC-Application: * two adjacent parallel waveguides, * conventional coupled mode theory model */ #include #include #include #include"metric.h" #define Pol TM #define MPnb 1.45 /* background refractive index */ #define MPng 3.40 /* guiding regions: refractive index */ #define Wavel 1.55 /* vacuum wavelength / mum */ double GPw = 0.20; /* width of the waveguides / mum */ double GPg = 0.30; /* gap between the waveguides / mum */ /* Scan over gap widths: parameter interval */ #define MinG 0.1 #define MaxG 0.9 #define NumG 100 /* display window, animation of field propagation */ #define DWx0 -1.5 #define DWx1 1.5 #define DWz0 0.0 #define DWz1 15.0 /* ---------------------------------------------------------- */ /* a single waveguide */ Waveguide iowgdef() { Waveguide g(1); g.hx(0) = -GPw/2.0; g.hx(1) = GPw/2.0; g.n(0) = MPnb; g.n(1) = MPng; g.n(2) = MPnb; g.lambda = Wavel; return g; } /* the two parallel waveguides */ Waveguide cpwgdef() { Waveguide g(3); g.hx(0) = -(GPg/2.0+GPw); g.hx(1) = -(GPg/2.0); g.hx(2) = GPg/2.0; g.hx(3) = GPg/2.0+GPw; g.n(0) = MPnb; g.n(1) = MPng; g.n(2) = MPnb; g.n(3) = MPng; g.n(4) = MPnb; g.lambda = Wavel; return g; } /* ---------------------------------------------------------- */ /* something went wrong: stop the program */ void cmtstop(const char *e) { fprintf(stderr, "CMT: ERROR: %s.\n", e); exit(1); } /* ---------------------------------------------------------- */ /* coupled waveguide simulation */ int main() { int j; char pc = polchr(Pol); Fcomp fcp = principalcomp(Pol); ModeArray ma, basis; Mode m_m, m_p; // the input/output waveguides and the coupling segment Waveguide wg, cp; fprintf(stderr, "\nTwo parallel waveguides:\n"); fprintf(stderr, "----------------------\n\n"); fprintf(stderr, "Device parameters:\n"); fprintf(stderr, "lambda = %g mum\n", Wavel); fprintf(stderr, "n_b = %g\n", MPnb); fprintf(stderr, "n_g = %g\n", MPng); fprintf(stderr, "w = %g mum\n", GPw); fprintf(stderr, "g = %g mum\n", GPg); fprintf(stderr, "\n"); // display window for mode plots Interval display(-GPg-2.0*GPw, GPg+2.0*GPw); /* modes of the port waveguides */ wg = iowgdef(); fprintf(stderr, "Port modes:\n"); modeanalysis(wg, Pol, ma); if(ma.num != 1) cmtstop("Input guide is not single mode"); m_m = ma(0); fprintf(stderr, "T%c0: beta = %g /mum, n_eff = %g\n", pc, m_m.beta, m_m.neff); m_m.translate(-GPg/2.0-GPw/2.0); m_m.plot(fcp, ORG, display, 500, 'p', 'm', 'V'); m_p = ma(0); m_p.translate(GPg/2.0+GPw/2.0); m_p.plot(fcp, ORG, display, 500, 'p', 'p', 'V'); /* the CMT mode basis */ basis.clear(); basis.add(m_m); basis.add(m_p); /* compute supermode propagation constants and amplitude vectors */ Dmatrix sigma; Dvector smpc; Dmatrix smav; cp = cpwgdef(); CMTsetup(basis, cp, sigma, smpc, smav); fprintf(stderr,"Sigma:\n"); sigma.write(stderr); fprintf(stderr,"Smpc:\n"); smpc.write(stderr); fprintf(stderr,"Smav:\n"); smav.write(stderr); fprintf(stderr, "\nCMT-supermodes:\n"); fprintf(stderr, " T%c0: beta = %g /mum, n_eff = %g\n", pc, smpc(0), smpc(0)/2.0/PI*Wavel); fprintf(stderr, " T%c1: beta = %g /mum, n_eff = %g\n", pc, smpc(1), smpc(1)/2.0/PI*Wavel); fprintf(stderr, " CMT: L_c = %g mum\n", PI/(smpc(0)-smpc(1))); /* The exact solution: */ fprintf(stderr, "\nExact coupling:\n"); cp = cpwgdef(); fprintf(stderr, "\nExact supermodes:\n"); modeanalysis(cp, Pol, ma); for(j=0; j<=ma.num-1; ++j) { fprintf(stderr, " T%c%d: beta = %g /mum, n_eff = %g\n",