/* * WMM - Wave-matching method for mode analysis of dielectric waveguides * https://wmm.computational-photonics.eu/ */ /* * mode.cpp * general mode definitions */ #include #include #include #include"inout.h" #include"complex.h" #include"matrix.h" #include"waveg.h" #include"maxweq.h" #include"wmmvis.h" #include"mode.h" /* error message */ void modeerror(const char *s) { fprintf(stderr, " ERR: Mode: %s.\n", s); exit(1); } /* store num values of component cp between (x0, y0) and (x1, y1) in a vector cp: EX - HZ, SZ */ Dvector Mode::fieldvec(Fcomp cp, int num, double x0, double y0, double x1, double y1) const { double x, y; double dx, dy; int j; if(num <= 0) modeerror("fieldvec: num <= 0"); Dvector f(num); if(num == 1) { x = 0.5*(x0+x1); y = 0.5*(x0+x1); f(0) = field(cp, x, y); return f; } dx = (x1-x0)/(num-1); dy = (y1-y0)/(num-1); for(j=0; j<=num-1; ++j) { x = x0+j*dx; y = y0+j*dy; f(j) = field(cp, x, y); } return f; } /* evaluate component cp on a rectangular npx x npy mesh on area disp cp: EX - HZ, SZ foa: ORG, MOD, SQR */ Dmatrix Mode::fieldmat(Rect disp, int npx, int npy, Fcomp cp, Afo foa) const { double x, y; double dx, dy; int i, j; if(npx <= 1) modeerror("fieldmat: npx <= 1"); if(npy <= 1) modeerror("fieldmat: npy <= 1"); Dmatrix f(npx, npy); double ft; 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(cp, x, y); switch(foa) { case MOD: ft = fabs(ft); break; case SQR: ft = ft*ft; break; case ORG: break; case REP: if(cp == EZ || cp == HZ) ft = 0.0; break; case IMP: if(cp != EZ && cp != HZ) ft = 0.0; break; } f(i, j) = ft; } } return f; } /* z component of the Poyntingvector, integrated over the entire x-y-plane */ /* TE part */ double Mode::tepower() const { return -0.5*recint(HX, EY, wg.xyplane); } /* TM part */ double Mode::tmpower() const { return 0.5*recint(EX, HY, wg.xyplane); } /* total power */ double Mode::totpower() const { return tepower()+tmpower(); } /* polarization character of a hybrid mode */ /* mod: 'E': (int E_y^2 dxdy)/(int E_x^2+E_y^2 dxdy) */ /* 'H': (int H_x^2 dxdy)/(int H_x^2+H_y^2 dxdy) */ /* 'S': (int S_z^TE dxdy)/(int S_z dxdy) */ double Mode::polcharacter(char mod) const { double te, tm; switch(mod) { case 'E': te = recint(EY, EY, wg.xyplane); tm = recint(EX, EX, wg.xyplane); break; case 'H': te = recint(HX, HX, wg.xyplane); tm = recint(HY, HY, wg.xyplane); break; default: te = tepower(); tm = tmpower(); break; } return te/(te+tm); } /* set maximum values for E, H, Sz (only a coarse approximation !) */ void Mode::setfieldmax() { int l, m; double x, y, w, s; int np = 10; Rect r; minE = 1.0e+300; maxE = -1.0e+300; minH = 1.0e+300; maxH = -1.0e+300; minS = 1.0e+300; maxS = -1.0e+300; for(l=1; l<=wg.nx; ++l) { for(m=1; m<=wg.ny;++m) { r = wg.rectbounds(l, m); s = (r.x1-r.x0)/np; y = (r.y1+r.y0)/2.0; for(x=r.x0; x<=r.x1; x+=s) { w = field(EX, x, y); if(w > maxE) maxE = w; if(w < minE) minE = w; w = field(EY, x, y); if(w > maxE) maxE = w; if(w < minE) minE = w; w = field(EZ, x, y); if(w > maxE) maxE = w; if(w < minE) minE = w; } s = (r.y1-r.y0)/np; x = (r.x1+r.x0)/2.0; for(y=r.y0; y<=r.y1; y+=s) { w = field(EX, x, y); if(w > maxE) maxE = w; if(w < minE) minE = w; w = field(EY, x, y); if(w > maxE) maxE = w; if(w < minE) minE = w; w = field(EZ, x, y); if(w > maxE) maxE = w; if(w < minE) minE = w; } s = (r.x1-r.x0)/np; y = (r.y1+r.y0)/2.0; for(x=r.x0; x<=r.x1; x+=s) { w = field(HX, x, y); if(w > maxH) maxH = w; if(w < minH) minH = w; w = field(HY, x, y); if(w > maxH) maxH = w; if(w < minH) minH = w; w = field(HZ, x, y); if(w > maxH) maxH = w; if(w < minH) minH = w; } s = (r.y1-r.y0)/np; x = (r.x1+r.x0)/2.0; for(y=r.y0; y<=r.y1; y+=s) { w = field(HX, x, y); if(w > maxH) maxH = w; if(w < minH) minH = w; w = field(HY, x, y); if(w > maxH) maxH = w; if(w < minH) minH = w; w = field(HZ, x, y); if(w > maxH) maxH = w; if(w < minH) minH = w; } s = (r.x1-r.x0)/np; y = (r.y1+r.y0)/2.0; for(x=r.x0; x<=r.x1; x+=s) { w = field(SZ, x, y); if(w > maxS) maxS = w; if(w < minS) minS = w; } s = (r.y1-r.y0)/np; x = (r.x1+r.x0)/2.0; for(y=r.y0; y<=r.y1; y+=s) { w = field(SZ, x, y); if(w > maxS) maxS = w; if(w < minS) minS = w; } } } return; } /* coordinates of the maximum mode intensity on rectangle r, average 1/e^2-radius in x- and y-direction only a coarse approximation ! */ double Mode::maxintensity(Rect r, double& xm, double& ym, double& rx, double& ry) const { double x, y, w, sx, sy; int np = 50; double max; double a, b; max = -1.0e+300; xm = 0.0; ym = 0.0; sx = (r.x1-r.x0)/np; sy = (r.y1-r.y0)/np; for(x=r.x0; x<=r.x1; x+=sx) { for(y=r.y0; y<=r.y1; y+=sy) { w = fabs(field(SZ, x, y)); if(w > max) { max = w; xm = x; ym = y; } } } x = xm; y = ym; w = fabs(field(SZ, x, y)); while(w > max*0.1353) { x -= sx; w = fabs(field(SZ, x, y)); } a = x; x = xm; y = ym; w = fabs(field(SZ, x, y)); while(w > max*0.1353) { x += sx; w = fabs(field(SZ, x, y)); } b = x; rx = 0.5*(b-a); x = xm; y = ym; w = fabs(field(SZ, x, y)); while(w > max*0.1353) { y -= sy; w = fabs(field(SZ, x, y)); } a = y; x = xm; y = ym; w = fabs(field(SZ, x, y)); while(w > max*0.1353) { y += sy; w = fabs(field(SZ, x, y)); } b = y; ry = 0.5*(b-a); return max; } /* normalize mode to totpower() == nrm */ void Mode::normalize(double nrm) { ampfac = 1.0; ampfac = sqrt(nrm)/sqrt(totpower()); setfieldmax(); return; } /* write field profile cp: EX - HZ, SZ foa: MOD, ORG, SQR fac: output factor (typically 1.0 or -1.0) disp: output region on the x-y-plane ext0, ext1: filename id characters npx, npy: number of points in output mesh adj: 1: adjust contour levels to cp adj: 0: adjust contour levels to minE, maxEm, minH ... */ void Mode::writeprofile(Fcomp cp, Afo foa, double fac, Rect disp, char ext0, char ext1, int npx, int npy, int adj) const { FILE *dat; int xi, yi, l, m; double x, y, w, s; double minz, maxz; double mi, ma; double maxf; char name[13] = "t_______.hvl"; Dmatrix f(npx+1,npy+1); switch(pol) { case QTE: name[1] = 'e'; break; case QTM: name[1] = 'm'; break; case VEC: name[0] = 'v'; name[1] = 'c'; break; } switch(sym) { case SYM: name[2] = 's'; break; case ASY: name[2] = 'a'; break; case NOS: name[2] = '-'; break; } 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; minz = 1.0e+300; maxz = -1.0e+300; maxf = -1.0e+300; for(xi=0; xi<=npx; ++xi) { x = disp.x0+(disp.x1-disp.x0)/npx*xi; for(yi=0; yi<=npy; ++yi) { y = disp.y0+(disp.y1-disp.y0)/npy*yi; w = field(cp, x, y); if(fabs(w) > maxf) maxf = fabs(w); w = fac*w; if(foa == MOD) w = fabs(w); if(foa == SQR) w = w*w; if(foa == IMP) w = 0.0; f(xi,yi) = w; if(w < minz) minz = w; if(w > maxz) maxz = w; } } fprintf(stderr, " max|%c%c| = %g; *(%g) >> %s\n", fldchr(cp), cpchr(cp), maxf, fac, name); if(adj == 0) { mi = 0.0; ma = 1.0; switch(fldchr(cp)) { case 'E': mi = minE; ma = maxE; break; case 'H': mi = minH; ma = maxH; break; case 'S': mi = minS; ma = maxS; break; } switch(foa) { case MOD: ma = fabs(ma); mi = fabs(ma); minz = 0.0; if(ma > mi) maxz = ma; else maxz = mi; break; case ORG: case REP: maxz = ma; minz = mi; break; case SQR: maxz = ma*ma; if(maxz < mi*mi) maxz = mi*mi; minz = 0.0; break; case IMP: maxz = 1.0; minz = 0.0; break; } } dat = fopen(name, "w+"); fprintf(dat, "c\n"); fprintf(dat, "z\n"); fprintf(dat, "%g\n", disp.y0); fprintf(dat, "%g\n", disp.x0); fprintf(dat, "%g\n", (disp.y1-disp.y0)/npy); fprintf(dat, "%g\n", (disp.x1-disp.x0)/npx); fprintf(dat, "%d\n", npy+1); fprintf(dat, "%d\n", npx+1); fprintf(dat, "v\n"); fprintf(dat, "19\n"); s = (maxz-minz)/20.0; for(xi=1; xi<=19; ++xi) fprintf(dat, "%g\n", minz+xi*s); for(xi=0; xi<=npx; ++xi) { x = disp.x0+(disp.x1-disp.x0)/npx*xi; for(yi=0; yi<=npy; ++yi) { y = disp.y0+(disp.y1-disp.y0)/npy*yi; w = f(xi, yi); fprintf(dat, "%g\n", w); } } fprintf(dat, "%g\n", minz); fprintf(dat, "%g\n", maxz); fprintf(dat, "%d\n", wg.nx); fprintf(dat, "%d\n", wg.ny); for(l=0; l<=wg.nx; ++l) fprintf(dat, "%g\n", wg.hx(l)); for(m=0; m<=wg.ny; ++m) fprintf(dat, "%g\n", wg.hy(m)); for(l=0; l<=wg.nx+1; ++l) { for(m=0; m<=wg.ny+1; ++m) fprintf(dat, "%g\n", wg.n(l, m)); } fprintf(dat, "%g\n", disp.x0); fprintf(dat, "%g\n", disp.x1); fprintf(dat, "%g\n", disp.y0); fprintf(dat, "%g\n", disp.y1); fclose(dat); return; } /* write profile section along line (x0,y0) -> (x1,y1) cp: EX - HZ, SZ ext0, ext1: filename id characters np: number of output points */ void Mode::writesec(Fcomp cp, char ext0, char ext1, int np, double x0, double y0, double x1, double y1) const { FILE *dat; int i; double x, y, dx, dy; char name[13] = "t______.xyf"; switch(pol) { case QTE: name[1] = 'e'; break; case QTM: name[1] = 'm'; break; case VEC: name[0] = 'v'; name[1] = 'c'; break; } switch(sym) { case SYM: name[2] = 's'; break; case ASY: name[2] = 'a'; break; case NOS: name[2] = '-'; break; } name[3] = fldchr(cp); name[4] = cpchr(cp); name[5] = ext0; name[6] = ext1; fprintf(stderr, ">> %s\n", name); dat = fopen(name, "w+"); dx = x1-x0; dy = y1-y0; for(i=0; i<=np; ++i) { x = x0+i*dx/np; y = y0+i*dy/np; if(dx == 0.0) fprintf(dat, "%g %g\n", y, field(cp, x, y)); else { if(dy == 0.0) fprintf(dat, "%g %g\n", x, field(cp, x, y)); else { fprintf(dat, "%g %g\n", i*sqrt(dx*dx+dy*dy)/np, field(cp, x, y)); } } } fclose(dat); return; } /* ------------------------------------------------------------------------ */ /* permittivity perturbation: propagation constant shift */ Complex Mode::phaseshift(Perturbation p) const { Complex pt; double nrm; Complex db; db = CC0; // xx pt = p.e(0, 0); if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EX, p.r); // yy pt = p.e(1, 1); if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EY, EY, p.r); // zz pt = p.e(2, 2); if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EZ, EZ, p.r); // xy, yx pt = p.e(0, 1) + p.e(1, 0); if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EY, p.r); // xz, zx pt = p.e(0, 2) - p.e(2, 0); pt = CCI*pt; if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EX, EZ, p.r); // yz, zy pt = p.e(1, 2) - p.e(2, 1); pt = CCI*pt; if(pt.re != 0.0 || pt.im != 0.0) db = db+pt*recint(EY, EZ, p.r); nrm = recint(HY, EX, wg.xyplane); nrm -= recint(HX, EY, wg.xyplane); nrm *= 2.0; db = db*(val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm); return db; } /* geometry variations: propagation constant shift due to moving the horizontal boundary between rectangles [l,m] and [l+1,m] ( moving hx(l) for hy(m-1) < y < hy(m) ) */ double Mode::horgeovar(int l, int m) const { double pt; double nrm; double y0, y1; if(l<0) return 0.0; if(m<0) return 0.0; if(l>wg.nx+1) return 0.0; if(m>wg.ny) return 0.0; if(m <= 0) y0 = wg.xyplane.y0; else y0 = wg.hy(m-1); if(m >= wg.ny+1) y1 = wg.xyplane.y1; else y1 = wg.hy(m); /* if(pol == QTM) { pt = 0.5*(horlineint(l+1, m, HY, HY, wg.hx(l), y0, y1) +horlineint(l, m, HY, HY, wg.hx(l), y0, y1)) /wg.eps(l,m)/wg.eps(l+1,m); pt += 0.5*(horlineint(l+1, m, EZ, EZ, wg.hx(l), y0, y1) +horlineint(l, m, EZ, EZ, wg.hx(l), y0, y1)) /val_invomep0(wg.lambda)/val_invomep0(wg.lambda) /beta/beta; pt *= (wg.eps(l,m)-wg.eps(l,m+1)); nrm = 0.0; for(int li=0; li<=wg.nx+1; ++li) { for(int mi=0; mi<=wg.ny+1; ++mi) nrm += recintlm(li, mi, HY, HY)/wg.eps(li,mi); } pt *= beta/2.0/nrm; return pt; } */ /* if(pol == QTE) { pt = 0.5*(horlineint(l+1, m, EY, EY, wg.hx(l), y0, y1) +horlineint(l, m, EY, EY, wg.hx(l), y0, y1)); pt *= (wg.eps(l,m)-wg.eps(l,m+1)); nrm = recint(EY, EY, wg.xyplane); pt *= val_k0(wg.lambda)*val_k0(wg.lambda)/2.0/beta/nrm; return pt; } */ pt = 0.5*horlineint(l+1, m, EX, EX, wg.hx(l), y0, y1) *wg.eps(l+1,m)/wg.eps(l,m); pt += 0.5*horlineint(l, m, EX, EX, wg.hx(l), y0, y1) *wg.eps(l,m)/wg.eps(l+1,m); pt += 0.5*(horlineint(l+1, m, EY, EY, wg.hx(l), y0, y1) +horlineint(l, m, EY, EY, wg.hx(l), y0, y1)); pt += 0.5*(horlineint(l+1, m, EZ, EZ, wg.hx(l), y0, y1) +horlineint(l, m, EZ, EZ, wg.hx(l), y0, y1)); pt *= (wg.eps(l,m)-wg.eps(l+1,m)); nrm = recint(HY, EX, wg.xyplane); nrm -= recint(HX, EY, wg.xyplane); nrm *= 2.0; pt *= val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm; return pt; } /* geometry variations: propagation constant shift due to moving the vertical boundary between rectangles [l,m] and [l,m+1] ( moving hy(m) for hx(l-1) < x < hx(l) ) */ double Mode::vergeovar(int l, int m) const { double pt; double nrm; double x0, x1; if(m<0) return 0.0; if(l<0) return 0.0; if(m>wg.ny) return 0.0; if(l>wg.nx+1) return 0.0; if(l == 0) x0 = wg.xyplane.x0; else x0 = wg.hx(l-1); if(l == wg.nx+1) x1 = wg.xyplane.x1; else x1 = wg.hx(l); if(pol == QTM) { pt = 0.5*(verlineint(l, m+1, HY, HY, x0, x1, wg.hy(m)) +verlineint(l, m, HY, HY, x0, x1, wg.hy(m))); pt *= (wg.eps(l,m)-wg.eps(l,m+1)); nrm = recint(HY, HY, wg.xyplane); pt *= val_k0(wg.lambda)*val_k0(wg.lambda)/2.0/beta/nrm; return pt; } pt = 0.5*(verlineint(l, m+1, EX, EX, x0, x1, wg.hy(m)) +verlineint(l, m, EX, EX, x0, x1, wg.hy(m))); pt += 0.5*verlineint(l, m+1, EY, EY, x0, x1, wg.hy(m)) *wg.eps(l,m+1)/wg.eps(l,m); pt += 0.5*verlineint(l, m, EY, EY, x0, x1, wg.hy(m)) *wg.eps(l,m)/wg.eps(l,m+1); pt += 0.5*(verlineint(l, m+1, EZ, EZ, x0, x1, wg.hy(m)) +verlineint(l, m, EZ, EZ, x0, x1, wg.hy(m))); pt *= (wg.eps(l,m)-wg.eps(l,m+1)); nrm = recint(HY, EX, wg.xyplane); nrm -= recint(HX, EY, wg.xyplane); nrm *= 2.0; pt *= val_k0(wg.lambda)*INVSQRTMU0/INVSQRTEP0/nrm; return pt; } /* ------------------------------------------------------------------------ */ /* output to file */ void Mode::write(FILE *dat) { writerdpid(dat); comment("MODE", dat); comment("pol", dat); fputint(pol, dat); comment("sym", dat); fputint(sym, dat); comment("beta", dat); fputdouble(beta, dat); comment("k0", dat); fputdouble(k0, dat); comment("neff", dat); fputdouble(neff, dat); comment("npcB", dat); fputdouble(npcB, dat); comment("minE", dat); fputdouble(minE, dat); comment("maxE", dat); fputdouble(maxE, dat); comment("minH", dat); fputdouble(minH, dat); comment("maxH", dat); fputdouble(maxH, dat); comment("minS", dat); fputdouble(minS, dat); comment("maxS", dat); fputdouble(maxS, dat); comment("ampfac", dat); fputdouble(ampfac, dat); wg.write(dat); writerdp(dat); return; } /* input from file */ void Mode::read(FILE *dat) { readrdpid(dat); pol = Polarization(fgetint(dat)); sym = Symmetry(fgetint(dat)); beta = fgetdouble(dat); k0 = fgetdouble(dat); neff = fgetdouble(dat); npcB = fgetdouble(dat); minE = fgetdouble(dat); maxE = fgetdouble(dat); minH = fgetdouble(dat); maxH = fgetdouble(dat); minS = fgetdouble(dat); maxS = fgetdouble(dat); ampfac = fgetdouble(dat); wg.read(dat); readrdp(dat); return; } /* output to default file */ void Mode::write_def(char ext0, char ext1) { char name[10]; FILE *dat; name[0] = 't'; switch(pol) { case QTE: name[1] = 'e'; break; case QTM: name[1] = 'm'; break; case VEC: name[0] = 'v'; name[1] = 'c'; break; } switch(sym) { case SYM: name[2] = 's'; break; case ASY: name[2] = 'a'; break; case NOS: name[2] = '-'; break; } name[3] = ext0; name[4] = ext1; name[5] = '.'; name[6] = 'm'; name[7] = 'o'; name[8] = 'd'; name[9] = 0; fprintf(stderr, "beta: %.10g, neff: %.10g, npcB: %.10g >> %s\n", beta, neff, npcB, name); dat = fopen(name, "w+"); write(dat); fclose(dat); return; } /* input from default file */ void Mode::read_def(Polarization p, Symmetry s, char ext0, char ext1) { char name[10]; FILE *dat; name[0] = 't'; switch(p) { case QTE: name[1] = 'e'; break; case QTM: name[1] = 'm'; break; case VEC: name[0] = 'v'; name[1] = 'c'; break; } switch(s) { case SYM: name[2] = 's'; break; case ASY: name[2] = 'a'; break; case NOS: name[2] = '-'; break; } name[3] = ext0; name[4] = ext1; name[5] = '.'; name[6] = 'm'; name[7] = 'o'; name[8] = 'd'; name[9] = 0; fprintf(stderr, "<< %s : ", name); dat = fopen(name, "r"); read(dat); fclose(dat); fprintf(stderr, "beta: %.10g, neff: %.10g, npcB: %.10g\n", beta, neff, npcB); return; } /* overlap 0.5*intint mode_EX*inif_HY - mode_EY*inif_HX dxdy with an initial field, with real HX and HY components, integrals are evaluated by gaussian quadrature formulas, separately on the fundamental rectangles, inside the domain rectangle r */ /* helper functions */ double Mode::ovleval(double (*inif)(Fcomp cp, double x, double y), double x, double y) const { double fex, fhy, fey, fhx; fhy = inif(HY, x, y); fex = 0; if(fhy != 0.0) fex = field(EX, x, y); fhx = inif(HX, x, y); fey = 0; if(fhx != 0.0) fey = field(EY, x, y); return 0.5*(fex*fhy-fey*fhx); } /* based on: Numerical Recipes in C --- The Art of Scientific Computing Press, Teukolsky, Vetterling, Flannery, Cambridge University Press, 1994 */ double Mode::ovlpix(double (*inif)(Fcomp cp, double x, double y), double x0, double x1, double y, int numx) const { int j; double xr,xm,dx,s; static double x[]={0.0,0.1488743389,0.4333953941, 0.6794095682,0.8650633666,0.9739065285}; static double w[]={0.0,0.2955242247,0.2692667193, 0.2190863625,0.1494513491,0.0666713443}; double ovl = 0.0; double a, b, step; int i; step = (x1-x0)/((double) numx); for(i=0; i<=numx-1; ++i) { a = x0+((double) i)*step; b = a+step; xm=0.5*(b+a); xr=0.5*(b-a); s=0.0; for (j=1;j<=5;j++) { dx=xr*x[j]; s += w[j]*(ovleval(inif, xm+dx, y) +ovleval(inif, xm-dx, y)); } ovl += s*xr; } return ovl; } double Mode::ovlpint(double (*inif)(Fcomp cp, double x, double y), Rect r, int numx, int numy) const { double ovl = 0.0; double a, b, step; int i; int j; double yr,ym,dy,s; static double y[]={0.0,0.1488743389,0.4333953941, 0.6794095682,0.8650633666,0.9739065285}; static double w[]={0.0,0.2955242247,0.2692667193, 0.2190863625,0.1494513491,0.0666713443}; step = (r.y1-r.y0)/((double) numy); for(i=0; i<=numy-1; ++i) { a = r.y0+((double) i)*step; b = a+step; ym=0.5*(b+a); yr=0.5*(b-a); s=0.0; for (j=1;j<=5;j++) { dy=yr*y[j]; s += w[j]*(ovlpix(inif, r.x0, r.x1, ym+dy, numx) +ovlpix(inif, r.x0, r.x1, ym-dy, numx)); } ovl += s*yr; } return ovl; } #define HDIST 1.0e-8 double Mode::ovlp(double (*inif)(Fcomp cp, double x, double y), Rect r, int numx, int numy) const { double ovl = 0.0; double xb, yl; double xt, yr; int l, m, dum; Rect rp; if(r.x0 > r.x1) { xb=r.x1; r.x1=r.x0; r.x0=xb; } if(r.y0 > r.y1) { yl=r.y1; r.y1=r.y0; r.y0=yl; } yl = r.y0; while(fabs(r.y1-yl) > HDIST) { wg.rectidx(0.0, yl+HDIST/2.0, dum, m); yr = wg.rectbounds(0,m).y1; if(r.y1 HDIST) { wg.rectidx(xb+HDIST/2.0, 0.0, l, dum); xt = wg.rectbounds(l,0).x1; if(r.x1> %s\n", fldchr(cp), cpchr(cp), name); dat = fopen(name, "w+"); mlout_title(dat, name, "WMM mode profile"); maxf = 1.0; minf = -1.0; switch(fldchr(cp)) { case 'E': switch(foa) { case ORG: case REP: minf = minE; maxf = maxE; break; case MOD: maxf = fabs(maxE); if(fabs(minE)>maxf) maxf=fabs(minE); minf = 0.0; break; case SQR: maxf = maxE*maxE; if(minE*minE>maxf) maxf=minE*minE; minf = 0.0; break; case IMP: // not reasonable maxf = 1.0; minf = 0.0; break; } break; case 'H': switch(foa) { case ORG: case REP: minf = minH; maxf = maxH; break; case MOD: maxf = fabs(maxH); if(fabs(minH)>maxf) maxf=fabs(minH); minf = 0.0; break; case SQR: maxf = maxH*maxH; if(minH*minH>maxf) maxf=minH*minH; minf = 0.0; break; case IMP: // not reasonable maxf = 1.0; minf = 0.0; break; } break; case 'S': switch(foa) { case ORG: case REP: case MOD: minf = minS; maxf = maxS; break; case SQR: // not reasonable maxf = maxS*maxS; minf = 0.0; break; case IMP: // not reasonable maxf = 1.0; minf = 0.0; break; } break; } mlout_geo(dat, wg, minf, maxf); mlout_meshxy(dat, disp, npx, npy); Dmatrix fld; fld = fieldmat(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); break; case 'S': mlout_surface(name, dat, cp, foa); break; case 'I': mlout_image(name, dat, cp, foa); break; default: break; } mlout_print(dat, name, 'e'); fclose(dat); return; } /* write all components to MATLAB .m file, surface plots ext0, ext1: filename id characters disp: output region on the x-y-plane npx, npy: number of points in output mesh */ void Mode::acmfile(char ext0, char ext1, Rect disp, int npx, int npy) const { FILE *dat; Dmatrix fld; char name[13] = "t____A.m"; switch(pol) { case QTE: name[1] = 'e'; break; case QTM: name[1] = 'm'; break; case VEC: name[0] = 'v'; name[1] = 'c'; break; } switch(sym) { case SYM: name[2] = 's'; break; case ASY: name[2] = 'a'; break; case NOS: name[2] = 'n'; break; } name[3] = ext0; name[4] = ext1; fprintf(stderr, "Ex - Hz >> %s\n", name); dat = fopen(name, "w+"); mlout_title(dat, name, "WMM mode profile"); name[6] = 0; mlout_geo(dat, wg, 0.0, 1.0); mlout_meshxy(dat, disp, npx, npy); fld = fieldmat(disp, npx, npy, EX, ORG); mlout_fld(dat, npx, npy, EX, fld); fld = fieldmat(disp, npx, npy, EY, ORG); mlout_fld(dat, npx, npy, EY, fld); fld = fieldmat(disp, npx, npy, EZ, ORG); mlout_fld(dat, npx, npy, EZ, fld); fld = fieldmat(disp, npx, npy, HX, ORG); mlout_fld(dat, npx, npy, HX, fld); fld = fieldmat(disp, npx, npy, HY, ORG); mlout_fld(dat, npx, npy, HY, fld); fld = fieldmat(disp, npx, npy, HZ, ORG); mlout_fld(dat, npx, npy, HZ, fld); mlout_acmfile(name, dat, minE, maxE, minH, maxH); mlout_print(dat, name, 'e'); fclose(dat); return; } /* write profile section along line (x0,y0) -> (x1,y1) to MATLAB .m file cp: EX - HZ, SZ ext0, ext1: filename id characters np: number of output points pltype: 'L': geometry information + plot commands (default) 'V': field, mesh, and plot command, to be included into *L.m --- at present implemented for horizontal or vertical lines only --- */ #define HDISTF 1.0e-8 void Mode::secmfile(Fcomp cp, char ext0, char ext1, int np, char pltype, double x0, double y0, double x1, double y1) const { FILE *dat; int i, j; double x, y, dx, dy, t; double xbd, ybd; double epsold, epsnew; double minf, maxf, f; int nsec, nbd; char name[13] = "t_______.m"; char ori; switch(pol) { case QTE: name[1] = 'e'; break; case QTM: name[1] = 'm'; break; case VEC: name[0] = 'v'; name[1] = 'c'; break; } switch(sym) { case SYM: name[2] = 's'; break; case ASY: name[2] = 'a'; break; case NOS: name[2] = '-'; break; } name[3] = fldchr(cp); name[4] = cpchr(cp); name[5] = ext0; name[6] = ext1; if(pltype != 'V') pltype = 'L'; name[7] = pltype; dx = fabs(x1-x0); dy = fabs(y1-y0); if(dx > 1.0e-10 && dy > 1.0e-10) { fprintf(stderr, "%c%c [%g, %g] -> [%g, %g] >> file \n", fldchr(cp), cpchr(cp), x0, y0, x1, y1); fprintf(stderr, " not implemented for tilted lines, use Mode.writesec().\n"); return; } if(dx < 1.0e-10 && dy < 1.0e-10) { fprintf(stderr, "Mode.secmfile(): Nothing to do !\n"); return; } if(dx < 1.0e-10) ori = 'h'; else ori = 'v'; if(ori == 'v') { if(x0 > x1) { t = x0; x0 = x1; x1 = t; }} else { if(y0 > y1) { t = y0; y0 = y1; y1 = t; }} fprintf(stderr, "%c%c [%g, %g] -> [%g, %g] >> %s\n", fldchr(cp), cpchr(cp), x0, y0, x1, y1, name); Dmatrix fld(wg.nx+2+wg.ny+2, np+1); Dmatrix pos(wg.nx+2+wg.ny+2, np+1); Ivector nump(wg.nx+2+wg.ny+2); Dvector of(np+1); Dvector op(np+1); nsec = 0; Dvector bd(wg.nx+1+wg.ny+1); Dvector ri(wg.nx+2+wg.ny+2); nbd = 0; dx = x1-x0; dy = y1-y0; nump(nsec) = 0; x = x0; y = y0; epsold = wg.eps(x, y); ri(nbd) = sqrt(epsold); if(ori=='h') pos(nsec, nump(nsec)) = y; else pos(nsec, nump(nsec)) = x; fld(nsec, nump(nsec)) = field(cp, x, y); ++nump(nsec); for(i=1; i<=np; ++i) { x = x0+i*dx/np; y = y0+i*dy/np; epsnew = wg.eps(x, y); if(fabs(epsnew-epsold)<1.0e-10) { if(ori=='h') pos(nsec, nump(nsec)) = y; else pos(nsec, nump(nsec)) = x; fld(nsec, nump(nsec)) = field(cp, x, y); ++nump(nsec); } else { if(ori=='h') { xbd = x0; ybd = (wg.rectbounds(x, y)).y0; } else { xbd = (wg.rectbounds(x, y)).x0; ybd = y0; } if(ori=='h') pos(nsec, nump(nsec)) = ybd; else pos(nsec, nump(nsec)) = xbd; fld(nsec, nump(nsec)) = field(cp, xbd-dx*HDISTF, ybd-dy*HDISTF); ++nump(nsec); if(ori=='h') bd(nbd) = ybd; else bd(nbd) = xbd; ++nbd; epsold = epsnew; ri(nbd) = sqrt(epsold); ++nsec; nump(nsec) = 0; if(ori=='h') pos(nsec, nump(nsec)) = ybd; else pos(nsec, nump(nsec)) = xbd; fld(nsec, nump(nsec)) = field(cp, xbd+dx*HDISTF, ybd+dy*HDISTF); ++nump(nsec); if(ori=='h') pos(nsec, nump(nsec)) = y; else pos(nsec, nump(nsec)) = x; fld(nsec, nump(nsec)) = field(cp, x, y); ++nump(nsec); } } ++nsec; dat = fopen(name, "w+"); mlout_title(dat, name, "WMM mode field section"); name[8] = 0; minf = -1.0; maxf = 1.0; if(pltype == 'L') { switch(fldchr(cp)) { case 'E': minf = minE; maxf = maxE; break; case 'H': minf = minH; maxf = maxH; break; case 'S': minf = minS; maxf = maxS; break; } mlout_geo(dat, wg, minf, maxf); mlout_Lsecgeo(dat, x0, y0, x1, y1, nbd, bd, ri); } minf = fld(0, 0); maxf = fld(0, 0); for(j=0; j<=nsec-1; ++j) { for(i=0; i<=nump(j)-1; ++i) { f = fld(j, i); if(f < minf) minf = f; if(f > maxf) maxf = f; of(i) = f; op(i) = pos(j, i); } if(pltype == 'L') mlout_sec1D(dat, cp, dig10(j), dig1(j), ' ', ' ', nump(j), of, op); else mlout_sec1D(dat, cp, dig10(j), dig1(j), ext0, ext1, nump(j), of, op); } if(pltype == 'L') { mlout_Lsecplot(name, dat, ori, minf, maxf, cp, nbd, nsec); mlout_print(dat, name, 'e'); } else mlout_Vsecplot(dat, cp, nsec, ext0, ext1); fclose(dat); return; } /* write single component to MATLAB .m file, fancy style :-) cp: EX - HZ, SZ disp: output region on the x-y-plane npx, npy: number of points in output mesh ext0, ext1: filename id characters */ void Mode::fancymfile(Fcomp cp, 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; char name[13] = "t______F.m"; switch(pol) { case QTE: name[1] = 'e'; break; case QTM: name[1] = 'm'; break; case VEC: name[0] = 'v'; name[1] = 'c'; break; } switch(sym) { case SYM: name[2] = 's'; break; case ASY: name[2] = 'a'; break; case NOS: name[2] = 'n'; break; } name[3] = fldchr(cp); name[4] = cpchr(cp); name[5] = ext0; name[6] = ext1; fprintf(stderr, "%c%c >> %s\n", fldchr(cp), cpchr(cp), name); dat = fopen(name, "w+"); mlout_title(dat, name, "WMM mode profile :-)"); name[8] = 0; switch(fldchr(cp)) { case 'E': mlout_geo(dat, wg, minE, maxE); break; case 'H': mlout_geo(dat, wg, minH, maxH); break; case 'S': mlout_geo(dat, wg, minS, maxS); break; } mlout_meshxy(dat, disp, npx, npy); Dmatrix fld; fld = fieldmat(disp, npx, npy, cp, ORG); mlout_fld(dat, npx, npy, cp, 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(cp, np, x0+HDIST, yp-HDIST, x1-HDIST, yp-HDIST); mlout_sec(dat, x0+HDIST, yp-HDIST, x1-HDIST, yp-HDIST, np, cp, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(cp, np, x0+HDIST, yp+HDIST, x1-HDIST, yp+HDIST); mlout_sec(dat, x0+HDIST, yp+HDIST, x1-HDIST, yp+HDIST, np, cp, 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(cp, np, xp-HDIST, y0+HDIST, xp-HDIST, y1-HDIST); mlout_sec(dat, xp-HDIST, y0+HDIST, xp-HDIST, y1-HDIST, np, cp, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(cp, np, xp+HDIST, y0+HDIST, xp+HDIST, y1-HDIST); mlout_sec(dat, xp+HDIST, y0+HDIST, xp+HDIST, y1-HDIST, np, cp, dig10(numc), dig1(numc), fv); ++numc; } } } } fv = fieldvec(cp, npx, disp.x0, disp.y0, disp.x1, disp.y0); mlout_sec(dat, disp.x0, disp.y0, disp.x1, disp.y0, npx, cp, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(cp, npy, disp.x1, disp.y0, disp.x1, disp.y1); mlout_sec(dat, disp.x1, disp.y0, disp.x1, disp.y1, npy, cp, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(cp, npx, disp.x1, disp.y1, disp.x0, disp.y1); mlout_sec(dat, disp.x1, disp.y1, disp.x0, disp.y1, npx, cp, dig10(numc), dig1(numc), fv); ++numc; fv = fieldvec(cp, npy, disp.x0, disp.y1, disp.x0, disp.y0); mlout_sec(dat, disp.x0, disp.y1, disp.x0, disp.y0, npy, cp, dig10(numc), dig1(numc), fv); ++numc; mlout_fancy(name, dat, cp, numc); mlout_print(dat, name, 'p'); fclose(dat); return; }