#ifndef BEPFLD_H #define BEPFLD_H /* * METRIC --- Mode expansion modeling in integrated optics / photonics * http://metric.computational-photonics.eu/ */ /* * bepfld.h * Bidirectional eigenmode expansion, sequences of arrays of modes */ class BepField { public: // the underlying sequence of waveguide segments SegWgStruct s; // initialize BepField(); BepField(SegWgStruct st); BepField(Circuit rc); BepField(Circuit rc, int opt); // opt = 1: optimize structure (default) // full BEP initialization: setup the container for the // bidirectional mode expansion & spectral discretization // rc: the structure under investigation, // p: fields of polarization p, // wx: vertical computational window, // Mx: number of expansion terms per vertical slice, // suppress log output, if quiet == 1 */ BepField(SegWgStruct rc, Polarization p, Interval wx, int Mx, int quiet); BepField(SegWgStruct rc, Polarization p, Interval wx, int Mx); BepField(Circuit rc, Polarization p, Interval wx, int Mx, int quiet); BepField(Circuit rc, Polarization p, Interval wx, int Mx); // opt = 1: optimize structure (default) BepField(SegWgStruct rc, int opt, Polarization p, Interval wx, int Mx, int quiet); BepField(Circuit rc, int opt, Polarization p, Interval wx, int Mx, int quiet); // destroy ~BepField(); // copy constructor BepField(const BepField& mef); // assignment BepField& operator=(const BepField& mef); // free allocated memory void clear(); // set all mode basis sets to ma void init(ModeArray ma); // input from FILE dat void read(FILE *dat); // input from file with default name void read(char ext0, char ext1); // output to FILE dat void write(FILE *dat) const; // output to file with default name void write(char ext0, char ext1) const; // subscripting: specify / read the mode basis sets on segment seg inline ModeArray& operator() (int seg) { return marr[seg]; } inline ModeArray operator() (int seg) const { return marr[seg]; } // expansion coefficients: specify amplitudes for segment seg // at left (pos = SBLEFT) or right (pos = SBRIGHT) segment boundary // ... for forward traveling fields void setf(int seg, SBorder pos, const Cvector& amp); // ... for backward traveling fields void setb(int seg, SBorder pos, const Cvector& amp); // expansion coefficients: retrieve amplitudes for segment seg // at left (pos = SBLEFT) or right (pos = SBRIGHT) segment boundary // ... for forward traveling fields Cvector f(int seg, SBorder pos) const; // ... for backward traveling fields Cvector b(int seg, SBorder pos) const; // --------------------------------------------------------------- // check the SegWgStruct for consistency with respect to // polarization, wavelength, array dimensions; exit, if inconsistent void check(); // ... & set global values for Polarization BF_pol; // polarization (obsolete, if BF_vM) double BF_lambda; // vacuum wavelength int BF_vM; // vTE and vTM modes Complex BF_ky; // with uniform wavenumber ky // --------------------------------------------------------------- // mode interference evaluation // all modes on a segment are assumed to belong to the same waveguide ! // field superposition at point x, z // cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ Complex field(Fcomp cp, double x, double z) const; // local field "strength", id: one of mE, mH, qE, qH, mW double lfs(FSid id, double x, double z) const; // local electromagnetic energy density at x, z double endens(double x, double z) const; // the six electric and magnetic components at x, z void emfield(double x, double z, Complex& ex, Complex& ey, Complex& ez, Complex& hx, Complex& hy, Complex& hz) const; // the total directional power of the mode superposition; // assumptions: // - mode orthorgonality, different modes of the same waveguide // - modes are normalized // - only forward or only backward propagating fields are present double dir_power(int seg, Propdir dir) const; // ... for vTE / vTM modes: power associated with polarization pol double dir_power(int seg, Polarization pol, Propdir dir) const; // the directional guided power of the mode superposition; // assumptions: // - mode orthorgonality, different modes of the same waveguide // - modes are normalized // - only forward or only backward propagating fields are present double dir_guided_power(int seg, Propdir dir) const; // ... for vTE / vTM modes: power associated with polarization pol double dir_guided_power(int seg, Polarization pol, Propdir dir) const; // --------------------------------------------------------------- // specification of input fields, additive: // mode m with amplitude A at border bdr void input(SBorder bdr, Complex A, const Mode& m); // mode of order o (number o) with amplitude A at border bdr void input(SBorder bdr, Complex A, int o); // tilted Gaussian beam gb, multiplied by amplitude A, at border bdr void input(SBorder bdr, const Gaussianbeam& gb, Complex A); // ... a normalized beam, unspecified phase void input(SBorder bdr, const Gaussianbeam& gb); // vM: amplitude of incoming mode of order o passing border bdr, // with polarization pol, evaluated at propagation coordinate p Complex Ain(SBorder bdr, Polarization pol, int o, double p) const; // --------------------------------------------------------------- // output amplitudes & powers for the full sequence, // amplitude of outgoing mode m at border bdr Complex Aout(SBorder bdr, const Mode& m) const; // amplitude of outgoing mode of order o (number o) at border bdr Complex Aout(SBorder bdr, int o) const; // ... for superpositions of vTE / vTM modes: // mode of polarzation pol, order o Complex Aout(SBorder bdr, Polarization pol, int o) const; // vM: amplitude of outgoing mode of order o passing border bdr, // with polarization pol, evaluated at propagation coordinate p Complex Aout(SBorder bdr, Polarization pol, int o, double p) const; // output power associated with mode m at border bdr double Pout(SBorder bdr, const Mode& m) const; // output power associated with mode of order o at border bdr double Pout(SBorder bdr, int o) const; // output power associated with mode // of Polarization pol, order o at border bdr (vM) double Pout(SBorder bdr, Polarization pol, int o) const; // guided output power crossing border bdr double Pgout(SBorder bdr) const; // guided output power associated with polarization pol (vM) // crossing border bdr double Pgout(SBorder bdr, Polarization pol) const; // total output power crossing border bdr (guided and nonguided fields) double Pout(SBorder bdr) const; // total output power associated with polarization pol // crossing border bdr (guided and nonguided fields) double Pout(SBorder bdr, Polarization pol) const; // --------------------------------------------------------------- // input power for the full sequence, // total input power crossing border bdr (guided and nonguided fields) double Pin(SBorder bdr) const; // --------------------------------------------------------------- // "modal probing": projection of a mode onto the local expansion // directional amplitude of probe mode m at position z Complex locModAmp(double z, const Mode& m, Propdir dir) const; // power associated with directional mode m at position z double locModPow(double z, const Mode& m, Propdir dir) const; // --------------------------------------------------------------- // compute the discretized bidirectional mode spectra: // modes of polarization p, Mx expansion terms per slice, // on a vertical computation window cw, // suppress log output, if quiet == 1 void discretize(Polarization p, Interval cw, int Mx, int quiet); // --------------------------------------------------------------- // bidirectional mode expansion solver: // given the structure with the accompanying mode expansion, // where input mode amplitudes are specified by setf(0, f), // find the remaining mode coefficients. // Assumption: no backwards traveling fields on the last segment void bepsim(); // bidirectional mode expansion solver: // given the structure with the accompanying mode expansion and // input mode amplitudes inp on segment 0, // find the amplitudes ref and trans of the reflected and transmitted // fields; fill the intermediate amplitudes. // Assumption: no backwards traveling fields on the last segment void bepsim(Cvector& inp, Cvector& ref, Cvector& trans); // bidirectional mode expansion solver: // given the structure with the accompanying mode expansion // find the mode reflectivity matrix ref // and transmittivity matrix trans // the mode amplitude vectors are not accessed // Assumption: no backwards traveling fields on the last segment void bepsim(Cmatrix& ref, Cmatrix& trans); // bidirectional mode expansion solver: // given the structure with the accompanying mode expansion, // where input mode amplitudes are specified by setf(0, f), // find the remaining mode coefficients. // Assumption: no backwards traveling fields on the last segment // tm_m, tm_p: already computed arrays of transfermatrices // for the juntions in this object, // as output from transfermats() void bepsim(Cmatrix *tm_m, Cmatrix *tm_p); // mode expansion solver: // given the structure with the accompanying mode expansion // find mode reflectivity and transmissivity matrices // kff, kfb, kbf, kbb, // that link the amplitudes b(0, SBRIGHT), f(s.nz+1, SBLEFT) // of outgoing modes to the incoming mode amplitudes // f(0, SBRIGHT), b(s.nz+1, SBLEFT) as // ( f(s.nz+1, SBLEFT) ) ( kff kfb ) ( f(0, SBRIGHT) ) // ( ) = ( ) ( ) // ( b(0, SBRIGHT) ) ( kbf kbb ) ( b(s.nz+1, SBLEFT) ) ; // the mode amplitude vectors are not accessed void bepsim(Cmatrix& kff, Cmatrix& kfb, Cmatrix& kbf, Cmatrix& kbb); // bidirectional mode expansion solver, for QUEP implementation // given the structure with the accompanying mode expansion, // initialize the BepField (amplitude arrays), // setup - transfer matrices for the input and output junction // - a matrix that connects amplitudes on the first and last // inner segment // - matrices that determine the amplitudes on all inner segments void quepsetup(Cmatrix& m0mbf, Cmatrix& m0mbb, Cmatrix& m0pff, Cmatrix& m0pfb, Cvector& t1, Cmatrix& d0, Cmatrix& d1, Cmatrix& mNmbf, Cmatrix& mNmbb, Cmatrix& mNpff, Cmatrix& mNpfb, Cvector& tN, Cmatrix& dN, Cmatrix& dNp, Cmatrix& kff, Cmatrix& kfb, Cmatrix& kbf, Cmatrix& kbb, Cmatrix *iff, Cmatrix *ifb, Cmatrix *ibf, Cmatrix *ibb, Cvector *tx); // QUEP implementation: overlaps of perpendiculary oriented mode // expansion fields Cmatrix crossoverlap(const ModeArray& mh, double zpos, Cvector& t1, Cvector& tN, Cmatrix& kff, Cmatrix& kfb, Cmatrix& kbf, Cmatrix& kbb, Cmatrix *iff, Cmatrix *ifb, Cmatrix *ibf, Cmatrix *ibb); // --------------------------------------------------------------- // vBEP / vQUEP implementation // prepare for the vTE/vTM mode representation, for given ky void vectorify(Complex ky); // compute the discretized bidirectional mode spectra: // vTE and vTM modes, uniform wavenumber ky, // 2*Mx expansion terms per vertical slice on a vertical // computational window cw, // type = 0: boundary conditions LIMD for TE and LIMN for TM modes, // type = 1: boundary conditions LIMN for TE and LIMD for TM modes void discretize(Complex ky, int type, Interval cw, int Mx, int quiet); // full vBEP initialization: setup the container for the // bidirectional mode expansion & spectral discretization // rc: the structure under investigation, // ky: fields of vTE and vTM polarizations, with y-wavenumber ky // wx: vertical computational window, // Mx: number of expansion terms per vertical slice, // suppress log output, if quiet == 1 */ BepField(SegWgStruct rc, Complex ky, Interval wx, int Mx, int quiet); BepField(SegWgStruct rc, Complex p, Interval wx, int Mx); BepField(Circuit rc, Complex ky, Interval wx, int Mx, int quiet); BepField(Circuit rc, Complex ky, Interval wx, int Mx); // ... with ky specified by the input field: // guided mode of polarization pol, order ord, // coming in at border bdr, with amplitude amp, // at angle of incidence theta (^o) BepField(SegWgStruct rc, SBorder bdr, Polarization pol, int ord, Complex amp, double theta, Interval wx, int Mx, int quiet); BepField(SegWgStruct rc, SBorder bdr, Polarization pol, int ord, Complex amp, double theta, Interval wx, int Mx); BepField(Circuit rc, SBorder bdr, Polarization pol, int ord, Complex amp, double theta, Interval wx, int Mx, int quiet); BepField(Circuit rc, SBorder bdr, Polarization pol, int ord, Complex amp, double theta, Interval wx, int Mx); // --------------------------------------------------------------- // preparation of physical field plots: // selection of the global phase / snapshot time 0 within one period // multiply the entire field by a phase factor exp(i phi), // the internal representation is modified accordingly void adjustphase(double phi); // adjust the global phase, such that a plot(cp, ORG) of the physical // basic field component (cp = EY for TE, cp = HY for TM) exhibits the // maximum field at position (x, z) void adjustphase(double x, double z); // adjust the global phase, such that a plot(cp, ORG) of the physical // basic field component (cp = EY for TE, cp = HY for TM) exhibits an // overall field maximum (-> time snapshot of an extremal state // in case of a configuration with partially standing waves) // ... determine the maximum on a regular mesh of numx * mumz points // on a reasonable computational window void adjustphase(int numx, int numz); // ... as before, with (coarse) defaults for numx, numz void adjustphase(); // --------------------------------------------------------------- // visualization: output to MATLAB .m files // evaluation of the mode interference pattern in the display window // xbeg <= x <= xend, zbeg <= z <= zend // npx, npz: number of plot points in the output mesh // ext0, ext1: filename id characters for the m.file // assumption: uniform polarization wavelength of all included modes // image plot corresponding to field component cp // cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ // foa: MOD, SQR, REP, IMP void plot(Fcomp cp, Afo foa, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const; // image plot of the relative phase // cp: one of EX, EY, EZ, HX, HY, HZ void phasemap(Fcomp cp, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const; // image plot of the field "strength", id: one of mE, mH, qE, qH, mW void plot(FSid id, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const; // electromagnetic energy density image void plot(double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const; // fancy style surface plot corresponding to component cp // cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ void fplot(Fcomp cp, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const; // animation of the light propagation // the frames show images of EY (TE) resp. HY (TM), at // ntfr equidistant times over one time period void movie(double xbeg, double xend, double zbeg, double zend, int npx, int npz, int ntfr, char ext0, char ext1) const; // animation of the light propagation, fancy style ;-) // the frames show surfaces of EY component (TE) resp. HY (TM), at // ntfr equidistant times over one time period void fmovie(double xbeg, double xend, double zbeg, double zend, int npx, int npz, int ntfr, char ext0, char ext1) const; // export full field data ("everything") into a viewer MATLAB file void viewer(double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const; private: ModeArray *marr; int nums; Cvector *afl; Cvector *afr; Cvector *abl; Cvector *abr; }; /* orthogonality matrix for a mode superposition, inverted */ Cmatrix orthomat_inv(ModeArray& ma); /* overlap matrix for two mode superpositions */ Cmatrix overlapmat(ModeArray& m1, ModeArray& m2); /* compute transfer matrices for a junction at z=0: m1: local modes for z<0 m2: local modes for z>0 m1ff - m2bb: output, overlap integrals between forward and backward propagating modes on both sides of the junction */ void transfermats(ModeArray& m1, ModeArray& m2, Cmatrix& m1ff, Cmatrix& m1fb, Cmatrix& m1bf, Cmatrix& m1bb, Cmatrix& m2ff, Cmatrix& m2fb, Cmatrix& m2bf, Cmatrix& m2bb); /* compute transfer matrices for a junction at z=0: m1: local modes for z<0 m2: local modes for z>0 mat1, mat2: output, overlap integrals between forward and backward propagating modes on both sides of the junction */ void transfermats(ModeArray& m1, ModeArray& m2, Cmatrix& mat1, Cmatrix& mat2); /* compute transfer matrices for a junction at z=0: m1: local modes for z<0 m2: local modes for z>0 d1, d2: given orthogonality matrices for mode sets m1, m2 mat1, mat2: output, overlap integrals between forward and backward propagating modes on both sides of the junction */ void transfermats(ModeArray& m1, ModeArray& m2, Cmatrix& d1, Cmatrix& d2, Cmatrix& mat1, Cmatrix& mat2); /* a long three segment coupler: m1: local modes for z<0 m2: local modes for 0