#ifndef WMMTRIF_H
#define WMMTRIF_H

/*
 * WMM
 * Wave-matching method for mode analysis of dielectric waveguides
 * Manfred Lohmeyer 
 * University of Osnabrueck, Department of Physics
 * Barbarastrasse 7, D-49069 Osnabrueck, Germany
 * (1999)
 * Manfred Hammer
 * University of Twente, Faculty of Mathematical Sciences
 * P.O. Box 217, 7500AE Enschede, The Netherlands
 * (2002)
 */

/*
 * wmmtrif.h
 * the basic solutions to the Helmholtz wave equation
 */


/* Fundamental solutions of the wave equation */
enum Fstype { FS_A = 0,     // cos cos
	      FS_B = 1,     // cos sin
	      FS_C = 2,     // sin cos
	      FS_D = 3,     // sin sin
	      FS_E = 4,     // exp exp
	      FS_F = 5,     // exp cos
	      FS_G = 6,     // exp sin
	      FS_H = 7,     // cos exp
	      FS_I = 8 };   // sin exp

/* trial functions:  parametrization */
class Trifun
{
public:
        int Amsi;         // start index in amplitude matrix
	double Gfcof;     // normalization
	Fstype GfTyp;     // type of trial function
	double AFp;       // wave vector paramter alpha
	char PfQf;        // p-q-modus
	double Pj;        // wavevector component p
	int Pvz;          // sign
	double Qj;        // wavevector component q
	int Qvz;          // sign 
	double Xr;        // coordinate offset x0
	double Yr;        // coordinate offset y0
	char Svz;         // sign due to symmetry

	double Amplit;    // amplitude in finally assembled mode 
		
	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat);

	// partial derivatives
	Fstype DGfTyp(Derlev d) const { return Fstype(dgftypes(int(d))); }; 
	// ... and corresponding coefficients 
	double DGfcof(Derlev d) const { return dgfcoeff(int(d)); }; 

	// initialize dependent parameters
	void initdeppars(double beta, double eps, double k0, Polarization p,
			 Rect r);

	// value of field/derivative d at position x, y
	double val(Derlev d, double x, double y) const;

	// for factorizing trial functions:
	// values of x- and y-part or derivative d at position x or y
	double xval(Derlev d, double x) const;
	double yval(Derlev d, double y) const;

	// 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);

	// translate the trial function by (dx, dy)  
	void translate(double dx, double dy);

private:
	Ivector dgftypes;
	Dvector dgfcoeff;
};

/* ------------------------------------------------------------------------ */

/* integration of trial functions along lines */
double tfintx(const Trifun& tf1, Derlev d1, const Trifun& tf2, Derlev d2,
	      double y, double xa, double xb);
double tfinty(const Trifun& tf1, Derlev d1, const Trifun& tf2, Derlev d2, 
	      double x, double ya, double yb);

/* integration of trial functions over rectangles */
double tfnrm(const Trifun& tf1, Derlev d1, const Trifun& tf2, Derlev d2,
             Rect r);

/* initialize tables for trial function integration */
void inittfinttables();

class Tfintegrate
{
public:
	// initialize
	Tfintegrate();
};

extern Tfintegrate TFINTINIT;

/* ------------------------------------------------------------------------ */

/* trial functions for a rectangular waveguide */
class Trifunfield 
{
public:
	// initialize
	Trifunfield();
	Trifunfield(int nx, int ny);
	// destroy
	~Trifunfield();
	// copy constructor
	Trifunfield(const Trifunfield& t);
	// assignment
	Trifunfield& operator=(const Trifunfield& t);

	// input
	void read(FILE *dat);
	// output
	void write(FILE *dat);

	// insert a trial function set
	void insert(Trifun *s, int nf, int cp, int l, int m); 

	// subscripting: -> Trifun* 
	Trifun* operator() 
		      (const int& cp, const int& l, const int& m, const int& j)
		      const
		      { return &(tmem[cplmtoid(cp, l, m)][j]); }
	// number of trial functions in cell l, m, component cp
	int ntf(const int& cp, int& l, int& m) const
	       { return tfn[cplmtoid(cp, l, m)]; }

	// remove trial functions with Amsi si
	void removetrifun(int si);

	// translate the trial function field by (dx, dy)  
	void translate(double dx, double dy);

private:
	int nlx;
	int nly;
	int *tfn;
	Trifun* (*tmem);
	int cplmtoid(const int& cp, const int& l, const int& m) const
		    { return cp+l*2+m*(nlx+2)*2; }
};

#endif  // WMMTRIF_H