/* * WMM - Wave-matching method for mode analysis of dielectric waveguides * https://wmm.computational-photonics.eu/ */ /* * wmmdisc.cpp * spectral discretisation procedures --- set up the trial function set */ #include #include #include #include"complex.h" #include"matrix.h" #include"wmmmat.h" #include"waveg.h" #include"maxweq.h" #include"wmmtrif.h" #include"mode.h" #include"wmmfld.h" #include"wmmarr.h" #include"wmm.h" #include"wmmgam.h" #include"wmmpol.h" #include"wmmglob.h" #include"wmmdisc.h" /* error message */ void wmmdiscerror(const char *msg) { fprintf(stderr, "WMMDISC: ERROR: %s\n", msg); exit(1); } /* ................................................................... */ /* set trial function parameters, depending on propagation constant */ void settfpars() { int mcp, cp, l, m, j; double e; Rect r; /* fprintf(stderr, "\n settfpars: "); */ mcp = 0; if(Pol == VEC) mcp = 1; for(m=0; m<=Wg.ny+1; ++m) { for(l=0; l<=Wg.nx+1; ++l) { e = Wg.eps(l,m); r = Wg.rectbounds(l,m); for(cp=0; cp<=mcp; ++cp) { for(j=0; j<=Phi.ntf(cp,l,m)-1; ++j) (*Phi(cp,l,m,j)).initdeppars( Beta, e, K0, Pol, r); } } } return; } /* ................................................................... */ /* reduce spectral discretization, enforce numerically positive normalization matrices */ void reducespecdisc(double b) { int l, m, j, dim, offs, eq; int maxi; int maxnaf = 0; double maxa; double *snmat; double *ev; int r; int mn; Beta = b; for(m=0; m<=Wg.ny+1; ++m) { for(l=0; l<=Wg.nx+1; ++l) { if(Pol != VEC) mn = Phi.ntf(0,l,m); else mn = Phi.ntf(0,l,m)+Phi.ntf(1,l,m); if(mn > maxnaf) maxnaf = mn; } } settfpars(); snmat = new double[maxnaf*maxnaf]; for(m=0; m<=Wg.ny+1; ++m) { for(l=0; l<=Wg.nx+1; ++l) { if(Swf(l,m) != 0.0) { matto0(snmat, maxnaf, maxnaf); dim = fillsnmat(l,m, snmat, 0); r = choldc(snmat, dim, 0); while(r != 0) { matto0(snmat, maxnaf, maxnaf); dim = fillsnmat(l, m, snmat, 0); ev = new double[dim]; smallesteigval(snmat, dim, 1, ev, 0, 0); maxa = fabs(ev[0]); maxi = 0; for(j=1; j<=dim-1; ++j) { if(fabs(ev[j]) > maxa) { maxa = fabs(ev[j]); maxi = j; } } delete[] ev; eq = 0; j = 0; offs = (*Phi(0,l,m,0)).Amsi; while(eq == 0 && j<=Phi.ntf(0,l,m)-1) { if((*Phi(0,l,m,j)).Amsi-offs == maxi) eq = 1; else ++j; } if(eq == 1) { Phi.removetrifun( (*Phi(0,l,m,j)).Amsi ); --AmDim; } else { if(Pol != VEC) wmmdiscerror("reducepqdisc 1"); eq = 0; j = 0; while(eq == 0 && j<=Phi.ntf(1,l,m)-1) { if((*Phi(1,l,m,j)).Amsi-offs == maxi) eq = 1; else ++j; } if(eq != 1) wmmdiscerror("reducepqdisc 2"); Phi.removetrifun( (*Phi(1,l,m,j)).Amsi ); --AmDim; } matto0(snmat, maxnaf, maxnaf); dim = fillsnmat(l, m, snmat, 0); r = choldc(snmat, dim, 0); } } } } delete[] snmat; return; } /* ................................................................... */ /* the spectral discretization procedure */ /* maximum number of trial functions per rectangle */ #define MNTFPR 1000 /* helper function */ int topqdc(int npq, Fstype typ, char pfqf, int pv, int qv, double *afv, Trifun* tf, int ntf) { int j; Trifun *fun; if(ntf+npq >= MNTFPR) wmmdiscerror("MNTFPR overflow !!!"); for(j=0; j<=npq-1; ++j) { fun = &tf[ntf+j]; (*fun).Amsi = AmDim; (*fun).AFp = afv[j]; (*fun).Pvz = pv; (*fun).Qvz = qv; (*fun).Svz = '+'; (*fun).PfQf = pfqf; (*fun).Xr = GPxr; (*fun).Yr = GPyr; (*fun).GfTyp = typ; (*fun).Gfcof = 1.0; (*fun).Amplit = 1.0; ++AmDim; } ntf += npq; return ntf; } /* helper function */ int addidgf( Fstype typ, char pfqf, int pv, int qv, double afp, int si, double xr, double yr, char sv, Trifun* tf, int ntf) { Trifun *fun; if(ntf+1 >= MNTFPR) wmmdiscerror("MNTFPR overflow !!!"); fun = &tf[ntf]; (*fun).Amsi = si; (*fun).AFp = afp; (*fun).Pvz = pv; (*fun).Qvz = qv; (*fun).Svz = sv; (*fun).PfQf = pfqf; (*fun).Xr = xr; (*fun).Yr = yr; (*fun).GfTyp = typ; (*fun).Gfcof = 1.0; (*fun).Amplit = 1.0; return ntf+1; } /* the spectral discretization procedure */ void specdisc(double sdm, double b0, double b1, double b2) { int l, m, mm, j; int cp, mcp; int v, npq, st; Fstype typ; char sv; double yr; double w; double afv[MTFPSET]; Trifun tf[MNTFPR]; Trifun *fun; int ntf, ontf; Rect r; /* fprintf(stderr, "specdisc(%g, %.16g, %.16g, %.16g)\n", sdm, b0, b1, b2); */ /* fprintf(stderr, "sdm: %g db: %.16g\n", sdm, b2-b0); */ GFDiff = sdm*Par.fin_d_alpha +(1.0-sdm)*Par.ini_d_alpha; MaxWvpA = sdm*Par.fin_alpha_max+(1.0-sdm)*Par.ini_alpha_max; Vlsg = int(sdm*Par.fin_N_alpha +(1.0-sdm)*Par.ini_N_alpha); Beta = b1; SDBeta = Beta; AmDim = 0; mcp = 0; if(Pol == VEC) mcp = 1; for(m=0; m<=Wg.ny+1; ++m) { for(l=0; l<=Wg.nx+1; ++l) { for(cp=0; cp<=mcp; ++cp) { ntf = 0; st = Styp(l,m); if(cp == 1) { switch(st) { case VSYM: st = VASY; break; case VASY: st = VSYM; break; case IDSY: st = IDAS; break; case IDAS: st = IDSY; break; } } if(st == CALC || st == VSYM || st == VASY) { w = K0*K0*Wg.eps(l,m)-Beta*Beta; if(w>=0.0) { v = 1; w = sqrt(w);} else { v = -1; w = sqrt(-w);} r = Wg.rectbounds(l,m); GPxa = r.x0; GPxb = r.x1; GPya = r.y0; GPyb = r.y1; GPw = w; switch(Rtyp(l,m)) { case 0: if(v == 1) { GPxr = (GPxa+GPxb)/2.0; GPyr = (GPya+GPyb)/2.0; if(st == CALC || st == VSYM) { npq = afdisc(FS_A, '-', 1, 1, afv); ntf = topqdc(npq, FS_A, '-', 1, 1, afv, tf, ntf); } if(st == CALC || st == VASY) { npq = afdisc(FS_B, '-', 1, 1, afv); ntf = topqdc(npq, FS_B, '-', 1, 1, afv, tf, ntf); } if(st == CALC || st == VSYM) { npq = afdisc(FS_C, '-', 1, 1, afv); ntf = topqdc(npq, FS_C, '-', 1, 1, afv, tf, ntf); } if(st == CALC || st == VASY) { npq = afdisc(FS_D, '-', 1, 1, afv); ntf = topqdc(npq, FS_D, '-', 1, 1, afv, tf, ntf); } GPyr = (GPya+GPyb)/2.0; GPxr = GPxb; if(st == CALC || st == VSYM) { npq = afdisc(FS_F, 'p', 1, 1, afv); ntf = topqdc(npq, FS_F, 'p', 1, 1, afv, tf, ntf); } if(st == CALC || st == VASY) { npq = afdisc(FS_G, 'p', 1,-1, afv); ntf = topqdc(npq, FS_G, 'p', 1,-1, afv, tf, ntf); } GPxr = GPxa; if(st == CALC || st == VSYM) { npq = afdisc(FS_F, 'p',-1, 1, afv); ntf = topqdc(npq, FS_F, 'p',-1, 1, afv, tf, ntf); } if(st == CALC || st == VASY) { npq = afdisc(FS_G, 'p',-1,-1, afv); ntf = topqdc(npq, FS_G, 'p',-1,-1, afv, tf, ntf); } GPxr = (GPxa+GPxb)/2.0; GPyr = GPyb; if(st == CALC || st == VASY || st == VSYM) { npq = afdisc(FS_H, 'q', 1, 1, afv); ntf = topqdc(npq, FS_H, 'q', 1, 1, afv, tf, ntf); } if(st == CALC || st == VASY || st == VSYM) { npq = afdisc(FS_I, 'q',-1, 1, afv); ntf = topqdc(npq, FS_I, 'q',-1, 1, afv, tf, ntf); } GPyr = GPya; if(st == CALC) { npq = afdisc(FS_H, 'q', 1,-1, afv); ntf = topqdc(npq, FS_H, 'q', 1,-1, afv, tf, ntf); } if(st == CALC) { npq = afdisc(FS_I, 'q',-1,-1, afv); ntf = topqdc(npq, FS_I, 'q',-1,-1, afv, tf, ntf); } } else { GPxr = GPxb; GPyr = GPyb; if(st == CALC || st == VSYM || st == VASY) { npq = afdisc(FS_E, '-', 1, 1, afv); ntf = topqdc(npq, FS_E, '-', 1, 1, afv, tf, ntf); } GPxr = GPxb; GPyr = GPya; if(st == CALC) { npq = afdisc(FS_E, '-', 1,-1, afv); ntf = topqdc(npq, FS_E, '-', 1,-1, afv, tf, ntf); } GPxr = GPxa; GPyr = GPyb; if(st == CALC || st == VSYM || st == VASY) { npq = afdisc(FS_E, '-',-1, 1, afv); ntf = topqdc(npq, FS_E, '-',-1, 1, afv, tf, ntf); } GPxr = GPxa; GPyr = GPya; if(st == CALC) { npq = afdisc(FS_E, '-',-1,-1, afv); ntf = topqdc(npq, FS_E, '-',-1,-1, afv, tf, ntf); } GPxr = (GPxa+GPxb)/2.0; GPyr = GPyb; if(st == CALC || st == VSYM || st == VASY) { npq = afdisc(FS_H, 'p', 1, 1, afv); ntf = topqdc(npq, FS_H, 'p', 1, 1, afv, tf, ntf); } GPyr = GPya; if(st == CALC) { npq = afdisc(FS_H, 'p', 1,-1, afv); ntf = topqdc(npq, FS_H, 'p', 1,-1, afv, tf, ntf); } GPyr = GPyb; if(st == CALC || st == VSYM || st == VASY) { npq = afdisc(FS_I, 'p',-1, 1, afv); ntf = topqdc(npq, FS_I, 'p',-1, 1, afv, tf, ntf); } GPyr = GPya; if(st == CALC) { npq = afdisc(FS_I, 'p',-1,-1, afv); ntf = topqdc(npq, FS_I, 'p',-1,-1, afv, tf, ntf); } GPyr = (GPya+GPyb)/2.0; GPxr = GPxb; if(st == CALC || st == VSYM) { npq = afdisc(FS_F, 'q', 1, 1, afv); ntf = topqdc(npq, FS_F, 'q', 1, 1, afv, tf, ntf); } GPxr = GPxa; if(st == CALC || st == VSYM) { npq = afdisc(FS_F, 'q',-1, 1, afv); ntf = topqdc(npq, FS_F, 'q',-1, 1, afv, tf, ntf); } GPxr = GPxb; if(st == CALC || st == VASY) { npq = afdisc(FS_G, 'q', 1,-1, afv); ntf = topqdc(npq, FS_G, 'q', 1,-1, afv, tf, ntf); } GPxr = GPxa; if(st == CALC || st == VASY) { npq = afdisc(FS_G, 'q',-1,-1, afv); ntf = topqdc(npq, FS_G, 'q',-1,-1, afv, tf, ntf); } } break; case 1: if(v==1) { GPyr = (GPya+GPyb)/2.0; GPxr = GPxa; if(st == CALC || st == VSYM) { npq = afdisc(FS_F, 'p',-1, 1, afv); ntf = topqdc(npq, FS_F, 'p',-1, 1, afv, tf, ntf); } if(st == CALC || st == VASY) { npq = afdisc(FS_G, 'p',-1, 1, afv); ntf = topqdc(npq, FS_G, 'p',-1, 1, afv, tf, ntf); } } else { GPxr = GPxa; GPyr = GPyb; if(st == CALC || st == VSYM || st == VASY) { npq = afdisc(FS_E, '-',-1, 1, afv); ntf = topqdc(npq, FS_E, '-',-1, 1, afv, tf, ntf); } GPxr = GPxa; GPyr = GPya; if(st == CALC) { npq = afdisc(FS_E, '-',-1,-1, afv); ntf = topqdc(npq, FS_E, '-',-1,-1, afv, tf, ntf); } GPyr = (GPya+GPyb)/2.0; GPxr = GPxa; if(st == CALC || st == VSYM) { npq = afdisc(FS_F, 'q',-1, 1, afv); ntf = topqdc(npq, FS_F, 'q',-1, 1, afv, tf, ntf); } if(st == CALC || st == VASY) { npq = afdisc(FS_G, 'q',-1, 1, afv); ntf = topqdc(npq, FS_G, 'q',-1, 1, afv, tf, ntf); } } break; case 2: GPxr = GPxa; GPyr = GPya; npq = afdisc(FS_E, '-',-1,-1, afv); ntf = topqdc(npq, FS_E, '-',-1,-1, afv, tf, ntf); break; case 3: if(v==1) { GPxr = (GPxa+GPxb)/2.0; GPyr = GPya; npq = afdisc(FS_H, 'q', 1,-1, afv); ntf = topqdc(npq, FS_H, 'q', 1,-1, afv, tf, ntf); npq = afdisc(FS_I, 'q', 1,-1, afv); ntf = topqdc(npq, FS_I, 'q', 1,-1, afv, tf, ntf); } else { GPxr = GPxb; GPyr = GPya; npq = afdisc(FS_E, '-', 1,-1, afv); ntf = topqdc(npq, FS_E, '-', 1,-1, afv, tf, ntf); GPxr = GPxa; GPyr = GPya; npq = afdisc(FS_E, '-',-1,-1, afv); ntf = topqdc(npq, FS_E, '-',-1,-1, afv, tf, ntf); GPxr = (GPxa+GPxb)/2.0; GPyr = GPya; npq = afdisc(FS_H, 'p', 1,-1, afv); ntf = topqdc(npq, FS_H, 'p', 1,-1, afv, tf, ntf); npq = afdisc(FS_I, 'p', 1,-1, afv); ntf = topqdc(npq, FS_I, 'p', 1,-1, afv, tf, ntf); } break; case 4: GPxr = GPxb; GPyr = GPya; npq = afdisc(FS_E, '-', 1,-1, afv); ntf = topqdc(npq, FS_E, '-', 1,-1, afv, tf, ntf); break; case 5: if(v==1) { GPyr = (GPya+GPyb)/2.0; GPxr = GPxb; if(st == CALC || st == VSYM) { npq = afdisc(FS_F, 'p', 1, 1, afv); ntf = topqdc(npq, FS_F, 'p', 1, 1, afv, tf, ntf); } if(st == CALC || st == VASY) { npq = afdisc(FS_G, 'p', 1, 1, afv); ntf = topqdc(npq, FS_G, 'p', 1, 1, afv, tf, ntf); } } else { GPxr = GPxb; GPyr = GPyb; if(st == CALC || st == VSYM || st == VASY) { npq = afdisc(FS_E, '-', 1, 1, afv); ntf = topqdc(npq, FS_E, '-', 1, 1, afv, tf, ntf); } GPxr = GPxb; GPyr = GPya; if(st == CALC) { npq = afdisc(FS_E, '-', 1,-1, afv); ntf = topqdc(npq, FS_E, '-', 1,-1, afv, tf, ntf); } GPyr = (GPya+GPyb)/2.0; GPxr = GPxb; if(st == CALC || st == VSYM) { npq = afdisc(FS_F, 'q', 1, 1, afv); ntf = topqdc(npq, FS_F, 'q', 1, 1, afv, tf, ntf); } if(st == CALC || st == VASY) { npq = afdisc(FS_G, 'q', 1, 1, afv); ntf = topqdc(npq, FS_G, 'q', 1, 1, afv, tf, ntf); } } break; case 6: GPxr = GPxb; GPyr = GPyb; npq = afdisc(FS_E, '-', 1, 1, afv); ntf = topqdc(npq, FS_E, '-', 1, 1, afv, tf, ntf); break; case 7: if(v==1) { GPxr = (GPxa+GPxb)/2.0; GPyr = GPyb; npq = afdisc(FS_H, 'q', 1, 1, afv); ntf = topqdc(npq, FS_H, 'q', 1, 1, afv, tf, ntf); npq = afdisc(FS_I, 'q', 1, 1, afv); ntf = topqdc(npq, FS_I, 'q', 1, 1, afv, tf, ntf); } else { GPxr = GPxb; GPyr = GPyb; npq = afdisc(FS_E, '-', 1, 1, afv); ntf = topqdc(npq, FS_E, '-', 1, 1, afv, tf, ntf); GPxr = GPxa; GPyr = GPyb; npq = afdisc(FS_E, '-',-1, 1, afv); ntf = topqdc(npq, FS_E, '-',-1, 1, afv, tf, ntf); GPxr = (GPxa+GPxb)/2.0; GPyr = GPyb; npq = afdisc(FS_H, 'p', 1, 1, afv); ntf = topqdc(npq, FS_H, 'p', 1, 1, afv, tf, ntf); npq = afdisc(FS_I, 'p', 1, 1, afv); ntf = topqdc(npq, FS_I, 'p', 1, 1, afv, tf, ntf); } break; case 8: GPxr = GPxa; GPyr = GPyb; npq = afdisc(FS_E, '-',-1, 1, afv); ntf = topqdc(npq, FS_E, '-',-1, 1, afv, tf, ntf); break; } } if(st == VSYM || st == VASY) { if(st == VSYM) sv = '+'; else sv = '-'; ontf = ntf; for(j=0; j<=ontf-1; ++j) { fun = &tf[j]; typ = (*fun).GfTyp; if(typ == FS_E || typ == FS_H || typ == FS_I) { ntf = addidgf(typ, (*fun).PfQf, (*fun).Pvz, -(*fun).Qvz, (*fun).AFp, (*fun).Amsi, (*fun).Xr, r.y0, sv, tf, ntf); } } } if(st == VSYM || st == VASY || st == CALC) { Phi.insert(tf, ntf, cp, l, m); } } } } for(m=0; m<=Wg.ny+1; ++m) { for(l=0; l<=Wg.nx+1; ++l) { for(cp=0; cp<=mcp; ++cp) { st = Styp(l,m); if(cp == 1) { switch(st) { case VSYM: st = VASY; break; case VASY: st = VSYM; break; case IDSY: st = IDAS; break; case IDAS: st = IDSY; break; } } if(st == IDSY || st == IDAS) { r = Wg.rectbounds(l, m); ntf = 0; if(st == IDSY) sv = '+'; else sv = '-'; mm = Wg.ny+1-m; npq = Phi.ntf(cp, l, mm); for(j=0; j<=npq-1; ++j) { fun = Phi(cp, l, mm, j); typ = (*fun).GfTyp; if(typ == FS_E || typ == FS_H || typ == FS_I) { if(-(*fun).Qvz > 0) yr = r.y1; else yr = r.y0; } else yr = 0.5*(r.y0+r.y1); ntf = addidgf(typ, (*fun).PfQf, (*fun).Pvz, -(*fun).Qvz, (*fun).AFp, (*fun).Amsi, (*fun).Xr, yr, sv, tf, ntf); } Phi.insert(tf, ntf, cp, l, m); } } } } settfpars(); j = AmDim; fprintf(stderr, " # %d ", AmDim); if(b0 == b2) reducespecdisc(b1); else { reducespecdisc(b1); reducespecdisc(b0); reducespecdisc(b2); } if(AmDim < j) fprintf(stderr, "> %d\n", AmDim); else fprintf(stderr, "\n"); return; } /* save spectral discretization for diagnosis purposes */ /* helper function: Fstype -> char */ char fstoch(Fstype t) { return 'A'+char(int(t)); } void diagspecdisc() { FILE *dat; int l, m, j; Trifun *fun; int cp, mcp; mcp = 0; if(Pol == VEC) mcp = 1; fprintf(stderr, ">> pjqj\n"); dat = fopen("pjqj", "w+"); fprintf(dat, "beta: %.15g\n", Beta); fprintf(dat,"# trial functions:\n"); for(m=0; m<=Wg.ny+1; ++m) { for(l=0; l<=Wg.nx+1; ++l) { if(Pol == VEC) fprintf(dat, " [%d,%d] - 0: %d 1: %d\n", l, m, Phi.ntf(0,l,m), Phi.ntf(1,l,m)); else fprintf(dat, " [%d,%d] - %d\n", l, m, Phi.ntf(0,l,m)); } } fprintf(dat,"amplitude matrix:\n"); fprintf(dat," dimension: %d \n", AmDim); fprintf(dat,"trial function parameters:\n"); for(m=0; m<=Wg.ny+1; ++m) { for(l=0; l<=Wg.nx+1; ++l) { for(cp=0; cp<=mcp; ++cp) { fprintf(dat, "[%d,%d] --- cp: %d\n", l, m, cp); fprintf(dat, " #GF: %d STYP: %d\n", Phi.ntf(cp,l,m), Styp(l,m)); for(j=0; j<=Phi.ntf(cp,l,m)-1; ++j) { fun = Phi(cp,l,m,j); fprintf(dat, "(%d)(ai: %d) p: %g q: %g", j, (*fun).Amsi, (*fun).Pj, (*fun).Qj); fprintf(dat, " Gtyp: %c*%g ", fstoch((*fun).GfTyp), (*fun).Gfcof); fprintf(dat, "#: (%g, %g)\n", (*fun).Xr, (*fun).Yr); } fprintf(dat, "\n"); } } } fclose(dat); return; }