#ifndef MODE_H
#define MODE_H

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

/*
 * mode.h
 * general mode definitions 
 */

/* general mode */
class Mode
{
public:
	// the corresponding waveguide  
	Waveguide wg;
	// Polarization
	Polarization pol;
	// Symmetry
	Symmetry sym;
	// propagation constant 
	double beta;
	// vacuum wavenumber
	double k0;
	// effective modal index 
	double neff;
	// normalized propagation constant 
	double npcB;

	// minimum / maximum field amplitudes
	double minE;
	double maxE;
	double minH;
	double maxH;
	double minS;
	double maxS;

	// global (normalization) factor
	double ampfac;

	// field values at point (x, y) 
	// 	Fcomp: EX - HZ, SZ
	virtual double field(Fcomp cp, double x, double y) const = 0;

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

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

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

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

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

	// longitudinal component of the Poyntingvector, 
	// integrated over the entire x-y-domain 
	// TE part 
	double tepower() const;
	// TM part 
	double tmpower() const;
	// total power 
	double totpower() const;

	// polarization character of a hybrid mode 
	// mod: 'E': (int E_y^2 dxdy)/(int E_x^2+E_y^2 dxdy)
	//      'H': (int H_x^2 dxdy)/(int H_x^2+H_y^2 dxdy)
	//      'S': (int S_z^TE dxdy)/(int S_z dxdy)
	double polcharacter(char mod) const;

	// set maximum values for E, H, Sz,  only a coarse approximation !
	void setfieldmax();

	// coordinates of the maximum mode intensity on rectangle r, 
	// average 1/e^2-radius in x- and y-direction 
	// only a coarse approximation ! 
	double maxintensity(Rect r, double& xm, double& ym, 
				    double& rx, double& ry) const;

	// normalize mode to totpower() == nrm
	void normalize(double nrm);

	// store num values of component cp 
	// between (x0, y0) and (x1, y1) in a vector 
	//   cp: EX - HZ, SZ 
	Dvector fieldvec(Fcomp cp, int num, 
                       double x0, double y0, double x1, double y1) const;

	// evaluate component cp on a rectangular npx x npy mesh on area disp 
	//    cp: EX - HZ, SZ 
	//    foa: ORG, MOD, SQR 
	Dmatrix fieldmat(Rect disp, int npx, int npy, Fcomp cp, Afo foa) const;

	// write field profile  
	// 	cp: EX - HZ, SZ
	// 	foa: MOD, ORG, SQR  
	//      fac: output factor (typically 1.0 or -1.0)
	//      disp: output region on the x-y-plane
	//      ext0, ext1: filename id characters
	//      npx, npy: number of points in output mesh
	//      adj: 1: adjust contour levels to cp
	//      adj: 0: adjust contour levels to minE, maxE, minH ...
	void writeprofile(Fcomp cp, Afo foa, double fac, Rect disp,
			  char ext0, char ext1, int npx, int npy, int adj) 
                          const;

	// write single component to MATLAB .m file  
	//	cp: EX - HZ, SZ
	//	foa: MOD, ORG, SQR  
	//	disp: output region on the x-y-plane
	//	npx, npy: number of points in output mesh
	//	ext0, ext1: filename id characters
	//	pltype: 'C': contour plot 
	//	        'S': surface plot
	//	        'I': intensity image 
	//	        'N': waveguide, field + mesh only, 
        //                   no plot commands (default)
	void mfile(Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
		   char ext0, char ext1, char pltype) const;

	// write all components to MATLAB .m file, surface plots  
	//	ext0, ext1: filename id characters
	//	disp: output region on the x-y-plane
	//	npx, npy: number of points in output mesh 
	void acmfile(char ext0, char ext1, Rect disp, int nx, int ny) const;

	// write profile section along line (x0,y0) -> (x1,y1)
	// 	cp: EX - HZ, SZ
	//      ext0, ext1: filename id characters
	//      np: number of output points 
	void writesec(Fcomp cp, char ext0, char ext1, int np,
		      double x0, double y0, double x1, double y1) const;

	// write profile section along line (x0,y0) -> (x1,y1) 
        // to MATLAB .m file
	// 	cp: EX - HZ, SZ
	//      ext0, ext1: filename id characters
	//      np: number of output points 
	//      pltype: 'L': geometry information + plot commands (default)
	//              'V': field, mesh, and plot comand, 
        //                   to be included into *L.m 
	// --- at present implemented for horizontal or vertical lines only ---
	void secmfile(Fcomp cp, char ext0, char ext1, int np, char pltype,
		      double x0, double y0, double x1, double y1) const;

	// write single component to MATLAB .m file, fancy style  :-)
	//	cp: EX - HZ, SZ
	//	disp: output region on the x-y-plane
	//	npx, npy: number of points in output mesh
	//	ext0, ext1: filename id characters
	void fancymfile(Fcomp cp, Rect disp, int npx, int npy, 
                        char ext0, char ext1) const;

	// permittivity perturbations: propagation constant shift  
	Complex phaseshift(Perturbation p) const;

	// geometry variations: propagation constant shift due to
	// moving the horizontal boundary between rectangles [l,m] and [l+1,m]
	// ( moving hx(l) for hy(m-1) < y < hy(m) )
	double horgeovar(int l, int m) const;

	// geometry variations: propagation constant shift due to
	// moving the vertical boundary between rectangles [l,m] and [l,m+1]
	// ( moving hy(m) for hx(l-1) < x < hx(l) )
	double vergeovar(int l, int m) const;

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

	// output to file with default name
	void write_def(char ext0, char ext1);
	// input from file with default name
	void read_def(Polarization p, Symmetry s, char ext0, char ext1);

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

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

        // overlap 0.5*intint mode_EX*inif_HY - mode_EY*inif_HX dxdy 
	// between mode and an initial field, with real HX and HY components, 
	// integrals are evaluated locally on wg.rectangles by gaussian
	// quadrature formulas, numx*numy times applied 
        double ovlp(double (*inif)(Fcomp, double, double),       
		    Rect r, int numx, int numy) const;

	// destroy
	virtual ~Mode(){}

private:
	double ovleval(double (*inif)(Fcomp cp, double x, double y), 
		       double x, double y) const;
        double ovlpix(double (*inif)(Fcomp cp, double x, double y), 
                      double x0, double x1, double y, int numx) const;
	double ovlpint(double (*inif)(Fcomp cp, double x, double y), Rect r,
		       int numx, int numy) const;
};

#endif  // MODE_H