/* * WMM - Wave-matching method for mode analysis of dielectric waveguides * https://wmm.computational-photonics.eu/ */ /* * wmmtrif.cpp * the basic solutions to the Helmholtz wave equation */ #include #include #include #include"inout.h" #include"complex.h" #include"matrix.h" #include"waveg.h" #include"maxweq.h" #include"wmmtrif.h" /* error message */ void wmmtriferror(const char *m) { fprintf(stderr, "\n ERR: %s\n", m); exit(1); } /* constructor */ Trifun::Trifun() { Amsi = 0; Gfcof = 1.0; GfTyp = FS_E; AFp = 1.0; PfQf = 'p'; Pj = 1.0; Pvz = 1; Qj = 1.0; Qvz = 1; Xr = 0.0; Yr = 0.0; Svz = '-'; Amplit = 1.0; dgftypes = Ivector(6); dgfcoeff = Dvector(6); dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; } /* constructor */ Trifun::Trifun( int amsi, double gfcof, Fstype gftyp, double afp, char pfqf, double pj, int pvz, double qj, int qvz, double xr, double yr, char svz, double amplit) { Amsi = amsi; Gfcof = gfcof; GfTyp = gftyp; AFp = afp; PfQf = pfqf; Pj = pj; Pvz = pvz; Qj = qj; Qvz = qvz; Xr = xr; Yr = yr; Svz = svz; Amplit = amplit; dgftypes = Ivector(6); dgfcoeff = Dvector(6); dgftypes(0) = int(GfTyp); dgfcoeff(0) = 1.0; } /* trial function output */ void Trifun::write(FILE *dat) { comment("TRIFUN", dat); fputint(Amsi, dat); fputdouble(Gfcof, dat); fputint(GfTyp, dat); fputdouble(AFp, dat); fputchar(PfQf, dat); fputdouble(Pj, dat); fputint(Pvz, dat); fputdouble(Qj, dat); fputint(Qvz, dat); fputdouble(Xr, dat); fputdouble(Yr, dat); fputchar(Svz, dat); fputdouble(Amplit, dat); dgftypes.write(dat); dgfcoeff.write(dat); return; } /* input trial function */ void Trifun::read(FILE *dat) { Amsi = fgetint(dat); Gfcof = fgetdouble(dat); GfTyp = Fstype(fgetint(dat)); AFp = fgetdouble(dat); PfQf = fgetchar(dat); Pj = fgetdouble(dat); Pvz = fgetint(dat); Qj = fgetdouble(dat); Qvz = fgetint(dat); Xr = fgetdouble(dat); Yr = fgetdouble(dat); Svz = fgetchar(dat); Amplit = fgetdouble(dat); dgftypes.read(dat); dgfcoeff.read(dat); return; } /* translate by (dx, dy) */ void Trifun::translate(double dx, double dy) { Xr += dx; Yr += dy; return; } /* initialize dependent parameters in Trifun */ void Trifun::initdeppars(double beta, double eps, double k0, Polarization p, Rect r) { double w, pf, nrmf; w = k0*k0*eps-beta*beta; if(w>=0.0) w = sqrt(w); else w = sqrt(-w); switch(GfTyp) { case FS_A: case FS_B: case FS_C: case FS_D: case FS_E: Pj = ((double) Pvz)*w*cos(AFp); Qj = ((double) Qvz)*w*sin(AFp); break; case FS_F: case FS_G: case FS_H: case FS_I: switch(PfQf) { case 'p': Pj = ((double) Pvz)*w*sinh(AFp); Qj = ((double) Qvz)*w*cosh(AFp); break; case 'q': Pj = ((double) Pvz)*w*cosh(AFp); Qj = ((double) Qvz)*w*sinh(AFp); break; } break; } pf = 1.0; if(p == QTM) pf = 1.0/eps; Gfcof = 1.0; Amplit = 1.0; dgftypes(0) = int(GfTyp); dgfcoeff(0) = 1.0; nrmf = 1.0/sqrt(pf*tfnrm(*this, FLD, *this, FLD, r)); if(Svz == '-') Gfcof = -nrmf; else Gfcof = nrmf; switch(GfTyp) { case FS_A: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_C); dgfcoeff(1) = -Pj*Gfcof; dgftypes(2) = int(FS_B); dgfcoeff(2) = -Qj*Gfcof; dgftypes(3) = int(FS_D); dgfcoeff(3) = Pj*Qj*Gfcof; dgftypes(4) = int(FS_A); dgfcoeff(4) = -Pj*Pj*Gfcof; dgftypes(5) = int(FS_A); dgfcoeff(5) = -Qj*Qj*Gfcof; break; case FS_B: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_D); dgfcoeff(1) = -Pj*Gfcof; dgftypes(2) = int(FS_A); dgfcoeff(2) = Qj*Gfcof; dgftypes(3) = int(FS_C); dgfcoeff(3) = -Pj*Qj*Gfcof; dgftypes(4) = int(FS_B); dgfcoeff(4) = -Pj*Pj*Gfcof; dgftypes(5) = int(FS_B); dgfcoeff(5) = -Qj*Qj*Gfcof; break; case FS_C: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_A); dgfcoeff(1) = Pj*Gfcof; dgftypes(2) = int(FS_D); dgfcoeff(2) = -Qj*Gfcof; dgftypes(3) = int(FS_B); dgfcoeff(3) = -Pj*Qj*Gfcof; dgftypes(4) = int(FS_C); dgfcoeff(4) = -Pj*Pj*Gfcof; dgftypes(5) = int(FS_C); dgfcoeff(5) = -Qj*Qj*Gfcof; break; case FS_D: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_B); dgfcoeff(1) = Pj*Gfcof; dgftypes(2) = int(FS_C); dgfcoeff(2) = Qj*Gfcof; dgftypes(3) = int(FS_A); dgfcoeff(3) = Pj*Qj*Gfcof; dgftypes(4) = int(FS_D); dgfcoeff(4) = -Pj*Pj*Gfcof; dgftypes(5) = int(FS_D); dgfcoeff(5) = -Qj*Qj*Gfcof; break; case FS_E: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_E); dgfcoeff(1) = Pj*Gfcof; dgftypes(2) = int(FS_E); dgfcoeff(2) = Qj*Gfcof; dgftypes(3) = int(FS_E); dgfcoeff(3) = Pj*Qj*Gfcof; dgftypes(4) = int(FS_E); dgfcoeff(4) = Pj*Pj*Gfcof; dgftypes(5) = int(FS_E); dgfcoeff(5) = Qj*Qj*Gfcof; break; case FS_F: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_F); dgfcoeff(1) = Pj*Gfcof; dgftypes(2) = int(FS_G); dgfcoeff(2) = -Qj*Gfcof; dgftypes(3) = int(FS_G); dgfcoeff(3) = -Pj*Qj*Gfcof; dgftypes(4) = int(FS_F); dgfcoeff(4) = Pj*Pj*Gfcof; dgftypes(5) = int(FS_F); dgfcoeff(5) = -Qj*Qj*Gfcof; break; case FS_G: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_G); dgfcoeff(1) = Pj*Gfcof; dgftypes(2) = int(FS_F); dgfcoeff(2) = Qj*Gfcof; dgftypes(3) = int(FS_F); dgfcoeff(3) = Pj*Qj*Gfcof; dgftypes(4) = int(FS_G); dgfcoeff(4) = Pj*Pj*Gfcof; dgftypes(5) = int(FS_G); dgfcoeff(5) = -Qj*Qj*Gfcof; break; case FS_H: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_I); dgfcoeff(1) = -Pj*Gfcof; dgftypes(2) = int(FS_H); dgfcoeff(2) = Qj*Gfcof; dgftypes(3) = int(FS_I); dgfcoeff(3) = -Pj*Qj*Gfcof; dgftypes(4) = int(FS_H); dgfcoeff(4) = -Pj*Pj*Gfcof; dgftypes(5) = int(FS_H); dgfcoeff(5) = Qj*Qj*Gfcof; break; case FS_I: dgftypes(0) = int(GfTyp); dgfcoeff(0) = Gfcof; dgftypes(1) = int(FS_H); dgfcoeff(1) = Pj*Gfcof; dgftypes(2) = int(FS_I); dgfcoeff(2) = Qj*Gfcof; dgftypes(3) = int(FS_H); dgfcoeff(3) = Pj*Qj*Gfcof; dgftypes(4) = int(FS_I); dgfcoeff(4) = -Pj*Pj*Gfcof; dgftypes(5) = int(FS_I); dgfcoeff(5) = Qj*Qj*Gfcof; break; } return; } /* ----------------------------------------------------------------- */ /* 1-D integration of products of harmonics/exponentials */ #define PQTolCC 1.0e-5 #define PQTolSS 1.0e-5 #define PQTolCS 1.0e-5 #define PQTolCE 1.0e-5 #define PQTolEE 1.0e-5 #define AWAY_2 5.0e+19 // smaller than [away] in "waveg.c" inline double intfCC(const double& p1, const double& r1, const double& p2, const double& r2, const double& x0, const double& x1) { double pr, a, b, s1, s2; if(fabs(p1-p2)AWAY_2) s1 = 0.0; else s1 = exp(p1*(x1-r1)+p2*(x1-r2)); if(fabs(x0)>AWAY_2) s2 = 0.0; else s2 = exp(p1*(x0-r1)+p2*(x0-r2)); return (s1-s2)/(p1+p2); } } inline double intfEC(const double& p1, const double& r1, const double& p2, const double& r2, const double& x0, const double& x1) { double da2, db2, s1, s2; da2 = p2*(x0-r2); db2 = p2*(x1-r2); if(fabs(x1)>AWAY_2) s1 = 0.0; else s1 = exp(p1*(x1-r1))*(p2*sin(db2)+p1*cos(db2)); if(fabs(x0)>AWAY_2) s2 = 0.0; else s2 = exp(p1*(x0-r1))*(p2*sin(da2)+p1*cos(da2)); return (s1-s2)/(p1*p1+p2*p2); } inline double intfES(const double& p1, const double& r1, const double& p2, const double& r2, const double& x0, const double& x1) { double da2, db2, s1, s2; da2 = p2*(x0-r2); db2 = p2*(x1-r2); if(fabs(x1)>AWAY_2) s1 = 0.0; else s1 = exp(p1*(x1-r1))*(p1*sin(db2)-p2*cos(db2)); if(fabs(x0)>AWAY_2) s2 = 0.0; else s2 = exp(p1*(x0-r1))*(p1*sin(da2)-p2*cos(da2)); return (s1-s2)/(p1*p1+p2*p2); } inline double intfCE(const double& p1, const double& r1, const double& p2, const double& r2, const double& x0, const double& x1) { double da1, db1, s1, s2; da1 = p1*(x0-r1); db1 = p1*(x1-r1); if(fabs(x1)>AWAY_2) s1 = 0.0; else s1 = exp(p2*(x1-r2))*(p1*sin(db1)+p2*cos(db1)); if(fabs(x0)>AWAY_2) s2 = 0.0; else s2 = exp(p2*(x0-r2))*(p1*sin(da1)+p2*cos(da1)); return (s1-s2)/(p1*p1+p2*p2); } inline double intfSE(const double& p1, const double& r1, const double& p2, const double& r2, const double& x0, const double& x1) { double da1, db1, s1, s2; da1 = p1*(x0-r1); db1 = p1*(x1-r1); if(fabs(x1)>AWAY_2) s1 = 0.0; else s1 = exp(p2*(x1-r2))*(p2*sin(db1)-p1*cos(db1)); if(fabs(x0)>AWAY_2) s2 = 0.0; else s2 = exp(p2*(x0-r2))*(p2*sin(da1)-p1*cos(da1)); return (s1-s2)/(p1*p1+p2*p2); } typedef double (*DtoDfunP) (double); typedef double (*IntfunP) (const double&, const double&, const double&, const double&, const double&, const double&); /* cos, sin, exp */ #define CSE_C 0 #define CSE_S 1 #define CSE_E 2 /* Fstype A-I -> cos, sin, exp id */ int xfstocse[9]; int yfstocse[9]; /* Fstype A-I -> cos, sin, exp - function */ DtoDfunP xfstocsefun[9]; DtoDfunP yfstocsefun[9]; /* cse, cse -> integral function */ IntfunP ccsseeint[3][3]; /* initialize tables for trial function integration */ void inittfinttables() { xfstocse[int(FS_A)] = CSE_C; xfstocse[int(FS_B)] = CSE_C; xfstocse[int(FS_C)] = CSE_S; xfstocse[int(FS_D)] = CSE_S; xfstocse[int(FS_E)] = CSE_E; xfstocse[int(FS_F)] = CSE_E; xfstocse[int(FS_G)] = CSE_E; xfstocse[int(FS_H)] = CSE_C; xfstocse[int(FS_I)] = CSE_S; yfstocse[int(FS_A)] = CSE_C; yfstocse[int(FS_B)] = CSE_S; yfstocse[int(FS_C)] = CSE_C; yfstocse[int(FS_D)] = CSE_S; yfstocse[int(FS_E)] = CSE_E; yfstocse[int(FS_F)] = CSE_C; yfstocse[int(FS_G)] = CSE_S; yfstocse[int(FS_H)] = CSE_E; yfstocse[int(FS_I)] = CSE_E; xfstocsefun[int(FS_A)] = &cos; xfstocsefun[int(FS_B)] = &cos; xfstocsefun[int(FS_C)] = &sin; xfstocsefun[int(FS_D)] = &sin; xfstocsefun[int(FS_E)] = &exp; xfstocsefun[int(FS_F)] = &exp; xfstocsefun[int(FS_G)] = &exp; xfstocsefun[int(FS_H)] = &cos; xfstocsefun[int(FS_I)] = &sin; yfstocsefun[int(FS_A)] = &cos; yfstocsefun[int(FS_B)] = &sin; yfstocsefun[int(FS_C)] = &cos; yfstocsefun[int(FS_D)] = &sin; yfstocsefun[int(FS_E)] = &exp; yfstocsefun[int(FS_F)] = &cos; yfstocsefun[int(FS_G)] = &sin; yfstocsefun[int(FS_H)] = &exp; yfstocsefun[int(FS_I)] = &exp; ccsseeint[CSE_C][CSE_C] = &intfCC; ccsseeint[CSE_C][CSE_S] = &intfCS; ccsseeint[CSE_S][CSE_C] = &intfSC; ccsseeint[CSE_S][CSE_S] = &intfSS; ccsseeint[CSE_E][CSE_E] = &intfEE; ccsseeint[CSE_E][CSE_C] = &intfEC; ccsseeint[CSE_E][CSE_S] = &intfES; ccsseeint[CSE_C][CSE_E] = &intfCE; ccsseeint[CSE_S][CSE_E] = &intfSE; // fprintf(stderr, "- Tfinit -\n"); } Tfintegrate::Tfintegrate() { inittfinttables(); } Tfintegrate TFINTINIT; /* ------------------------------------------------------------------ */ /* value of field/derivative d at position x, y */ double Trifun::val(Derlev d, double x, double y) const { return Amplit*DGfcof(d) *xfstocsefun[int(DGfTyp(d))](Pj*(x-Xr)) *yfstocsefun[int(DGfTyp(d))](Qj*(y-Yr)); } /* for factorizing trial functions: values of x- and y-part or derivative d at position x or y */ double Trifun::xval(Derlev d, double x) const { return xfstocsefun[int(DGfTyp(d))](Pj*(x-Xr)); } double Trifun::yval(Derlev d, double y) const { return yfstocsefun[int(DGfTyp(d))](Qj*(y-Yr)); } /* ------------------------------------------------------------------ */ /* integration of trial functions along vertical lines */ double tfintx(const Trifun& tf1, Derlev d1, const Trifun& tf2, Derlev d2, double y, double xa, double xb) { double f1, f2, xi; f1 = tf1.yval(d1, y); f2 = tf2.yval(d2, y); xi = ccsseeint[xfstocse[int(tf1.DGfTyp(d1))]] [xfstocse[int(tf2.DGfTyp(d2))]] (tf1.Pj, tf1.Xr, tf2.Pj, tf2.Xr, xa, xb); return f1*f2*xi*tf1.Amplit*tf2.Amplit*tf1.DGfcof(d1)*tf2.DGfcof(d2); } /* integration of trial functions along horizontal lines */ double tfinty(const Trifun& tf1, Derlev d1, const Trifun& tf2, Derlev d2, double x, double ya, double yb) { double f1, f2, yi; f1 = tf1.xval(d1, x); f2 = tf2.xval(d2, x); yi = ccsseeint[yfstocse[int(tf1.DGfTyp(d1))]] [yfstocse[int(tf2.DGfTyp(d2))]] (tf1.Qj, tf1.Yr, tf2.Qj, tf2.Yr, ya, yb); return f1*f2*yi*tf1.Amplit*tf2.Amplit*tf1.DGfcof(d1)*tf2.DGfcof(d2); } /* integration of trial functions over rectangles */ double tfnrm(const Trifun& tf1, Derlev d1, const Trifun& tf2, Derlev d2, Rect r) { double xi, yi; /* fprintf(stderr, "n: %c *%g * %c *%g, [%g, %g | %g, %g] -> ", 'A'+int(tf1.DGfTyp(d1)), tf1.DGfcof(d1), 'A'+int(tf2.DGfTyp(d2)), tf2.DGfcof(d2), r.x0, r.y0, r.x1, r.y1); */ xi = ccsseeint[xfstocse[int(tf1.DGfTyp(d1))]] [xfstocse[int(tf2.DGfTyp(d2))]] (tf1.Pj, tf1.Xr, tf2.Pj, tf2.Xr, r.x0, r.x1); yi = ccsseeint[yfstocse[int(tf1.DGfTyp(d1))]] [yfstocse[int(tf2.DGfTyp(d2))]] (tf1.Qj, tf1.Yr, tf2.Qj, tf2.Yr, r.y0, r.y1); /* fprintf(stderr, "%g\n", xi*yi*tf1.Amplit*tf2.Amplit*tf1.DGfcof(d1)*tf2.DGfcof(d2)); */ return xi*yi*tf1.Amplit*tf2.Amplit*tf1.DGfcof(d1)*tf2.DGfcof(d2); } /* ----------------------------------------------------------------- */ /* initialize trial function handling */ Trifunfield::Trifunfield() { nlx = 0; nly = 0; tfn = NULL; tmem = NULL; } Trifunfield::Trifunfield(int nx, int ny) { nlx = nx; nly = ny; tfn = new int[2*(nlx+2)*(nly+2)]; tmem = new Trifun*[2*(nlx+2)*(nly+2)]; int j; for(j=0; j<=2*(nlx+2)*(nly+2)-1; ++j) { tfn[j] = 0; tmem[j] = NULL; } } /* destroy */ Trifunfield::~Trifunfield() { int j; if(tmem != NULL) { for(j=0; j<=2*(nlx+2)*(nly+2)-1; ++j) { if(tmem[j] != NULL) delete[] tmem[j]; } delete[] tfn; delete[] tmem; } } /* copy constructor */ Trifunfield::Trifunfield(const Trifunfield& t) { int d; nlx = t.nlx; nly = t.nly; tfn = new int[2*(nlx+2)*(nly+2)]; tmem = new Trifun*[2*(nlx+2)*(nly+2)]; for(int j=0; j<=2*(nlx+2)*(nly+2)-1; ++j) { d = t.tfn[j]; tfn[j] = d; tmem[j] = new Trifun[d]; for(int k=0; k<=d-1; ++k) { tmem[j][k] = t.tmem[j][k]; } } } /* assignment */ Trifunfield& Trifunfield::operator=(const Trifunfield& t) { int d; if(this != &t) { if(tmem != NULL) { for(int j=0; j<=2*(nlx+2)*(nly+2)-1; ++j) { if(tmem[j] != NULL) delete[] tmem[j]; } delete[] tfn; delete[] tmem; } nlx = t.nlx; nly = t.nly; tfn = new int[2*(nlx+2)*(nly+2)]; tmem = new Trifun*[2*(nlx+2)*(nly+2)]; for(int j=0; j<=2*(nlx+2)*(nly+2)-1; ++j) { d = t.tfn[j]; tfn[j] = d; tmem[j] = new Trifun[d]; for(int k=0; k<=d-1; ++k) { tmem[j][k] = t.tmem[j][k]; } } } return *this; } /* insert a trial function set */ void Trifunfield::insert(Trifun *s, int nf, int cp, int l, int m) { int j, k; j = cplmtoid(cp, l, m); tfn[j] = nf; if(tmem[j] != NULL) delete[] tmem[j]; tmem[j] = new Trifun[nf]; for(k=0; k<=nf-1; ++k) tmem[j][k] = s[k]; return; } /* remove trial functions with Amsi si */ void Trifunfield::removetrifun(int si) { int l, m, j, s, amsi, k, i; int sfl[10]; int sfm[10]; int sfj[10]; int sfc[10]; int numsf = -1; int cp, cplm; if(tfn==NULL) return; for(m=0; m<=nly+1; ++m) { for(l=0; l<=nlx+1; ++l) { for(cp=0; cp<=1; ++cp) { for(j=0; j<=tfn[cplmtoid(cp, l, m)]-1; ++j) { amsi = (tmem[cplmtoid(cp, l, m)][j]).Amsi; if(amsi == si) { ++numsf; sfl[numsf] = l; sfm[numsf] = m; sfc[numsf] = cp; sfj[numsf] = j; } if(amsi > si) (tmem[cplmtoid(cp, l, m)][j]).Amsi = amsi-1; } } } } /* fprintf(stderr, "TF removed: \n"); */ for(s=0; s<=numsf; ++s) { l = sfl[s]; m = sfm[s]; cp = sfc[s]; j = sfj[s]; cplm = cplmtoid(cp,l,m); /* fprintf(stderr, " *[cp,l,m,j] = [%d, %d, %d, %d]\n", cp,l,m,j); fprintf(stderr, " Amsi = %d\n", (tmem[cplm][j]).Amsi); fprintf(stderr, " AFp = %g\n", (tmem[cplm][j]).AFp); fprintf(stderr, " Pvz = %d\n", (tmem[cplm][j]).Pvz); fprintf(stderr, " Qvz = %d\n", (tmem[cplm][j]).Qvz); fprintf(stderr, " Svz = %c\n", (tmem[cplm][j]).Svz); fprintf(stderr, " PfQf = %c\n", (tmem[cplm][j]).PfQf); fprintf(stderr, " Xr = %g\n", (tmem[cplm][j]).Xr); fprintf(stderr, " Yr = %g\n", (tmem[cplm][j]).Yr); fprintf(stderr, " GfTyp= %c\n", (tmem[cplm][j]).GfTyp); fprintf(stderr, " Pj = %g\n", (tmem[cplm][j]).Pj); fprintf(stderr, " Qj = %g\n", (tmem[cplm][j]).Qj); fprintf(stderr, " Gfcof= %g\n", (tmem[cplm][j]).Gfcof); */ for(k=j; k<=tfn[cplm]-2; ++k) { tmem[cplm][k] = tmem[cplm][k+1]; } --tfn[cplm]; for(i=s+1; i<=numsf; ++i) { if(sfl[i] == l && sfm[i] == m && sfc[i] == cp && sfj[i] > j) --sfj[i]; } } return; } /* translate by (dx, dy) */ void Trifunfield::translate(double dx, double dy) { int cp, l, m, j, e; for(cp=0; cp<=1; ++cp) { for(l=0; l<=nlx+1; ++l) { for(m=0; m<=nly+1; ++m) { e = cplmtoid(cp, l, m); for(j=0; j<=tfn[e]-1; ++j) tmem[e][j].translate(dx, dy); } } } return; } /* output to file */ void Trifunfield::write(FILE *dat) { int cp, l, m, j, e; comment("TRIFUNFIELD", dat); fputint(nlx, dat); fputint(nly, dat); for(cp=0; cp<=1; ++cp) { for(l=0; l<=nlx+1; ++l) { for(m=0; m<=nly+1; ++m) { e = cplmtoid(cp, l, m); comment("CELL", dat); fputint(cp, dat); fputint(l, dat); fputint(m, dat); comment("ntf", dat); fputint(tfn[e], dat); for(j=0; j<=tfn[e]-1; ++j) { comment("tf", dat); fputint(j, dat); tmem[e][j].write(dat); } } } } return; } /* input from file */ void Trifunfield::read(FILE *dat) { int cp, l, m, e; if(dat == stderr || dat == stdout) { wmmtriferror("Trifunfield: read"); } if(tmem != NULL) { for(int i=0; i<=2*(nlx+2)*(nly+2)-1; ++i) { if(tmem[i] != NULL) delete[] tmem[i]; } delete[] tfn; delete[] tmem; } nlx = fgetint(dat); nly = fgetint(dat); tfn = new int[2*(nlx+2)*(nly+2)]; tmem = new Trifun*[2*(nlx+2)*(nly+2)]; for(int i=0; i<=2*(nlx+2)*(nly+2)-1; ++i) { cp = fgetint(dat); l = fgetint(dat); m = fgetint(dat); e = cplmtoid(cp, l, m); tfn[e] = fgetint(dat); tmem[e] = new Trifun[tfn[e]]; for(int j=0; j<=tfn[e]-1; ++j) { fgetint(dat); tmem[e][j].read(dat); } } return; }