/* * WMM - Wave-matching method for mode analysis of dielectric waveguides * https://wmm.computational-photonics.eu/ */ /* * wmmarr.cpp * arrays of WMM modes, mode interference evaluation */ #include #include #include #include"inout.h" #include"complex.h" #include"matrix.h" #include"waveg.h" #include"maxweq.h" #include"wmmmat.h" #include"wmmtrif.h" #include"wmmvis.h" #include"mode.h" #include"wmmfld.h" #include"wmmarr.h" /* WMM_ModeArray: Array of WMM_Modes */ /* error message */ void wmmarrerror(const char *msg) { fprintf(stderr, "ERR: %s.\n", msg); exit(1); } /* initialize */ WMM_ModeArray::WMM_ModeArray() : num(0) { int i; for(i=0; i<=MAXNUMMODES-1; ++i) mvec[i] = NULL; } /* destroy */ WMM_ModeArray::~WMM_ModeArray() { int i; for(i=0; i<=MAXNUMMODES-1; ++i) { if(mvec[i] != NULL) { delete mvec[i]; mvec[i] = NULL; } } num = 0; } /* copy constructor */ WMM_ModeArray::WMM_ModeArray(const WMM_ModeArray& ma) : num(ma.num) { int i; for(i=0; i<=num-1; ++i) { mvec[i] = new WMM_Mode; *mvec[i] = *(ma.mvec[i]); } } /* assignment */ WMM_ModeArray& WMM_ModeArray::operator=(const WMM_ModeArray& ma) { int i; if(this != &ma) { for(i=0; i<=MAXNUMMODES-1; ++i) { if(mvec[i] != NULL) { delete mvec[i]; mvec[i] = NULL; } } num = ma.num; for(i=0; i<=num-1; ++i) { mvec[i] = new WMM_Mode; *mvec[i] = *(ma.mvec[i]); } } return *this; } /* delete all Mode entries */ void WMM_ModeArray::clear() { int i; for(i=0; i<=MAXNUMMODES-1; ++i) { if(mvec[i] != NULL) { delete mvec[i]; mvec[i] = NULL; } } num = 0; } /* input from FILE dat */ void WMM_ModeArray::read(FILE *dat) { int i; if(dat == stderr || dat == stdout) { wmmarrerror("WMM_ModeArray: read"); } else { for(i=0; i<=MAXNUMMODES-1; ++i) { if(mvec[i] != NULL) { delete mvec[i]; mvec[i] = NULL; } } num = 0; num = fgetint(dat); if(num >= MAXNUMMODES) { fprintf(stderr, "WMM_ModeArray: read: NUM reduced to %d.\n", MAXNUMMODES); num = MAXNUMMODES; } for(i=0; i<=num-1; ++i) { mvec[i] = new WMM_Mode; (*mvec[i]).read(dat); } } return; } /* output to FILE dat */ void WMM_ModeArray::write(FILE *dat) { int i; comment("WMM_ModeArray", dat); comment("num", dat); fputint(num, dat); for(i=0; i<=num-1; ++i) { comment("*** -------- ***", dat); (*mvec[i]).write(dat); } return; } /* input from default file */ void WMM_ModeArray::read_def(char ext0, char ext1) { char name[10]; FILE *dat; int i; name[0] = 'w'; name[1] = 'm'; name[2] = 'm'; name[3] = ext0; name[4] = ext1; name[5] = '.'; name[6] = 'm'; name[7] = 'a'; name[8] = 'r'; 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) %c%c%c, beta: %.10g, neff: %.10g, npcB: %.10g\n", i, poltochar1((*mvec[i]).pol), poltochar2((*mvec[i]).pol), symtochar((*mvec[i]).sym), (*mvec[i]).beta, (*mvec[i]).neff, (*mvec[i]).npcB); return; } /* output to default file */ void WMM_ModeArray::write_def(char ext0, char ext1) { char name[10]; FILE *dat; int i; name[0] = 'w'; name[1] = 'm'; name[2] = 'm'; name[3] = ext0; name[4] = ext1; name[5] = '.'; name[6] = 'm'; name[7] = 'a'; name[8] = 'r'; name[9] = 0; fprintf(stderr, ">> %s #%d\n", name, num); for(i=0; i<=num-1; ++i) fprintf(stderr, " (%d) %c%c%c, beta: %.10g, neff: %.10g, npcB: %.10g\n", i, poltochar1((*mvec[i]).pol), poltochar2((*mvec[i]).pol), symtochar((*mvec[i]).sym), (*mvec[i]).beta, (*mvec[i]).neff, (*mvec[i]).npcB); dat = fopen(name, "w+"); write(dat); fclose(dat); return; } /* subscripting */ WMM_Mode& WMM_ModeArray::operator() (int i) { if(i<0 || i>=num) wmmarrerror("WMM_ModeArray: () out of range"); return *mvec[i]; } WMM_Mode WMM_ModeArray::operator() (int i) const { if(i<0 || i>=num) wmmarrerror("WMM_ModeArray: () out of range"); return *mvec[i]; } /* add a mode */ void WMM_ModeArray::add(WMM_Mode m) { if(num >= MAXNUMMODES) wmmarrerror("WMM_ModeArray: add: MAXNUMMODES exceeded"); if(mvec[num] != NULL) wmmarrerror("WMM_ModeArray: add: mvec corrupted"); mvec[num] = new WMM_Mode; *mvec[num] = m; ++num; return; } /* delete a mode entry */ void WMM_ModeArray::remove(int i) { int j; if(i<0 || i>=num) wmmarrerror("WMM_ModeArray: remove: argument out of range"); if(mvec[i] == NULL) wmmarrerror("WMM_ModeArray: remove: mvec corrupted"); delete mvec[i]; mvec[i] = NULL; if(i==num-1) { --num; return; } for(j=i; j<=num-2; ++j) mvec[j] = mvec[j+1]; mvec[num-1] = NULL; --num; return; } /* add an entire WMM_ModeArray nma, nma is cleared ! */ void WMM_ModeArray::merge(WMM_ModeArray& ma) { int j; if(num+ma.num >= MAXNUMMODES) wmmarrerror("WMM_ModeArray: merge: MAXNUMMODES exceeded"); if(mvec[num] != NULL) wmmarrerror("WMM_ModeArray: merge: mvec corrupted"); for(j=0; j<=ma.num-1; ++j) { if(mvec[num] != NULL) wmmarrerror("WMM_ModeArray: merge: mvec corrupted"); mvec[num] = ma.mvec[j]; ++num; ma.mvec[j] = NULL; } ma.num = 0; return; } /* sort the array by propagation constants, highest first */ void WMM_ModeArray::sort() { int j, k, maxi; double maxb; WMM_ModeP t; if(num<=1) return; for(j=0; j<=num-2; ++j) { maxi = j; maxb = (*mvec[j]).beta; for(k=j+1; k<=num-1; ++k) { if(maxb<(*mvec[k]).beta) { maxb = (*mvec[k]).beta; maxi = k; } } t = mvec[j]; mvec[j] = mvec[maxi]; mvec[maxi] = t; } return; } /* ................................................................... */ /* * mode interference evaluation and visualization * all modes are assumed to belong to the same waveguide ! */ /* field superposition at point (x, y, z), amp: complex amplitudes at z=0 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z) pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) cp: EX, EY, EZ, HX, HY, HZ, SZ */ Complex WMM_ModeArray::field(Cvector amp, Cvector pert, Fcomp cp, double x, double y, double z) const { Complex s, iphase; double f, sum; Cvector a; int i, j, k; a = amp; if(amp.nel < num) { a = Cvector(num); a.init(CC0); for(i=0; i<=amp.nel-1; ++i) a(i) = amp(i); } Cvector cbet(num); for(i=0; i<=num-1; ++i) cbet(i) = pert(i)+(*mvec[i]).beta; if(z != 0.0) { for(i=0; i<=num-1; ++i) { iphase = CCI*cbet(i)*(-z); a(i) = a(i)*exp(iphase); } } s = CC0; if(cp != SZ) { for(i=0; i<=num-1; ++i) { f = (*mvec[i]).field(cp, x, y); if(cp != EZ && cp != HZ) s = s+a(i)*f; else { a(i) = CCI*a(i); s = s+a(i)*f; } } return s; } sum = 0.0; Dvector fEx(num); Dvector fEy(num); Dvector fHx(num); Dvector fHy(num); for(j=0; j<=num-1; ++j) fEx(j) = (*mvec[j]).field(EX, x, y); for(j=0; j<=num-1; ++j) fEy(j) = (*mvec[j]).field(EY, x, y); for(j=0; j<=num-1; ++j) fHx(j) = (*mvec[j]).field(HX, x, y); for(j=0; j<=num-1; ++j) fHy(j) = (*mvec[j]).field(HY, x, y); for(j=0; j<=num-1; ++j) { for(k=0; k<=num-1; ++k) { f = (a(j).re*a(k).re + a(j).im*a(k).im); sum += f*(fEx(j)*fHy(k)-fEy(j)*fHx(k)); } } s = Complex(0.5*sum); return s; } /* evaluate component cp on a rectangular npx x npy mesh on area disp at position z amp: complex amplitudes at z=0 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z) pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) cp: EX - HZ, SZ foa: MOD, SQR, REP, IMP */ Dmatrix WMM_ModeArray::fieldmat(Cvector amp, Cvector pert, double z, Rect disp, int npx, int npy, Fcomp cp, Afo foa) const { double x, y; double dx, dy; int i, j; if(npx <= 1) wmmarrerror("fieldmat: npx <= 1"); if(npy <= 1) wmmarrerror("fieldmat: npy <= 1"); Dmatrix f(npx, npy); Complex ft; double ftd; dx = (disp.x1-disp.x0)/(npx-1); dy = (disp.y1-disp.y0)/(npy-1); for(i=0; i<=npx-1; ++i) { x = disp.x0+i*dx; for(j=0; j<=npy-1; ++j) { y = disp.y0+j*dy; ft = field(amp, pert, cp, x, y, z); ftd = 0.0; switch(foa) { case MOD: case ORG: ftd = ft.abs(); break; case SQR: ftd = ft.sqabs(); break; case REP: ftd = ft.re; break; case IMP: ftd = ft.im; break; } f(i, j) = ftd; } } return f; } /* store nump values of component cp between (x0, y0) and (x1, y1) in a vector at position z amp: complex amplitudes at z=0 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z) pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) cp: EX - HZ, SZ foa: MOD, SQR, REP, IMP */ Dvector WMM_ModeArray::fieldvec(Cvector amp, Cvector pert, double z, Fcomp cp, Afo foa, int nump, double x0, double y0, double x1, double y1) const { double x, y; double dx, dy; int j; if(nump <= 1) wmmarrerror("fieldvec: nump <= 1"); Dvector f(nump); Complex ft; double ftd; dx = (x1-x0)/(nump-1); dy = (y1-y0)/(nump-1); for(j=0; j<=nump-1; ++j) { x = x0+j*dx; y = y0+j*dy; ft = field(amp, pert, cp, x, y, z); ftd = 0.0; switch(foa) { case MOD: case ORG: ftd = ft.abs(); break; case SQR: ftd = ft.sqabs(); break; case REP: ftd = ft.re; break; case IMP: ftd = ft.im; break; } f(j) = ftd; } return f; } /* - 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 WMM_ModeArray::pweight(const WMM_Mode& imode, int m, const WMM_Mode& omode) { return lscalprod(omode, (*mvec[m]))*lscalprod((*mvec[m]), imode); } /* single relative power level for device length l */ double WMM_ModeArray::iopower(const WMM_Mode& imode, const WMM_Mode& omode, 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 WMM_ModeArray::iopower(const WMM_Mode& imode, const WMM_Mode& omode, 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) wmmarrerror("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 WMM_ModeArray::writeiopower(const WMM_Mode& imode, const WMM_Mode& omode, Cvector pert, int numl, double lbeg, double lend, char ext0, char ext1) { Dvector p; int i; double l, dl; char name[13] = "__pow__.xyf"; FILE *dat; if(numl <= 0) wmmarrerror("iopower: numl <= 0"); switch((*mvec[0]).pol) { case QTE: name[0] = 't'; name[1] = 'e'; break; case QTM: name[0] = 't'; name[1] = 'm'; break; case VEC: name[0] = 'v'; name[1] = 'c'; break; } 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 ---------------------------------------- */ /* write single component of the interference field at position z to MATLAB .m file amp: complex amplitudes at z=0 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z) pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) cp: EX - HZ, SZ foa: MOD, ORG, SQR, REP, IMP disp: output region on the x-y-plane npx, npy: number of points in output mesh ext0, ext1: filename id characters pltype: 'C': contour plot 'S': surface plot 'I': intensity image 'N': field + mesh only, no plot commands (default) */ void WMM_ModeArray::mfile(Cvector amp, Cvector pert, double z, Fcomp cp, Afo foa, Rect disp, int npx, int npy, char ext0, char ext1, char pltype) const { FILE *dat; char name[13] = "int______.m"; double rz; double max, f, maxsum, minf, maxf; int i; rz = floor(z*10.0+0.5)/10.0; switch(foa) { case MOD: name[3] = 'm'; break; case ORG: name[3] = 'f'; break; case SQR: name[3] = 'q'; break; case REP: name[3] = 'r'; break; case IMP: name[3] = 'i'; break; } name[4] = fldchr(cp); name[5] = cpchr(cp); name[6] = ext0; name[7] = ext1; name[8] = pltype; fprintf(stderr, "%c%c(z=%g) >> %s\n", fldchr(cp), cpchr(cp), z, name); dat = fopen(name, "w+"); mlout_title(dat, name, "WMM interference field"); maxf = 1.0; minf = -1.0; switch(fldchr(cp)) { case 'E': maxsum = 0.0; for(i=0; i<=num-1; ++i) { max = fabs((*mvec[i]).maxE)*amp(i).abs(); f = fabs((*mvec[i]).minE)*amp(i).abs(); if(max < f) max = f; maxsum += max; } switch(foa) { case MOD: case ORG: maxf = maxsum; minf = 0.0; break; case SQR: maxf = maxsum*maxsum; minf = 0.0; break; case REP: case IMP: maxf = maxsum; minf = -maxsum; break; } break; case 'H': maxsum = 0.0; for(i=0; i<=num-1; ++i) { max = fabs((*mvec[i]).maxH)*amp(i).abs(); f = fabs((*mvec[i]).minH)*amp(i).abs(); if(max < f) max = f; maxsum += max; } switch(foa) { case MOD: case ORG: maxf = maxsum; minf = 0.0; break; case SQR: maxf = maxsum*maxsum; minf = 0.0; break; case REP: case IMP: maxf = maxsum; minf = -maxsum; break; } break; case 'S': maxsum = 0.0; for(i=0; i<=num-1; ++i) { max = fabs((*mvec[i]).maxS)*amp(i).abs()*amp(i).abs(); f = fabs((*mvec[i]).minS)*amp(i).abs()*amp(i).abs(); if(max < f) max = f; maxsum += sqrt(max); } switch(foa) { case MOD: case ORG: case REP: maxf = maxsum*maxsum; minf = 0.0; break; case SQR: // not reasonable maxf = maxsum*maxsum*maxsum*maxsum; minf = 0.0; break; case IMP: // not reasonable maxf = 1.0; minf = 0.0; } break; } mlout_geo(dat, (*mvec[0]).wg, minf, maxf); mlout_meshxy(dat, disp, npx, npy); Dmatrix fld; fld = fieldmat(amp, pert, z, disp, npx, npy, cp, foa); mlout_fld(dat, npx, npy, cp, fld); name[9] = 0; switch(pltype) { case 'C': mlout_contour(name, dat, cp, foa); mlout_annotatezpos(dat, rz); break; case 'S': mlout_surface(name, dat, cp, foa); break; case 'I': mlout_image(name, dat, cp, foa); mlout_annotatezpos(dat, rz); break; default: break; } mlout_print(dat, name, 'e'); fclose(dat); return; } /* write interference field at position z to MATLAB .m file, fancy style :-) amp: complex amplitudes at z=0 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z) pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) disp: output region on the x-y-plane npx, npy: number of points in output mesh ext0, ext1: filename id characters */ #define HDIST 1.0e-8 void WMM_ModeArray::fancymfile(Cvector amp, Cvector pert, double z, Rect disp, int npx, int npy, char ext0, char ext1) const { FILE *dat; int np, l, m; double x0, x1, xp, y0, y1, yp; int numc; double max, f, maxsum, minf, maxf; int i; Waveguide wg; char name[13] = "intfSz__F.m"; name[6] = ext0; name[7] = ext1; fprintf(stderr, "Sz(z=%g) >> %s\n", z, name); dat = fopen(name, "w+"); mlout_title(dat, name, "WMM interference field :-)"); name[9] = 0; wg = (*mvec[0]).wg; maxsum = 0.0; for(i=0; i<=num-1; ++i) { max = fabs((*mvec[i]).maxS)*amp(i).abs()*amp(i).abs(); f = fabs((*mvec[i]).minS)*amp(i).abs()*amp(i).abs(); if(max < f) max = f; maxsum += sqrt(max); } maxf = maxsum*maxsum; minf = 0.0; mlout_geo(dat, wg, minf, maxf); mlout_meshxy(dat, disp, npx, npy); Dmatrix fld; fld = fieldmat(amp, pert, z, disp, npx, npy, SZ, ORG); mlout_fld(dat, npx, npy, SZ, fld); Dvector fv; numc = 0; for(l=0; l<=wg.nx+1; ++l) { if(l==0) x0 = disp.x0; else x0 = wg.hx(l-1); if(l==wg.nx+1) x1 = disp.x1; else x1 = wg.hx(l); for(m=0; m<=wg.ny; ++m) { if(fabs(wg.n(l,m)-wg.n(l,m+1)) > 1.0e-10) { yp = wg.hy(m); np = (int)(((double) npx) *(x1-x0)/(disp.x1-disp.x0)); if(np >= 2) { fv = fieldvec(amp, pert, z, SZ, ORG, np, x0+HDIST, yp-HDIST, x1-HDIST, yp-HDIST); mlout_sec(dat, x0+HDIST, yp-HDIST, x1-HDIST, yp-HDIST, np, SZ, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(amp, pert, z, SZ, ORG, np, x0+HDIST, yp+HDIST, x1-HDIST, yp+HDIST); mlout_sec(dat, x0+HDIST, yp+HDIST, x1-HDIST, yp+HDIST, np, SZ, dig10(numc), dig1(numc), fv); ++numc; } } } } for(m=0; m<=wg.ny+1; ++m) { if(m==0) y0 = disp.y0; else y0 = wg.hy(m-1); if(m==wg.ny+1) y1 = disp.y1; else y1 = wg.hy(m); for(l=0; l<=wg.nx; ++l) { if(fabs(wg.n(l,m)-wg.n(l+1,m)) > 1.0e-10) { xp = wg.hx(l); np = (int)(((double)npy) *(y1-y0)/(disp.y1-disp.y0)); if(np >= 2) { fv = fieldvec(amp, pert, z, SZ, ORG, np, xp-HDIST, y0+HDIST, xp-HDIST, y1-HDIST); mlout_sec(dat, xp-HDIST, y0+HDIST, xp-HDIST, y1-HDIST, np, SZ, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(amp, pert, z, SZ, ORG, np, xp+HDIST, y0+HDIST, xp+HDIST, y1-HDIST); mlout_sec(dat, xp+HDIST, y0+HDIST, xp+HDIST, y1-HDIST, np, SZ, dig10(numc), dig1(numc), fv); ++numc; } } } } fv = fieldvec(amp, pert, z, SZ, ORG, npx, disp.x0, disp.y0, disp.x1, disp.y0); mlout_sec(dat, disp.x0, disp.y0, disp.x1, disp.y0, npx, SZ, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(amp, pert, z, SZ, ORG, npy, disp.x1, disp.y0, disp.x1, disp.y1); mlout_sec(dat, disp.x1, disp.y0, disp.x1, disp.y1, npy, SZ, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(amp, pert, z, SZ, ORG, npx, disp.x1, disp.y1, disp.x0, disp.y1); mlout_sec(dat, disp.x1, disp.y1, disp.x0, disp.y1, npx, SZ, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(amp, pert, z, SZ, ORG, npy, disp.x0, disp.y1, disp.x0, disp.y0); mlout_sec(dat, disp.x0, disp.y1, disp.x0, disp.y0, npy, SZ, dig10(numc), dig1(numc), fv); ++numc; mlout_fancy(name, dat, SZ, numc); mlout_print(dat, name, 'p'); fclose(dat); return; } /* animate the interference field amp: complex amplitudes at z=0 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z) pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) cp: EX - HZ, SZ foa: MOD, ORG, SQR, REP, IMP disp: output region on the x-y-plane npx, npy: number of points in output mesh pltype: 'C': contour plot 'S': surface plot 'I': intensity image 'F': fancy plot, cp and foa are set to SZ, MOD (default) z0, z1, numz: .m files are generated for z-positions z0+i*(z1-z0)/numz, i=0, ..., numz-1 */ void WMM_ModeArray::movie(Cvector amp, Cvector pert, Fcomp cp, Afo foa, Rect disp, int npx, int npy, char pltype, int numz, double z0, double z1) const { int i; double z; FILE *dat; char name[13] = "int___plM.m"; if(numz >= 100) numz = 100; if(pltype == 'F') { cp = SZ; foa = ORG; } switch(foa) { case MOD: name[3] = 'm'; break; case ORG: name[3] = 'f'; break; case SQR: name[3] = 'q'; break; case REP: name[3] = 'r'; break; case IMP: name[3] = 'i'; break; } name[4] = fldchr(cp); name[5] = cpchr(cp); for(i=0; i<=numz-1; ++i) { z = z0+i*(z1-z0)/numz; switch(pltype) { case 'C': case 'S': case 'I': mfile(amp, pert, z, cp, foa, disp, npx, npy, dig10(i), dig1(i), pltype); break; default: fancymfile(amp, pert, z, disp, npx, npy, dig10(i), dig1(i)); break; } } fprintf(stderr, ">> %s\n", name); dat = fopen(name, "w+"); mlout_title(dat, name, "WMM interference field animation"); name[8] = pltype; name[9] = 0; mlout_play(dat, name, numz); fclose(dat); return; } /* interference pattern on the horizontal plane ybeg <= y <= yend, zbeg <= z <= zend at position x image plot corresponding to the squareroot of the local intensity (SZ) amp: complex amplitudes at z=0 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z) pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) npy, npz: number of points in output mesh ext0, ext1: filename id characters */ void WMM_ModeArray::prop(Cvector amp, Cvector pert, double x, double ybeg, double yend, int npy, double zbeg, double zend, int npz, char ext0, char ext1) const { FILE *dat; char name[13] = "propSz__I.m"; double max, f, maxsum, minf, maxf; double y, dy, z, dz; Complex ph; int i, j, k, m; name[6] = ext0; name[7] = ext1; fprintf(stderr, "prop(z=%g -> %g) >> %s\n", zbeg, zend, name); dat = fopen(name, "w+"); mlout_title(dat, name, "WMM interference field"); maxsum = 0.0; for(i=0; i<=num-1; ++i) { max = fabs((*mvec[i]).maxS)*amp(i).abs()*amp(i).abs(); f = fabs((*mvec[i]).minS)*amp(i).abs()*amp(i).abs(); if(max < f) max = f; maxsum += sqrt(max); } maxf = maxsum; minf = 0.0; mlout_geoyz(dat, (*mvec[0]).wg, minf, maxf); mlout_meshyz(dat, ybeg, yend, npy, zbeg, zend, npz); Dmatrix fld(npy, npz); dy = (yend-ybeg)/(npy-1); dz = (zend-zbeg)/(npz-1); Dmatrix pmprofEx(num, npy); Dmatrix pmprofEy(num, npy); Dmatrix pmprofHx(num, npy); Dmatrix pmprofHy(num, npy); for(i=0; i<=npy-1; ++i) { y = ybeg+i*dy; for(m=0; m<=num-1; ++m) { pmprofEx(m, i) = (*mvec[m]).field(EX, x, y); pmprofEy(m, i) = (*mvec[m]).field(EY, x, y); pmprofHx(m, i) = (*mvec[m]).field(HX, x, y); pmprofHy(m, i) = (*mvec[m]).field(HY, x, y); } } Cvector cbet(num); Cvector ma(num); for(m=0; m<=num-1; ++m) cbet(m) = pert(m)+(*mvec[m]).beta; for(j=0; j<=npz-1; ++j) { z = zbeg+j*dz; for(m=0; m<=num-1; ++m) { ph = CCI*cbet(m)*(-z); ma(m) = amp(m)*exp(ph); } for(i=0; i<=npy-1; ++i) { y = ybeg+i*dy; f = 0.0; for(m=0; m<=num-1; ++m) { for(k=0; k<=num-1; ++k) { f += (ma(m).re*ma(k).re+ma(m).im*ma(k).im) *(pmprofEx(m,i)*pmprofHy(k,i)-pmprofEy(m,i)*pmprofHx(k,i)); } } fld(i, j) = sqrt(0.5*fabs(f)); } } mlout_fld(dat, npy, npz, SZ, fld); name[9] = 0; mlout_propimage(name, dat); mlout_print(dat, name, 'e'); fclose(dat); return; } /* Three segment coupler, interference pattern on the horizontal plane ybeg <= y <= yend, zbeg <= z <= zend at elevation x image plot corresponding to the squareroot of the local intensity (SZ) z < 0: intensity profile of the input modes 0 < z < l: mode interference pattern l < z: intensity profiles of the output modes, scaled by the output amplitudes (!) this(j) are assumed to be normalized. iamp: complex input mode amplitudes imodes: input modes, excite the modes this(j) at z=0, belonging to well separated port waveguides iwg: the geometry of the combined input waveguides omodes: output modes, belonging to well separated port waveguides owg: the geometry of the combined output waveguides iwg and owg are for displaying purposes only pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) npy, npz: number of points in output mesh ext0, ext1: filename id characters */ #define ZDIST 1.0e-8 void WMM_ModeArray::ioprop(Cvector iamp, const WMM_ModeArray& imodes, Waveguide iwg, const WMM_ModeArray& omodes, Waveguide owg, Cvector pert, double x, double l, double ybeg, double yend, int npy, double zbeg, double zend, int npz, char ext0, char ext1) const { FILE *dat; char name[13] = "propSz__I.m"; double max, f, maxsum, minf, maxf; int i, j, m, k; double y, z, dy, dz; Complex a, aj; Complex ph; double s; name[6] = ext0; name[7] = ext1; if(zbeg >= zend) wmmarrerror("ioprop: invalid z-range"); fprintf(stderr, "ioprop(z=%g -> %g) >> %s\n", zbeg, zend, name); dat = fopen(name, "w+"); mlout_title(dat, name, "WMM interference field"); maxsum = 0.0; for(i=0; i<=imodes.num-1; ++i) { max = fabs(imodes(i).maxS)*iamp(i).abs()*iamp(i).abs(); f = fabs(imodes(i).minS)*iamp(i).abs()*iamp(i).abs(); if(max < f) max = f; maxsum += sqrt(max); } maxf = maxsum; minf = 0.0; mlout_iopropgeo(dat, iwg, (*mvec[0]).wg, owg, l, minf, maxf); mlout_meshyz(dat, ybeg, yend, npy, zbeg, zend, npz); s = 0.0; for(i=0; i<=imodes.num-1; ++i) { fprintf(stderr, " A^I_%d = %g + i %g, |A|^2: %g\n", i, iamp(i).re, iamp(i).im, iamp(i).sqabs()); s += iamp(i).sqabs(); } fprintf(stderr, " total input power: %g\n", s); Cvector pamp(num); for(i=0; i<=num-1; ++i) { a = CC0; for(j=0; j<=imodes.num-1; ++j) { aj = iamp(j)*lscalprod((*mvec[i]), imodes(j)); a = a+aj; } pamp(i) = a; } s=0.0; for(i=0; i<=num-1; ++i) { fprintf(stderr, " A^II_%d = %g + i %g, |A|^2: %g\n", i, pamp(i).re, pamp(i).im, pamp(i).sqabs()); s += pamp(i).sqabs(); } fprintf(stderr, " propagating power: %g\n", s); Cvector oamp(omodes.num); for(i=0; i<=omodes.num-1; ++i) { a = CC0; for(j=0; j<=num-1; ++j) { aj = pamp(j)*lscalprod(omodes(i), (*mvec[j])); ph = pert(j)+(*mvec[j]).beta; ph = CCI*ph*(-l); aj = aj*exp(ph); a = a+aj; } oamp(i) = a; } s=0.0; for(i=0; i<=omodes.num-1; ++i) { fprintf(stderr, " A^III_%d = %g + i %g, |A|^2: %g\n", i, oamp(i).re, oamp(i).im, oamp(i).sqabs()); s += oamp(i).sqabs(); } fprintf(stderr, " total output power: %g\n", s); Dmatrix fld(npy, npz); dy = (yend-ybeg)/(npy-1); dz = (zend-zbeg)/(npz-1); Dvector improf(npy); Dvector omprof(npy); for(i=0; i<=npy-1; ++i) { y = ybeg+i*dy; improf(i) = 0.0; for(m=0; m<=imodes.num-1; ++m) { improf(i) += imodes(m).field(SZ, x, y) *iamp(m).sqabs(); } omprof(i) = 0.0; for(m=0; m<=omodes.num-1; ++m) { omprof(i) += omodes(m).field(SZ, x, y) *oamp(m).sqabs(); } } Dmatrix pmprofEx(num, npy); Dmatrix pmprofEy(num, npy); Dmatrix pmprofHx(num, npy); Dmatrix pmprofHy(num, npy); for(i=0; i<=npy-1; ++i) { y = ybeg+i*dy; for(m=0; m<=num-1; ++m) { pmprofEx(m, i) = (*mvec[m]).field(EX, x, y); pmprofEy(m, i) = (*mvec[m]).field(EY, x, y); pmprofHx(m, i) = (*mvec[m]).field(HX, x, y); pmprofHy(m, i) = (*mvec[m]).field(HY, x, y); } } Cvector cbet(num); Cvector ma(num); for(m=0; m<=num-1; ++m) cbet(m) = pert(m)+(*mvec[m]).beta; for(j=0; j<=npz-1; ++j) { z = zbeg+j*dz; if(z<= 0.0) { for(i=0; i<=npy-1; ++i) fld(i, j) = sqrt(improf(i)); } else { if(z>=l) { for(i=0; i<=npy-1; ++i) fld(i, j) = sqrt(omprof(i)); } else { for(m=0; m<=num-1; ++m) { ph = CCI*cbet(m)*(-z); ma(m) = pamp(m)*exp(ph); } for(i=0; i<=npy-1; ++i) { y = ybeg+i*dy; f = 0.0; for(m=0; m<=num-1; ++m) { for(k=0; k<=num-1; ++k) { f += (ma(m).re*ma(k).re+ma(m).im*ma(k).im) *(pmprofEx(m,i)*pmprofHy(k,i)-pmprofEy(m,i)*pmprofHx(k,i)); } } fld(i, j) = sqrt(0.5*fabs(f)); } } } } mlout_fld(dat, npy, npz, SZ, fld); name[9] = 0; mlout_iopropimage(name, dat); mlout_print(dat, name, 'e'); fclose(dat); return; } /* ... as before, but with a single input mode imode only */ void WMM_ModeArray::ioprop(const WMM_Mode& imode, const WMM_ModeArray& omodes, Waveguide owg, Cvector pert, double x, double l, double ybeg, double yend, int npy, double zbeg, double zend, int npz, char ext0, char ext1) const { WMM_ModeArray imodes; imodes.add(imode); Cvector iamp(1); iamp(0) = CC1; ioprop(iamp, imodes, imode.wg, omodes, owg, pert, x, l, ybeg, yend, npy, zbeg, zend, npz, ext0, ext1); return; } /* ... as before, with single input and output modes imode and omode */ void WMM_ModeArray::ioprop(const WMM_Mode& imode, const WMM_Mode& omode, Cvector pert, double x, double l, double ybeg, double yend, int npy, double zbeg, double zend, int npz, char ext0, char ext1) const { WMM_ModeArray omodes; omodes.add(omode); ioprop(imode, omodes, omode.wg, pert, x, l, ybeg, yend, npy, zbeg, zend, npz, ext0, ext1); return; }