/* * METRIC --- Mode expansion modeling in integrated optics / photonics * http://metric.computational-photonics.eu/ */ /* * bepfld.cpp * Bidirectional eigenmode expansion, sequences of arrays of modes */ #include #include #include #include"inout.h" #include"complex.h" #include"matrix.h" #include"structure.h" #include"gengwed.h" #include"slamode.h" #include"slamarr.h" #include"slams.h" #include"bepfld.h" #include"quepfld.h" #include"matlvis.h" /* error message */ void bepfielderror(const char *m) { fprintf(stderr, "\nbepfld: %s.\n", m); exit(1); } /* ------------------------------------------------------------------- */ /* initialize */ BepField::BepField() { s = SegWgStruct(-1); marr = new ModeArray[s.nz+2]; afl = new Cvector[s.nz+2]; afr = new Cvector[s.nz+2]; abl = new Cvector[s.nz+2]; abr = new Cvector[s.nz+2]; nums = s.nz+2; for(int i=0; i<=s.nz+1; ++i) marr[i].clear(); BF_vM = 0; } BepField::BepField(SegWgStruct st) { s = st; marr = new ModeArray[s.nz+2]; afl = new Cvector[s.nz+2]; afr = new Cvector[s.nz+2]; abl = new Cvector[s.nz+2]; abr = new Cvector[s.nz+2]; nums = s.nz+2; for(int i=0; i<=s.nz+1; ++i) marr[i].clear(); BF_vM = 0; } BepField::BepField(Circuit rc) { marr = NULL; afl = afr = abl = abr = NULL; SegWgStruct st(rc); (*this) = BepField(st); } BepField::BepField(Circuit rc, int opt) { marr = NULL; afl = afr = abl = abr = NULL; SegWgStruct st(rc, opt); (*this) = BepField(st); } /* full BEP initialization: setup the container for the quadridirectional mode expansion, compute the 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 optimize structure, if opt == 1 */ BepField::BepField(Circuit rc, int opt, Polarization p, Interval wx, int Mx, int quiet) { marr = NULL; afl = afr = abl = abr = NULL; if(quiet != 1) { fprintf(stderr, "\n------------- BEP -------------------------------------------------- '04 ---\n"); fprintf(stderr, "T%c ", polCHR(p)); fprintf(stderr, "x in (%.10g, %.10g), Mx: %d\n", wx.x0, wx.x1, Mx); fprintf(stderr, "----------------------------------------------------------------------------\n"); fprintf(stderr, " Nx: %d\n", rc.nx); fprintf(stderr, " hx: "); for(int j=0; j<=rc.nx; ++j) fprintf(stderr, "%6.10g ", rc.hx(j)); fprintf(stderr, "\n"); fprintf(stderr, " Nz: %d\n", rc.nz); fprintf(stderr, " hz: "); for(int j=0; j<=rc.nz; ++j) fprintf(stderr, "%6.10g ", rc.hz(j)); fprintf(stderr, "\n"); if(rc.special) fprintf(stderr, " eps:\n"); else fprintf(stderr, " n:\n"); for(int i=rc.nx+1; i>=0; --i) { fprintf(stderr,"layer %d: ", i); for(int j=0; j<=rc.nz+1; ++j) fprintf(stderr, "%6.10g ", rc.n(i, j)); fprintf(stderr,"\n"); } fprintf(stderr, "lambda: %.10g k0: %g", rc.lambda, val_k0(rc.lambda)); if(rc.special) fprintf(stderr, " (!)\n"); else fprintf(stderr, "\n"); fprintf(stderr, "----------------------------------------------------------------------------\n"); } (*this) = BepField(rc, opt); discretize(p, wx, Mx, quiet); } BepField::BepField(Circuit rc, Polarization p, Interval wx, int Mx, int quiet) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc, 1, p, wx, Mx, 0); } BepField::BepField(Circuit rc, Polarization p, Interval wx, int Mx) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc, p, wx, Mx, 0); } BepField::BepField(SegWgStruct rc, Polarization p, Interval wx, int Mx, int quiet) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc.circuit(), p, wx, Mx, quiet); } BepField::BepField(SegWgStruct rc, Polarization p, Interval wx, int Mx) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc.circuit(), p, wx, Mx, 0); } BepField::BepField(SegWgStruct rc, int opt, Polarization p, Interval wx, int Mx, int quiet) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc.circuit(), opt, p, wx, Mx, quiet); } /* destroy */ BepField::~BepField() { if(marr != NULL) delete[] marr; if(afl != NULL) delete[] afl; if(afr != NULL) delete[] afr; if(abl != NULL) delete[] abl; if(abr != NULL) delete[] abr; marr = NULL; afl = NULL; afr = NULL; abl = NULL; abr = NULL; nums = 0; s.clear(); } /* copy constructor */ BepField::BepField(const BepField& mef) { s = mef.s; marr = new ModeArray[s.nz+2]; nums = s.nz+2; for(int i=0; i<=s.nz+1; ++i) marr[i] = mef.marr[i]; afl = new Cvector[s.nz+2]; afr = new Cvector[s.nz+2]; abl = new Cvector[s.nz+2]; abr = new Cvector[s.nz+2]; for(int i=0; i<=s.nz+1; ++i) afl[i] = mef.afl[i]; for(int i=0; i<=s.nz+1; ++i) afr[i] = mef.afr[i]; for(int i=0; i<=s.nz+1; ++i) abl[i] = mef.abl[i]; for(int i=0; i<=s.nz+1; ++i) abr[i] = mef.abr[i]; BF_pol = mef.BF_pol; BF_lambda = mef.BF_lambda; BF_vM = mef.BF_vM; BF_ky = mef.BF_ky; } /* assignment */ BepField& BepField::operator=(const BepField& mef) { if(this != &mef) { if(marr != NULL) delete[] marr; if(afl != NULL) delete[] afl; if(afr != NULL) delete[] afr; if(abl != NULL) delete[] abl; if(abr != NULL) delete[] abr; marr = NULL; afl = NULL; afr = NULL; abl = NULL; abr = NULL; nums = 0; s = mef.s; marr = new ModeArray[s.nz+2]; nums = s.nz+2; for(int i=0; i<=s.nz+1; ++i) marr[i] = mef.marr[i]; afl = new Cvector[s.nz+2]; afr = new Cvector[s.nz+2]; abl = new Cvector[s.nz+2]; abr = new Cvector[s.nz+2]; for(int i=0; i<=s.nz+1; ++i) afl[i] = mef.afl[i]; for(int i=0; i<=s.nz+1; ++i) afr[i] = mef.afr[i]; for(int i=0; i<=s.nz+1; ++i) abl[i] = mef.abl[i]; for(int i=0; i<=s.nz+1; ++i) abr[i] = mef.abr[i]; BF_pol = mef.BF_pol; BF_lambda = mef.BF_lambda; BF_vM = mef.BF_vM; BF_ky = mef.BF_ky; } return *this; } /* delete all segment entries */ void BepField::clear() { s.clear(); if(marr != NULL) delete[] marr; if(afl != NULL) delete[] afl; if(afr != NULL) delete[] afr; if(abl != NULL) delete[] abl; if(abr != NULL) delete[] abr; marr = NULL; afl = NULL; afr = NULL; abl = NULL; abr = NULL; nums = 0; } /* set mode fields on all segments to ma */ void BepField::init(ModeArray ma) { if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) bepfielderror("init: SegWgStruct cleared" ); if(nums < s.nz+2) bepfielderror("init: structure mismatch"); for(int i=0; i<=s.nz+1; ++i) marr[i] = ma; } /* subscripting ModeArray& BepField::operator() (int seg) { if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) bepfielderror("access: SegWgStruct cleared" ); if(seg < 0) bepfielderror("access ( < 0)"); if(seg > s.nz+1) bepfielderror("access ( > s.nz+1)"); if(seg >= nums) bepfielderror("access: struct mismatch"); return marr[seg]; } ModeArray BepField::operator() (int seg) const { if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) bepfielderror("access: SegWgStruct cleared" ); if(seg < 0) bepfielderror("access ( < 0)"); if(seg > s.nz+1) bepfielderror("access ( > s.nz+1)"); if(seg >= nums) bepfielderror("access: struct mismatch"); return marr[seg]; } */ /* input from FILE dat */ void BepField::read(FILE *dat) { if(dat == stderr || dat == stdout) { bepfielderror("BepField: read"); } else { if(marr != NULL) delete[] marr; if(afl != NULL) delete[] afl; if(afr != NULL) delete[] afr; if(abl != NULL) delete[] abl; if(abr != NULL) delete[] abr; marr = NULL; afl = NULL; afr = NULL; abl = NULL; abr = NULL; nums = 0; s.read(dat); if(s.nz<-1) bepfielderror("read: s.nz < -1"); marr = new ModeArray[s.nz+2]; nums = s.nz+2; for(int i=0; i<=s.nz+1; ++i) (marr[i]).read(dat); afl = new Cvector[s.nz+2]; afr = new Cvector[s.nz+2]; abl = new Cvector[s.nz+2]; abr = new Cvector[s.nz+2]; for(int i=0; i<=s.nz+1; ++i) (afl[i]).read(dat); for(int i=0; i<=s.nz+1; ++i) (afr[i]).read(dat); for(int i=0; i<=s.nz+1; ++i) (abl[i]).read(dat); for(int i=0; i<=s.nz+1; ++i) (abr[i]).read(dat); } check(); return; } /* output to FILE dat */ void BepField::write(FILE *dat) const { if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) bepfielderror("write: SegWgStruct cleared" ); comment("BepField", dat); s.write(dat); for(int i=0; i<=s.nz+1; ++i) (marr[i]).write(dat); comment("fl's", dat); for(int i=0; i<=s.nz+1; ++i) (afl[i]).write(dat); comment("fr's", dat); for(int i=0; i<=s.nz+1; ++i) (afr[i]).write(dat); comment("bl's", dat); for(int i=0; i<=s.nz+1; ++i) (abl[i]).write(dat); comment("br's", dat); for(int i=0; i<=s.nz+1; ++i) (abr[i]).write(dat); return; } /* input from default file */ void BepField::read(char ext0, char ext1) { char name[10]; FILE *dat; name[0] = 's'; name[1] = 'l'; name[2] = 'a'; name[3] = ext0; name[4] = ext1; name[5] = '.'; name[6] = 'm'; name[7] = 'e'; name[8] = 'f'; name[9] = 0; fprintf(stderr, "<< %s\n", name); dat = fopen(name, "r"); read(dat); fclose(dat); return; } /* output to default file */ void BepField::write(char ext0, char ext1) const { char name[10]; FILE *dat; name[0] = 's'; name[1] = 'l'; name[2] = 'a'; name[3] = ext0; name[4] = ext1; name[5] = '.'; name[6] = 'm'; name[7] = 'e'; name[8] = 'f'; name[9] = 0; fprintf(stderr, ">> %s\n", name); dat = fopen(name, "w+"); write(dat); fclose(dat); return; } /* expansion coefficients: specify amplitudes of forward traveling fields on segment seg */ void BepField::setf(int seg, SBorder pos, const Cvector& amp) { if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) bepfielderror("setf: SegWgStruct cleared" ); if(seg < 0) bepfielderror("setf ( < 0)"); if(seg > s.nz+1) bepfielderror("setf ( > s.nz+1)"); if(marr[seg].num != amp.nel) bepfielderror("setf: dimensions"); if(pos != SBLEFT && pos != SBRIGHT) bepfielderror("setf: pos"); if(seg >= nums) bepfielderror("setf: struct mismatch"); if(pos == SBLEFT) afl[seg] = amp; else afr[seg] = amp; return; } /* expansion coefficients: specify amplitudes of backward traveling fields on segment seg */ void BepField::setb(int seg, SBorder pos, const Cvector& amp) { if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) bepfielderror("setb: SegWgStruct cleared" ); if(seg < 0) bepfielderror("setb ( < 0)"); if(seg > s.nz+1) bepfielderror("setb ( > s.nz+1)"); if(marr[seg].num != amp.nel) bepfielderror("setb: dimensions "); if(pos != SBLEFT && pos != SBRIGHT) bepfielderror("setb: pos"); if(seg >= nums) bepfielderror("setb: struct mismatch"); if(pos == SBLEFT) abl[seg] = amp; else abr[seg] = amp; return; } /* expansion coefficients: retrieve amplitudes of forward traveling fields on segment seg */ Cvector BepField::f(int seg, SBorder pos) const { if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) bepfielderror("f: SegWgStruct cleared" ); if(seg < 0) bepfielderror("f ( < 0)"); if(seg > s.nz+1) bepfielderror("f ( > s.nz+1)"); if(pos != SBLEFT && pos != SBRIGHT) bepfielderror("f: pos"); if(seg >= nums) bepfielderror("f: struct mismatch"); if(pos == SBLEFT) return afl[seg]; else return afr[seg]; } /* expansion coefficients: retrieve amplitudes of backward traveling fields on segment seg */ Cvector BepField::b(int seg, SBorder pos) const { if(marr==NULL || afl==NULL || afr==NULL || abl==NULL || abr==NULL) bepfielderror("b: SegWgStruct cleared" ); if(seg < 0) bepfielderror("b ( < 0)"); if(seg > s.nz+1) bepfielderror("b ( > s.nz+1)"); if(pos != SBLEFT && pos != SBRIGHT) bepfielderror("b: pos"); if(seg >= nums) bepfielderror("b: struct mismatch"); if(pos == SBLEFT) return abl[seg]; else return abr[seg]; } /* ------------------------------------------------------------------- */ /* check the SegWgStruct for consistency with respect to polarization, wavelength, array dimensions; exit, if inconsistent */ void BepField::check() { if(nums != s.nz+2) bepfielderror("check: general array dimensions"); /* if(marr == NULL) bepfielderror("check: SegWgStruct cleared (marr)"); if(afl == NULL) bepfielderror("check: SegWgStruct cleared (afl)" ); if(afr == NULL) bepfielderror("check: SegWgStruct cleared (afr)" ); if(abl == NULL) bepfielderror("check: SegWgStruct cleared (abl)" ); if(abr == NULL) bepfielderror("check: SegWgStruct cleared (abr)" ); for(int i=0; i<=s.nz+1; ++i) { if(afl[i].nel != marr[i].num) bepfielderror("check: array dimensions (afl)"); if(afr[i].nel != marr[i].num) bepfielderror("check: array dimensions (afr)"); if(abl[i].nel != marr[i].num) bepfielderror("check: array dimensions (abl)"); if(abr[i].nel != marr[i].num) bepfielderror("check: array dimensions (abr)"); if(afl[i].nel <= 0) bepfielderror("check: afl.nel == 0"); if(afr[i].nel <= 0) bepfielderror("check: afr.nel == 0"); if(abl[i].nel <= 0) bepfielderror("check: abl.nel == 0"); if(abr[i].nel <= 0) bepfielderror("check: abr.nel == 0"); } */ BF_pol = (marr[0])(0).pol; BF_lambda = (marr[0])(0).wg.lambda; BF_vM = (marr[0])(0).vM; BF_ky = (marr[0])(0).ky; for(int i=0; i<=s.nz+1; ++i) { int mni = marr[i].num; Mode mij; for(int j=0; j<=mni-1; ++j) { mij = (marr[i])(j); if(mij.wg.lambda != BF_lambda) bepfielderror("check: wavelengths mixed"); if(BF_vM) { if(mij.vM != BF_vM) bepfielderror("check: nonuniform vM"); if((mij.ky-BF_ky).sqabs() > SLAMS_BQTOL) bepfielderror("check: nonuniform ky"); } else { if(mij.pol != BF_pol) bepfielderror("check: polarizations mixed"); } } } return; } /* ------------------------------------------------------------------- */ /* 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 BepField::field(Fcomp cp, double x, double z) const { if(s.nz == -1) return marr[0].field(afl[0], abl[0], cp, x, z); int zi = s.segidx(z); if(zi <= 0) return marr[0].field(afr[0], abr[0], cp, x, z-s.hz(0)); if(zi >= s.nz+1) return marr[s.nz+1].field(afl[s.nz+1], abl[s.nz+1], cp, x, z-s.hz(s.nz)); return marr[zi].field(afl[zi], abr[zi], cp, x, z, s.hz(zi-1), s.hz(zi)); } /* local field "strength", id: one of mE, mH, qE, qH, mW */ double BepField::lfs(FSid id, double x, double z) const { int zi = s.segidx(z); if(zi <= 0) return marr[0].lfs(afr[0], abr[0], id, x, z-s.hz(0)); if(zi >= s.nz+1) return marr[s.nz+1].lfs(afl[s.nz+1], abl[s.nz+1], id, x, z-s.hz(s.nz)); return marr[zi].lfs(afl[zi], abr[zi], id, x, z, s.hz(zi-1), s.hz(zi)); } /* local electromagnetic energy density at x, z */ double BepField::endens(double x, double z) const { return lfs(mW, x, z); } /* the six electric and magnetic components at x, z */ void BepField::emfield(double x, double z, Complex& ex, Complex& ey, Complex& ez, Complex& hx, Complex& hy, Complex& hz) const { if(s.nz == -1) { marr[0].emfield(afl[0], abl[0], x, z, ex, ey, ez, hx, hy, hz); return; } int zi = s.segidx(z); if(zi <= 0) { marr[0].emfield(afr[0], abr[0], x, z-s.hz(0), ex, ey, ez, hx, hy, hz); return; } if(zi >= s.nz+1) { marr[s.nz+1].emfield(afl[s.nz+1], abl[s.nz+1], x, z-s.hz(s.nz), ex, ey, ez, hx, hy, hz); return; } marr[zi].emfield(afl[zi], abr[zi], x, z, s.hz(zi-1), s.hz(zi), ex, ey, ez, hx, hy, hz); return; } /* the total directional power of the mode superposition on segment seg; assumptions: - mode orthorgonality, different modes of the same waveguide - modes are normalized - only forward or only backward propagating fields are present */ double BepField::dir_power(int seg, Propdir dir) const { if(BF_vM) { double a, p = 0.0; if(seg < 0) seg = 0; if(seg > s.nz+1) seg = s.nz+1; ModeArray ma = marr[seg]; int n = ma.num; switch(dir) { case FORW: for(int i=0; i<=n-1; ++i) { if(ma(i).cVppt == PROPAG) { if(seg == 0) a = afr[0](i).sqabs(); else a = afl[seg](i).sqabs(); p += a; } } return p; break; case BACK: for(int i=0; i<=n-1; ++i) { if(ma(i).cVppt == PROPAG) { if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs(); else a = abr[seg](i).sqabs(); p += a; } } return p; break; default: bepfielderror("dir_power, vM: dir"); break; } return 0.0; } double a, p = 0.0; if(seg < 0) seg = 0; if(seg > s.nz+1) seg = s.nz+1; ModeArray ma = marr[seg]; int n = ma.num; switch(dir) { case FORW: for(int i=0; i<=n-1; ++i) { if(ma(i).ppt == PROPAG) { if(seg == 0) a = afr[0](i).sqabs(); else a = afl[seg](i).sqabs(); p += a; } } return p; break; case BACK: for(int i=0; i<=n-1; ++i) { if(ma(i).ppt == PROPAG) { if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs(); else a = abr[seg](i).sqabs(); p += a; } } return p; break; default: bepfielderror("dir_power: dir"); break; } return 0.0; } /* ... modes of polarization pol (vM) */ double BepField::dir_power(int seg, Polarization pol, Propdir dir) const { if(BF_vM == 0) return dir_power(seg, dir); double a, p = 0.0; if(seg < 0) seg = 0; if(seg > s.nz+1) seg = s.nz+1; ModeArray ma = marr[seg]; int n = ma.num; switch(dir) { case FORW: for(int i=0; i<=n-1; ++i) { if(ma(i).pol == pol && ma(i).cVppt == PROPAG) { if(seg == 0) a = afr[0](i).sqabs(); else a = afl[seg](i).sqabs(); p += a; } } return p; break; case BACK: for(int i=0; i<=n-1; ++i) { if(ma(i).pol == pol && ma(i).cVppt == PROPAG) { if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs(); else a = abr[seg](i).sqabs(); p += a; } } return p; break; default: bepfielderror("dir_power, pol: dir"); break; } return 0.0; } /* the guided directional power of the mode superposition on segment seg; assumptions: - mode orthorgonality, different modes of the same waveguide - modes are normalized - only forward or only backward propagating fields are present */ double BepField::dir_guided_power(int seg, Propdir dir) const { if(BF_vM) { double a, p = 0.0; if(seg < 0) seg = 0; if(seg > s.nz+1) seg = s.nz+1; ModeArray ma = marr[seg]; int n = ma.num; switch(dir) { case FORW: for(int i=0; i<=n-1; ++i) { if(ma(i).npcB > 0.0 && ma(i).cVppt == PROPAG) { if(seg == 0) a = afr[0](i).sqabs(); else a = afl[seg](i).sqabs(); p += a; } } return p; break; case BACK: for(int i=0; i<=n-1; ++i) { if(ma(i).npcB > 0.0 && ma(i).cVppt == PROPAG) { if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs(); else a = abr[seg](i).sqabs(); p += a; } } return p; break; default: bepfielderror("cVdir_guided_power: dir"); break; } return 0.0; } double a, p = 0.0; if(seg < 0) seg = 0; if(seg > s.nz+1) seg = s.nz+1; ModeArray ma = marr[seg]; int n = ma.num; switch(dir) { case FORW: for(int i=0; i<=n-1; ++i) { if(ma(i).npcB > 0.0 && ma(i).ppt == PROPAG) { if(seg == 0) a = afr[0](i).sqabs(); else a = afl[seg](i).sqabs(); p += a; } } return p; break; case BACK: for(int i=0; i<=n-1; ++i) { if(ma(i).npcB > 0.0 && ma(i).ppt == PROPAG) { if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs(); else a = abr[seg](i).sqabs(); p += a; } } return p; break; default: bepfielderror("dir_guided_power: dir"); break; } return 0.0; } /* ... of polarization pol (vM) */ double BepField::dir_guided_power(int seg, Polarization pol, Propdir dir) const { if(BF_vM == 0) return dir_guided_power(seg, dir); double a, p = 0.0; if(seg < 0) seg = 0; if(seg > s.nz+1) seg = s.nz+1; ModeArray ma = marr[seg]; int n = ma.num; switch(dir) { case FORW: for(int i=0; i<=n-1; ++i) { if(ma(i).pol == pol && ma(i).npcB > 0.0 && ma(i).cVppt == PROPAG) { if(seg == 0) a = afr[0](i).sqabs(); else a = afl[seg](i).sqabs(); p += a; } } return p; break; case BACK: for(int i=0; i<=n-1; ++i) { if(ma(i).pol == pol && ma(i).npcB > 0.0 && ma(i).cVppt == PROPAG) { if(seg == s.nz+1) a = abl[s.nz+1](i).sqabs(); else a = abr[seg](i).sqabs(); p += a; } } return p; break; default: bepfielderror("dir_guided_power, pol: dir"); break; } return 0.0; } /* --- specification of input fields, additive ----------------------- */ /* mode m with amplitude A at border bdr */ void BepField::input(SBorder bdr, Complex A, const Mode& m) { if(m.ppt != PROPAG) bepfielderror("input, EVANESC mode"); Cvector amp; ModeArray ma; int n; switch(bdr) { case SBRIGHT: amp = b(s.nz+1, SBLEFT); ma = marr[s.nz+1]; n = ma.num; for(int j=0; j<=n-1; ++j) { if(ma(j).ppt == PROPAG) amp(j) -= overlap(ma(j), BACK, m, BACK); } setb(s.nz+1, SBLEFT, amp); return; case SBLEFT: amp = f(0, SBRIGHT); ma = marr[0]; n = ma.num; for(int j=0; j<=n-1; ++j) { if(ma(j).ppt == PROPAG) amp(j) += overlap(ma(j), FORW, m, FORW); } setf(0, SBRIGHT, amp); return; } return; } /* mode of order o with amplitude A at border bdr */ void BepField::input(SBorder bdr, Complex A, int o) { Cvector amp; switch(bdr) { case SBRIGHT: amp = b(s.nz+1, SBLEFT); if(o<0 || o>amp.nel-1) bepfielderror("input: order, SBRIGHT"); if((marr[s.nz+1])(o).ppt == EVANESC) bepfielderror("input, EVANESC mode, SBRIGHT"); amp(o) += A; setb(s.nz+1, SBLEFT, amp); return; case SBLEFT: amp = f(0, SBRIGHT); if(o<0 || o>amp.nel-1) bepfielderror("input: order, SBLEFT"); if((marr[0])(o).ppt == EVANESC) bepfielderror("input, EVANESC mode, SBLEFT"); amp(o) += A; setf(0, SBRIGHT, amp); return; } return; } /* Gaussian beam gb, multiplied by amplitude A, at border bdr */ void BepField::input(SBorder bdr, const Gaussianbeam& gb, Complex A) { Cvector amp; ModeArray ma; int n; switch(bdr) { case SBRIGHT: amp = b(s.nz+1, SBLEFT); ma = marr[s.nz+1]; n = ma.num; for(int j=0; j<=n-1; ++j) { if(ma(j).ppt == PROPAG) amp(j) -= A*(ma(j)).overlap(gb, BACK); } setb(s.nz+1, SBLEFT, amp); return; case SBLEFT: amp = f(0, SBRIGHT); ma = marr[0]; n = ma.num; for(int j=0; j<=n-1; ++j) { if(ma(j).ppt == PROPAG) amp(j) += A*(ma(j)).overlap(gb, FORW); } setf(0, SBRIGHT, amp); return; } return; } /* ... a normalized beam, unspecified phase */ void BepField::input(SBorder bdr, const Gaussianbeam& gb) { input(bdr, gb, CC1); return; } /* vM: amplitude of incoming mode of order o passing border bdr, with polarization pol, evaluated at propagation coordinate p */ Complex BepField::Ain(SBorder bdr, Polarization pol, int o, double p) const { if(BF_vM == 0) bepfielderror("Ain, vM = 0, pol, pos: not implemented"); int num, mx, mi; double zw; Complex a, ef; switch(bdr) { case SBRIGHT: num = marr[s.nz+1].num; if(num/2 != ((int) ((double) num)/2.0)) bepfielderror("Ain, pol, o, p; confused"); mx = num/2; if(o<0 || o>mx-1) bepfielderror("Ain, pol, o, p; RIGHT: val o"); if(pol == TE) mi = o; else mi = mx+o; if(marr[s.nz+1](mi).pol != pol) bepfielderror("Ain, pol, o, p; RIGHT: pol"); if(marr[s.nz+1](mi).ord != o) bepfielderror("Ain, pol, o, p; RIGHT: ord"); a = (b(s.nz+1, SBLEFT))(mi); zw = s.hz(s.nz); ef = (marr[s.nz+1](mi).efac(BACK, p-zw)); return a*ef; case SBLEFT: num = marr[0].num; if(num/2 != ((int) ((double) num)/2.0)) bepfielderror("Ain, pol, o, p; confused"); mx = num/2; if(o<0 || o>mx-1) bepfielderror("Ain, pol, o, p; LEFT: val o"); if(pol == TE) mi = o; else mi = mx+o; if(marr[0](mi).pol != pol) bepfielderror("Ain, pol, o, p; LEFT: pol"); if(marr[0](mi).ord != o) bepfielderror("Ain, pol, o, p; LEFT: ord"); a = (f(0, SBRIGHT))(mi); zw = s.hz(0); ef = (marr[0](mi).efac(FORW, p-zw)); return a*ef; } return CC0; } /* --- output amplitudes & powers, assuming normalized modes --- */ /* amplitude of outgoing mode m at border bdr */ Complex BepField::Aout(SBorder bdr, const Mode& m) const { if(BF_vM) { if(m.vM == 0) bepfielderror("Aout, vM, m: vM = 0"); if(m.cVppt != PROPAG) bepfielderror("Aout, vM, m: EVANESC mode"); Complex a = CC0; Cvector amp; ModeArray ma; int n; switch(bdr) { case SBRIGHT: amp = f(s.nz+1, SBLEFT); ma = marr[s.nz+1]; n = ma.num; for(int j=0; j<=n-1; ++j) { if(ma(j).cVppt == PROPAG) a += amp(j)*overlap(m, FORW, ma(j), FORW); } return a/overlap(m, FORW, m, FORW); case SBLEFT: amp = b(0, SBRIGHT); ma = marr[0]; n = ma.num; for(int j=0; j<=n-1; ++j) { if(ma(j).cVppt == PROPAG) a -= amp(j)*overlap(m, BACK, ma(j), BACK); } return a/overlap(m, BACK, m, BACK); } return CC0; } if(m.ppt != PROPAG) bepfielderror("Aout, EVANESC mode"); Complex a = CC0; Cvector amp; ModeArray ma; int n; switch(bdr) { case SBRIGHT: amp = f(s.nz+1, SBLEFT); ma = marr[s.nz+1]; n = ma.num; for(int j=0; j<=n-1; ++j) { if(ma(j).ppt == PROPAG) a += amp(j)*overlap(m, FORW, ma(j), FORW); } return a/m.power(); case SBLEFT: amp = b(0, SBRIGHT); ma = marr[0]; n = ma.num; for(int j=0; j<=n-1; ++j) { if(ma(j).ppt == PROPAG) a -= amp(j)*overlap(m, BACK, ma(j), BACK); } return a/m.power(); } return CC0; } /* amplitude of outgoing mode of order o at border bdr */ Complex BepField::Aout(SBorder bdr, int o) const { if(BF_vM) bepfielderror("Aout, o: BF_vM"); switch(bdr) { case SBRIGHT: if(o<0 || o>marr[s.nz+1].num-1) bepfielderror("Aout: order, RIGHT"); return (f(s.nz+1, SBLEFT))(o); case SBLEFT: if(o<0 || o>marr[0].num-1) bepfielderror("Aout: order, LEFT"); return (b(0, SBRIGHT))(o); } return CC0; } /* amplitude of outgoing mode of polarzation pol, order o at border bdr (vM) */ Complex BepField::Aout(SBorder bdr, Polarization pol, int o) const { if(BF_vM == 0) bepfielderror("Aout, vM = 0, pol: not implemented"); switch(bdr) { case SBRIGHT: return Aout(bdr, pol, o, s.hz(s.nz)); case SBLEFT: return Aout(bdr, pol, o, s.hz(0)); } return CC0; } /* vM: amplitude of outgoing mode of order o passing border bdr, with polarization pol, evaluated at propagation coordinate p */ Complex BepField::Aout(SBorder bdr, Polarization pol, int o, double p) const { if(BF_vM == 0) bepfielderror("Aout, vM = 0, pol, pos: not implemented"); int num, mx, mi; double zw; Complex a, ef; switch(bdr) { case SBRIGHT: num = marr[s.nz+1].num; if(num/2 != ((int) ((double) num)/2.0)) bepfielderror("Aout, pol, o, p; confused"); mx = num/2; if(o<0 || o>mx-1) bepfielderror("Aout, pol, o, p; RIGHT: val o"); if(pol == TE) mi = o; else mi = mx+o; if(marr[s.nz+1](mi).pol != pol) bepfielderror("Aout, pol, o, p; RIGHT: pol"); if(marr[s.nz+1](mi).ord != o) bepfielderror("Aout, pol, o, p; RIGHT: ord"); if(marr[s.nz+1](mi).cVppt != PROPAG) return CC0; a = (f(s.nz+1, SBLEFT))(mi); zw = s.hz(s.nz); ef = (marr[s.nz+1](mi).efac(FORW, p-zw)); return a*ef; case SBLEFT: num = marr[0].num; if(num/2 != ((int) ((double) num)/2.0)) bepfielderror("Aout, pol, o, p; confused"); mx = num/2; if(o<0 || o>mx-1) bepfielderror("Aout, pol, o, p; LEFT: val o"); if(pol == TE) mi = o; else mi = mx+o; if(marr[0](mi).pol != pol) bepfielderror("Aout, pol, o, p; LEFT: pol"); if(marr[0](mi).ord != o) bepfielderror("Aout, pol, o, p; LEFT: ord"); if(marr[0](mi).cVppt != PROPAG) return CC0; a = (b(0, SBRIGHT))(mi); zw = s.hz(0); ef = (marr[0](mi).efac(BACK, p-zw)); return a*ef; } return CC0; } /* output power associated with mode m at border bdr */ double BepField::Pout(SBorder bdr, const Mode& m) const { return Aout(bdr, m).sqabs(); } /* output power associated with mode of order o at border bdr */ double BepField::Pout(SBorder bdr, int o) const { return Aout(bdr, o).sqabs(); } /* output power associated with mode of Polarization pol, order o at border bdr (vM) */ double BepField::Pout(SBorder bdr, Polarization pol, int o) const { return Aout(bdr, pol, o).sqabs(); } /* guided output power crossing border bdr */ double BepField::Pgout(SBorder bdr) const { switch(bdr) { case SBRIGHT: return dir_guided_power(s.nz+1, FORW); case SBLEFT: return dir_guided_power(0, BACK); } return 0.0; } /* guided output power associated with polarization pol (vM) crossing border bdr */ double BepField::Pgout(SBorder bdr, Polarization pol) const { switch(bdr) { case SBRIGHT: return dir_guided_power(s.nz+1, pol, FORW); case SBLEFT: return dir_guided_power(0, pol, BACK); } return 0.0; } /* total output power crossing border bdr (guided and nonguided fields) */ double BepField::Pout(SBorder bdr) const { switch(bdr) { case SBRIGHT: return dir_power(s.nz+1, FORW); case SBLEFT: return dir_power(0, BACK); } return 0.0; } /* total output power associated with polarization pol crossing border bdr (guided and nonguided fields) */ double BepField::Pout(SBorder bdr, Polarization pol) const { switch(bdr) { case SBRIGHT: return dir_power(s.nz+1, pol, FORW); case SBLEFT: return dir_power(0, pol, BACK); } return 0.0; } /* ------------------------------------------------------------------- */ /* total input power crossing border bdr (guided and nonguided fields) */ double BepField::Pin(SBorder bdr) const { switch(bdr) { case SBLEFT: return dir_power(0, FORW); case SBRIGHT: return dir_power(s.nz+1, BACK); } return 0.0; } /* --- "modal probing": projection of a mode onto the local expansion ---- */ /* directional amplitude of probe mode m at position z */ Complex BepField::locModAmp(double z, const Mode& m, Propdir dir) const { if(BF_vM) bepfielderror("locModAmp: vM not implemented"); Complex a = CC0; Cvector famp, bamp; ModeArray ma; int n; int si = s.segidx(z); Interval seg = s.segment(z); if(si <= 0) { famp = f(0, SBRIGHT); bamp = b(0, SBRIGHT); ma = marr[0]; n = ma.num; for(int j=0; j<=n-1; ++j) { a += famp(j)*overlap(m, dir, ma(j), FORW) *ma(j).efac(FORW, z-seg.x1); a += bamp(j)*overlap(m, dir, ma(j), BACK) *ma(j).efac(BACK, z-seg.x1); } return a/m.power(); } if(si >= s.nz+1) { famp = f(s.nz+1, SBLEFT); bamp = b(s.nz+1, SBLEFT); ma = marr[s.nz+1]; n = ma.num; for(int j=0; j<=n-1; ++j) { a += famp(j)*overlap(m, dir, ma(j), FORW) *ma(j).efac(FORW, z-seg.x0); a += bamp(j)*overlap(m, dir, ma(j), BACK) *ma(j).efac(BACK, z-seg.x0); } return a/m.power(); } famp = f(si, SBLEFT); bamp = b(si, SBRIGHT); ma = marr[si]; n = ma.num; for(int j=0; j<=n-1; ++j) { a += famp(j)*overlap(m, dir, ma(j), FORW) *ma(j).efac(FORW, z-seg.x0); a += bamp(j)*overlap(m, dir, ma(j), BACK) *ma(j).efac(BACK, z-seg.x1); } return a/m.power(); } /* power associated with directional mode m at position z */ double BepField::locModPow(double z, const Mode& m, Propdir dir) const { if(BF_vM) bepfielderror("locModPow: vM not implemented"); return locModAmp(z, m, dir).sqabs()*m.power(); } /* ------------------------------------------------------------------- */ /* compute the discretized bidirectional mode spectra: modes of polarization p, Mx expansion terms per vertical slice on a vertical computational window cw */ void BepField::discretize(Polarization p, Interval cw, int Mx, int quiet) { Waveguide wg; ModeArray ma; if(quiet != 1) fprintf(stderr, "SpecDisc: "); for(int i=0; i<=s.nz+1; ++i) { int done = 0; int j=0; while(j marr[i].num) rfmin = marr[i].num; if(rfmax < marr[i].num) rfmax = marr[i].num; } if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin); else fprintf(stderr, " (%d-%d)) ", rfmin, rfmax); if(s.nz == -1) { Cvector zero(marr[0].num); zero.init(CC0); afl[0] = afr[0]; afr[0] = zero; abl[0] = zero; abr[0] = zero; return; } Dvector *p; p = new Dvector[s.nz+2]; for(int i=0; i<=s.nz+1; ++i) { int done = 0; int o = 0; while(done == 0 && o<=i-1) { if(marr[o].equal(marr[i]) == 1) { p[i] = p[o]; done = 1; } ++o; } if(done == 0) { ModeArray ma = marr[i]; int n = ma.num; Dvector tp(n); for(int m=0; m<=n-1; ++m) tp(m) = sqrt(fabs(ma(m).power())); p[i] = tp; } } Cmatrix *dmat; Ivector filld(s.nz+2); Ivector idd(s.nz+2); dmat = new Cmatrix[s.nz+2]; filld.init(0); for(int j=0; j<=s.nz+1; ++j) { int done = 0; int o = 0; while(done == 0 && o<=j-1) { if(marr[o].equal(marr[j]) == 1 && filld(o) == 1) { idd(j) = o; done = 1; } ++o; } if(done == 0) { dmat[j] = orthomat_inv(marr[j]); fprintf(stderr, "."); idd(j) = j; filld(j) = 1; } } filld.strip(); Cmatrix m_m, m_p; Cmatrix *mmff, *mmfb, *mmbf, *mmbb; Cmatrix *mpff, *mpfb, *mpbf, *mpbb; Ivector fill(s.nz+1); Ivector id(s.nz+1); mmff = new Cmatrix[s.nz+1]; mmfb = new Cmatrix[s.nz+1]; mmbf = new Cmatrix[s.nz+1]; mmbb = new Cmatrix[s.nz+1]; mpff = new Cmatrix[s.nz+1]; mpfb = new Cmatrix[s.nz+1]; mpbf = new Cmatrix[s.nz+1]; mpbb = new Cmatrix[s.nz+1]; fill.init(0); for(int j=0; j<=s.nz; ++j) { int done = 0; int o = 0; while(done == 0 && o<=j-1) { if(marr[o ].equal(marr[j ]) == 1 && marr[o+1].equal(marr[j+1]) == 1 && fill(o) == 1) { id(j) = o+1; done = 1; } if(marr[o ].equal(marr[j+1]) == 1 && marr[o+1].equal(marr[j ]) == 1 && fill(o) == 1) { id(j) = -1*(o+1); done = 1; } ++o; } if(done == 0) { int nl, nr; nl = marr[j ].num; nr = marr[j+1].num; if(tm_m != NULL && tm_p != NULL) { m_m = tm_m[j]; m_p = tm_p[j]; } else transfermats(marr[j], marr[j+1], dmat[idd(j)], dmat[idd(j+1)], m_m, m_p); fprintf(stderr, "."); mmff[j] = m_m.submatrix(nl, nr, 0, 0 ); mmfb[j] = m_m.submatrix(nl, nr, 0, nr); mmbf[j] = m_m.submatrix(nl, nr, nl, 0 ); mmbb[j] = m_m.submatrix(nl, nr, nl, nr); m_m.strip(); mpff[j] = m_p.submatrix(nr, nl, 0, 0 ); mpfb[j] = m_p.submatrix(nr, nl, 0, nl); mpbf[j] = m_p.submatrix(nr, nl, nr, 0 ); mpbb[j] = m_p.submatrix(nr, nl, nr, nl); m_p.strip(); id(j) = j+1; fill(j) = 1; } } fill.strip(); idd.strip(); delete[] dmat; fprintf(stderr, ":"); Cvector *t = NULL; if(s.nz >= 1) { t = new Cvector[s.nz]; for(int i=0; i<=s.nz-1; ++i) { ModeArray ma = marr[i+1]; int ni = ma.num; double l = s.hz(i+1)-s.hz(i); Cvector tmpt(ni); for(int j=0; j<=ni-1; ++j) { /* if(ma(j).ppt == PROPAG) tmpt(j) = expi(-l*ma(j).beta); else tmpt(j) = exp(-l*ma(j).beta); */ tmpt(j) = ma(j).efac(FORW, l); } t[i] = tmpt; } fprintf(stderr, "."); } Cmatrix u; Cmatrix ri, rr; Cmatrix *r; r = new Cmatrix[s.nz+1]; u.unit(marr[s.nz].num); if(id(s.nz) > 0) ri = sub(u, mult(mmbf[ id(s.nz)-1], mpfb[ id(s.nz)-1])); else ri = sub(u, mult(mpbf[-id(s.nz)-1], mmfb[-id(s.nz)-1])); ri.inverse(); if(id(s.nz) > 0) r[s.nz] = mult(ri, mult(mmbf[ id(s.nz)-1], mpff[ id(s.nz)-1])); else r[s.nz] = mult(ri, mult(mpbf[-id(s.nz)-1], mmff[-id(s.nz)-1])); fprintf(stderr, "."); for(int j=s.nz-1; j>=0; --j) { u.unit(marr[j].num); if(id(j) > 0) rr = add(mmbf[ id(j)-1], mult(mmbb[ id(j)-1], diagmult(t[j], diagmult(r[j+1],t[j])))); else rr = add(mpbf[-id(j)-1], mult(mpbb[-id(j)-1], diagmult(t[j], diagmult(r[j+1],t[j])))); if(id(j) > 0) ri = sub(u, mult(rr, mpfb[ id(j)-1])); else ri = sub(u, mult(rr, mmfb[-id(j)-1])); ri.inverse(); if(id(j) > 0) r[j] = mult(ri, mult(rr, mpff[ id(j)-1])); else r[j] = mult(ri, mult(rr, mmff[-id(j)-1])); rr.strip(); ri.strip(); fprintf(stderr, "."); } u.strip(); delete[] mpbf; delete[] mpbb; delete[] mmbf; delete[] mmbb; for(int m=0; m<=marr[0].num-1; ++m) (afr[0])(m)/=(p[0])(m); abr[0] = mult(r[0], afr[0]); r[0].strip(); for(int j=1; j<=s.nz; ++j) { if(id(j-1) > 0) afl[j] = add(mult(mpff[ id(j-1)-1], afr[j-1]), mult(mpfb[ id(j-1)-1], abr[j-1])); else afl[j] = add(mult(mmff[-id(j-1)-1], afr[j-1]), mult(mmfb[-id(j-1)-1], abr[j-1])); afr[j] = diagmult(t[j-1], afl[j]); abr[j] = mult(r[j], afr[j]); r[j].strip(); abl[j] = diagmult(t[j-1], abr[j]); fprintf(stderr, "."); } if(id(s.nz) > 0) afl[s.nz+1] = add(mult(mpff[ id(s.nz)-1], afr[s.nz]), mult(mpfb[ id(s.nz)-1], abr[s.nz])); else afl[s.nz+1] = add(mult(mmff[-id(s.nz)-1], afr[s.nz]), mult(mmfb[-id(s.nz)-1], abr[s.nz])); delete[] mpff; delete[] mpfb; delete[] mmff; delete[] mmfb; delete[] r; for(int m=0; m<=marr[0].num-1; ++m) { (afr[0])(m)*=(p[0])(m); (abr[0])(m)*=(p[0])(m); } for(int j=1; j<=s.nz; ++j) { for(int m=0; m<=marr[j].num-1; ++m) { (afl[j])(m)*=(p[j])(m); (afr[j])(m)*=(p[j])(m); (abl[j])(m)*=(p[j])(m); (abr[j])(m)*=(p[j])(m); } } for(int m=0; m<=marr[s.nz+1].num-1; ++m) (afl[s.nz+1])(m)*=(p[s.nz+1])(m); fprintf(stderr, " Ok.\n"); Cvector zero; zero = Cvector(marr[0].num); zero.init(CC0); abl[0] = zero; afl[0] = zero; zero = Cvector(marr[s.nz+1].num); zero.init(CC0); abl[s.nz+1] = zero; afr[s.nz+1] = zero; abr[s.nz+1] = zero; delete[] p; if(t != NULL) delete[] t; return; } /* 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 BepField::bepsim() { Cmatrix *tm_m; Cmatrix *tm_p; tm_m = NULL; tm_p = NULL; bepsim(tm_m, tm_p); return; } /* 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 BepField::bepsim(Cvector& inp, Cvector& ref, Cvector& trans) { setf(0, SBRIGHT, inp); bepsim(); ref = b(0, SBRIGHT); trans = f(s.nz+1, SBLEFT); return; } /* 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 BepField::bepsim(Cmatrix& ref, Cmatrix& trans) { fprintf(stderr, "BEP ("); if(s.nz < -1 ) bepfielderror("bepsim: nz < -1"); if(s.nz+2 != nums) bepfielderror("bepsim: structure mismatch"); fprintf(stderr, "%d x", s.nz+2); int rfmin, rfmax; rfmin = rfmax = marr[0].num; for(int i=1; i<=s.nz+1; ++i) { if(rfmin > marr[i].num) rfmin = marr[i].num; if(rfmax < marr[i].num) rfmax = marr[i].num; } if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin); else fprintf(stderr, " (%d-%d)) ", rfmin, rfmax); if(s.nz == -1) { ref = Cmatrix(marr[0].num, marr[0].num); ref.init(CC0); trans.unit(marr[0].num); return; } Cmatrix *dmat; Ivector filld(s.nz+2); Ivector idd(s.nz+2); dmat = new Cmatrix[s.nz+2]; filld.init(0); for(int j=0; j<=s.nz+1; ++j) { int done = 0; int o = 0; while(done == 0 && o<=j-1) { if(marr[o].equal(marr[j]) == 1 && filld(o) == 1) { idd(j) = o; done = 1; } ++o; } if(done == 0) { dmat[j] = orthomat_inv(marr[j]); fprintf(stderr, "."); idd(j) = j; filld(j) = 1; } } filld.strip(); Cmatrix m_m, m_p; Cmatrix *mmff, *mmfb, *mmbf, *mmbb; Cmatrix *mpff, *mpfb, *mpbf, *mpbb; Ivector fill(s.nz+1); Ivector id(s.nz+1); mmff = new Cmatrix[s.nz+1]; mmfb = new Cmatrix[s.nz+1]; mmbf = new Cmatrix[s.nz+1]; mmbb = new Cmatrix[s.nz+1]; mpff = new Cmatrix[s.nz+1]; mpfb = new Cmatrix[s.nz+1]; mpbf = new Cmatrix[s.nz+1]; mpbb = new Cmatrix[s.nz+1]; fill.init(0); for(int j=0; j<=s.nz; ++j) { int done = 0; int o = 0; while(done == 0 && o<=j-1) { if(marr[o ].equal(marr[j ]) == 1 && marr[o+1].equal(marr[j+1]) == 1 && fill(o) == 1) { id(j) = o+1; done = 1; } if(marr[o ].equal(marr[j+1]) == 1 && marr[o+1].equal(marr[j ]) == 1 && fill(o) == 1) { id(j) = -1*(o+1); done = 1; } ++o; } if(done == 0) { int nl, nr; nl = marr[j ].num; nr = marr[j+1].num; transfermats(marr[j], marr[j+1], dmat[idd(j)], dmat[idd(j+1)], m_m, m_p); fprintf(stderr, "."); mmff[j] = m_m.submatrix(nl, nr, 0, 0 ); mmfb[j] = m_m.submatrix(nl, nr, 0, nr); mmbf[j] = m_m.submatrix(nl, nr, nl, 0 ); mmbb[j] = m_m.submatrix(nl, nr, nl, nr); m_m.strip(); mpff[j] = m_p.submatrix(nr, nl, 0, 0 ); mpfb[j] = m_p.submatrix(nr, nl, 0, nl); mpbf[j] = m_p.submatrix(nr, nl, nr, 0 ); mpbb[j] = m_p.submatrix(nr, nl, nr, nl); m_p.strip(); id(j) = j+1; fill(j) = 1; } } fill.strip(); idd.strip(); delete[] dmat; fprintf(stderr, ":"); Cvector *t = NULL; if(s.nz >= 1) { t = new Cvector[s.nz]; for(int i=0; i<=s.nz-1; ++i) { ModeArray ma = marr[i+1]; int ni = ma.num; double l = s.hz(i+1)-s.hz(i); Cvector tmpt(ni); for(int j=0; j<=ni-1; ++j) { /* if(ma(j).ppt == PROPAG) tmpt(j) = expi(-l*ma(j).beta); else tmpt(j) = exp(-l*ma(j).beta); */ tmpt(j) = ma(j).efac(FORW, l); } t[i] = tmpt; } fprintf(stderr, "."); } Cmatrix u; Cmatrix ri, rr; Cmatrix *r; r = new Cmatrix[s.nz+1]; u.unit(marr[s.nz].num); if(id(s.nz) > 0) ri = sub(u, mult(mmbf[ id(s.nz)-1], mpfb[ id(s.nz)-1])); else ri = sub(u, mult(mpbf[-id(s.nz)-1], mmfb[-id(s.nz)-1])); ri.inverse(); if(id(s.nz) > 0) r[s.nz] = mult(ri, mult(mmbf[ id(s.nz)-1], mpff[ id(s.nz)-1])); else r[s.nz] = mult(ri, mult(mpbf[-id(s.nz)-1], mmff[-id(s.nz)-1])); fprintf(stderr, "."); for(int j=s.nz-1; j>=0; --j) { u.unit(marr[j].num); if(id(j) > 0) rr = add(mmbf[ id(j)-1], mult(mmbb[ id(j)-1], diagmult(t[j], diagmult(r[j+1],t[j])))); else rr = add(mpbf[-id(j)-1], mult(mpbb[-id(j)-1], diagmult(t[j], diagmult(r[j+1],t[j])))); if(id(j) > 0) ri = sub(u, mult(rr, mpfb[ id(j)-1])); else ri = sub(u, mult(rr, mmfb[-id(j)-1])); ri.inverse(); if(id(j) > 0) r[j] = mult(ri, mult(rr, mpff[ id(j)-1])); else r[j] = mult(ri, mult(rr, mmff[-id(j)-1])); rr.strip(); ri.strip(); fprintf(stderr, "."); } u.strip(); delete[] mpbf; delete[] mpbb; delete[] mmbf; delete[] mmbb; ref = r[0]; Cmatrix tf, tb, tmp; tf.unit(marr[0].num); tb = r[0]; r[0].strip(); for(int j=1; j<=s.nz; ++j) { if(id(j-1) > 0) tmp = add(mult(mpff[ id(j-1)-1], tf), mult(mpfb[ id(j-1)-1], tb)); else tmp = add(mult(mmff[-id(j-1)-1], tf), mult(mmfb[-id(j-1)-1], tb)); tf = diagmult(t[j-1], tmp); tb = mult(r[j], tf); r[j].strip(); fprintf(stderr, "."); } if(id(s.nz) > 0) trans = add(mult(mpff[ id(s.nz)-1], tf), mult(mpfb[ id(s.nz)-1], tb)); else trans = add(mult(mmff[-id(s.nz)-1], tf), mult(mmfb[-id(s.nz)-1], tb)); tf.strip(); tb.strip(); delete[] mpff; delete[] mpfb; delete[] mmff; delete[] mmfb; delete[] r; fprintf(stderr, " Ok.\n"); if(t != NULL) delete[] t; return; } /* 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 BepField::bepsim(Cmatrix& kff, Cmatrix& kfb, Cmatrix& kbf, Cmatrix& kbb) { fprintf(stderr, "BEP ("); if(s.nz+2 != nums) bepfielderror("bepsim: structure mismatch"); if(s.nz < 1 ) bepfielderror("bepsim: nz < 1"); fprintf(stderr, "%d x", s.nz+2); int rfmin, rfmax; rfmin = rfmax = marr[0].num; for(int i=1; i<=s.nz; ++i) { if(rfmin > marr[i].num) rfmin = marr[i].num; if(rfmax < marr[i].num) rfmax = marr[i].num; } if(rfmin > marr[s.nz+1].num) rfmin = marr[s.nz+1].num; if(rfmax < marr[s.nz+1].num) rfmax = marr[s.nz+1].num; if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin); else fprintf(stderr, " (%d - %d)) ", rfmin, rfmax); Cmatrix *dmat; Ivector filld(s.nz+2); Ivector idd(s.nz+2); dmat = new Cmatrix[s.nz+2]; filld.init(0); for(int j=0; j<=s.nz+1; ++j) { int done = 0; int o = 0; while(done == 0 && o<=j-1) { if(marr[o].equal(marr[j]) == 1 && filld(o) == 1) { idd(j) = o; done = 1; } ++o; } if(done == 0) { dmat[j] = orthomat_inv(marr[j]); fprintf(stderr, "."); idd(j) = j; filld(j) = 1; } } filld.strip(); Cmatrix m_m, m_p; Cmatrix *mmff, *mmfb, *mmbf, *mmbb; Cmatrix *mpff, *mpfb, *mpbf, *mpbb; Ivector fill(s.nz+1); Ivector id(s.nz+1); mmff = new Cmatrix[s.nz+1]; mmfb = new Cmatrix[s.nz+1]; mmbf = new Cmatrix[s.nz+1]; mmbb = new Cmatrix[s.nz+1]; mpff = new Cmatrix[s.nz+1]; mpfb = new Cmatrix[s.nz+1]; mpbf = new Cmatrix[s.nz+1]; mpbb = new Cmatrix[s.nz+1]; fill.init(0); for(int j=0; j<=s.nz; ++j) { int done = 0; int o = 0; while(done == 0 && o<=j-1) { if(marr[o ].equal(marr[j ]) == 1 && marr[o+1].equal(marr[j+1]) == 1 && fill(o) == 1) { id(j) = o+1; done = 1; } if(marr[o ].equal(marr[j+1]) == 1 && marr[o+1].equal(marr[j ]) == 1 && fill(o) == 1) { id(j) = -1*(o+1); done = 1; } ++o; } if(done == 0) { int nl, nr; nl = marr[j ].num; nr = marr[j+1].num; transfermats(marr[j], marr[j+1], dmat[idd(j)], dmat[idd(j+1)], m_m, m_p); fprintf(stderr, "."); mmff[j] = m_m.submatrix(nl, nr, 0, 0 ); mmfb[j] = m_m.submatrix(nl, nr, 0, nr); mmbf[j] = m_m.submatrix(nl, nr, nl, 0 ); mmbb[j] = m_m.submatrix(nl, nr, nl, nr); m_m.strip(); mpff[j] = m_p.submatrix(nr, nl, 0, 0 ); mpfb[j] = m_p.submatrix(nr, nl, 0, nl); mpbf[j] = m_p.submatrix(nr, nl, nr, 0 ); mpbb[j] = m_p.submatrix(nr, nl, nr, nl); m_p.strip(); id(j) = j+1; fill(j) = 1; } } fill.strip(); idd.strip(); delete[] dmat; fprintf(stderr, ":"); Cvector *t = NULL; t = new Cvector[s.nz]; for(int i=0; i<=s.nz-1; ++i) { ModeArray ma = marr[i+1]; int ni = ma.num; double l = s.hz(i+1)-s.hz(i); Cvector tmpt(ni); for(int j=0; j<=ni-1; ++j) { /* if(ma(j).ppt == PROPAG) tmpt(j) = expi(-l*ma(j).beta); else tmpt(j) = exp(-l*ma(j).beta); */ tmpt(j) = ma(j).efac(FORW, l); } t[i] = tmpt; } fprintf(stderr, "."); Cmatrix u, ri; Cmatrix *jff, *jfb, *jbf, *jbb; jff = new Cmatrix[s.nz+1]; jfb = new Cmatrix[s.nz+1]; jbf = new Cmatrix[s.nz+1]; jbb = new Cmatrix[s.nz+1]; for(int j=0; j<=s.nz; ++j) { u.unit(marr[j+1].num); if(id(j) > 0) ri = sub(u, mult(mpfb[ id(j)-1], mmbf[ id(j)-1])); else ri = sub(u, mult(mmfb[-id(j)-1], mpbf[-id(j)-1])); ri.inverse(); if(id(j) > 0) jff[j] = mult(ri, mpff[ id(j)-1]); else jff[j] = mult(ri, mmff[-id(j)-1]); if(id(j) > 0) jfb[j] = mult(ri, mult(mpfb[ id(j)-1], mmbb[ id(j)-1])); else jfb[j] = mult(ri, mult(mmfb[-id(j)-1], mpbb[-id(j)-1])); u.unit(marr[j].num); if(id(j) > 0) ri = sub(u, mult(mmbf[ id(j)-1], mpfb[ id(j)-1])); else ri = sub(u, mult(mpbf[-id(j)-1], mmfb[-id(j)-1])); ri.inverse(); if(id(j) > 0) jbf[j] = mult(ri, mult(mmbf[ id(j)-1], mpff[ id(j)-1])); else jbf[j] = mult(ri, mult(mpbf[-id(j)-1], mmff[-id(j)-1])); if(id(j) > 0) jbb[j] = mult(ri, mmbb[ id(j)-1]); else jbb[j] = mult(ri, mpbb[-id(j)-1]); fprintf(stderr, "."); } delete[] mmff; delete[] mmfb; delete[] mmbf; delete[] mmbb; delete[] mpff; delete[] mpfb; delete[] mpbf; delete[] mpbb; Cmatrix *kkff, *kkfb, *kkbf, *kkbb; kkff = new Cmatrix[s.nz+1]; kkfb = new Cmatrix[s.nz+1]; kkbf = new Cmatrix[s.nz+1]; kkbb = new Cmatrix[s.nz+1]; Cmatrix *aff, *afb, *abf, *abb; aff = new Cmatrix[s.nz]; afb = new Cmatrix[s.nz]; abf = new Cmatrix[s.nz]; abb = new Cmatrix[s.nz]; kkff[0] = jff[0]; kkfb[0] = jfb[0]; kkbf[0] = jbf[0]; kkbb[0] = jbb[0]; for(int j=0; j<=s.nz-1; ++j) { u.unit(marr[j+1].num); ri = sub(u, mult(diagmult(kkfb[j ], t[j+1-1]), diagmult(jbf[j+1], t[j+1-1]))); ri.inverse(); aff[j] = mult(ri, kkff[j]); afb[j] = mult(ri, mult(kkfb[j], diagmult(t[j+1-1], jbb[j+1]))); ri = sub(u, mult(diagmult(jbf[j+1], t[j+1-1]), diagmult(kkfb[j ], t[j+1-1]))); ri.inverse(); abf[j] = mult(ri, mult(jbf[j+1], diagmult(t[j+1-1], kkff[j]))); abb[j] = mult(ri, jbb[j+1]); ri = diagmult(jff[j+1], t[j+1-1]); kkff[j+1] = mult(ri, aff[j]); kkfb[j+1] = add(mult(ri, afb[j]), jfb[j+1]); ri = diagmult(kkbb[j], t[j+1-1]); kkbf[j+1] = add(mult(ri, abf[j]), kkbf[j]); kkbb[j+1] = mult(ri, abb[j]); fprintf(stderr, "."); } kff = kkff[s.nz]; kfb = kkfb[s.nz]; kbf = kkbf[s.nz]; kbb = kkbb[s.nz]; delete[] kkff; delete[] kkfb; delete[] kkbf; delete[] kkbb; delete[] jff; delete[] jfb; delete[] jbf; delete[] jbb; delete[] aff; delete[] afb; delete[] abf; delete[] abb; delete[] t; fprintf(stderr, " Ok.\n"); return; } /* ------------------------------------------------------------------- */ /* 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 BepField::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) { fprintf(stderr, "BEP ("); if(s.nz < 3) bepfielderror("quepsetup: nz < 3"); if(s.nz+2 != nums) bepfielderror("quepsetup: structure mismatch"); fprintf(stderr, "%d x", s.nz+2); int rfmin, rfmax; if(afr[0].nel != marr[0].num) bepfielderror("quepsetup: dimensions 0"); abr[0] = Cvector(marr[0].num); afl[0] = Cvector(marr[0].num); abl[0] = Cvector(marr[0].num); rfmin = rfmax = marr[0].num; for(int i=1; i<=s.nz; ++i) { afl[i] = Cvector(marr[i].num); afr[i] = Cvector(marr[i].num); abl[i] = Cvector(marr[i].num); abr[i] = Cvector(marr[i].num); if(rfmin > marr[i].num) rfmin = marr[i].num; if(rfmax < marr[i].num) rfmax = marr[i].num; } if(abl[s.nz+1].nel != marr[s.nz+1].num) bepfielderror("quepsetup: dimensions nz+1"); afl[s.nz+1] = Cvector(marr[s.nz+1].num); afr[s.nz+1] = Cvector(marr[s.nz+1].num); abr[s.nz+1] = Cvector(marr[s.nz+1].num); if(rfmin > marr[s.nz+1].num) rfmin = marr[s.nz+1].num; if(rfmax < marr[s.nz+1].num) rfmax = marr[s.nz+1].num; if(rfmin == rfmax) fprintf(stderr, " %d) ", rfmin); else fprintf(stderr, " (%d - %d)) ", rfmin, rfmax); Cmatrix *dmat; Ivector filld(s.nz+2); Ivector idd(s.nz+2); dmat = new Cmatrix[s.nz+2]; filld.init(0); for(int j=0; j<=s.nz+1; ++j) { int done = 0; int o = 0; while(done == 0 && o<=j-1) { if(marr[o].equal(marr[j]) == 1 && filld(o) == 1) { idd(j) = o; done = 1; } ++o; } if(done == 0) { dmat[j] = orthomat_inv(marr[j]); fprintf(stderr, "."); idd(j) = j; filld(j) = 1; } } filld.strip(); d0 = dmat[idd(0)]; d1 = dmat[idd(1)]; dN = dmat[idd(s.nz)]; dNp = dmat[idd(s.nz+1)]; Cmatrix m_m, m_p; Cmatrix *mmff, *mmfb, *mmbf, *mmbb; Cmatrix *mpff, *mpfb, *mpbf, *mpbb; Ivector fill(s.nz+1); Ivector id(s.nz+1); mmff = new Cmatrix[s.nz+1]; mmfb = new Cmatrix[s.nz+1]; mmbf = new Cmatrix[s.nz+1]; mmbb = new Cmatrix[s.nz+1]; mpff = new Cmatrix[s.nz+1]; mpfb = new Cmatrix[s.nz+1]; mpbf = new Cmatrix[s.nz+1]; mpbb = new Cmatrix[s.nz+1]; fill.init(0); for(int j=0; j<=s.nz; ++j) { int done = 0; int o = 0; while(done == 0 && o<=j-1) { if(marr[o ].equal(marr[j ]) == 1 && marr[o+1].equal(marr[j+1]) == 1 && fill(o) == 1) { id(j) = o+1; done = 1; } if(marr[o ].equal(marr[j+1]) == 1 && marr[o+1].equal(marr[j ]) == 1 && fill(o) == 1) { id(j) = -1*(o+1); done = 1; } ++o; } if(done == 0) { int nl, nr; nl = marr[j ].num; nr = marr[j+1].num; transfermats(marr[j], marr[j+1], dmat[idd(j)], dmat[idd(j+1)], m_m, m_p); fprintf(stderr, "."); mmff[j] = m_m.submatrix(nl, nr, 0, 0 ); mmfb[j] = m_m.submatrix(nl, nr, 0, nr); mmbf[j] = m_m.submatrix(nl, nr, nl, 0 ); mmbb[j] = m_m.submatrix(nl, nr, nl, nr); m_m.strip(); mpff[j] = m_p.submatrix(nr, nl, 0, 0 ); mpfb[j] = m_p.submatrix(nr, nl, 0, nl); mpbf[j] = m_p.submatrix(nr, nl, nr, 0 ); mpbb[j] = m_p.submatrix(nr, nl, nr, nl); m_p.strip(); id(j) = j+1; fill(j) = 1; } } fill.strip(); idd.strip(); delete[] dmat; if(id(0) > 0) { m0mbf = mmbf[id(0)-1]; m0mbb = mmbb[id(0)-1]; m0pff = mpff[id(0)-1]; m0pfb = mpfb[id(0)-1]; } else { m0mbf = mpbf[-id(0)-1]; m0mbb = mpbb[-id(0)-1]; m0pff = mmff[-id(0)-1]; m0pfb = mmfb[-id(0)-1]; } if(id(s.nz) > 0) { mNmbf = mmbf[id(s.nz)-1]; mNmbb = mmbb[id(s.nz)-1]; mNpff = mpff[id(s.nz)-1]; mNpfb = mpfb[id(s.nz)-1]; } else { mNmbf = mpbf[-id(s.nz)-1]; mNmbb = mpbb[-id(s.nz)-1]; mNpff = mmff[-id(s.nz)-1]; mNpfb = mmfb[-id(s.nz)-1]; } fprintf(stderr, ":"); Cvector *t = NULL; t = new Cvector[s.nz]; for(int i=0; i<=s.nz-1; ++i) { ModeArray ma = marr[i+1]; int ni = ma.num; double l = s.hz(i+1)-s.hz(i); Cvector tmpt(ni); for(int j=0; j<=ni-1; ++j) { /* if(ma(j).ppt == PROPAG) tmpt(j) = expi(-l*ma(j).beta); else tmpt(j) = exp(-l*ma(j).beta); */ tmpt(j) = ma(j).efac(FORW, l); } t[i] = tmpt; } fprintf(stderr, "."); t1 = t[0]; for(int i=1; i<=s.nz-2; ++i) tx[i+1] = t[i]; tN = t[s.nz-1]; Cmatrix u, ri; Cmatrix *jff, *jfb, *jbf, *jbb; jff = new Cmatrix[s.nz]; jfb = new Cmatrix[s.nz]; jbf = new Cmatrix[s.nz]; jbb = new Cmatrix[s.nz]; for(int j=1; j<=s.nz-1; ++j) { u.unit(marr[j+1].num); if(id(j) > 0) ri = sub(u, mult(mpfb[ id(j)-1], mmbf[ id(j)-1])); else ri = sub(u, mult(mmfb[-id(j)-1], mpbf[-id(j)-1])); ri.inverse(); if(id(j) > 0) jff[j] = mult(ri, mpff[ id(j)-1]); else jff[j] = mult(ri, mmff[-id(j)-1]); if(id(j) > 0) jfb[j] = mult(ri, mult(mpfb[ id(j)-1], mmbb[ id(j)-1])); else jfb[j] = mult(ri, mult(mmfb[-id(j)-1], mpbb[-id(j)-1])); u.unit(marr[j].num); if(id(j) > 0) ri = sub(u, mult(mmbf[ id(j)-1], mpfb[ id(j)-1])); else ri = sub(u, mult(mpbf[-id(j)-1], mmfb[-id(j)-1])); ri.inverse(); if(id(j) > 0) jbf[j] = mult(ri, mult(mmbf[ id(j)-1], mpff[ id(j)-1])); else jbf[j] = mult(ri, mult(mpbf[-id(j)-1], mmff[-id(j)-1])); if(id(j) > 0) jbb[j] = mult(ri, mmbb[ id(j)-1]); else jbb[j] = mult(ri, mpbb[-id(j)-1]); fprintf(stderr, "."); } delete[] mmff; delete[] mmfb; delete[] mmbf; delete[] mmbb; delete[] mpff; delete[] mpfb; delete[] mpbf; delete[] mpbb; Cmatrix *kkff, *kkfb, *kkbf, *kkbb; kkff = new Cmatrix[s.nz]; kkfb = new Cmatrix[s.nz]; kkbf = new Cmatrix[s.nz]; kkbb = new Cmatrix[s.nz]; Cmatrix *aff, *afb, *abf, *abb; aff = new Cmatrix[s.nz-1]; afb = new Cmatrix[s.nz-1]; abf = new Cmatrix[s.nz-1]; abb = new Cmatrix[s.nz-1]; kkff[1] = jff[1]; kkfb[1] = jfb[1]; kkbf[1] = jbf[1]; kkbb[1] = jbb[1]; for(int j=1; j<=s.nz-2; ++j) { u.unit(marr[j+1].num); ri = sub(u, mult(diagmult(kkfb[j ], t[j+1-1]), diagmult(jbf[j+1], t[j+1-1]))); ri.inverse(); aff[j] = mult(ri, kkff[j]); afb[j] = mult(ri, mult(kkfb[j], diagmult(t[j+1-1], jbb[j+1]))); ri = sub(u, mult(diagmult(jbf[j+1], t[j+1-1]), diagmult(kkfb[j ], t[j+1-1]))); ri.inverse(); abf[j] = mult(ri, mult(jbf[j+1], diagmult(t[j+1-1], kkff[j]))); abb[j] = mult(ri, jbb[j+1]); ri = diagmult(jff[j+1], t[j+1-1]); kkff[j+1] = mult(ri, aff[j]); kkfb[j+1] = add(mult(ri, afb[j]), jfb[j+1]); ri = diagmult(kkbb[j], t[j+1-1]); kkbf[j+1] = add(mult(ri, abf[j]), kkbf[j]); kkbb[j+1] = mult(ri, abb[j]); fprintf(stderr, "."); } kff = kkff[s.nz-1]; kfb = kkfb[s.nz-1]; kbf = kkbf[s.nz-1]; kbb = kkbb[s.nz-1]; delete[] kkff; delete[] kkfb; delete[] kkbf; delete[] kkbb; delete[] jff; delete[] jfb; delete[] jbf; delete[] jbb; iff[s.nz-1] = aff[s.nz-2]; ifb[s.nz-1] = afb[s.nz-2]; ibf[s.nz-1] = abf[s.nz-2]; ibb[s.nz-1] = abb[s.nz-2]; for(int j=s.nz-1; j>=3; --j) { ri = diagmult(afb[j-2], t[j-1]); iff[j-1] = add(aff[j-2], mult(ri, ibf[j])); ifb[j-1] = mult(ri, ibb[j]); ri = diagmult(abb[j-2], t[j-1]); ibf[j-1] = add(abf[j-2], mult(ri, ibf[j])); ibb[j-1] = mult(ri, ibb[j]); } delete[] aff; delete[] afb; delete[] abf; delete[] abb; delete[] t; fprintf(stderr, " Ok.\n"); return; } /* for the QUEP implementation: overlaps of perpendiculary oriented mode expansion fields */ Cmatrix BepField::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) { Cmatrix g, vt; double p = zpos; Cmatrix vff, vfb, vbf, vbb; Cmatrix *gff, *gfb, *gbf, *gbb; gff = new Cmatrix[s.nz+1]; gfb = new Cmatrix[s.nz+1]; gbf = new Cmatrix[s.nz+1]; gbb = new Cmatrix[s.nz+1]; for(int l=1; l<=s.nz; ++l) { mh.crossoverlap(marr[l], p, s.segment(l), gff[l], gfb[l], gbf[l], gbb[l]); } g = gff[1]; vt = g; g = gfb[1]; vt.addeq(mult(g, diagmult(kbf, t1))); for(int l=2; l<=s.nz-1; ++l) { g = gff[l]; vt.addeq(mult(g, diagmult(iff[l], t1))); g = gfb[l]; vt.addeq(mult(g, diagmult(ibf[l], t1))); } g = gff[s.nz]; vt.addeq(mult(g, diagmult(kff, t1))); vff = vt; g = gbf[1]; vt = g; g = gbb[1]; vt.addeq(mult(g, diagmult(kbf, t1))); for(int l=2; l<=s.nz-1; ++l) { g = gbf[l]; vt.addeq(mult(g, diagmult(iff[l], t1))); g = gbb[l]; vt.addeq(mult(g, diagmult(ibf[l], t1))); } g = gbf[s.nz]; vt.addeq(mult(g, diagmult(kff, t1))); vbf = vt; g = gfb[1]; vt = mult(g, diagmult(kbb, tN)); for(int l=2; l<=s.nz-1; ++l) { g = gff[l]; vt.addeq(mult(g, diagmult(ifb[l], tN))); g = gfb[l]; vt.addeq(mult(g, diagmult(ibb[l], tN))); } g = gff[s.nz]; vt.addeq(mult(g, diagmult(kfb, tN))); g = gfb[s.nz]; vt.addeq(g); vfb = vt; g = gbb[1]; vt = mult(g, diagmult(kbb, tN)); for(int l=2; l<=s.nz-1; ++l) { g = gbf[l]; vt.addeq(mult(g, diagmult(ifb[l], tN))); g = gbb[l]; vt.addeq(mult(g, diagmult(ibb[l], tN))); } g = gbf[s.nz]; vt.addeq(mult(g, diagmult(kfb, tN))); g = gbb[s.nz]; vt.addeq(g); vbb = vt; vt.strip(); g.strip(); for(int l=1; l<=s.nz; ++l) { gff[l].strip(); gfb[l].strip(); gbf[l].strip(); gbb[l].strip(); } delete[] gff; delete[] gfb; delete[] gbf; delete[] gbb; Cmatrix vmat(mh.num+mh.num, marr[1].num+marr[s.nz].num); vmat.setsubmatrix(vff, 0, 0 ); vmat.setsubmatrix(vfb, 0, marr[1].num); vmat.setsubmatrix(vbf, mh.num, 0 ); vmat.setsubmatrix(vbb, mh.num, marr[1].num); vff.strip(); vfb.strip(); vbf.strip(); vbb.strip(); return vmat; } /* ------------------------------------------------------------------- */ /* for the vBEP / vQUEP implementaion */ /* prepare for the vTE/vTM mode representation, for given ky */ void BepField::vectorify(Complex ky) { for(int i=0; i<=s.nz+1; ++i) marr[i].vectorify(ky); return; } /* 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 BepField::discretize(Complex ky, int type, Interval cw, int Mx, int quiet) { Waveguide wg; ModeArray ma; if(quiet != 1) fprintf(stderr, "SpecDisc (vM):"); for(int i=0; i<=s.nz+1; ++i) { int done = 0; int j=0; while(j=0; --i) { fprintf(stderr,"layer %d: ", i); for(int j=0; j<=rc.nz+1; ++j) fprintf(stderr, "%6.10g ", rc.n(i, j)); fprintf(stderr,"\n"); } fprintf(stderr, "lambda: %.10g k0: %g", rc.lambda, val_k0(rc.lambda)); if(rc.special) fprintf(stderr, " (!)\n"); else fprintf(stderr, "\n"); fprintf(stderr, "----------------------------------------------------------------------------\n"); } (*this) = BepField(rc, 1); discretize(ky, 0, wx, Mx, quiet); } BepField::BepField(Circuit rc, Complex ky, Interval wx, int Mx) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc, ky, wx, Mx, 0); } BepField::BepField(SegWgStruct rc, Complex ky, Interval wx, int Mx, int quiet) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc.circuit(), ky, wx, Mx, quiet); } BepField::BepField(SegWgStruct rc, Complex ky, Interval wx, int Mx) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc.circuit(), ky, wx, Mx, 0); } /* ... 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::BepField(Circuit rc, SBorder bdr, Polarization pol, int ord, Complex amp, double theta, Interval wx, int Mx, int quiet) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc, 1); Waveguide ewg; switch(bdr) { case SBLEFT: ewg = s(0); break; case SBRIGHT: ewg = s(s.nz+1); break; } ModeArray ma; modeanalysis(ewg, pol, ma, 1); if(ord < 0 || ord >= ma.num) bepfielderror("vBEP, setup, mode order"); double N = ma(ord).neff; double b = ma(ord).beta; if(theta <= -90.0 || theta >= 90.0) bepfielderror("vBEP, setup, theta"); Complex ky = Complex(b*sin(theta/180.0*PI), 0.0); if(quiet != 1) { fprintf(stderr, "\n----------------------------------------------------------------------------\n"); fprintf(stderr, "in: %s, neff = %g, beta = %g; amp = %g+i%g; angle = %g^o\n", ma(ord).ids, N, b, amp.re, amp.im, theta); } (*this) = BepField(rc, ky, wx, Mx, quiet); if(pol == TE) input(bdr, amp, ord); else input(bdr, amp, Mx+ord); } BepField::BepField(Circuit rc, SBorder bdr, Polarization pol, int ord, Complex amp, double theta, Interval wx, int Mx) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc, bdr, pol, ord, amp, theta, wx, Mx, 0); } BepField::BepField(SegWgStruct rc, SBorder bdr, Polarization pol, int ord, Complex amp, double theta, Interval wx, int Mx, int quiet) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc.circuit(), bdr, pol, ord, amp, theta, wx, Mx, quiet); } BepField::BepField(SegWgStruct rc, SBorder bdr, Polarization pol, int ord, Complex amp, double theta, Interval wx, int Mx) { marr = NULL; afl = afr = abl = abr = NULL; (*this) = BepField(rc.circuit(), bdr, pol, ord, amp, theta, wx, Mx, 0); } /* ------------------------------------------------------------------- */ /* preparation of physical field plots: selection of the global phase / snapshot time 0 within one period, a uniform polarization is assumed */ /* multiply the entire field by a phase factor exp(i phi), the internal representation is modified accordingly */ void BepField::adjustphase(double phi) { Complex e = expi(phi); for(int j=0; j<=nums-1; ++ j) { afl[j].multeq(e); afr[j].multeq(e); abl[j].multeq(e); abr[j].multeq(e); } return; } /* 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 BepField::adjustphase(double x, double z) { int seg = s.segidx(z); if(marr[seg].num <= 0) return; Polarization p = (marr[seg])(0).pol; Fcomp cp = principalcomp(p); adjustphase(-field(cp, x, z).arg()); return; } /* 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 BepField::adjustphase(int numx, int numz) { if(marr[0].num <= 0) return; Polarization p = (marr[0])(0).pol; Fcomp cp = principalcomp(p); Interval wx, wz; s.defaultwindow(wx, wz); double x, z; double xs, zs; xs = wx.len()/((double) numx); zs = wz.len()/((double) numz); x = wx.x0; double xm = 0.0, zm = 0.0, maxf = 0.0; double f; for(int xi=0; xi <=numx; ++xi) { z = wz.x0; for(int zi=0; zi <=numz; ++zi) { f = field(cp, x, z).sqabs(); if(f > maxf) { xm=x; zm=z; maxf=f; } z += zs; } x += xs; } adjustphase(-field(cp, xm, zm).arg()); return; } /* ... as before, with (coarse) defaults for numx, numz */ void BepField::adjustphase() { Interval wx, wz; s.defaultwindow(wx, wz); int numx, numz; numx = (int) ceil(sqrt(2500.0*wx.len()/wz.len())); numz = numx*((int) ceil(wz.len()/wx.len())); adjustphase(numx, numz); return; } /* ------------------------------------------------------------------- */ /* visualization of the light propagation, output: MATLAB scripts */ /* mode interference plots, image of component cp cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ foa: MOD, SQR, REP, IMP, ORG = REP xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ext0, ext1: filename id characters for the m.file */ void BepField::plot(Fcomp cp, Afo foa, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { QuepField qfld; qfld.h = (*this); qfld.QF_pol = BF_pol; qfld.QF_lambda = BF_lambda; qfld.QF_vM = BF_vM; qfld.QF_ky = BF_ky; qfld.plot(cp, foa, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* mode interference plots, fancy style surface ;-) cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ext0, ext1: filename id characters for the m.file */ void BepField::fplot(Fcomp cp, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { QuepField qfld; qfld.h = (*this); qfld.QF_pol = BF_pol; qfld.QF_lambda = BF_lambda; qfld.QF_vM = BF_vM; qfld.QF_ky = BF_ky; qfld.fplot(cp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* image of the relative phase cp: one of EX, EY, EZ, HX, HY, HZ xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ext0, ext1: filename id characters for the m.file */ void BepField::phasemap(Fcomp cp, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { QuepField qfld; qfld.h = (*this); qfld.QF_pol = BF_pol; qfld.QF_lambda = BF_lambda; qfld.QF_vM = BF_vM; qfld.QF_ky = BF_ky; qfld.phasemap(cp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* image plot of the field "strength", id: one of mE, mH, qE, qH, mW xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ext0, ext1: filename id characters for the m.file */ void BepField::plot(FSid id, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { QuepField qfld; qfld.h = (*this); qfld.QF_pol = BF_pol; qfld.QF_lambda = BF_lambda; qfld.QF_vM = BF_vM; qfld.QF_ky = BF_ky; qfld.plot(id, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* electromagnetic energy density image xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ext0, ext1: filename id characters for the m.file */ void BepField::plot(double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { QuepField qfld; qfld.h = (*this); qfld.QF_pol = BF_pol; qfld.QF_lambda = BF_lambda; qfld.QF_vM = BF_vM; qfld.QF_ky = BF_ky; qfld.plot(xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* animation of the light propagation the frames show the EY component (TE) resp. HY component (TM), at ntfr equidistant times over one time period xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ext0, ext1: filename id characters for the m.file */ void BepField::movie(double xbeg, double xend, double zbeg, double zend, int npx, int npz, int ntfr, char ext0, char ext1) const { QuepField qfld; qfld.h = (*this); qfld.QF_pol = BF_pol; qfld.QF_lambda = BF_lambda; qfld.QF_vM = BF_vM; qfld.QF_ky = BF_ky; qfld.movie(xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1); return; } /* animation of the light propagation, fancy style ;-) the frames show the EY component (TE) resp. HY component (TM), at ntfr equidistant times over one time period xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ext0, ext1: filename id characters for the m.file */ void BepField::fmovie(double xbeg, double xend, double zbeg, double zend, int npx, int npz, int ntfr, char ext0, char ext1) const { QuepField qfld; qfld.h = (*this); qfld.QF_pol = BF_pol; qfld.QF_lambda = BF_lambda; qfld.QF_vM = BF_vM; qfld.QF_ky = BF_ky; qfld.fmovie(xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1); return; } /* export full field data ("all") into a viewer m-file xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ext0, ext1: filename id characters */ void BepField::viewer(double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { QuepField qfld; qfld.h = (*this); qfld.QF_pol = BF_pol; qfld.QF_lambda = BF_lambda; qfld.QF_vM = BF_vM; qfld.QF_ky = BF_ky; qfld.viewer(xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* ------------------------------------------------------------------------ */ /* orthogonality matrix for a mode superposition, inverted */ Cmatrix orthomat_inv(ModeArray& ma) { Cmatrix dff, dfb, dbf, dbb; Cmatrix d; int n = ma.num; dff = Cmatrix(n, n); dfb = Cmatrix(n, n); dbf = Cmatrix(n, n); dbb = Cmatrix(n, n); for(int l=0; l<=n-1; ++l) { for(int m=0; m<=n-1; ++m) { bidir_overlaps(ma(l), ma(m), dff(l, m), dfb(l, m), dbf(l, m), dbb(l, m)); } } d = Cmatrix(n+n, n+n); d.setsubmatrix(dff, 0, 0); d.setsubmatrix(dfb, 0, n); d.setsubmatrix(dbf, n, 0); d.setsubmatrix(dbb, n, n); dff.strip(); dfb.strip(); dbf.strip(); dbb.strip(); d.inverse(); return d; } /* overlap matrix for two mode superpositions */ Cmatrix overlapmat(ModeArray& m2, ModeArray& m1) { Cmatrix off, ofb, obf, obb; Cmatrix o; off = Cmatrix(m2.num, m1.num); ofb = Cmatrix(m2.num, m1.num); obf = Cmatrix(m2.num, m1.num); obb = Cmatrix(m2.num, m1.num); for(int l=0; l<=m2.num-1; ++l) { for(int m=0; m<=m1.num-1; ++m) { bidir_overlaps(m2(l), m1(m), off(l, m), ofb(l, m), obf(l, m), obb(l, m)); } } o = Cmatrix(m2.num+m2.num, m1.num+m1.num); o.setsubmatrix(off, 0, 0 ); o.setsubmatrix(ofb, 0, m1.num); o.setsubmatrix(obf, m2.num, 0 ); o.setsubmatrix(obb, m2.num, m1.num); off.strip(); ofb.strip(); obf.strip(); obb.strip(); return o; } /* 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) { if(m1.num <= 0 || m2.num <= 0) return; Cmatrix o; Cmatrix d; o = overlapmat(m2, m1); d = orthomat_inv(m2); mat2 = mult(d, o); d.strip(); d = orthomat_inv(m1); o.adjoint(); mat1 = mult(d, o); d.strip(); o.strip(); return; } /* 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) { if(m1.num <= 0 || m2.num <= 0) return; Cmatrix o; o = overlapmat(m2, m1); mat2 = mult(d2, o); o.adjoint(); mat1 = mult(d1, o); o.strip(); return; } /* 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) { if(m1.num <= 0 || m2.num <= 0) return; Cmatrix mat1, mat2; transfermats(m1, m2, mat1, mat2); m2ff = mat2.submatrix(m2.num, m1.num, 0, 0 ); m2fb = mat2.submatrix(m2.num, m1.num, 0, m1.num); m2bf = mat2.submatrix(m2.num, m1.num, m2.num, 0 ); m2bb = mat2.submatrix(m2.num, m1.num, m2.num, m1.num); mat2.strip(); m1ff = mat1.submatrix(m1.num, m2.num, 0, 0 ); m1fb = mat1.submatrix(m1.num, m2.num, 0, m2.num); m1bf = mat1.submatrix(m1.num, m2.num, m1.num, 0 ); m1bb = mat1.submatrix(m1.num, m2.num, m1.num, m2.num); mat1.strip(); return; } /* - specific driver routines --------------------------------------------- */ /* a long three segment coupler: m1: local modes for z<0 m2: local modes for 0 0.0) ++ng; nr = m2.num-ng; fprintf(stderr, "Coupling segment: %d guided modes,\n", ng); fprintf(stderr, " %d modes in the radiation field.\n", nr); Cmatrix m1ffg, m1ffr, m1fbg, m1fbr, m1bfg, m1bfr, m1bbg, m1bbr; fprintf(stderr, " >>> "); m1ffg = m_1.submatrix(m1.num, ng, 0, 0 ); m1ffr = m_1.submatrix(m1.num, nr, 0, ng ); m1fbg = m_1.submatrix(m1.num, ng, 0, m2.num ); m1fbr = m_1.submatrix(m1.num, nr, 0, m2.num+ng); m1bfg = m_1.submatrix(m1.num, ng, m1.num, 0 ); m1bfr = m_1.submatrix(m1.num, nr, m1.num, ng ); m1bbg = m_1.submatrix(m1.num, ng, m1.num, m2.num ); m1bbr = m_1.submatrix(m1.num, nr, m1.num, m2.num+ng); m_1.strip(); Cmatrix m2fgf, m2fgb, m2frf, m2frb, m2bgf, m2bgb, m2brf, m2brb; m2fgf = m_2.submatrix(ng, m1.num, 0, 0 ); m2fgb = m_2.submatrix(ng, m1.num, 0, m1.num); m2frf = m_2.submatrix(nr, m1.num, ng, 0 ); m2frb = m_2.submatrix(nr, m1.num, ng, m1.num); m2bgf = m_2.submatrix(ng, m1.num, m2.num, 0 ); m2bgb = m_2.submatrix(ng, m1.num, m2.num, m1.num); m2brf = m_2.submatrix(nr, m1.num, m2.num+ng, 0 ); m2brb = m_2.submatrix(nr, m1.num, m2.num+ng, m1.num); m_2.strip(); Cmatrix q2fgf, q2fgb, q2frf, q2frb, q2bgf, q2bgb, q2brf, q2brb; q2fgf = q_2.submatrix(ng, m3.num, 0, 0 ); q2fgb = q_2.submatrix(ng, m3.num, 0, m3.num); q2frf = q_2.submatrix(nr, m3.num, ng, 0 ); q2frb = q_2.submatrix(nr, m3.num, ng, m3.num); q2bgf = q_2.submatrix(ng, m3.num, m2.num, 0 ); q2bgb = q_2.submatrix(ng, m3.num, m2.num, m3.num); q2brf = q_2.submatrix(nr, m3.num, m2.num+ng, 0 ); q2brb = q_2.submatrix(nr, m3.num, m2.num+ng, m3.num); q_2.strip(); Cmatrix q3ffg, q3ffr, q3fbg, q3fbr, q3bfg, q3bfr, q3bbg, q3bbr; q3ffg = q_3.submatrix(m3.num, ng, 0, 0 ); q3ffr = q_3.submatrix(m3.num, nr, 0, ng ); q3fbg = q_3.submatrix(m3.num, ng, 0, m2.num ); q3fbr = q_3.submatrix(m3.num, nr, 0, m2.num+ng); q3bfg = q_3.submatrix(m3.num, ng, m3.num, 0 ); q3bfr = q_3.submatrix(m3.num, nr, m3.num, ng ); q3bbg = q_3.submatrix(m3.num, ng, m3.num, m2.num ); q3bbr = q_3.submatrix(m3.num, nr, m3.num, m2.num+ng); q_3.strip(); fprintf(stderr, "."); Cmatrix tfp(ng, ng); Cmatrix tbp(ng, ng); Cmatrix tfm(ng, ng); Cmatrix tbm(ng, ng); Complex e; tfp.init(CC0); tbp.init(CC0); tfm.init(CC0); tbm.init(CC0); for(int j=0; j<=ng-1; ++j) { e = CCI*(Complex(m2(j).beta, 0.0)+fpert(j))*l*(-1.0); tfp(j, j) = exp(e); e = CCI*(Complex(m2(j).beta, 0.0)+bpert(j))*l; tbp(j, j) = exp(e); e = CCI*(Complex(m2(j).beta, 0.0)+fpert(j))*l; tfm(j, j) = exp(e); e = CCI*(Complex(m2(j).beta, 0.0)+bpert(j))*l*(-1.0); tbm(j, j) = exp(e); } fprintf(stderr, "."); Cmatrix u; Cmatrix ra, rb, rc, rd, re, rf; u.unit(m1.num); ra = add(u, mult(-1.0, mult(m1bfr, m2frb))); rb = add(mult(m1bfg, mult(tfm, q2fgf)), mult(m1bbg, mult(tbm, q2bgf))); rc = mult(m1bfr, m2frf); u.unit(m3.num); rd = add(u, mult(-1.0, mult(q3fbr, q2brf))); re = add(mult(q3ffg, mult(tfp, m2fgf)), mult(q3fbg, mult(tbp, m2bgf))); rf = add(mult(q3ffg, mult(tfp, m2fgb)), mult(q3fbg, mult(tbp, m2bgb))); fprintf(stderr, "."); ra.inverse(); Cmatrix r1, r2, r; fprintf(stderr, "."); r1 = add(rd, mult(-1.0, mult(rf, mult(ra, rb)))); fprintf(stderr, "."); r1.inverse(); fprintf(stderr, "."); r2 = add(re, mult(rf, mult(ra, rc))); r = mult(r1, r2); r1.strip(); r2.strip(); fprintf(stderr, "."); Dvector p1(m1.num); Dvector p2(m2.num); Dvector p3(m3.num); for(int i=0; i<=m1.num-1; ++i) p1(i) = sqrt(fabs(m1(i).power())); for(int i=0; i<=m2.num-1; ++i) p2(i) = sqrt(fabs(m2(i).power())); for(int i=0; i<=m3.num-1; ++i) p3(i) = sqrt(fabs(m3(i).power())); Cvector f1(m1.num); for(int i=0; i<=m1.num-1; ++i) f1(i) = amp1_f(i)/p1(i); fprintf(stderr, "."); amp3_f = mult(r, f1); amp1_b = add(mult(mult(ra, rb), amp3_f), mult(mult(ra, rc), f1)); fprintf(stderr, "."); Cvector v2(m2.num); v2.init(CC0); v2.setsubvector(add(mult(m2fgf, amp1_f), mult(m2fgb, amp1_b)), 0); v2.setsubvector(add(mult(m2frf, amp1_f), mult(m2frb, amp1_b)), ng); amp2_f_a = v2; v2.init(CC0); v2.setsubvector(add(mult(m2bgf, amp1_f), mult(m2bgb, amp1_b)), 0); amp2_b_a = v2; v2.init(CC0); v2.setsubvector(mult(q2fgf, amp3_f), 0); amp2_f_b = v2; v2.init(CC0); v2.setsubvector(mult(q2bgf, amp3_f), 0); v2.setsubvector(mult(q2brf, amp3_f), ng); amp2_b_b = v2; fprintf(stderr, "."); for(int i=0; i<=m3.num-1; ++i) amp3_f(i) = amp3_f(i)*p3(i); for(int i=0; i<=m1.num-1; ++i) amp1_b(i) = amp1_b(i)*p1(i); for(int i=0; i<=m2.num-1; ++i) { amp2_f_a(i) = amp2_f_a(i)*p2(i); amp2_b_a(i) = amp2_b_a(i)*p2(i); amp2_f_b(i) = amp2_f_b(i)*p2(i); amp2_b_b(i) = amp2_b_b(i)*p2(i); } fprintf(stderr, " <<<\n"); /* Cvector t, q; Cvector vf1, vb1; Cvector vf2ag, vf2ar, vb2ag, vb2ar; Cvector vf2bg, vf2br, vb2bg, vb2br; Cvector vf3, vb3; vf1 = amp1_f; vb1 = amp1_b; vf3 = amp3_f; vb3 = Cvector(m3.num); vb3.init(CC0); vf2ag = amp2_f_a.subvector(ng, 0); vf2ar = amp2_f_a.subvector(nr, ng); vb2ag = amp2_b_a.subvector(ng, 0); vb2ar = amp2_b_a.subvector(nr, ng); vf2bg = amp2_f_b.subvector(ng, 0); vf2br = amp2_f_b.subvector(nr, ng); vb2bg = amp2_b_b.subvector(ng, 0); vb2br = amp2_b_b.subvector(nr, ng); t = add( add( mult(m1bfg, vf2ag), mult(m1bfr, vf2ar)), mult(m1bbg, vb2ag)); q = add( vb1, mult(-1.0, t)); fprintf(stderr, "(1): %g\n", q.norm()); t = add( mult(m2fgf, vf1), mult(m2fgb, vb1)); q = add( vf2ag, mult(-1.0, t)); fprintf(stderr, "(2): %g\n", q.norm()); t = add( mult(m2frf, vf1), mult(m2frb, vb1)); q = add( vf2ar, mult(-1.0, t)); fprintf(stderr, "(3): %g\n", q.norm()); t = add( mult(m2bgf, vf1), mult(m2bgb, vb1)); q = add( vb2ag, mult(-1.0, t)); fprintf(stderr, "(4): %g\n", q.norm()); t = mult(q2fgf, vf3); q = add( vf2bg, mult(-1.0, t)); fprintf(stderr, "(5): %g\n", q.norm()); t = mult(q2bgf, vf3); q = add( vb2bg, mult(-1.0, t)); fprintf(stderr, "(6): %g\n", q.norm()); t = mult(q2brf, vf3); q = add( vb2br, mult(-1.0, t)); fprintf(stderr, "(7): %g\n", q.norm()); t = add( add( mult(q3ffg, vf2bg), mult(q3fbg, vb2bg)), mult(q3fbr, vb2br)); q = add( vf3, mult(-1.0, t)); fprintf(stderr, "(8): %g\n", q.norm()); t = add( add( mult(m1ffg, vf2ag), mult(m1ffr, vf2ar)), mult(m1fbg, vb2ag)); q = add( vf1, mult(-1.0, t)); fprintf(stderr, "(9): %g\n", q.norm()); t = add( add( mult(q3bfg, vf2bg), mult(q3bbg, vb2bg)), mult(q3bbr, vb2br)); q = t; fprintf(stderr, "(10): %g\n", q.norm()); t = add( mult(m2brf, vf1), mult(m2brb, vb1)); q = t; fprintf(stderr, "(11): %g\n", q.norm()); t = add( mult(q2frf, vf3), mult(q2frb, vb3)); q = t; fprintf(stderr, "(12): %g\n", q.norm()); */ fprintf(stderr, "Ok.\n"); return; } /* image plot of the light propagation across the long three segment coupler, output: MATLAB script m1: local modes for z<0 m2: local modes for 0 %g) >> %s\n", zbeg, zend, name); dat = fopen(name, "w+"); mlout_title(dat, name, "Interference field"); mlout_l3scgeoxz_2(dat, m1(m1.num-1).wg, m2(m2.num-1).wg, m3(m3.num-1).wg, l, xbeg, xend); mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz); Dmatrix fld(npx, npz); dx = (xend-xbeg)/(npx-1); dz = (zend-zbeg)/(npz-1); Cvector amp3_b(amp3_f.nel); amp3_b.init(CC0); for(j=0; j<=npz-1; ++j) { z = zbeg+j*dz; for(i=0; i<=npx-1; ++i) { x = xbeg+i*dx; if(z<=0.0) f = m1.field(amp1_f, amp1_b, SZ, x, z).abs(); else { if(z <= l) { // if(z <= l/2.0) // f = m2.field(amp2_f_a, amp2_b_a, // SZ, x, z).abs(); // else // f = m2.field(amp2_f_b, amp2_b_b, // SZ, x, z-l).abs(); f = (l-z)/l*m2.field(amp2_f_a, amp2_b_a, SZ, x, z).abs() + z/l*m2.field(amp2_f_b, amp2_b_b, SZ, x, z-l).abs(); } else f = m3.field(amp3_f, amp3_b, SZ, x, z-l).abs(); } fld(i, j) = sqrt(f); } } mlout_fld(dat, npx, npz, SZ, fld); name[7] = 0; mlout_1Dl3scimage(name, dat); mlout_print(dat, name, 'e'); fclose(dat); return; }