/* * METRIC --- Mode expansion modeling in integrated optics / photonics * http://metric.computational-photonics.eu/ */ /* * structure.cpp * Multilayer slab waveguides, waveguide sequences, * rectangular optical circuits. */ #include #include #include #include"inout.h" #include"complex.h" #include"matrix.h" #include"structure.h" #include"gengwed.h" #include"matlvis.h" /* error message */ void structerror(const char *m) { fprintf(stderr, "\nstructure: %s.\n", m); exit(1); } /* ----------------------------------------------------------------------- */ /* a normalized parameter that specifies whether localized fields related to pieces of a waveguide are used to establish approximate modal basis sets */ double DECOUPLED = 1.0e-6; /* ------------------------------------------------------------------- */ /* interval from the waveguide cross section */ /* initialize */ Interval::Interval() { x0 = 0.0; x1 = 1.0; } Interval::Interval(double xa, double xb) { if(xa <= xb) { x0 = xa; x1 = xb; } else { x0 = xb; x1 = xa; } } /* output to FILE dat */ void Interval::write(FILE *dat) const { if(dat == stderr || dat == stdout) { fprintf(dat, "[%g, %g]\n", x0, x1); } else { fputdouble(x0, dat); fputdouble(x1, dat); } return; } /* input from FILE dat */ void Interval::read(FILE *dat) { if(dat == stderr || dat == stdout) { structerror("Interval: read"); } else { x0 = fgetdouble(dat); x1 = fgetdouble(dat); } return; } /* compare the present object to i, return 1 if equal, otherwise return 0 */ int Interval::equal(const Interval& i) const { if(this == &i) return 1; if(fabs(x0-i.x0) > COMPTOL_IV) return 0; if(fabs(x1-i.x1) > COMPTOL_IV) return 0; return 1; } /* restrict iv to the part inside the present object, return value: 0: no overlap, 1: something from iv remains */ int Interval::clip(Interval& iv) const { if(x0 >= x1) return 0; if(iv.x0 >= iv.x1) return 0; if(iv.x1 <= x0) return 0; if(iv.x0 >= x1) return 0; if(x0 > iv.x0) iv.x0 = x0; if(x1 < iv.x1) iv.x1 = x1; return 1; } /* ------------------------------------------------------------------- */ /* waveguide geometry + permittivity profile; vacuum wavelength */ /* initialize */ Waveguide::Waveguide() { nx = 0; hx = Dvector(); lambda = 1.0; n = Dvector(); special = 0; } Waveguide::Waveguide(int vnx) { nx = vnx; hx = Dvector(nx+1); lambda = 1.0; n = Dvector(nx+2); special = 0; } Waveguide::Waveguide(int vnx, Dvector vhx, double l, Dvector vn) { nx = vnx; hx = vhx; lambda = l; n = vn; special = 0; } /* free allocated memory */ void Waveguide::strip() { nx = 0; hx.strip(); n.strip(); } /* get Interval index corresponding to position x */ int Waveguide::layeridx(double x) const { int l=0; while(l<=nx && x>hx(l)) ++l; return l; } /* get layer boundaries corresponding to index l */ Interval Waveguide::layer(int l) const { Interval i; if(l <= 0) i.x0 = -AWAY; else i.x0 = hx(l-1); if(l >= nx+1) i.x1 = AWAY; else i.x1 = hx(l); return i; } /* get layer boundaries corresponding to position x */ Interval Waveguide::layer(double x) const { return layer(layeridx(x)); } /* permittivity on layer l */ double Waveguide::eps(int l) const { double e; e = n(l); if(special) return e; else return e*e; } /* permittivity at position x */ double Waveguide::eps(double x) const { double e; int l; l = layeridx(x); e = n(l); if(special) return e; else return e*e; } /* upper bound for effective mode permittivities, default value */ double Waveguide::defaultepseffmax() const { int l; double eeffmax = -MAXREFIND*MAXREFIND; for(l=1; l<=nx; ++l) { if(eeffmax < eps(l)) eeffmax = eps(l); } return eeffmax; } /* lower bound for effective mode permittivities, default value */ double Waveguide::defaultepseffmin() const { double eeffmin; eeffmin = eps(0); if(eeffmin hx+dx */ void Waveguide::translate(double dx) { int l; for(l=0; l<=nx; ++l) hx(l) += dx; return; } /* output to FILE dat */ void Waveguide::write(FILE *dat) const { int i; if(dat == stderr || dat == stdout) { comment("Waveguide", dat); fprintf(dat, "Nx: %d, Hx: ", nx); for(i=0; i<=nx; ++i) fprintf(dat, "%g ", hx(i)); fprintf(dat, "\n"); fprintf(dat, "lambda: %g\n", lambda); fprintf(dat, "n:\n"); for(i=0; i<=nx+1; ++i) fprintf(dat, "%g ", n(i)); fprintf(dat, "\n"); fprintf(dat, "eps:\n"); for(i=0; i<=nx+1; ++i) fprintf(dat, "%g ", eps(i)); fprintf(dat, "\n"); fprintf(dat, "special: %d\n", special); } else { comment("Waveguide", dat); comment("nx", dat); fputint(nx, dat); comment("hx[0..nx]", dat); hx.write(dat); comment("lambda", dat); fputdouble(lambda, dat); comment("n[0..nx+1]", dat); n.write(dat); comment("special", dat); fputint(special, dat); } return; } /* input from FILE dat */ void Waveguide::read(FILE *dat) { if(dat == stderr || dat == stdout) { structerror("Waveguide: read"); } nx = fgetint(dat); hx.read(dat); lambda = fgetdouble(dat); n.read(dat); special = fgetint(dat); return; } /* test matching of boundary positions with another waveguide, all boundaries of this object should appear in w */ int Waveguide::bdmatch(Waveguide w) const { int j, k, e; for(j=0; j<=w.nx; ++j) { e = 0; for(k=0; k<=nx; ++k) { if(fabs(hx(k)-w.hx(j))<=HDIST) e = 1; } if(e != 1) return 1; } return 0; } /* compare the present object to wg, return 1 if equal, otherwise return 0 */ int Waveguide::equal(const Waveguide& wg) const { if(this == &wg) return 1; if(nx != wg.nx) return 0; if(fabs(lambda-wg.lambda) > COMPTOL_LAMBDA) return 0; for(int l=0; l<=nx+1; ++l) if(fabs(n(l)-wg.n(l)) > COMPTOL_N) return 0; for(int l=0; l<=nx; ++l) if(fabs(hx(l)-wg.hx(l)) > COMPTOL_HX) return 0; if(special != wg.special) return 0; return 1; } /* test for constant refractive index, return 1, if constant, 0 otherwise */ int Waveguide::constmed() const { double tn = n(0); for(int l=1; l<=nx+1; ++l) if(fabs(tn - n(l))> COMPTOL_N) return 0; return 1; } /* remove unnecessary boundaries from waveguide, down to a minimum of minnx inner layers */ void Waveguide::optimize(int minnx) { int done = 0; int rl = 0; while(nx >= (minnx+1) && done != 1) { done = 1; for(int l=0; l<=nx; ++l) { if(done == 1 && fabs(n(l)-n(l+1)) < COMPTOL_N) { done = 0; rl = l; } } if(done != 1) { for(int l=rl; l<=nx; ++l) n(l) = n(l+1); for(int l=rl; l<=nx-1; ++l) hx(l) = hx(l+1); --nx; } } Dvector n_n = n.subvector(nx+2, 0); Dvector n_hx = hx.subvector(nx+1, 0); n = n_n; hx = n_hx; return; } /* ... down to a minimum of one inner layer */ void Waveguide::optimize() { optimize(1); return; } /* check consistency of the waveguide object: error&exit, if not ok */ void Waveguide::consistency() const { if(nx <= 0 || nx >= 1000) structerror("Waveguide: consistency: nx"); if(hx.nel != nx+1) structerror("Waveguide: consistency, hx: #elements"); for(int b=0; b<=nx-1; ++b) { if(hx(b) > hx(b+1) - HDIST) structerror("Waveguide: consistency, hx: order/distance"); } if(lambda < 0.001 && lambda > 100.0) structerror("Waveguide: consistency, lambda: magnitude"); if(n.nel != nx+2) structerror("Waveguide: consistency, n: #elements"); if(special) { for(int l=0; l<=nx+1; ++l) { if( eps(l)<-MAXREFIND*MAXREFIND || eps(l)> MAXREFIND*MAXREFIND) structerror("Waveguide: consistency(special), eps: magnitude"); } } else { for(int l=0; l<=nx+1; ++l) { if(n(l)MAXREFIND) structerror("Waveguide: consistency, n: magnitude"); } } return; } /* check the symmetry of the waveguide, return value l l == 0: no symmetry l >= 1: symmetric structure with respect to layer l */ int Waveguide::checksymmetry() const { if(nx % 2 == 0) return 0; double c = (hx(nx)+hx(0))/2.0; int cl = layeridx(c); for(int l=0; l<=nx; ++l) { if(fabs((c-hx(l)) - (hx(nx-l)-c))> COMPTOL_HX) return 0; } for(int l=0; l<=nx+1; ++l) { if(fabs(n(l)-n(nx-l+1)) > COMPTOL_N) return 0; } return cl; } /* split the present waveguide into two, at layer l the pieces are stored in wg1 and wg2 */ void Waveguide::split(int sl, Waveguide& wg1, Waveguide& wg2) const { if(sl == 1) { wg1 = Waveguide(1); wg1.hx(0) = hx(0)-LIMBDSEP/2.0; wg1.hx(1) = hx(0); wg1.n(0) = n(0); wg1.n(1) = n(0); wg1.n(2) = n(1); wg1.lambda = lambda; wg1.special = special; } else { int nnx = sl-1; wg1 = Waveguide(nnx); for(int l=0; l<=nnx; ++l) wg1.hx(l) = hx(l); for(int l=0; l<=nnx+1; ++l) wg1.n(l) = n(l); wg1.lambda = lambda; wg1.special = special; } if(sl == nx) { wg2 = Waveguide(1); wg2.hx(0) = hx(nx); wg2.hx(1) = hx(nx)+LIMBDSEP/2.0; wg2.n(0) = n(nx); wg2.n(1) = n(nx+1); wg2.n(2) = n(nx+1); wg2.lambda = lambda; wg2.special = special; } else { int nnx = nx-sl; wg2 = Waveguide(nnx); for(int l=0; l<=nnx; ++l) wg2.hx(l) = hx(sl+l); for(int l=0; l<=nnx+1; ++l) wg2.n(l) = n(sl+l); wg2.lambda = lambda; wg2.special = special; } return; } /* expand the waveguide towards a structure that is symmetric with respect to xs; xs > hx(nx) */ Waveguide Waveguide::expand(double xs) { if(xs < hx(nx)+HDIST) structerror("Waveguide: expand: xs"); int nnx = 2*nx+1; Waveguide wg(nnx); for(int b=0; b<=nx; ++b) { wg.hx(b) = hx(b); wg.hx(nnx-b) = xs+(xs-hx(b)); } for(int l=0; l<=nx; ++l) { wg.n(l) = n(l); wg.n(nnx+1-l) = n(l); } wg.n(nx+1) = n(nx+1); wg.lambda = lambda; wg.special = special; return wg; } /* mirror the waveguide at xs; xs > hx(nx) */ Waveguide Waveguide::mirror(double xs) { if(xs < hx(nx)+HDIST) structerror("Waveguide: mirror: xs"); Waveguide wg(nx); for(int b=0; b<=nx; ++b) { wg.hx(nx-b) = xs+(xs-hx(b)); } for(int l=0; l<=nx+1; ++l) { wg.n(nx+1-l) = n(l); } wg.lambda = lambda; wg.special = special; return wg; } /* determine a limit for the effective permittivity, beyond which the inner layer l optically decouples the neighbouring structures */ double Waveguide::decoupepseff(int l) { if(l<=0 || l>=nx+1) structerror("Waveguide: decoupepseff: l"); double d = hx(l)-hx(l-1); double lnD = log(DECOUPLED); lnD = lnD/2/PI; return eps(l)+lnD*lnD*lambda*lambda/d/d; } /* determine a limit for the effective permittivity, beyond which the lower boundary at position p becomes irrelevant */ double Waveguide::lbdecepseff(double p) { if(p >= hx(0)-HDIST) structerror("Waveguide: lbdecepseff: p"); double d = hx(0)-p; d *= 2.0; double lnD = log(DECOUPLED); lnD = lnD/2/PI; return eps(0)+lnD*lnD*lambda*lambda/d/d; } /* determine a limit for the effective permittivity, beyond which the upper boundary at position p becomes irrelevant */ double Waveguide::ubdecepseff(double p) { if(p <= hx(nx)+HDIST) structerror("Waveguide: ubdecepseff: p"); double d = p-hx(nx); d *= 2.0; double lnD = log(DECOUPLED); lnD = lnD/2/PI; return eps(nx+1)+lnD*lnD*lambda*lambda/d/d; } /* refractive index profile -> MATLAB .m file disp: output region on the x axis ext0, ext1: filename id characters */ void Waveguide::plot(Interval disp, char ext0, char ext1) const { FILE *dat; char name[13] = "n__.m"; name[1] = ext0; name[2] = ext1; fprintf(stderr, "n [%g, %g] >> %s\n", disp.x0, disp.x1, name); if(special) fprintf(stderr, " (special: eps<->n !)\n"); dat = fopen(name, "w+"); mlout_title(dat, name, "Waveguide profile"); name[3] = 0; double minn, maxn; minn = n.min(); maxn = n.max(); if(special) { minn = sgn_sqrt(minn); maxn = sgn_sqrt(maxn); } mlout_geo(dat, (*this), minn, maxn); mlout_refindplot(name, dat, disp); mlout_print(dat, name, 'e'); fclose(dat); return; } void Waveguide::plot(char ext0, char ext1) const { Interval disp; double w = (hx(nx)-hx(0))/((double) nx); disp.x0 = hx(0)-w; disp.x1 = hx(nx)+w; plot(disp, ext0, ext1); return; } /* ------------------------------------------------------------------- */ /* segmented planar waveguide: geometry + permittivity profile; vacuum wavelength ... */ /* initialize */ SegWgStruct::SegWgStruct() { nz = -1; wgvec = new Waveguide[nz+2]; Waveguide wg(1); wg.n(0) = 1.0; wg.n(1) = 1.0; wg.n(2) = 1.0; wg.hx(0) = -1.0; wg.hx(1) = 1.0; wg.lambda = 1.0; wg.special = 0; for(int i=0; i<=nz+1; ++i) wgvec[i] = wg; hz = Dvector(nz+1); for(int i=0; i<=nz; ++i) hz(i) = ((double) i); } SegWgStruct::SegWgStruct(int vnz) { if(vnz<-1) structerror("SegWgStruct: init < -1"); nz = vnz; wgvec = new Waveguide[nz+2]; Waveguide wg(1); wg.n(0) = 1.0; wg.n(1) = 1.0; wg.n(2) = 1.0; wg.hx(0) = -1.0; wg.hx(1) = 1.0; wg.lambda = 1.0; wg.special = 0; for(int i=0; i<=nz+1; ++i) wgvec[i] = wg; hz = Dvector(nz+1); for(int i=0; i<=nz; ++i) hz(i) = ((double) i); } /* convert a Circuit to a SegWgStruct opt = 1: merge equal segments & layers (default) */ SegWgStruct::SegWgStruct(Circuit rc, int opt) { nz = rc.nz; hz = rc.hz; wgvec = new Waveguide[nz+2]; for(int s=0; s<=nz+1; ++s) { Waveguide wg(rc.nx); wg.hx = rc.hx; wg.n = rc.n.col(s); wg.lambda = rc.lambda; wg.special = rc.special; if(opt == 1) wg.optimize(); wgvec[s] = wg; } if(opt == 1) optimize(); } SegWgStruct::SegWgStruct(Circuit rc) { wgvec = NULL; (*this) = SegWgStruct(rc, 1); } /* destroy */ SegWgStruct::~SegWgStruct() { if(wgvec != NULL) delete[] wgvec; wgvec = NULL; nz = 0; hz.strip(); } /* copy constructor */ SegWgStruct::SegWgStruct(const SegWgStruct& swg) { nz = swg.nz; hz = swg.hz; wgvec = new Waveguide[nz+2]; for(int i=0; i<=nz+1; ++i) wgvec[i] = swg.wgvec[i]; } /* assignment */ SegWgStruct& SegWgStruct::operator=(const SegWgStruct& swg) { if(this != &swg) { if(wgvec != NULL) delete[] wgvec; wgvec = NULL; nz = swg.nz; hz = swg.hz; wgvec = new Waveguide[nz+2]; for(int i=0; i<=nz+1; ++i) wgvec[i] = swg.wgvec[i]; } return *this; } /* delete all segment entries */ void SegWgStruct::clear() { nz = -1; hz.strip(); if(wgvec != NULL) delete[] wgvec; wgvec = NULL; } /* set all segments to wg */ void SegWgStruct::init(Waveguide wg) { for(int i=0; i<=nz+1; ++i) wgvec[i] = wg; } /* subscripting Waveguide& SegWgStruct::operator() (int s) { if(s < 0) structerror("SegWgStruct: access ( < 0)"); if(s > nz+1) structerror("SegWgStruct: access ( > nz+1)"); return wgvec[s]; } Waveguide SegWgStruct::operator() (int s) const { if(s < 0) structerror("SegWgStruct: access ( < 0)"); if(s > nz+1) structerror("SegWgStruct: access ( > nz+1)"); return wgvec[s]; } */ /* input from FILE dat */ void SegWgStruct::read(FILE *dat) { if(dat == stderr || dat == stdout) { structerror("SegWgStruct: read"); } else { if(wgvec != NULL) delete[] wgvec; wgvec = NULL; nz = fgetint(dat); if(nz<-1) structerror("SegWgStruct: read: nz < -1"); hz.read(dat); wgvec = new Waveguide[nz+2]; for(int i=0; i<=nz+1; ++i) (wgvec[i]).read(dat); } return; } /* output to FILE dat */ void SegWgStruct::write(FILE *dat) const { comment("SegWgStruct", dat); comment("nz", dat); fputint(nz, dat); comment("hz", dat); hz.write(dat); for(int i=0; i<=nz+1; ++i) (wgvec[i]).write(dat); return; } /* vacuum wavelength: should be unique in all segments */ double SegWgStruct::lambda() { return wgvec[0].lambda; } /* check consistency of the SegWgStruct object: error&exit, if not ok */ void SegWgStruct::consistency() const { double l; l = wgvec[0].lambda; for(int i=1; i<=nz+1; ++i) { if(wgvec[i].lambda !=l) structerror("SegWgStruct: lambda not unique"); } int sp; sp = wgvec[0].special; for(int i=1; i<=nz+1; ++i) { if(wgvec[i].special !=sp) structerror("SegWgStruct: wg.special not unique"); } for(int i=0; i<=nz+1; ++i) wgvec[i].consistency(); return; } /* get segment index corresponding to position z */ int SegWgStruct::segidx(double z) const { int s=0; while(s<=nz && z>hz(s)) ++s; return s; } /* get segment boundaries corresponding to index s */ Interval SegWgStruct::segment(int s) const { Interval i; // if(s < 0) structerror("SegWgStruct: segment( < 0)"); // if(s > nz+1) structerror("SegWgStruct: segment( > nz+1)"); if(s <= 0) i.x0 = -AWAY; else i.x0 = hz(s-1); if(s >= nz+1) i.x1 = AWAY; else i.x1 = hz(s); return i; } /* get segment boundaries corresponding to position z */ Interval SegWgStruct::segment(double z) const { return segment(segidx(z)); } /* get refractive index at position x, z */ double SegWgStruct::n(double x, double z) const { int s = segidx(z); Waveguide wg = wgvec[s]; double nv = wg.n(wg.layeridx(x)); if(wg.special) return sgn_sqrt(nv); return nv; } /* get permittivity at position x, z */ double SegWgStruct::eps(double x, double z) const { int s = segidx(z); Waveguide wg = wgvec[s]; double nv = wg.n(wg.layeridx(x)); if(wg.special) return nv; return nv*nv; } /* translate: hz -> hz+dz */ void SegWgStruct::ztranslate(double dz) { for(int s=0; s<=nz; ++s) hz(s) += dz; return; } /* translate: hx -> hx+dx for all segments */ void SegWgStruct::xtranslate(double dx) { for(int s=0; s<=nz+1; ++s) (wgvec[s]).translate(dx); return; } /* combine equal adjacent segments */ void SegWgStruct::optimize() { int done = 0; int rs = 0; while(nz >= 2 && done != 1) { done = 1; for(int s=0; s<=nz; ++s) { if(done == 1 && (wgvec[s]).equal((wgvec[s+1])) != 0) { done = 0; rs = s; } } if(done != 1) { for(int s=rs; s<=nz; ++s) wgvec[s] = wgvec[s+1]; for(int s=rs; s<=nz-1; ++s) hz(s) = hz(s+1); --nz; } } Waveguide *n_wgvec; n_wgvec = new Waveguide[nz+2]; for(int s=0; s<=nz+1; ++s) n_wgvec[s] = wgvec[s]; delete[] wgvec; wgvec = n_wgvec; Dvector n_hz = hz.subvector(nz+1, 0); hz = n_hz; return; } /* QUEP: construct the corresponding vertical segmented waveguide sequence */ SegWgStruct SegWgStruct::rotate() const { consistency(); int spl = wgvec[0].special; int b=0; for(int s=0; s<=nz+1; ++s) b += wgvec[s].nx+1; Dvector tbd(b); Dvector bd; b=0; for(int s=0; s<=nz+1; ++s) { tbd.setsubvector(wgvec[s].hx, b); b += wgvec[s].nx+1; } tbd.sort(); bd = tbd; if(bd.nel<2) structerror("SegWgStruct: rotated: (1) #x_bounds < 2."); tbd(0) = bd(0); b=1; for(int i=1; i<=bd.nel-1; ++i) { if(fabs(tbd(b-1)-bd(i)) > HDIST) { tbd(b) = bd(i); ++b; } } bd = tbd.subvector(b, 0); if(bd.nel<2) structerror("SegWgStruct: rotated: (2) #x_bounds < 2."); int nnz = bd.nel-1; SegWgStruct nst(nnz); nst.hz = bd; Waveguide twg(nz); twg.lambda = wgvec[0].lambda; twg.hx = hz; twg.n.init(1.0); nst.init(twg); Interval itv; double x, z, r; for(int s=0; s<=nst.nz+1; ++s) { itv = nst.segment(s); x = 0.5*(itv.x0+itv.x1); for(int l=0; l<=twg.nx+1; ++l) { itv = twg.layer(l); z = 0.5*(itv.x0+itv.x1); int idx = wgvec[segidx(z)].layeridx(x); r = wgvec[segidx(z)].n(idx); nst(s).n(l) = r; } } for(int s=0; s<=nst.nz+1; ++s) nst(s).optimize(); nst.optimize(); for(int s=0; s<=nst.nz+1; ++s) nst(s).special = spl; return nst; } /* QUEP: enlarge the sequence by identical input and output segments */ SegWgStruct SegWgStruct::expand(Interval w) const { int nnz = nz+2; SegWgStruct nst(nnz); if(w.x0 >= hz(0)) structerror("SegWgStruct: expand: illegal window (1)"); if(w.x1 <= hz(nz)) structerror("SegWgStruct: expand: illegal window (2)"); nst.hz(0) = w.x0; for(int b=0; b<=nz; ++b) nst.hz(b+1) = hz(b); nst.hz(nnz) = w.x1; nst(0) = wgvec[0]; for(int b=0; b<=nz+1; ++b) nst(b+1) = wgvec[b]; nst(nnz+1) = wgvec[nz+1]; return nst; } /* QUEP: remove the additional outer segments */ SegWgStruct SegWgStruct::shrink() const { if(nz <= 2) structerror("SegWgStruct: shrink: nz<=2"); int nnz = nz-2; SegWgStruct nst(nnz); for(int b=0; b<=nnz; ++b) nst.hz(b) = hz(b+1); for(int b=0; b<=nnz+1; ++b) nst(b) = wgvec[b+1]; return nst; } /* determine a default window around the "inner" structure, e.g. for purposes of plotting */ void SegWgStruct::defaultwindow(Interval& wx, Interval& wz) const { wx = wz = Interval(-1.0, 1.0); if(nz < -1) return; Waveguide wg = wgvec[0]; double x0 = wg.hx(0); double x1 = wg.hx(wg.nx); double wl = wg.lambda; double emin = wg.defaultepseffmin(); double xt; for(int j=1; j<=nz+1; ++j) { wg = wgvec[j]; xt = wg.hx(0); if(xt < x0) x0 = xt; xt = wg.hx(wg.nx); if(xt > x1) x1 = xt; if(wg.eps(0) < emin) emin = wg.eps(0); if(wg.eps(wg.nx+1) < emin) emin = wg.eps(wg.nx+1); } if(emin < 1.0) emin=1.0; double d = 1.5*wl/sqrt(emin); x0 = x0-d; x1 = x1+d; double z0 = -d; double z1 = d; if(nz >= 0) { z0 = hz(0)-d; z1 = hz(nz)+d; } wx = Interval(x0, x1); wz = Interval(z0, z1); return; } /* convert this to a Circuit object */ Circuit SegWgStruct::circuit() const { SegWgStruct h; if(nz >= 1) h = (*this); else { h = SegWgStruct(1); if(nz == -1) { Waveguide wg = wgvec[0]; h(0) = wg; h(1) = wg; h(2) = wg; double d; d = wg.hx(wg.nx) - wg.hx(0); h.hz(0) = -d/2.0; h.hz(1) = d/2.0; } else { if(nz == 0) { Waveguide wg = wgvec[0]; h(0) = wg; h(1) = wg; h(2) = wgvec[1]; double d; d = wg.hx(wg.nx) - wg.hx(0); h.hz(0) = hz(0)-d; h.hz(1) = hz(0); } else structerror("SegWgStruct: circuit: nz"); } } SegWgStruct v; v = h.rotate(); Circuit circ(v.nz, h.nz); circ.lambda = h.wgvec[0].lambda; circ.special = h.wgvec[0].special; circ.hx = v.hz; circ.hz = h.hz; Waveguide wg; double x; for(int s=0; s<=h.nz+1; ++s) { wg = h.wgvec[s]; for(int l=0; l<=v.nz+1; ++l) { if(l>=1 && l<=v.nz) x = 0.5*(v.hz(l-1)+v.hz(l)); else { if(l==0) x = v.hz(0)-1.0; else x = v.hz(v.nz)+1.0; } circ.n(l, s) = wg.n(wg.layeridx(x)); } } return circ; } /* modelling of periodically corrugated waveguides: extend this periodically over the display window zbeg, zend */ SegWgStruct SegWgStruct::unwrap(double zbeg, double zend) const { if(zbeg >= zend) structerror("SegWgStruct: unwrap: zbeg, zend"); double L = hz(nz)-hz(0); double z0 = hz(0); while(zbeg <= z0) z0 -= L; double zn = hz(nz); while(zend >= zn) zn += L; int nnz = (int) floor((zn-z0)/L+0.5)*nz; SegWgStruct st(nnz); int j=0; st(j) = (*this)(0); st.hz(j) = z0; ++j; int k = 1; while(j<=nnz) { st(j) = (*this)(k); st.hz(j) = st.hz(j-1)+(hz(k)-hz(k-1)); ++j; ++k; if(k > nz) k = 1; } st(nnz+1) = (*this)(nz+1); st.optimize(); return st; } /* modelling of periodically corrugated waveguides: remove outermost equal slices */ SegWgStruct SegWgStruct::reduce() const { if(nz < 3) structerror("SegWgStruct: reduce: nz<3"); if((*this)(0).equal((*this)(1)) != 1) structerror("SegWgStruct: reduce: (0)!=(1)"); if((*this)(nz+1).equal((*this)(nz)) != 1) structerror("SegWgStruct: reduce: (nz)!=(nz+1)"); SegWgStruct st(nz-2); for(int l=1; l<=nz; ++l) st(l-1) = (*this)(l); for(int l=1; l<=nz-1; ++l) st.hz(l-1) = (*this).hz(l); return st; } /* refractive index profile -> MATLAB .m file ext0, ext1: filename id characters xbeg, xend, zbeg, zend: display window in the x-z-plane */ void SegWgStruct::plot(double xbeg, double xend, double zbeg, double zend, char ext0, char ext1) const { FILE *dat; char name[13] = "n__.m"; name[1] = ext0; name[2] = ext1; fprintf(stderr, "n [%g, %g] x [%g, %g] >> %s\n", xbeg, xend, zbeg, zend, name); if(special()) fprintf(stderr, " (special: eps<->n !)\n"); dat = fopen(name, "w+"); mlout_title(dat, name, "2D refractive index profile"); name[3] = 0; mlout_gengeoxz(dat, (*this), xbeg, xend, zbeg, zend); mlout_windowxz(dat, xbeg, xend, zbeg, zend); mlout_gengeorefindpatches(dat, (*this)); mlout_refindoverview(name, dat); mlout_print(dat, name, 'e'); fclose(dat); return; } void SegWgStruct::plot(char ext0, char ext1) const { Interval wx, wz; defaultwindow(wx, wz); plot(wx.x0, wx.x1, wz.x0, wz.x1, ext0, ext1); return; } /* special() != 0: n stores permittivity values in place of refractive indices, consistency checks for ref.ind. and perm. are disabled; a misuse of the Circuit container for problems with negative permittivity */ int SegWgStruct::special() const { int sp = wgvec[0].special; for(int i=1; i<=nz+1; ++i) { if(wgvec[i].special !=sp) structerror("SegWgStruct: special() not unique"); } return sp; } /* ------------------------------------------------------------------- */ /* rectangular 2D optical circuit, geometry & permittivity profile, vacuum wavelength; (so far) meant only as a help to define a SegWgStruct */ /* initialize */ Circuit::Circuit() { (*this) = Circuit(1, 1); } Circuit::Circuit(int vnx, int vnz) { if(vnx < 1) structerror("Circuit: vnx < 1"); if(vnz < 1) structerror("Circuit: vnz < 1"); nx = vnx; hx = Dvector(nx+1); nz = vnz; hz = Dvector(nz+1); n = Dmatrix(nx+2, nz+2); special = 0; } /* refractive index profile -> MATLAB .m file ext0, ext1: filename id characters xbeg, xend, zbeg, zend: display window in the x-z-plane */ void Circuit::plot(double xbeg, double xend, double zbeg, double zend, char ext0, char ext1) const { SegWgStruct s(*this); s.plot(xbeg, xend, zbeg, zend, ext0, ext1); return; } void Circuit::plot(char ext0, char ext1) const { SegWgStruct s(*this); s.plot(ext0, ext1); return; } /* set Rectangle index (l, m) corresponding to position (x,z) */ void Circuit::rectidx(double x, double z, int& l, int& m) const { l = 0; while(l<=nx && x>hx(l)) ++l; m = 0; while(m<=nz && z>hz(m)) ++m; return; } /* ------------------------------------------------------------------- */ /* for the VEIMS solver: cross section (2D) of a 3D rectangular waveguide, defined as a sequence of (1D) Waveguides */ /* initialize */ SegWgCrs::SegWgCrs() { ny = -1; wgvec = new Waveguide[ny+2]; Waveguide wg(1); wg.n(0) = 1.0; wg.n(1) = 1.0; wg.n(2) = 1.0; wg.hx(0) = -1.0; wg.hx(1) = 1.0; wg.lambda = 1.0; wg.special = 0; for(int i=0; i<=ny+1; ++i) wgvec[i] = wg; hy = Dvector(ny+1); for(int i=0; i<=ny; ++i) hy(i) = ((double) i); } SegWgCrs::SegWgCrs(int vny) { if(vny<-1) structerror("SeWgCrs: init < -1"); ny = vny; wgvec = new Waveguide[ny+2]; Waveguide wg(1); wg.n(0) = 1.0; wg.n(1) = 1.0; wg.n(2) = 1.0; wg.hx(0) = -1.0; wg.hx(1) = 1.0; wg.lambda = 1.0; wg.special = 0; for(int i=0; i<=ny+1; ++i) wgvec[i] = wg; hy = Dvector(ny+1); for(int i=0; i<=ny; ++i) hy(i) = ((double) i); } SegWgCrs::SegWgCrs(WgCrs rc) { ny = rc.ny; hy = rc.hy; wgvec = new Waveguide[ny+2]; for(int s=0; s<=ny+1; ++s) { Waveguide wg(rc.nx); wg.hx = rc.hx; wg.n = rc.n.col(s); wg.lambda = rc.lambda; wg.special = 0; wg.optimize(); wgvec[s] = wg; } } /* destroy */ SegWgCrs::~SegWgCrs() { if(wgvec != NULL) delete[] wgvec; wgvec = NULL; ny = 0; hy.strip(); } /* copy constructor */ SegWgCrs::SegWgCrs(const SegWgCrs& swg) { ny = swg.ny; hy = swg.hy; wgvec = new Waveguide[ny+2]; for(int i=0; i<=ny+1; ++i) wgvec[i] = swg.wgvec[i]; } /* assignment */ SegWgCrs& SegWgCrs::operator=(const SegWgCrs& swg) { if(this != &swg) { if(wgvec != NULL) delete[] wgvec; wgvec = NULL; ny = swg.ny; hy = swg.hy; wgvec = new Waveguide[ny+2]; for(int i=0; i<=ny+1; ++i) wgvec[i] = swg.wgvec[i]; } return *this; } /* delete all segment entries */ void SegWgCrs::clear() { ny = -1; hy.strip(); if(wgvec != NULL) delete[] wgvec; wgvec = NULL; } /* set all segments to wg */ void SegWgCrs::init(Waveguide wg) { for(int i=0; i<=ny+1; ++i) wgvec[i] = wg; } /* subscripting Waveguide& SegWgCrs::operator() (int s) { if(s < 0) structerror("SegWgCrs: access ( < 0)"); if(s > ny+1) structerror("SegWgCrs: access ( > ny+1)"); return wgvec[s]; } Waveguide SegWgCrs::operator() (int s) const { if(s < 0) structerror("SegWgCrs: access ( < 0)"); if(s > ny+1) structerror("SegWgCrs: access ( > ny+1)"); return wgvec[s]; } */ /* input from FILE dat */ void SegWgCrs::read(FILE *dat) { if(dat == stderr || dat == stdout) { structerror("SegWgCrs: read"); } else { if(wgvec != NULL) delete[] wgvec; wgvec = NULL; ny = fgetint(dat); if(ny<-1) structerror("read: ny < -1"); hy.read(dat); wgvec = new Waveguide[ny+2]; for(int i=0; i<=ny+1; ++i) (wgvec[i]).read(dat); } return; } /* output to FILE dat */ void SegWgCrs::write(FILE *dat) const { comment("SegWgCrs", dat); comment("ny", dat); fputint(ny, dat); comment("hy", dat); hy.write(dat); for(int i=0; i<=ny+1; ++i) (wgvec[i]).write(dat); return; } /* vacuum wavelength: should be unique in all segments */ double SegWgCrs::lambda() { return wgvec[0].lambda; } /* check consistency of the SegWgCrs object: error&exit, if not ok */ void SegWgCrs::consistency() const { double l; l = wgvec[0].lambda; for(int i=1; i<=ny+1; ++i) { if(wgvec[i].lambda !=l) structerror("SegWgCrs: lambda not unique"); } int sp; sp = wgvec[0].special; for(int i=1; i<=ny+1; ++i) { if(wgvec[i].special !=sp) structerror("SegWgCrs: wg.special not unique"); } for(int i=0; i<=ny+1; ++i) wgvec[i].consistency(); return; } /* get segment index corresponding to position y */ int SegWgCrs::segidx(double y) const { int s=0; while(s<=ny && y>hy(s)) ++s; return s; } /* get segment boundaries corresponding to index s */ Interval SegWgCrs::segment(int s) const { Interval i; // if(s < 0) structerror("SegWgCrs: segment( < 0)"); // if(s > ny+1) structerror("SegWgCrs: segment( > ny+1)"); if(s <= 0) i.x0 = -AWAY; else i.x0 = hy(s-1); if(s >= ny+1) i.x1 = AWAY; else i.x1 = hy(s); return i; } /* get segment boundaries corresponding to position y */ Interval SegWgCrs::segment(double y) const { return segment(segidx(y)); } /* get refractive index at position x, y */ double SegWgCrs::n(double x, double y) const { int s = segidx(y); Waveguide wg = wgvec[s]; double nv = wg.n(wg.layeridx(x)); if(wg.special) return sgn_sqrt(nv); return nv; } /* get permittivity at position x, y */ double SegWgCrs::eps(double x, double y) const { int s = segidx(y); Waveguide wg = wgvec[s]; double nv = wg.n(wg.layeridx(x)); if(wg.special) return nv; return nv*nv; } /* translate: hy -> hy+dy */ void SegWgCrs::ytranslate(double dy) { for(int s=0; s<=ny; ++s) hy(s) += dy; return; } /* translate: hx -> hx+dx for all segments */ void SegWgCrs::xtranslate(double dx) { for(int s=0; s<=ny+1; ++s) (wgvec[s]).translate(dx); return; } /* construct the corresponding vertical segmented waveguide sequence */ SegWgCrs SegWgCrs::rotate() const { consistency(); int b=0; for(int s=0; s<=ny+1; ++s) b += wgvec[s].nx+1; Dvector tbd(b); Dvector bd; b=0; for(int s=0; s<=ny+1; ++s) { tbd.setsubvector(wgvec[s].hx, b); b += wgvec[s].nx+1; } tbd.sort(); bd = tbd; if(bd.nel<2) structerror("SegWgCrs: rotated: (1) #x_bounds < 2."); tbd(0) = bd(0); b=1; for(int i=1; i<=bd.nel-1; ++i) { if(fabs(tbd(b-1)-bd(i)) > HDIST) { tbd(b) = bd(i); ++b; } } bd = tbd.subvector(b, 0); if(bd.nel<2) structerror("SegWgCrs: rotated: (2) #x_bounds < 2."); int nny = bd.nel-1; SegWgCrs nst(nny); nst.hy = bd; Waveguide twg(ny); twg.lambda = wgvec[0].lambda; twg.hx = hy; twg.n.init(1.0); nst.init(twg); Interval itv; double x, y, r; for(int s=0; s<=nst.ny+1; ++s) { itv = nst.segment(s); x = 0.5*(itv.x0+itv.x1); for(int l=0; l<=twg.nx+1; ++l) { itv = twg.layer(l); y = 0.5*(itv.x0+itv.x1); int idx = wgvec[segidx(y)].layeridx(x); r = wgvec[segidx(y)].n(idx); nst(s).n(l) = r; } } for(int s=0; s<=nst.ny+1; ++s) nst(s).optimize(); return nst; } /* determine a default window around the "inner" structure, e.g. for purposes of plotting */ void SegWgCrs::defaultwindow(Interval& wx, Interval& wy) const { wx = wy = Interval(-1.0, 1.0); if(ny < -1) return; Waveguide wg = wgvec[0]; double x0 = wg.hx(0); double x1 = wg.hx(wg.nx); double wl = wg.lambda; double emin = wg.defaultepseffmin(); double xt; for(int j=1; j<=ny+1; ++j) { wg = wgvec[j]; xt = wg.hx(0); if(xt < x0) x0 = xt; xt = wg.hx(wg.nx); if(xt > x1) x1 = xt; if(wg.eps(0) < emin) emin = wg.eps(0); if(wg.eps(wg.nx+1) < emin) emin = wg.eps(wg.nx+1); } if(emin < 1.0) emin=1.0; double d = 1.5*wl/sqrt(emin); x0 = x0-d; x1 = x1+d; double y0 = -d; double y1 = d; if(ny >= 0) { y0 = hy(0)-d; y1 = hy(ny)+d; } wx = Interval(x0, x1); wy = Interval(y0, y1); return; } /* convert this to a WgCrs object */ WgCrs SegWgCrs::wgcrs() const { SegWgCrs h; if(ny >= 1) h = (*this); else { h = SegWgCrs(1); if(ny == -1) { Waveguide wg = wgvec[0]; h(0) = wg; h(1) = wg; h(2) = wg; double d; d = wg.hx(wg.nx) - wg.hx(0); h.hy(0) = -d/2.0; h.hy(1) = d/2.0; } else { if(ny == 0) { Waveguide wg = wgvec[0]; h(0) = wg; h(1) = wg; h(2) = wgvec[1]; double d; d = wg.hx(wg.nx) - wg.hx(0); h.hy(0) = hy(0)-d; h.hy(1) = hy(0); } else structerror("SegWgCrs: circuit: ny"); } } SegWgCrs v; v = h.rotate(); WgCrs circ(v.ny, h.ny); circ.lambda = h.wgvec[0].lambda; circ.hx = v.hy; circ.hy = h.hy; Waveguide wg; double x; for(int s=0; s<=h.ny+1; ++s) { wg = h.wgvec[s]; for(int l=0; l<=v.ny+1; ++l) { if(l>=1 && l<=v.ny) x = 0.5*(v.hy(l-1)+v.hy(l)); else { if(l==0) x = v.hy(0)-1.0; else x = v.hy(v.ny)+1.0; } circ.n(l, s) = wg.n(wg.layeridx(x)); } } return circ; } /* refractive index profile -> MATLAB .m file ext0, ext1: filename id characters xbeg, xend, ybeg, yend: display window in the x-y-plane */ void SegWgCrs::plot(double xbeg, double xend, double ybeg, double yend, char ext0, char ext1) const { FILE *dat; char name[13] = "n__.m"; name[1] = ext0; name[2] = ext1; fprintf(stderr, "n [%g, %g] x [%g, %g] >> %s\n", xbeg, xend, ybeg, yend, name); dat = fopen(name, "w+"); mlout_title(dat, name, "2D waveguide cross section"); name[3] = 0; mlout_gengeoxy(dat, (*this), Rect(xbeg, ybeg, xend, yend)); mlout_windowxy(dat, xbeg, xend, ybeg, yend); mlout_gengeorefindpatches(dat, (*this)); mlout_wgcrsoverview(name, dat); mlout_print(dat, name, 'e'); fclose(dat); return; } void SegWgCrs::plot(Rect disp, char ext0, char ext1) const { plot(disp.x0, disp.x1, disp.y0, disp.y1, ext0, ext1); return; } void SegWgCrs::plot(char ext0, char ext1) const { Interval wx, wy; defaultwindow(wx, wy); plot(wx.x0, wx.x1, wy.x0, wy.x1, ext0, ext1); return; } /* ------------------------------------------------------------------- */ /* for the VEIMS solver: cross section (2D) of a 3D rectangular waveguide, geometry & permittivity profile, vacuum wavelength; (so far) meant as a help to define a SegWgCrs */ /* initialize */ WgCrs::WgCrs() { (*this) = WgCrs(1, 1); } WgCrs::WgCrs(int vnx, int vny) { if(vnx < 1) structerror("WgCrs: vnx < 1"); if(vny < 1) structerror("WgCrs: vny < 1"); nx = vnx; hx = Dvector(nx+1); ny = vny; hy = Dvector(ny+1); n = Dmatrix(nx+2, ny+2); } WgCrs::WgCrs(int vnx, Dvector vhx, int vny, Dvector vhy, double l, Dmatrix vn) { nx = vnx; ny = vny; hx = vhx; hy = vhy; lambda = l; n = vn; } /* free allocated memory */ void WgCrs::strip() { nx = 0; ny = 0; hx.strip(); hy.strip(); n.strip(); } /* set Rectangle index (l, m) corresponding to position (x,y) */ void WgCrs::rectidx(double x, double y, int& l, int& m) const { l = 0; while(l<=nx && x>hx(l)) ++l; m = 0; while(m<=ny && y>hy(m)) ++m; return; } /* get rectangle boundaries corresponding to index (l,m) */ Rect WgCrs::rectbounds(int l, int m) const { Rect r; if(l <= 0) r.x0 = -AWAY; else r.x0 = hx(l-1); if(l >= nx+1) r.x1 = AWAY; else r.x1 = hx(l); if(m <= 0) r.y0 = -AWAY; else r.y0 = hy(m-1); if(m >= ny+1) r.y1 = AWAY; else r.y1 = hy(m); return r; } /* get rectangle boundaries corresponding to position (x,y) */ Rect WgCrs::rectbounds(double x, double y) const { int l, m; rectidx(x, y, l, m); return rectbounds(l, m); } /* test neighbourhood of two rectangles */ int WgCrs::testconnect(int l0, int m0, int l1, int m1) const { int ld, md; ld = l1-l0; if(ld<0) ld = -ld; md = m1-m0; if(md<0) md = -md; if(ld==1 && md==0) return 1; if(ld==0 && md==1) return 1; if(ld==0 && md==0) return 1; return 0; } /* test matching of boundary positions with another waveguide */ int WgCrs::bdmatch(WgCrs w) const { int j, k, e; for(j=0; j<=w.nx; ++j) { e = 0; for(k=0; k<=nx; ++k) { if(fabs(hx(k)-w.hx(j))<=HDIST) e = 1; } if(e != 1) return 1; } for(j=0; j<=w.ny; ++j) { e = 0; for(k=0; k<=ny; ++k) { if(fabs(hy(k)-w.hy(j))<=HDIST) e = 1; } if(e != 1) return 1; } return 0; } /* permittivity on rectangle l, m */ double WgCrs::eps(int l, int m) const { double e; e = n(l,m); return e*e; } /* permittivity at position x, y */ double WgCrs::eps(double x, double y) const { double e; int l, m; rectidx(x, y, l, m); e = n(l,m); return e*e; } /* translate: hx -> hx+dx, hy -> hy+dy */ void WgCrs::translate(double dx, double dy) { int l, m; for(l=0; l<=nx; ++l) hx(l) += dx; for(m=0; m<=ny; ++m) hy(m) += dy; return; } /* refractive index profile -> MATLAB .m file ext0, ext1: filename id characters xbeg, xend, ybeg, yend: display window in the x-y-plane */ void WgCrs::plot(double xbeg, double xend, double ybeg, double yend, char ext0, char ext1) const { SegWgCrs s(*this); s.plot(xbeg, xend, ybeg, yend, ext0, ext1); return; } void WgCrs::plot(Rect disp, char ext0, char ext1) const { SegWgCrs s(*this); s.plot(disp, ext0, ext1); return; } void WgCrs::plot(char ext0, char ext1) const { SegWgCrs s(*this); s.plot(ext0, ext1); return; } /* Rect: rectangular section from the waveguide cross section */ /* initialize */ Rect::Rect() { x0 = 0.0; y0 = 0.0; x1 = 1.0; y1 = 1.0; } Rect::Rect(double xa, double ya, double xb, double yb) { if(xa <= xb) { x0 = xa; x1 = xb; } else { x0 = xb; x1 = xa; } if(ya <= yb) { y0 = ya; y1 = yb; } else { y0 = yb; y1 = ya; } } /* output to FILE dat */ void Rect::write(FILE *dat) const { if(dat == stderr || dat == stdout) { fprintf(dat, "[%g,%g|%g,%g]", x0, y0, x1, y1); } else { fputdouble(x0, dat); fputdouble(y0, dat); fputdouble(x1, dat); fputdouble(y1, dat); } return; } /* input from FILE dat */ void Rect::read(FILE *dat) { if(dat == stderr || dat == stdout) { structerror("Rect: read"); } else { x0 = fgetdouble(dat); y0 = fgetdouble(dat); x1 = fgetdouble(dat); y1 = fgetdouble(dat); } return; } Rect XYplane(-AWAY, -AWAY, AWAY, AWAY); /* ----------------------------------------------------------------------- */ /* more general specification of a dielectric structure as a stack of possibly overlapping geometric shapes with constant permittivity, individual shapes */ /* initialize: generic */ Overlay::Overlay() { type = HLAY; eps = 1.0; n = 1.0; ubd = 1.0; lbd = -1.0; xo = 0.0; zo = 0.0; } /* initialize: horizontal (HLAY) or vertical layer (VLAY), refr.ind. nc, thickness d, center at x=ofs (HLAY) or z=ofs (VLAY) */ Overlay::Overlay(OvlType t, double nc, double d, double ofs) { if(t != HLAY && t != VLAY) structerror("Overlay: type"); type = t; eps = nc*nc; n = nc; ubd = ofs+0.5*d; lbd = ofs-0.5*d; xo = ofs; zo = ofs; } /* initialize: RING cavity, refr.ind. nc, radius r, core width d, origin at ox, oz */ Overlay::Overlay(OvlType t, double nc, double r, double d, double ox, double oz) { if(t != RING) structerror("Overlay: type"); type = t; eps = nc*nc; n = nc; ubd = r; lbd = r-d; xo = ox; zo = oz; } /* initialize: DISK cavity, refr.ind. nc, radius R, origin at ox, oz */ Overlay::Overlay(OvlType t, double nc, double r, double ox, double oz) { if(t != DISK) structerror("Overlay: type"); type = t; eps = nc*nc; n = nc; ubd = r; xo = ox; zo = oz; } /* initialize: rectangular region RECT, defined by Intervals rx, rz, refr.ind. nc */ Overlay::Overlay(OvlType t, double nc, Interval rx, Interval rz) { if(t != RECT) structerror("Overlay: type"); type = t; eps = nc*nc; n = nc; if(rx.x0 <= rx.x1) { lbd = rx.x0; ubd = rx.x1; } else { lbd = rx.x1; ubd = rx.x0; } if(rz.x0 <= rz.x1) { xo = rz.x0; zo = rz.x1; } else { xo = rz.x1; zo = rz.x0; } } /* decide whether a position falls inside the shape (1) or not (0) */ int Overlay::inside(double x, double z) { double dx, dz, r; switch(type) { case HLAY: if(x > ubd) return 0; if(x < lbd) return 0; return 1; break; case VLAY: if(z > ubd) return 0; if(z < lbd) return 0; return 1; break; case RING: dx = x-xo; dz = z-zo; r = sqrt(dx*dx+dz*dz); if(r > ubd) return 0; if(r < lbd) return 0; return 1; case DISK: dx = x-xo; dz = z-zo; r = sqrt(dx*dx+dz*dz); if(r > ubd) return 0; return 1; case RECT: if(x > ubd) return 0; if(x < lbd) return 0; if(z > zo) return 0; if(z < xo) return 0; return 1; break; } return 0; } /* read from FILE dat */ void Overlay::read(FILE *dat) { if(dat == stderr || dat == stdout) { structerror("Overlay: read"); } comment("Overlay", dat); type = (OvlType)fgetint(dat); eps = fgetdouble(dat); n = fgetdouble(dat); xo = fgetdouble(dat); zo = fgetdouble(dat); lbd = fgetdouble(dat); ubd = fgetdouble(dat); return; } /* output to FILE dat */ void Overlay::write(FILE *dat) const { comment("Overlay", dat); comment("type", dat); fputint(((int)type), dat); comment("eps", dat); fputdouble(eps, dat); comment("n", dat); fputdouble(n, dat); comment("xo", dat); fputdouble(xo, dat); comment("zo", dat); fputdouble(zo, dat); comment("lbd", dat); fputdouble(lbd, dat); comment("ubd", dat); fputdouble(ubd, dat); return; } /* the entire stack */ /* initialize */ OvlStruct::OvlStruct() { no = 0; lambda = 1.0; bgeps = 1.0; bgn = 1.0; ovlvec = NULL; } /* initialize */ OvlStruct::OvlStruct(double b, double l) { no = 0; bgeps = b; bgn = sqrt(fabs(b)); lambda = l; ovlvec = NULL; } /* destroy */ OvlStruct::~OvlStruct() { if(ovlvec != NULL) delete[] ovlvec; ovlvec = NULL; no = 0; } /* copy constructor */ OvlStruct::OvlStruct(const OvlStruct& o) { no = o.no; bgeps = o.bgeps; bgn = o.bgn; lambda = o.lambda; if(no >= 1) { ovlvec = new Overlay[no]; for(int i=0; i<=no-1; ++i) ovlvec[i] = o.ovlvec[i]; } } /* assignment */ OvlStruct& OvlStruct::operator=(const OvlStruct& o) { if(this != &o) { if(ovlvec != NULL) delete[] ovlvec; ovlvec = NULL; no = o.no; bgeps = o.bgeps; bgn = o.bgn; lambda = o.lambda; if(no >= 1) { ovlvec = new Overlay[no]; for(int i=0; i<=no-1; ++i) ovlvec[i] = o.ovlvec[i]; } } return *this; } /* delete all entries */ void OvlStruct::clear() { no = 0; if(ovlvec != NULL) delete[] ovlvec; ovlvec = NULL; } /* place an Overlay on the top of the stack */ void OvlStruct::place(Overlay o) { if(no<0) no = 0; Overlay *nov; nov = new Overlay[no+1]; Overlay* np = nov; if(no >= 1) { Overlay* mp = ovlvec; for(int i=0; i<=no-1; ++i) *np++ = *mp++; } *np = o; if(ovlvec != NULL) delete[] ovlvec; ovlvec = nov; ++no; return; } /* input from FILE dat */ void OvlStruct::read(FILE *dat) { if(dat == stderr || dat == stdout) { structerror("OvlStruct: read"); } else { if(ovlvec != NULL) delete[] ovlvec; ovlvec = NULL; bgeps = fgetdouble(dat); bgn = fgetdouble(dat); lambda = fgetdouble(dat); no = fgetint(dat); if(no<0) structerror("OvlStruct: read: no < 0"); if(no >= 1) { ovlvec = new Overlay[no]; for(int i=0; i<=no-1; ++i) (ovlvec[i]).read(dat); } } return; } /* output to FILE dat */ void OvlStruct::write(FILE *dat) const { if(dat==stdout || dat==stderr) { fprintf(dat, "OvlStruct, shapes: %d, n = %g, eps = %g, lambda = %g\n", no, bgn, bgeps, lambda); for(int i=0; i<=no-1; ++i) { fprintf(dat, "[%d] type %d, eps = %g, n = %g\n", i, ovlvec[i].type, ovlvec[i].eps, ovlvec[i].n); fprintf(dat, " xo = %g, zo = %g, lbd = %g, ubd = %g\n", ovlvec[i].xo, ovlvec[i].zo, ovlvec[i].lbd, ovlvec[i].ubd); } return; } comment("OvlStruct", dat); comment("bgeps", dat); fputdouble(bgeps, dat); comment("bgn", dat); fputdouble(bgn, dat); comment("lambda", dat); fputdouble(lambda, dat); comment("no", dat); fputint(no, dat); if(no >= 1) { for(int i=0; i<=no-1; ++i) (ovlvec[i]).write(dat); } return; } /* permittivity at position x, z */ double OvlStruct::eps(double x, double z) const { if(no >= 1) { for(int i=no-1; i>=0; --i) { if((ovlvec[i]).inside(x, z) != 0) return ovlvec[i].eps; } } return bgeps; } /* refractive index at position x, z */ double OvlStruct::n(double x, double z) const { if(no >= 1) { for(int i=no-1; i>=0; --i) { if(ovlvec[i].inside(x, z) != 0) return ovlvec[i].n; } } return bgn; } /* for purely rectangular structures: convert this to a Circuit object */ Circuit OvlStruct::circuit() const { int nx, nz; Dvector hx, hz; Circuit c; if(no == 0) { nx = 1; hx = Dvector(nx+1); hx(0) = -0.5; hx(1) = 0.5; nz = 1; hz = Dvector(nz+1); hz(0) = -0.5; hz(1) = 0.5; c = Circuit(nx, nz); c.hx = hx; c.hz = hz; c.lambda = lambda; c.n = Dmatrix(nx+2, nz+2); c.n.init(bgn); c.special = 0; return c; } nx = 0; nz = 0; hx.strip(); hz.strip(); for(int i=0; i<=no-1; ++i) { switch(ovlvec[i].type) { case HLAY: hx.append(ovlvec[i].ubd); hx.append(ovlvec[i].lbd); break; case VLAY: hz.append(ovlvec[i].ubd); hz.append(ovlvec[i].lbd); break; case RECT: hx.append(ovlvec[i].ubd); hx.append(ovlvec[i].lbd); hz.append(ovlvec[i].zo); hz.append(ovlvec[i].xo); break; default: structerror( "OvlStruct: circuit: non rectangular Overlay"); return c; break; } } Dvector thx = hx; hx.strip(); thx.sort(); for(int j=0; j<=thx.nel-1; ++j) { double x = thx(j); if(fabs(x) < AWAY/10.0) { if(hx.nel <= 0) hx.append(x); else if(fabs(x-hx(hx.nel-1)) > HDIST) hx.append(x); } } if(hx.nel <= 0) { hx = Dvector(2); hx(0) = -0.5; hx(1) = 0.5; } if(hx.nel == 1) { thx = hx; hx = Dvector(2); hx(0) = thx(0); hx(1) = thx(0)+1.0; } Dvector thz = hz; hz.strip(); thz.sort(); for(int j=0; j<=thz.nel-1; ++j) { double z = thz(j); if(fabs(z) < AWAY/10.0) { if(hz.nel <= 0) hz.append(z); else if(fabs(z-hz(hz.nel-1)) > HDIST) hz.append(z); } } if(hz.nel <= 0) { hz = Dvector(2); hz(0) = -0.5; hz(1) = 0.5; } if(hz.nel == 1) { thz = hz; hz = Dvector(2); hz(0) = thz(0); hz(1) = thz(0)+1.0; } nx = hx.nel-1; nz = hz.nel-1; c = Circuit(nx, nz); c.hx = hx; c.hz = hz; for(int i=0; i<=nx+1; ++i) { double x0, x1; if(i == 0 ) x0 = -AWAY; else x0 = hx(i-1); if(i == nx+1) x1 = AWAY; else x1 = hx(i); for(int j=0; j<=nz+1; ++j) { double z0, z1; if(j == 0 ) z0 = -AWAY; else z0 = hz(j-1); if(j == nz+1) z1 = AWAY; else z1 = hz(j); c.n(i, j) = n(0.5*(x0+x1), 0.5*(z0+z1)); } } c.lambda = lambda; c.special = 0; return c; } /* for purely rectangular structures: convert this to a SegWgStruct object */ SegWgStruct OvlStruct::segwgstruct() const { Circuit c = circuit(); SegWgStruct s(c); return s; }