/*
 * 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<stdio.h>
#include<stdlib.h>
#include<math.h>
#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)<PQTolCC)
	{
		if(p1 == p2)
			s1 = 0.25*(x1+x1-x0-x0)*cos(0.5*(p2+p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2+p1)*(r1-r2); 
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s1 = 0.25*(b-a)*cos(pr)
			     -0.0625*(p2-p1)*(b*b-a*a)*sin(pr);
		}
	}
	else 
	{
		a = 0.5*(p2*(x0-r2)-p1*(x0-r1));
		b = 0.5*(p2*(x1-r2)-p1*(x1-r1));
		s1 = cos(b+a)*sin(b-a)/(p2-p1);
	}
	if(fabs(p1+p2)<PQTolCC)
	{
		if(p1 == -p2) 
			s2 = 0.25*(x1+x1-x0-x0)*cos(0.5*(p2-p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2-p1)*(r1-r2); 
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s2 = 0.25*(b-a)*cos(pr)
			     -0.0625*(p2+p1)*(b*b-a*a)*sin(pr);
		}
	}
	else 
	{
		a = 0.5*(p1*(x0-r1)+p2*(x0-r2));
		b = 0.5*(p1*(x1-r1)+p2*(x1-r2));
		s2 = cos(b+a)*sin(b-a)/(p2+p1);
	}
	return s1+s2;
}

inline double intfCS(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)<PQTolCS)
	{
		if(p1 == p2)
			s1 = 0.25*(x1+x1-x0-x0)*sin(0.5*(p1+p2)*(r1-r2));
		else
		{
			pr = 0.5*(p1+p2)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s1 = 0.25*(b-a)*sin(pr)
			     +0.0625*(p2-p1)*(b*b-a*a)*cos(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)-p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)-p1*(x0-r1));
		s1 = -sin(a+b)*sin(a-b)/(p2-p1);
	}
	if(fabs(p1+p2)<PQTolCS)
	{
		if(p1 == -p2)
			s2 = 0.25*(x1+x1-x0-x0)*sin(0.5*(p2-p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2-p1)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s2 = 0.25*(b-a)*sin(pr)
			     +0.0625*(p2+p1)*(b*b-a*a)*cos(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)+p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)+p1*(x0-r1));
		s2 = -sin(a+b)*sin(a-b)/(p2+p1);
	}
	return s1+s2;
}

inline double intfSC(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)<PQTolCS)
	{
		if(p1 == p2)
			s1 = -0.25*(x1+x1-x0-x0)*sin(0.5*(p1+p2)*(r1-r2));
		else
		{
			pr = 0.5*(p1+p2)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s1 = -0.25*(b-a)*sin(pr)
			     +0.0625*(p2-p1)*(b*b-a*a)*cos(pr);
		}	
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)-p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)-p1*(x0-r1));
		s1 = sin(a+b)*sin(a-b)/(p2-p1);
	}
	if(fabs(p1+p2)<PQTolCS)
	{
		if(p1 == -p2)
			s2 = 0.25*(x1+x1-x0-x0)*sin(0.5*(p2-p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2-p1)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s2 = 0.25*(b-a)*sin(pr)
			     +0.0625*(p2+p1)*(b*b-a*a)*cos(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)+p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)+p1*(x0-r1));
		s2 = -sin(a+b)*sin(a-b)/(p2+p1);
	}
	return s1+s2;
}

inline double intfSS(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)<PQTolSS)
	{
		if(p1 == p2)
			s1 = 0.25*(x1+x1-x0-x0)*cos(0.5*(p1+p2)*(r1-r2));
		else
		{
			pr = 0.5*(p1+p2)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s1 = 0.25*(b-a)*cos(pr)
			     -0.0625*(p2-p1)*(b*b-a*a)*sin(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)-p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)-p1*(x0-r1));
		s1 = cos(b+a)*sin(b-a)/(p2-p1);
	}
	if(fabs(p1+p2)<PQTolSS)
	{
		if(p1 == -p2)
			s2 = -0.25*(x1+x1-x0-x0)*cos(0.5*(p2-p1)*(r1-r2));
		else
		{
			pr = 0.5*(p2-p1)*(r1-r2);
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			s2 = -0.25*(b-a)*cos(pr)
			     -0.0625*(p2+p1)*(b*b-a*a)*sin(pr);
		}
	}
	else 
	{
		b = 0.5*(p2*(x1-r2)+p1*(x1-r1));
		a = 0.5*(p2*(x0-r2)+p1*(x0-r1));
		s2 = -cos(b+a)*sin(b-a)/(p2+p1);
	}
	return s1+s2;
}

inline double intfEE(const double& p1, const double& r1,  
		     const double& p2, const double& r2,
		     const double& x0, const double& x1)
{
	double a, b, s1, s2;
	if(fabs(p1+p2)<PQTolEE)
	{
		if(p1 == -p2)
			return 0.5*(x1+x1-x0-x0)*exp(0.5*(p2-p1)*(r1-r2)); 
		else
		{
			a = x0+x0-r1-r2;
			b = x1+x1-r1-r2;
			return 0.5*(b-a+0.25*(p2+p1)*(b*b-a*a))
				  *exp(0.5*(p2-p1)*(r1-r2)); 
		}
	}
	else 
	{
		if(fabs(x1)>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;
}