#ifndef VEIMS_H
#define VEIMS_H

/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * http://metric.computational-photonics.eu/ 
 */

/*
 * veims.h
 * Variational Effective Index Mode Solver VEIMS 
 * mode analysis for waveguides with 2D rectangular cross section  
 * --- the simplest version, separabale mode profile approximations
 */

/* a single guided mode */
class EIMode 
{
public:
	// the supporting waveguide structure, alternative representations
	WgCrs wg;
	SegWgCrs st;
	// polarization: TE, TM
	Polarization pol;
	// propagation constant 
	double beta;
	// vacuum wavenumber
	double k0;
	// effective modal index 
	double neff;
	
	// minimum / maximum field amplitudes
	double minE;
	double maxE;
	double minH;
	double maxH;
	double minS;
	double maxS;
	
	// mode profile, field values at position (x, y) 
	// cp: EX - HZ, SZ
	double field(Fcomp cp, 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 rectangle r
 	// cp1, cp2 != Sz
	double recint(Fcomp cp1, Fcomp cp2, Rect r) const;

	// 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 power() 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 field(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 field(Rect disp, int npx, int npy, Fcomp cp, Afo foa) 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 plot(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 acplot(Rect disp, int nx, int ny, char ext0, char ext1) 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 secplot(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 fplot(Fcomp cp, Rect disp, int npx, int npy, 
                   char ext0, char ext1) const;

	// export full mode profile data ("all") into a viewer m-file
	// 	disp: the plot window
	// 	npx, npy: number of plot points in the x and y directions
	// 	ext0, ext1: filename id characters 
	void viewer(Rect disp, int npx, int npy, char ext0, char ext1) const;

	// permittivity perturbations: propagation constant shift  
	Complex phaseshift(ChWgPert 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)
	void translate(double dx, double dy);

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

	// the VEIM mode profile approximation:
	// index of the reference segment 
	int rsegi;
	// vertical mode profile 
	Mode vmp;
	// horizontal effective refractive index profile 
	Waveguide hwg;
	Dvector hwg_q;
	Dvector hwg_ec;
	// horizontal mode profile 
	Mode hmp;

	// VEIMS mode profile representation: separable fields 
	void fldrep(Fcomp cp, double &fac, FldorDer &vfod, FldorDer &hfod, double x, double y) const;
	void fldrep(Fcomp cp, double &fac, FldorDer &vfod, FldorDer &hfod, int l, int m) const;

	// indication of lateral leakage: 
	// 0: guided mode
	// 1: mode must be suspected to be leaky 
	int leaky;

private:

	// evaluation of mode profile integrals: helper functions
	double prodint(int l, int m, Fcomp cp1, Fcomp cp2, Rect r) const;
	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;

};

/* an array of guided modes */
class EIModeArray 
{
public:
	// number of modes included
	int num;
	// initialize
	EIModeArray();
	// destroy
	~EIModeArray();
	// copy constructor 
	EIModeArray(const EIModeArray& ma);
	// assignment 
	EIModeArray& operator=(const EIModeArray& s);

	// free allocated memory
	void clear();

	// subscripting
	EIMode& operator() (int i); 
	EIMode operator() (int i) const; 

	// add a mode
	void add(EIMode m);
	// delete a mode entry
	void remove(int i);
	// add an entire EIModeArray nma 
	void merge(EIModeArray& nma);

	// sort the array by propagation constants, highest first
	void sort();
	
	// -----------------------------------------------------------
	// mode interference evaluation 
	// all modes are assumed to belong to the same waveguide ! 

	// field superposition at point (x, y, z), 
	//    amp:  complex amplitudes at z=0
	//          amplitudes evolve according to 
        //          amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
	//    pert: propagation constant perturbations, 
	//          cbet(j) = this(j).beta+pert(j)
	//    cp:   EX, EY, EZ, HX, HY, HZ, SZ 
	Complex field(Cvector amp, Cvector pert, Fcomp cp, 
                      double x, double y, double z) const;

	// evaluate component cp on a rectangular npx x npy mesh on area disp 
	// at position z
	//    amp: complex amplitudes at z=0
	//         amplitudes evolve according to 
	//         amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
	//    pert: propagation constant perturbations, 
	//          cbet(j) = this(j).beta+pert(j)
	//    cp: EX - HZ, SZ 
	//    foa: MOD, SQR, REP, IMP 
	Dmatrix field(Cvector amp, Cvector pert, double z,
                      Rect disp, int npx, int npy, Fcomp cp, Afo foa) const;

	// store num values of component cp at position z
	// between (x0, y0) and (x1, y1) in a vector 
	//    amp: complex amplitudes at z=0
	// 	   amplitudes evolve according to 
	//         amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
	//    pert: propagation constant perturbations, 
	//          cbet(j) = this(j).beta+pert(j)
	//    cp: EX - HZ, SZ 
	//    foa: MOD, SQR, REP, IMP */
	Dvector field(Cvector amp, Cvector pert, double z, 
	              Fcomp cp, Afo foa, int nump, 
	              double x0, double y0, double x1, double y1) const;

	// ----------------------------------------------------------- 
	// Three segment coupler, power transfer evaluation 
	// imode: input mode, excites the array's modes 
	//        with rel. power 1 at z=0 
	// omode: output mode, relative power is returned 
	// pert: propagation constant perturbations, 
	//       cbet(j) = this(j).beta+pert(j) 
	// --- this(j) are assumed to be normalized ! --- 

	// weights for the power evaluation:
	// w = ( omode | this(m) ) ( this(m) | imode ) 
	double pweight(const EIMode& imode, int m, const EIMode& omode);
   
	// single relative power level for device length l    
	double iopower(const EIMode& imode, const EIMode& omode, 
                       Cvector pert, double l);

	// output power for numl devices of lengths between lmin and lmax     
	Dvector iopower(const EIMode& imode, const EIMode& omode, 
                        Cvector pert, int numl, double lbeg, double lend);

	// ... write to file     
	void writeiopower(const EIMode& imode, const EIMode& omode, 
                          Cvector pert, int numl, double lbeg, double lend,
                          char ext0, char ext1);

	// ----------------------------------------------------------- 
	// Visualization: output to MATLAB .m files  

	// write single component of the interference field at position z 
	// to MATLAB .m file  
	//    amp: complex amplitudes at z=0
	//         amplitudes evolve according to 
	//         amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
	//         (adjust maxf in the output file for |amp| >= 1 )
	//    pert: propagation constant perturbations, 
	//          cbet(j) = this(j).beta+pert(j)
	//    cp: EX - HZ, SZ
	//    foa: MOD, SQR, REP, IMP  
	//    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': field + mesh only, no plot commands (default) 
	void plot(Cvector amp, Cvector pert, double z, 
                  Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
	          char ext0, char ext1, char pltype) const;

	// write interference field at position z to MATLAB .m file, 
	// fancy style :-)  
	//    amp: complex amplitudes at z=0
	//         amplitudes evolve according to 
	//         amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
	//         (adjust maxf in the output file for |amp| >= 1 )
	//    pert: propagation constant perturbations, 
        //          cbet(j) = this(j).beta+pert(j)
	//    disp: output region on the x-y-plane
	//    npx, npy: number of points in output mesh
	//    ext0, ext1: filename id characters */
	void fplot(Cvector amp, Cvector pert, double z,
                   Rect disp, int npx, int npy, 
                   char ext0, char ext1) const; 

	// animate the interference field  
	//    amp: complex amplitudes at z=0
	//         amplitudes evolve according to 
        //         amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
	//    pert: propagation constant perturbations, 
	//          cbet(j) = this(j).beta+pert(j)
	//    cp: EX - HZ, SZ
	//    foa: MOD, ORG, SQR, REP, IMP 
	//    disp: output region on the x-y-plane
	//    npx, npy: number of points in output mesh
	//    pltype: 'C': contour plot 
	//            'S': surface plot
	//            'I': intensity image 
	//            'F': fancy plot, cp and foa are set to SZ, MOD (default)
	//    z0, z1, numz: .m files are generated for z-positions 
	//                   z0+i*(z1-z0)/numz, i=0, ..., numz-1 
	void movie(Cvector amp, Cvector pert,  
                   Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
		   char pltype,
                   int numz, double z0, double z1) const; 

	// interference pattern on the horizontal plane 
	// ybeg <= y <= yend, zbeg <= z <= zend at elevation x
	// image plot corresponding to the squareroot 
        // of the local intensity (SZ)
	// amp: complex amplitudes at z=0
	//      amplitudes evolve according to 
	//      amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
	// pert: propagation constant perturbations, 
	//       cbet(j) = this(j).beta+pert(j)
	// npy, npz: number of points in output mesh
	// ext0, ext1: filename id characters 
	void prop(Cvector amp, Cvector pert, double x,
                  double ybeg, double yend, int npy,
                  double zbeg, double zend, int npz,
                  char ext0, char ext1) const;

	// Three segment coupler, interference pattern on the horizontal plane 
	// ybeg <= y <= yend, zbeg <= z <= zend at elevation x
	// image plot corresponding to the squareroot 
	// of the local intensity (SZ)
	//     z < 0: intensity profile of the input modes
	// 0 < z < l: mode interference pattern
	//     l < z: intensity profiles of the output modes, 
	//            scaled by the output amplitudes
	// (!) this(j) are assumed to be normalized.
	//    iamp: complex input mode amplitudes
	//    imodes: input modes, excite the modes this(j) at z=0, 
       	//            belonging to well separated port waveguides 
	//    iwg: the geometry of the combined input waveguides
	//    omodes: output modes, belonging to well separated port waveguides 
	//    owg: the geometry of the combined output waveguides
	//         iwg and owg are for displaying purposes only
	//    pert: propagation constant perturbations, 
	//          cbet(j) = this(j).beta+pert(j)
	//    npy, npz: number of points in output mesh
	//    ext0, ext1: filename id characters 
	void ioprop(Cvector iamp, 
                    const EIModeArray& imodes, WgCrs iwg,
                    const EIModeArray& omodes, WgCrs owg,
                    Cvector pert, double x, double l,
                    double ybeg, double yend, int npy,
                    double zbeg, double zend, int npz,
                    char ext0, char ext1) const;

	// ... as before, but with a single input mode imode only 
	void ioprop(const EIMode& imode, 
                    const EIModeArray& omodes, WgCrs owg,
                    Cvector pert, double x, double l,
                    double ybeg, double yend, int npy,
                    double zbeg, double zend, int npz,
                    char ext0, char ext1) const;

	// ... as before, with single input and output modes imode and omode
	void ioprop(const EIMode& imode, 
                    const EIMode& omode, 
                    Cvector pert, double x, double l,
                    double ybeg, double yend, int npy,
                    double zbeg, double zend, int npz,
                    char ext0, char ext1) const;

private:
	EIMode *mvec;
};

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

/* integration of a fieldproduct between two modes
   over an arbitrary rectangle */
double twomrecint(const EIMode& mode1, Fcomp cp1,
                  const EIMode& 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 EIMode& m1, const EIMode& 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 EIMode& m1, const EIMode& m2);

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

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

/* the VEIMS mode solver, complete guided mode analysis, 
   returns the number of found modes
   wg, st: the structure under consideration
   rsn: reference slice number; -1: determ. automatically 
   pol: polarization type, TE or TM
   ma: the modes found (output)
   quiet == 1: suppress log output */
int veims(SegWgCrs wg, Polarization pol, int rsn, EIModeArray& mam, int quiet);
// ... with automatically determined reference slice
int veims(SegWgCrs st, Polarization pol, EIModeArray& ma, int quiet);
int veims(SegWgCrs st, Polarization pol, EIModeArray& ma);
int veims(WgCrs wg, Polarization pol, EIModeArray& ma, int quiet);
int veims(WgCrs wg, Polarization pol, EIModeArray& ma);

/* VEIMS algorithm: solve the effective 1D problem */
int veims_hsolve(Waveguide hwg, Dvector hwg_ec, ModeArray& ma);

#endif // VEIMS_H