/* * METRIC --- Mode expansion modeling in integrated optics / photonics * http://metric.computational-photonics.eu/ */ /* * slamarr.cpp * Arrays of Modes, mode interference evaluation */ #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"bepfld.h" #include"matlvis.h" /* error message */ void slamarrerror(const char *msg) { fprintf(stderr, "slamarr: %s.\n", msg); exit(1); } /* ModeArray: Array of Modes */ /* initialize */ ModeArray::ModeArray() : num(0) { mvec = NULL; } /* destroy */ ModeArray::~ModeArray() { if(mvec != NULL) delete[] mvec; mvec = NULL; num = 0; } /* copy constructor */ ModeArray::ModeArray(const ModeArray& ma) : num(ma.num) { mvec = new Mode[num]; Mode* ap = ma.mvec; Mode* mp = mvec; for(int i=0; i<=num-1; ++i) *mp++ = *ap++; } /* assignment */ ModeArray& ModeArray::operator=(const ModeArray& ma) { if(this != &ma) { if(mvec != NULL) delete[] mvec; num = ma.num; mvec = new Mode[num]; Mode *ap = ma.mvec; Mode *mp = mvec; for(int i=0; i<=num-1; ++i) *mp++ = *ap++; } return *this; } /* delete all Mode entries */ void ModeArray::clear() { if(mvec != NULL) delete[] mvec; mvec = NULL; num = 0; } /* input from FILE dat */ void ModeArray::read(FILE *dat) { if(dat == stderr || dat == stdout) slamarrerror("ModeArray: read"); if(mvec != NULL) delete[] mvec; mvec = NULL; num = 0; num = fgetint(dat); if(num <= 0) { num = 0; return; } mvec = new Mode[num]; for(int i=0; i<=num-1; ++i) { mvec[i].read(dat); } return; } /* output to FILE dat */ void ModeArray::write(FILE *dat) { comment("ModeArray", dat); comment("num", dat); fputint(num, dat); for(int i=0; i<=num-1; ++i) { comment("*** -------- ***", dat); mvec[i].write(dat); } return; } /* input from default file */ void ModeArray::read(char ext0, char ext1) { char name[10]; FILE *dat; int i; name[0] = 's'; name[1] = 'l'; name[2] = 'a'; name[3] = ext0; name[4] = ext1; name[5] = '.'; name[6] = 's'; name[7] = 'm'; name[8] = 'a'; name[9] = 0; fprintf(stderr, "<< %s ", name); dat = fopen(name, "r"); read(dat); fclose(dat); fprintf(stderr, "#%d\n", num); for(i=0; i<=num-1; ++i) fprintf(stderr, " (%d) %s, beta: %.10g, neff: %.10g, npcB: %.10g\n", i, mvec[i].ids, mvec[i].beta, mvec[i].neff, mvec[i].npcB); return; } /* output to default file */ void ModeArray::write(char ext0, char ext1) { char name[10]; FILE *dat; int i; name[0] = 's'; name[1] = 'l'; name[2] = 'a'; name[3] = ext0; name[4] = ext1; name[5] = '.'; name[6] = 's'; name[7] = 'm'; name[8] = 'a'; name[9] = 0; fprintf(stderr, ">> %s #%d\n", name, num); for(i=0; i<=num-1; ++i) fprintf(stderr, " (%d) %s, beta: %.10g, neff: %.10g, npcB: %.10g\n", i, mvec[i].ids, mvec[i].beta, mvec[i].neff, mvec[i].npcB); dat = fopen(name, "w+"); write(dat); fclose(dat); return; } /* subscripting Mode& ModeArray::operator() (int i) { if(i<0 || i>=num) slamarrerror("ModeArray: () out of range"); return mvec[i]; } Mode ModeArray::operator() (int i) const { if(i<0 || i>=num) slamarrerror("ModeArray: () out of range"); return mvec[i]; } */ /* add a mode */ void ModeArray::add(Mode m) { Mode *avec; avec = new Mode[num+1]; Mode* ap = avec; Mode* mp = mvec; for(int i=0; i<=num-1; ++i) *ap++ = *mp++; *ap = m; if(mvec != NULL) delete[] mvec; mvec = avec; ++num; return; } /* delete a mode entry */ void ModeArray::remove(int i) { if(i<0 || i>=num) slamarrerror("ModeArray: remove: argument out of range"); if(num == 1) { delete[] mvec; mvec = NULL; num = 0; return; } Mode *avec; avec = new Mode[num-1]; Mode* ap = avec; Mode* mp = mvec; for(int j=0; j<=i-1; ++j) *ap++ = *mp++; mp++; for(int j=i+1; j<=num-1; ++j) *ap++ = *mp++; if(mvec != NULL) delete[] mvec; mvec = avec; --num; return; } /* add an entire ModeArray ma */ void ModeArray::merge(ModeArray& ma) { Mode *avec; avec = new Mode[num+ma.num]; Mode* ap = avec; Mode* mp = mvec; for(int i=0; i<=num-1; ++i) *ap++ = *mp++; Mode* np = ma.mvec; for(int i=0; i<=ma.num-1; ++i) *ap++ = *np++; if(mvec != NULL) delete[] mvec; mvec = avec; num += ma.num; return; } /* sort the array by squared propagation constants, highest first */ void ModeArray::sort() { int j, k, maxi; double maxb; int minn; Mode t; if(num<=1) return; for(j=0; j<=num-2; ++j) { maxi = j; maxb = mvec[j].betaq; minn = mvec[j].nodes(); for(k=j+1; k<=num-1; ++k) { if(maxb mvec[k].nodes())) { maxb = mvec[k].betaq; minn = mvec[k].nodes(); maxi = k; } } t = mvec[j]; mvec[j] = mvec[maxi]; mvec[maxi] = t; } return; } /* compare the present object to ma, return 1 if equal, otherwise return 0 */ int ModeArray::equal(ModeArray& ma) { if(this == &ma) return 1; if(num != ma.num) return 0; for(int i=0; i<=num-1; ++i) if(mvec[i].equal(ma(i)) != 1) return 0; return 1; } /* ................................................................... */ /* * mode interference evaluation and visualization * all modes are assumed to belong to the same waveguide ! */ /* bidirectional field superposition at point x, z famp: complex amplitudes of the forward traveling modes at z=0 bamp: complex amplitudes of the backward traveling modes at z=0 (z may be negative) cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */ Complex ModeArray::field(const Cvector& famp, const Cvector& bamp, Fcomp cp, double x, double z) const { return field(famp, bamp, cp, x, z, 0.0, 0.0); } /* bidirectional field superposition at point x, z in [z0, z1], famp: complex amplitudes of the forward traveling modes at z=z0 bamp: complex amplitudes of the backward traveling modes at z=z1 cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */ Complex ModeArray::field(const Cvector& famp, const Cvector& bamp, Fcomp cp, double x, double z, double z0, double z1) const { Complex s, ef, eb, f, b; Complex ex, ey, ez, hx, hy, hz; if(famp.nel < num) slamarrerror("field: array length mismatch"); if(bamp.nel < num) slamarrerror("field: array length mismatch"); s = CC0; double dzf = z-z0; double dzb = z-z1; switch(cp) { case EX: case EY: case EZ: case HX: case HY: case HZ: for(int i=0; i<=num-1; ++i) { mvec[i].cfield(cp, x, f, b); mvec[i].efac(dzf, ef, dzb, eb); // s += famp(i)*f*ef; // s += bamp(i)*b*eb; // if(isfinite(ef.re) && isfinite(ef.im)) s += famp(i)*f*ef; // if(isfinite(eb.re) && isfinite(eb.im)) s += bamp(i)*b*eb; if(fabs(famp(i).re) > NEGLECT_FLD_CONTRIB || fabs(famp(i).im) > NEGLECT_FLD_CONTRIB) if(std::isfinite(ef.re) && std::isfinite(ef.im)) s += famp(i)*f*ef; if(fabs(bamp(i).re) > NEGLECT_FLD_CONTRIB || fabs(bamp(i).im) > NEGLECT_FLD_CONTRIB) if(std::isfinite(eb.re) && std::isfinite(eb.im)) s += bamp(i)*b*eb; } return s; case SX: ey = field(famp, bamp, EY, x, z, z0, z1); ez = field(famp, bamp, EZ, x, z, z0, z1); hy = field(famp, bamp, HY, x, z, z0, z1); hz = field(famp, bamp, HZ, x, z, z0, z1); s = ey*hz.conj()-ez*hy.conj(); return Complex(0.5*s.re, 0.0); case SY: ex = field(famp, bamp, EX, x, z, z0, z1); ez = field(famp, bamp, EZ, x, z, z0, z1); hx = field(famp, bamp, HX, x, z, z0, z1); hz = field(famp, bamp, HZ, x, z, z0, z1); s = ez*hx.conj()-ex*hz.conj(); return Complex(0.5*s.re, 0.0); case SZ: ex = field(famp, bamp, EX, x, z, z0, z1); ey = field(famp, bamp, EY, x, z, z0, z1); hx = field(famp, bamp, HX, x, z, z0, z1); hy = field(famp, bamp, HY, x, z, z0, z1); s = ex*hy.conj()-ey*hx.conj(); return Complex(0.5*s.re, 0.0); default: slamarrerror("field: mysterious cp"); } return CC0; } /* bidirectional field superposition, local field strength at x, z id: one of mE, mH, qE, qH, mW famp: complex amplitudes of the forward traveling modes at z=0 bamp: complex amplitudes of the backward traveling modes at z=0 (z may be negative) */ double ModeArray::lfs(const Cvector& famp, const Cvector& bamp, FSid id, double x, double z) const { Complex ex, ey, ez, hx, hy, hz; double eps_x; switch(id) { case mE: ex = field(famp, bamp, EX, x, z); ey = field(famp, bamp, EY, x, z); ez = field(famp, bamp, EZ, x, z); return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs()); case mH: hx = field(famp, bamp, HX, x, z); hy = field(famp, bamp, HY, x, z); hz = field(famp, bamp, HZ, x, z); return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs()); case qE: ex = field(famp, bamp, EX, x, z); ey = field(famp, bamp, EY, x, z); ez = field(famp, bamp, EZ, x, z); return ex.sqabs()+ey.sqabs()+ez.sqabs(); case qH: hx = field(famp, bamp, HX, x, z); hy = field(famp, bamp, HY, x, z); hz = field(famp, bamp, HZ, x, z); return hx.sqabs()+hy.sqabs()+hz.sqabs(); case mW: ex = field(famp, bamp, EX, x, z); ey = field(famp, bamp, EY, x, z); ez = field(famp, bamp, EZ, x, z); hx = field(famp, bamp, HX, x, z); hy = field(famp, bamp, HY, x, z); hz = field(famp, bamp, HZ, x, z); eps_x = mvec[num-1].wg.eps(x); return 0.25*( MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs()) +EP0*eps_x*(ex.sqabs()+ey.sqabs()+ez.sqabs())); break; default: slamarrerror("lfs(), FSid?"); break; } return 0.0; } /* bidirectional field superposition, local field strength at x, z at point x, z in [z0, z1], id: one of mE, mH, qE, qH, mW famp: complex amplitudes of the forward traveling modes at z=z0 bamp: complex amplitudes of the backward traveling modes at z=z1 */ double ModeArray::lfs(const Cvector& famp, const Cvector& bamp, FSid id, double x, double z, double z0, double z1) const { Complex ex, ey, ez, hx, hy, hz; double eps_x; switch(id) { case mE: ex = field(famp, bamp, EX, x, z, z0, z1); ey = field(famp, bamp, EY, x, z, z0, z1); ez = field(famp, bamp, EZ, x, z, z0, z1); return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs()); case mH: hx = field(famp, bamp, HX, x, z, z0, z1); hy = field(famp, bamp, HY, x, z, z0, z1); hz = field(famp, bamp, HZ, x, z, z0, z1); return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs()); case qE: ex = field(famp, bamp, EX, x, z, z0, z1); ey = field(famp, bamp, EY, x, z, z0, z1); ez = field(famp, bamp, EZ, x, z, z0, z1); return ex.sqabs()+ey.sqabs()+ez.sqabs(); case qH: hx = field(famp, bamp, HX, x, z, z0, z1); hy = field(famp, bamp, HY, x, z, z0, z1); hz = field(famp, bamp, HZ, x, z, z0, z1); return hx.sqabs()+hy.sqabs()+hz.sqabs(); case mW: ex = field(famp, bamp, EX, x, z, z0, z1); ey = field(famp, bamp, EY, x, z, z0, z1); ez = field(famp, bamp, EZ, x, z, z0, z1); hx = field(famp, bamp, HX, x, z, z0, z1); hy = field(famp, bamp, HY, x, z, z0, z1); hz = field(famp, bamp, HZ, x, z, z0, z1); eps_x = mvec[num-1].wg.eps(x); return 0.25*( MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs()) +EP0*eps_x*(ex.sqabs()+ey.sqabs()+ez.sqabs())); break; default: slamarrerror("lfs(z0, z1), FSid?"); break; } return 0.0; } /* bidirectional field superposition, local energy density at point x, z famp: complex amplitudes of the forward traveling modes at z=0 bamp: complex amplitudes of the backward traveling modes at z=0 (z may be negative) */ double ModeArray::endens(const Cvector& famp, const Cvector& bamp, double x, double z) const { return lfs(famp, bamp, mW, x, z); } /* bidirectional field superposition, local energy density at point x, z in [z0, z1], famp: complex amplitudes of the forward traveling modes at z=z0 bamp: complex amplitudes of the backward traveling modes at z=z1 */ double ModeArray::endens(const Cvector& famp, const Cvector& bamp, double x, double z, double z0, double z1) const { return lfs(famp, bamp, mW, x, z, z0, z1); } /* the six electric and magnetic components at x, z */ void ModeArray::emfield(const Cvector& famp, const Cvector& bamp, double x, double z, double z0, double z1, Complex& ex, Complex& ey, Complex& ez, Complex& hx, Complex& hy, Complex& hz) const { Complex s, ef, eb, f, b; if(famp.nel < num) slamarrerror("field: array length mismatch"); if(bamp.nel < num) slamarrerror("field: array length mismatch"); ex = ey = ez = hx = hy = hz = CC0; double dzf = z-z0; double dzb = z-z1; Complex af, ab; bool icf, icb; for(int i=0; i<=num-1; ++i) { af = famp(i); ab = bamp(i); mvec[i].efac(dzf, ef, dzb, eb); icf = true; if(fabs(af.re) < NEGLECT_FLD_CONTRIB && fabs(af.im) < NEGLECT_FLD_CONTRIB) icf = false; if(!std::isfinite(ef.re)) icf = false; if(!std::isfinite(ef.im)) icf = false; icb = true; if(fabs(ab.re) < NEGLECT_FLD_CONTRIB && fabs(ab.im) < NEGLECT_FLD_CONTRIB) icb = false; if(!std::isfinite(eb.re)) icb = false; if(!std::isfinite(eb.im)) icb = false; if(icf==false && icb==false) return; af *= ef; ab *= eb; mvec[i].cfield(EX, x, f, b); if(icf) ex += af*f; if(icb) ex += ab*b; mvec[i].cfield(EY, x, f, b); if(icf) ey += af*f; if(icb) ey += ab*b; mvec[i].cfield(EZ, x, f, b); if(icf) ez += af*f; if(icb) ez += ab*b; mvec[i].cfield(HX, x, f, b); if(icf) hx += af*f; if(icb) hx += ab*b; mvec[i].cfield(HY, x, f, b); if(icf) hy += af*f; if(icb) hy += ab*b; mvec[i].cfield(HZ, x, f, b); if(icf) hz += af*f; if(icb) hz += ab*b; } return; } void ModeArray::emfield(const Cvector& famp, const Cvector& bamp, double x, double z, Complex& ex, Complex& ey, Complex& ez, Complex& hx, Complex& hy, Complex& hz) const { emfield(famp, bamp, x, z, 0.0, 0.0, ex, ey, ez, hx, hy, hz); return; } /* field superposition at point (x, z), amp: complex amplitudes of forward traveling waves at z=0, amplitudes evolve according to amp_j(z) = amp_j(0)*exp(-i cbet_j*z) pert: propagation constant perturbations, with cbet_j = this(j).cbeta(FORW)+pert(j) cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */ Complex ModeArray::pfield(const Cvector& amp, const Cvector& pert, Fcomp cp, double x, double z) const { Complex s, e, f; Complex ex, ey, ez, hx, hy, hz; if( amp.nel < num) slamarrerror("pfield: array length mismatch"); if(pert.nel < num) slamarrerror("pfield: array length mismatch"); s = CC0; switch(cp) { case EX: case EY: case EZ: case HX: case HY: case HZ: for(int i=0; i<=num-1; ++i) { e = mvec[i].efac(FORW, pert(i), z); f = mvec[i].cfield(cp, FORW, x); s += amp(i)*f*e; } return s; case SX: ey = pfield(amp, pert, EY, x, z); ez = pfield(amp, pert, EZ, x, z); hy = pfield(amp, pert, HY, x, z); hz = pfield(amp, pert, HZ, x, z); s = ey*hz.conj()-ez*hy.conj(); return Complex(0.5*s.re, 0.0); case SY: ex = pfield(amp, pert, EX, x, z); ez = pfield(amp, pert, EZ, x, z); hx = pfield(amp, pert, HX, x, z); hz = pfield(amp, pert, HZ, x, z); s = ez*hx.conj()-ex*hz.conj(); return Complex(0.5*s.re, 0.0); case SZ: ex = pfield(amp, pert, EX, x, z); ey = pfield(amp, pert, EY, x, z); hx = pfield(amp, pert, HX, x, z); hy = pfield(amp, pert, HY, x, z); s = ex*hy.conj()-ey*hx.conj(); return Complex(0.5*s.re, 0.0); default: slamarrerror("pfield: mysterious cp"); } return CC0; } /* assuming that all modes belong to the same mode spectrum: check orthogonality & normalization */ double ModeArray::orthonormality() { double devn = 0.0; double devo = 0.0; double o; Mode mj; for(int j=0; j<=num-1; ++j) { mj = mvec[j]; if(mj.ppt == PROPAG) { o = (mj.overlap(FORW, mj, FORW)).re; if(fabs(1-o) > devn) devn = fabs(1-o); } else { o = (mj.overlap(FORW, mj, BACK)).abs(); if(fabs(1-o) > devn) devn = fabs(1-o); } } for(int j=0; j<=num-1; ++j) { mj = mvec[j]; for(int k=0; k<=num-1; ++k) { o = (mj.overlap(FORW, mvec[k], FORW)).abs(); if(j != k) { if(o > devo) devo = o; } } } fprintf(stderr, "DevNorm: %g, DevOrthog: %g\n", devn, devo); if(devo > devn) return devo; return devn; } /* related to the classification of FB modes: remove all non propagating modes from the array */ void ModeArray::restrict2prop() { int done = 0; while(done == 0 && num >= 1) { done=1; int r=-1; for(int j=0; j<=num-1; ++j) { if((mvec[j]).ppt==EVANESC) { done = 0; r = j; } } if(done == 0) remove(r); } return; } /* - Three segment coupler, power transfer evaluation -------------------- imode: input mode, excites the modes this(j) with rel. power 1 at z=0 omode: output mode, relative power is returned pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) (!) this(j) are assumed to be normalized. */ /* weights for the power evaluation: w = ( omode | this(m) ) ( this(m) | imode ) */ double ModeArray::pweight(const Mode& imode, int m, const Mode& omode) { Interval i(-AWAY, AWAY); double om, mi; om = 0.25*(omode.integrate(EX, (mvec[m]), HY, i) -omode.integrate(EY, (mvec[m]), HX, i) +(mvec[m]).integrate(EX, omode, HY, i) -(mvec[m]).integrate(EY, omode, HX, i)); mi = 0.25*((mvec[m]).integrate(EX, imode, HY, i) -(mvec[m]).integrate(EY, imode, HX, i) +imode.integrate(EX, (mvec[m]), HY, i) -imode.integrate(EY, (mvec[m]), HX, i)); return om*mi; } /* single relative power level for device length l */ double ModeArray::iopower(const Mode& imode, const Mode& omode, const Cvector& pert, double l) { Complex a, am, cbet; int m; a=CC0; for(m=0; m<=num-1; ++m) { cbet = pert(m)+mvec[m].beta; cbet = CCI*cbet*(-l); am = exp(cbet)*pweight(imode, m, omode); a = a+am; } return a.sqabs(); } /* output power for numl devices of lengths between lmin and lmax */ Dvector ModeArray::iopower(const Mode& imode, const Mode& omode, const Cvector& pert, int numl, double lbeg, double lend) { Complex a, am, cbet; a=CC0; Dvector w(num); double l, dl; int li, m; if(numl <= 0) slamarrerror("iopower: numl <= 0"); Dvector p(numl); if(numl == 1) { p(0) = iopower(imode, omode, pert, 0.5*(lbeg+lend)); return p; } for(m=0; m<=num-1; ++m) { w(m) = pweight(imode, m, omode); } dl = (lend-lbeg)/(numl-1.0); for(li=0; li<=numl-1; ++li) { a=CC0; l = lbeg+li*dl; for(m=0; m<=num-1; ++m) { cbet = pert(m)+mvec[m].beta; cbet = CCI*cbet*(-l); am = exp(cbet)*w(m); a = a+am; } p(li) = a.sqabs(); } return p; } /* ... write to file */ void ModeArray::writeiopower(const Mode& imode, const Mode& omode, const Cvector& pert, int numl, double lbeg, double lend, char ext0, char ext1) { Dvector p; int i; double l, dl; char name[13] = "t_pow__.xyf"; FILE *dat; if(numl <= 0) slamarrerror("iopower: numl <= 0"); name[1] = polchr(mvec[0].pol); name[5] = ext0; name[6] = ext1; fprintf(stderr, "P(L = %g -> %g) >> %s\n", lbeg, lend, name); p = iopower(imode, omode, pert, numl, lbeg, lend); dat = fopen(name, "w+"); if(p.nel == 1) { fprintf(dat, "%g %.10g\n", 0.5*(lbeg+lend), p(0)); fclose(dat); return; } dl = (lend-lbeg)/(p.nel-1.0); for(i=0; i<=p.nel-1; ++i) { l = lbeg+i*dl; fprintf(dat, "%g %.10g\n", l, p(i)); } fclose(dat); return; } /* - Output to MATLAB .m files ---------------------------------------- */ /* mode interference plots, image of component cp amp_f: amplitudes of the forward propagating modes at z=0 amp_b: amplitudes of the backward propagating modes at z=0 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 ModeArray::plot(const Cvector& amp_f, const Cvector& amp_b, Fcomp cp, Afo foa, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { SegWgStruct st(-1); st(0) = mvec[num-1].wg; BepField fld(st); fld(0) = *this; fld.setf(0, SBRIGHT, amp_f); fld.setf(0, SBLEFT, amp_f); fld.setb(0, SBRIGHT, amp_b); fld.setb(0, SBLEFT, amp_b); fld.check(); fld.plot(cp, foa, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* image of the relative phase amp_f: amplitudes of the forward propagating modes at z=0 amp_b: amplitudes of the backward propagating modes at z=0 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 ModeArray::phasemap(const Cvector& amp_f, const Cvector& amp_b, Fcomp cp, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { SegWgStruct st(-1); st(0) = mvec[num-1].wg; BepField fld(st); fld(0) = *this; fld.setf(0, SBRIGHT, amp_f); fld.setf(0, SBLEFT, amp_f); fld.setb(0, SBRIGHT, amp_b); fld.setb(0, SBLEFT, amp_b); fld.check(); fld.phasemap(cp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* electromagnetic energy density image amp_f: amplitudes of the forward propagating modes at z=0 amp_b: amplitudes of the backward propagating modes at z=0 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 ModeArray::plot(const Cvector& amp_f, const Cvector& amp_b, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { SegWgStruct st(-1); st(0) = mvec[num-1].wg; BepField fld(st); fld(0) = *this; fld.setf(0, SBRIGHT, amp_f); fld.setf(0, SBLEFT, amp_f); fld.setb(0, SBRIGHT, amp_b); fld.setb(0, SBLEFT, amp_b); fld.check(); fld.plot(xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* mode interference plots, fancy style surface ;-) amp_f: amplitudes of the forward propagating modes at z=0 amp_b: amplitudes of the backward propagating modes at z=0 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 ModeArray::fplot(const Cvector& amp_f, const Cvector& amp_b, Fcomp cp, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { SegWgStruct st(-1); st(0) = mvec[num-1].wg; BepField fld(st); fld(0) = *this; fld.setf(0, SBRIGHT, amp_f); fld.setf(0, SBLEFT, amp_f); fld.setb(0, SBRIGHT, amp_b); fld.setb(0, SBLEFT, amp_b); fld.check(); fld.fplot(cp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* animation of the light propagation along the waveguide, the frames show images of EY (TE) resp. HY (TM), at ntfr equidistant times over one time period amp_f: amplitudes of the forward propagating modes at z=0 amp_b: amplitudes of the backward propagating modes at z=0 xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ntfr: number of frames per time period ext0, ext1: filename id characters for the m.file */ void ModeArray::movie(const Cvector& amp_f, const Cvector& amp_b, double xbeg, double xend, double zbeg, double zend, int npx, int npz, int ntfr, char ext0, char ext1) { SegWgStruct st(-1); st(0) = mvec[num-1].wg; BepField fld(st); fld(0) = *this; fld.setf(0, SBRIGHT, amp_f); fld.setf(0, SBLEFT, amp_f); fld.setb(0, SBRIGHT, amp_b); fld.setb(0, SBLEFT, amp_b); fld.check(); fld.movie(xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1); return; } /* animation of the light propagation along the waveguide, fancy style ;-) the frames show surfaces of EY (TE) resp. HY (TM), at ntfr equidistant times over one time period amp_f: amplitudes of the forward propagating modes at z=0 amp_b: amplitudes of the backward propagating modes at z=0 xbeg, xend, zbeg, zend: x resp. z extensions of the plot window npx, npz: number of plot points in the x and z directions ntfr: number of frames per time period ext0, ext1: filename id characters for the m.file */ void ModeArray::fmovie(const Cvector& amp_f, const Cvector& amp_b, double xbeg, double xend, double zbeg, double zend, int npx, int npz, int ntfr, char ext0, char ext1) { SegWgStruct st(-1); st(0) = mvec[num-1].wg; BepField fld(st); fld(0) = *this; fld.setf(0, SBRIGHT, amp_f); fld.setf(0, SBLEFT, amp_f); fld.setb(0, SBRIGHT, amp_b); fld.setb(0, SBLEFT, amp_b); fld.check(); fld.fmovie(xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1); return; } /* export full field data ("all") into a viewer m-file amp_f: amplitudes of the forward propagating modes at z=0 amp_b: amplitudes of the backward propagating modes at z=0 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 ModeArray::viewer(const Cvector& amp_f, const Cvector& amp_b, double xbeg, double xend, double zbeg, double zend, int npx, int npz, char ext0, char ext1) const { SegWgStruct st(-1); st(0) = mvec[num-1].wg; BepField fld(st); fld(0) = *this; fld.setf(0, SBRIGHT, amp_f); fld.setf(0, SBLEFT, amp_f); fld.setb(0, SBRIGHT, amp_b); fld.setb(0, SBLEFT, amp_b); fld.check(); fld.viewer(xbeg, xend, zbeg, zend, npx, npz, ext0, ext1); return; } /* - for the QUEP implementation ---------------------------------------- */ /* overlaps of perpendiculary propagating modes */ Cmatrix ModeArray::crossoverlap(Propdir hdir, const ModeArray& mv, Propdir vdir, double zpos, Interval i) const { fprintf(stderr, "marr: cro: used.\n"); Cmatrix g(num, mv.num); for(int r=0; r<=num-1; ++r) { for(int m=0; m<=mv.num-1; ++m) g(r, m) = mvec[r].crossoverlap(hdir, mv(m), vdir, zpos, i); } return g; } void ModeArray::crossoverlap(const ModeArray& mv, double zpos, Interval i, Cmatrix& cff, Cmatrix& cfb, Cmatrix& cbf, Cmatrix& cbb) const { cff = Cmatrix(num, mv.num); cfb = Cmatrix(num, mv.num); cbf = Cmatrix(num, mv.num); cbb = Cmatrix(num, mv.num); for(int r=0; r<=num-1; ++r) { for(int m=0; m<=mv.num-1; ++m) { mvec[r].crossoverlap(mv(m), zpos, i, cff(r, m), cfb(r, m), cbf(r, m), cbb(r, m)); } } return; } /* - for the vQUEP implementation --------------------------------------- */ /* prepare vTE / vTM mode representation, for given ky, for all modes */ void ModeArray::vectorify(Complex ky) { for(int r=0; r<=num-1; ++r) mvec[r].vectorify(ky); return; }