#ifndef WMMFLD_H
#define WMMFLD_H

/*
 * WMM - Wave-matching method for mode analysis of dielectric waveguides
 * https://wmm.computational-photonics.eu/
 */

/*
 * wmmfld.h
 * WMM_Mode - procedures for evaluating the mode profile
 */

/* WMM-computed waveguide mode */ 
class WMM_Mode : public Mode
{
public:
	// field values at point (x, y) 
	// 	Fcomp: EX - HZ, SZ
	double field(Fcomp cp, double x, double y) const;

	// product of field values at point (x,y);  cp1, cp2 != Sz
	double fieldprod(Fcomp cp1, Fcomp cp2, double x, double y) const; 

	// integration of a fieldproduct along the line (x,ya)->(x,yb),
	// fields on rectangle l,m are evaluated;  cp1, cp2 != Sz
	double horlineint(int l, int m, Fcomp cp1, Fcomp cp2, 
				  double x, double ya, double yb) const;

	// integration of a fieldproduct along the line (xa,y)->(xb,y),
	// fields on rectangle l,m are evaluated;  cp1, cp2 != Sz
	double verlineint(int l, int m, Fcomp cp1, Fcomp cp2, 
				  double xa, double xb, double y) const;

	// integration of a fieldproduct over the entire rectangle l, m
 	// cp1, cp2 != Sz
	double recintlm(int l, int m, Fcomp cp1, Fcomp cp2) const;

	// integration of a fieldproduct over an arbitrary rectangle
 	// cp1, cp2 != Sz
	double recint(Fcomp cp1, Fcomp cp2, Rect r) const;

	// translate mode by (dx,dy)
	void translate(double dx, double dy);

	// output to file dat, representation dependent part
	void writerdpid(FILE *dat);
	void writerdp(FILE *dat);
	// input from file dat, representation dependent part
	void readrdpid(FILE *dat);
	void readrdp(FILE *dat);

	// vectorial formulation 
	Vecform vform;

	// trial functions 
	Trifunfield phi;

	// destroy
        ~WMM_Mode(){}

private:
	double fieldsum(int cp, Derlev abl, int l, int m, double x, double y)
			const ;
	double lint(int dir, Derlev fod0, int cp0, Derlev fod1, int cp1, 
		    int l, int m, double a, double b, double c) const;
	double lineint(int dir, int l, int m, Fcomp cp1, Fcomp cp2,
		       double a, double b, double c) const;
	double pint(Derlev fod0, int cp0, Derlev fod1, int cp1, int l, int m,
		    Rect r) const; 
	double prodint(int l, int m, Fcomp cp1, Fcomp cp2, Rect r) const;
};

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

/* integration of a fieldproduct between two modes 
   over an arbitrary rectangle */
double twomrecint(const WMM_Mode& mode1, Fcomp cp1, 
                  const WMM_Mode& mode2, Fcomp cp2, Rect r);

/* product 0.5\int\int(E_1x H_2y - E_1y H_2x)dxdy of two modes */
double scalprod(const WMM_Mode& m1, const WMM_Mode& m2); 

/* symmetrical product 
   0.25\int\int(E_1x H_2y - E_1y H_2x + E_2x H_1y - E_2y H_1x)dxdy 
   of two modes */
double lscalprod(const WMM_Mode& m1, const WMM_Mode& m2);

/* for a superposition of two WMM_Modes m1, m2 with coefficients c1, c2:
   intensity behind a polarizer at an angle alpha versus TE-polarization */
double polarizeroutput(Complex c1, const WMM_Mode& m1, 
                       Complex c2, const WMM_Mode& m2, 
		       double alpha);

#endif // WMMFLD_H