#ifndef SLAMODE_H
#define SLAMODE_H

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

/*
 * slamode.h
 * Modes of multilayer slab waveguides, mode overlaps
 */

/* field evaluation:
   neglect contributions with amplitudes below NEGLECT_FLD_CONTRIB */
extern double NEGLECT_FLD_CONTRIB;

enum FldBhv {HYP, OSC}; 

/* internal mode representation: field on a homogeneous layer */
class SmPiece
{
public:
	// Polarization
	Polarization pol;
	// the layer interfaces 
	double bif; // bottom
	double tif; // top
	// vacuum wavenumber 
	double k0;
	// two elementary amplitudes
	double A;
	double B;
	// "transverse" wave number
	double g;
	// local coordinate offsets
	double xa;
	double xb;
	// field behaviour: HYP, OSC
	enum FldBhv bhv; 
	// the local permittivity times the vacuum wavenumber 
	double k0nq;
	// 1/(n*n) 
	double inve;
	
	// initiate: set polarization, interfaces, permittivity
	void init(Polarization p, Interval i, double e, double lambda);
	// ... VEIMS nonstandard effective TM profile 
	void init(Interval i, double e, double ec, double lambda);

	// mode solver, transverse resonance evaluation:
	// start,
	void start_tre(double bq, double& phi, double &phis, 
	               Boundary_type bdt);
	// continue,
	void cont_tre(double bq, double& phi, double &phis);
	// finish.
	double finish_tre(double bq, double& phi, double &phis, 
	                  Boundary_type bdt);
	// set parameters in an upper segment for a modal solution 
	void polish_ub(Boundary_type bdt);

	// mode solver, reverse transverse resonance evaluation:
	// start,
	void start_revtre(double bq, double& phi, double &phis, 
	                  Boundary_type bdt);
	// continue,
	void cont_revtre(double bq, double& phi, double &phis);
	// finish.
	double finish_revtre(double bq, double& phi, double &phis, 
	                     Boundary_type bdt);
	// set parameters in a lower segment for a modal solution 
	void polish_lb(Boundary_type bdt);

	// number of nodes in the principal field component
	int nodes(double x0, double x1) const;
	// number of nodes in the derivative of the principal field component
	int dnodes(double x0, double x1) const;
	// maximum absolute level of the principal field component
	double maxf(Interval i) const;

	// field evaluation within a specific layer:
	// principal component
	double bfld(double x) const;
	// the derivative of the principal component with respect to x
	double bder(double x) const;

	// integrate the squared principal field component over the layer 
	double intbfq() const;
	// integrate the squared principal field component over interval x0, x1 
	double intbfq(double x0, double x1) const;

	// integrate the present principal field component 
	// times the principal component of modepiece mp over interval iv
	double intff(const SmPiece& mp, const Interval& iv) const;
	// integrate the present principal field component times the derivative 
	// of the principal component of modepiece mp over interval iv
	double intfd(const SmPiece& mp, const Interval& iv) const;
	// integrate the derivative of the present principal field component 
	// times the principal component of modepiece mp over interval iv
	double intdf(const SmPiece& mp, const Interval& iv) const;
	// integrate the derivative of the present principal field component 
	// times the deivative of the b. c. of modepiece mp over interval iv
	double intdd(const SmPiece& mp, const Interval& iv) const;

	// integrate the principal field component times a phase term
	Complex crintfCiS(double ce, double xr, 
	                          const Interval& iv) const;
	// integrate the derivative of the b. f. c. times a phase term
	Complex crintdCiS(double ce, double xr, 
	                          const Interval& iv) const;
	// integrate the principal field component times an exponential
	Complex crintfExp(double ce, double xr, 
	                          const Interval& iv) const;
	// integrate the derivative of the b. f. c. times an exponential
	Complex crintdExp(double ce, double xr, 
	                          const Interval& iv) const;

	// integrate the principal field component times a complex exponential 
	// (exp -i kr (x-xr) over interval iv 
	Complex crintfcE(Complex kr, double xr, const Interval& iv) const;

	// integrate the derivative o. t. pr. f. c. times a complex exponential 
	// (exp -i kr (x-xr) over interval iv 
	Complex crintdcE(Complex kr, double xr, const Interval& iv) const;

	// translate layer representation by dx
	void translate(double dx);

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

/* a slab mode */
class Mode
{
	friend class ModeArray;
public:
	// initialize
	Mode();
	Mode(const Waveguide& g,
             Boundary_type bdtb, double bpb, Boundary_type bdtt, double bpt,
             Polarization p);
	// VEIMS: nonstandard TM effective horizontal profile 
	Mode(const Waveguide& g, const Dvector& ec, 
             Boundary_type bdtb, double bpb, Boundary_type bdtt, double bpt);
	// destroy 
	~Mode();
	// copy constructor
	Mode(const Mode& m);
	// assignment
	Mode& operator=(const Mode& m);

	// the corresponding waveguide  
	Waveguide wg;
	// Polarization
	Polarization pol;
	// vacuum wavenumber
	double k0;
	// 1/(omega*mu0)
	double invommu0;
	// 1/(omega*epsilon0)
	double invomep0;

	// propagation constant
	double beta;
	// square of the propagation constant
	double betaq;
	// effective mode index 
	double neff;
	// normalized propagation constant 
	double npcB;

	// propagation type: PROPAG, EVANESC
	Propagation_type ppt;	

	// boundary type: OPEN, LIMD, LIMN
	Boundary_type bdt_b; // bottom
	Boundary_type bdt_t; // top
	// for bdt == LIMD, LIMN: boundary positions 
	double bp_b; // bottom
	double bp_t; // top 

	// mode order 
	int ord;
	// determine the number of nodes in the principal component 
	int nodes() const;
	// determine the number of nodes in the derivative of the principal comp. 
	int dnodes() const;
	// classify: determine and set mode order if possible,
        //           otherwise: set the special switch 
	void classify();
	// the token of the mode: 'TE0', 'TM5', 'TE110x' 
	char ids[8];
	// set this string
	void setids();
	// ... for given order 
	void setids(int o);

	// principal field component, maximum absolute level 
	double maxF;
	// determine that value 
	void setmaxF();

	// for guided modes only: extent of the mode profile; 
	// outside the returned interval the field strength 
	// (absolute value, principal component) is lower than fac * maxF; 
	// typically fac << 1.0
	Interval extent(double fac) const;

	// complex propagation constant:
	//    beta (real)      for ppt == PROPAG 
	// -i*beta (imaginary) for ppt == EVANESC
	// *   1.0  for d == FORW
	// * (-1.0) for d == BACK
	// invoke as cbeta(Propdir d)
	Cvector cbeta;

	// mode profile, 
	// the principal field component at transverse coordinate x
	double field(double x) const;
	// and its derivative
	double d_field(double x) const;

	// the mode profile, 
	// values of component cp: EX - HZ at transverse coordinate x
	// Note the convention for real mode profile components! 
	// ... for the forward traveling mode
	double field(Fcomp cp, double x) const;
	// ... for forward or backward traveling modes
	double field(Fcomp cp, Propdir d, double x) const;

	// the complex mode profile,
	// values of component cp: EX - HZ at transverse coordinate x,
	// for propagation direction d: FORW, BACK
	Complex cfield(Fcomp cp, Propdir d, double x) const;
	// forward and backward field
	void cfield(Fcomp cp, double x, Complex& ff, Complex& bf) const;

	// if viewed as a mode propagating along x: complex vert. mode profiles 
	Complex vertfield(Fcomp cp, Propdir d, double z) const;      

	// relevant values for the HCMT procedures
	// void cfield(Propdir d, double x, Cvector &fv) const;
	// void vertfield(Propdir d, double z, Cvector &fv) const;

	// factor of phase/amplitude change along a propagation distance z
	// d: FORW, BACK selects the forward / backward traveling field
	Complex efac(Propdir d, double z) const;
	// forward and backward values, for distances dzf and dzb
	void efac(double dzf, Complex& ef, double dzb, Complex& eb) const;
	// with the propagation constant shifted by pert
	Complex efac(Propdir d, Complex pert, double z) const;

	// mode field evolution along z, given a unit amplitude at z=0
	// d: FORW, BACK
	Complex cfield(Fcomp cp, Propdir d, double x, double z) const;

	// check the validity of the Maxwell-equations at a number of points 
	void checkmode() const;
	// for vTE / vTM variants
	void cVcheckmode() const;
	// check the continuity of the principal field components 
	// returns a maximum error quantity
	double checkcontinuity() const;

	// principal component, integrate a product of field or 
	// derivatives over interval i 
	double integrate(FldorDer fod1, FldorDer fod2, Interval i) const;

	// integration of a fieldproduct over interval i, 
        // cp1, cp2 != Sz
	double integrate(Fcomp cp1, Fcomp cp2, Interval i) const;
	// for bidirectional mode versions, 
        // d1, d2 = FORW, BACK,  
	// cc1, cc2 == 1: take the complex conjugate of the component 
	Complex integrate(Fcomp cp1, Propdir d1, int cc1,  
                          Fcomp cp2, Propdir d2, int cc2, Interval i) const;  

	// longitudinal component of the Poyntingvector, 
	// integrated over the entire cross section, ppt EVANESC: normalization 
	double power() const;
	
	// normalize mode to power() == 1 
	void normalize();
	// ... to power() == p 
	void normalize(double p);

	// fieldproduct with another mode, integral over interval i, 
	// principal components or their derivatives 
	double integrate(FldorDer fod1, const Mode& sm2, FldorDer fod2, 
	                 Interval i) const;

	// integrate a fieldproduct with another mode over interval i, 
	// real version 
	double integrate(Fcomp cp1, const Mode& sm2, Fcomp cp2, 
	                 Interval i) const;

	// ... integrate directional versions of a product with another mode; 
        // cc1, cc2 == 1: take the complex conjugate of the component 
	Complex integrate(Propdir d1, Fcomp c1, int cc1,
                          const Mode& m2, Propdir d2, Fcomp c2, int cc2,
                          Interval i) const;
	// the product  
        // 0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
  	// of the present object with another mode, directional variants
	Complex overlap(Propdir d1, const Mode& m2, Propdir d2) const; 

	// the product  
        // 0.25\int(E_1x H_2y - E_1y H_2x + E_2x H_1y - E_2y H_1x)dx 
	// (non complex conjugate fields)
  	// of the present object with another mode, directional variants
	Complex overlap_nc(Propdir d1, const Mode& m2, Propdir d2) const; 

	// the product  
        // 0.5\int(E_1x H_2y - E_1y H_2x)dx 
	// (non complex conjugate fields)
  	// of the present object with another mode, directional variants
	Complex overlap_ncns(Propdir d1, const Mode& m2, Propdir d2) const; 

	// for the QUEP implementation:
	// overlaps of perpendiculary propagating modes & helper functions 
	Complex cross_int(Fcomp cp, int l,
	                  Complex ceexp, double xoffs, Interval i) const;
	Complex cross_int(Fcomp cp, 
	                  Complex ceexp, double xoffs, Interval i) const;

	// for the QUEP implementation: 
	// overlaps of perpendiculary propagating modes 
	Complex crossoverlap(Propdir hdir,
	                     const Mode& mv, Propdir vdir, 
	                     double zpos, Interval i);
	void crossoverlap(const Mode& mv, double zpos, Interval i,
	                  Complex& cff, Complex& cfb,
	                  Complex& cbf, Complex& cbb);

	// product
        // 0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx
	// of a mode with a Gaussian beam gb, numerical evaluation of the 
        // integrals, alternatively for forward or backward propagating fields
	Complex overlap(const Gaussianbeam& gb, Propdir pd) const;

        // ... helper function: numerical integration by Gaussian quadrature
	Complex cnintegrate(Fcomp cpm, const Gaussianbeam& gb, Fcomp cpb,
                            Interval lyr, int numx) const;

	// field profile -> .xyf data file  
	// 	cp: EX - HZ
	// 	foa: MOD, ORG, SQR, REP, IMP  
	//      disp: output window on the x-axis
	//      ext0, ext1: filename id characters
	//      np: number of output points
	void writeprofile(Fcomp cp, Afo foa, 
	                  Interval disp, int np, char ext0, char ext1) const;
	// as before, with default parameters
	void writeprofile(char ext0, char ext1) const;

	// field profile -> MATLAB .m file  
	//	cp: EX - HZ
	//	foa: MOD, ORG, SQR, REP, IMP 
	//	disp: output window on the x axis
	//	np: number of output points 
	//	ext0, ext1: filename id characters
	//      pltype: 'L': geometry information & plot commands (default)
	//              'V': field, mesh, and plot comand, 
        //                   to be included into *L.m 
	void plot(Fcomp cp, Afo foa, Interval disp, int np, 
		  char ext0, char ext1, char pltype) const;
	// as before, with default parameters
	void plot(char ext0, char ext1, char pltype) const;
	void plot(char ext0, char ext1) const;

	// small permittivity perturbation: first order propagation constant
	// shift due to perturbation p, for propagation direction dir 
	Complex phaseshift(Perturbation p, Propdir dir);

	// small interface displacement: first order propagation constant shift
	// due to a shift of the dielectric interface hx(b) in the modes 
	// waveguide by a distance dx, for propagation direction dir 
	Complex phaseshift(int b, double dx, Propdir dir);

	// small displacement of a boundary position:
	// first order propagation constant shift due to a shift dx of the
	// boundary position at the bottom (p=='b') or at the top (p == 't')
	// of the computational window, for propagation direction dir, 
	// so far evaluated as a finite difference expression only ...
	Complex phaseshift_bdp(char p, double dx, Propdir dir);

	// small alteration of the vacuum wavelength: first order propagation
	// constant shift due to a shift of the vacuum wavelength by dl, 
	// for the version of the mode with propagation direction dir 
	Complex phaseshift(double dl, Propdir dir);
	//... the corresponding expression for the alteration of the 
	// effective mode index due to a vaccum wavelength shift dl
	Complex neffshift(double dl, Propdir dir);

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

	// output to file dat
	void write(FILE *dat) const;
	// output to file with default name
	void write(char ext0, char ext1);
	// input from file dat
	void read(FILE *dat);
	// input from file with default name
	void read(Polarization p, char ext0, char ext1);


	// compare the present object to m, return 1 if equal, 0 otherwise
	// both modes are assumed to be properly normalized eigenfunctions
	int equal(const Mode& m) const;

	/* - mode solver routines ---------------------------------- */

	// transverse resonance evaluation 
        double travres(double tbq); 

	// reverse transverse resonance evaluation 
        double revtravres(double tbq); 

	// node counting: determine the number of modes (returned) with squared
	// propagation constants larger than the present betaq, the return 
	// value is equal to the order+1 of the highest order mode above the 
	// present betaq; meant to be applied to an unfinished mode object
	int nummodesabove(double tbq);

	// ... inspect these quantities & exit; for diagnosis purposes. 
	void travresinspect(double bqmin, double bqmax);

	// localize a zero in travres between bql and bqu, up to accuracy lim 
	double localize(double bql, double bqu, double lim); 
	
	// set the relevant variables for an actual modal solution
	void polish(double bq);

	// expand the present mode object with LIMN or LIMD condition 
	// at the upper boundary bp_t towards a symmetric (LIMN) or 
	// antisymmetric (LIMD) mode on the doubled waveguide
	void expand();

	// mirror the present mode object with OPEN condition at the
	// upper boundary with respect to the boundary position p 
	void mirror(double p);

	// if nonzero: mode corresponds to a part of a decoupled structure,
	//             or a proper classification is not possible 
	int special;
	
	/* ---------------------------------------------------------- */

	// the internal mode representation
	int N;	
	SmPiece *piece;

private:
	// field evaluation, map: Fcomp -> field
	Ivector zero;
	Dvector fac;
	Ivector sfd;
	Ivector dnq;
	Ivector img;
	Ivector brv;
	Ivector evi;

public:
	// determine these quantities 
	void setfieldmap();

	/* ---------------------------------------------------------- */
	// for the vQUEP routines: extension to vTE, vTM modes 

	// vM = 1: (*this) is meant as a vTE, vTM mode
	int vM;
	
	// wave vector components in y- and z-direction
	Complex ky;
	Complex kz;

	// prepare vTE / vTM mode representation, for given ky
	void vectorify(Complex val_ky); 

	// complex directional wavenumber, z-axis
	// kz *   1.0  for d == FORW
	//    * (-1.0) for d == BACK
	// invoke as ckz(Propdir d)
	Cvector ckz;

	// vTE/vTM modes, complex mode profile
	// values of component cp: EX - HZ at transverse coordinate x,
	// for propagation direction d: FORW, BACK
	Complex cVfield(Fcomp cp, Propdir d, double x) const;

	// if viewed as a mode propagating along x: complex vert. mode profiles 
	Complex vertVfield(Fcomp cp, Propdir d, double z) const;     

	// integrate a fieldproduct of two vTE/ vTM modes (FORW) 
	// over Interval i, complex conjugation of the first factors 
	Complex cVintegrate(Fcomp cp, const Mode& m, Fcomp cpm, 
                            Interval i) const;
	// check: numerical evaluation of the integral 
	Complex num_cVintegrate(Fcomp cp, const Mode& m, Fcomp cpm, 
	                        Interval i) const;

	// the product  
        // 0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
  	// of the present vM object with another mode, directional variants
	Complex cVoverlap(Propdir d1, const Mode& m2, Propdir d2) const; 

	// overlaps of perpendiculary propagating vTE and vTM modes 
	// & helper functions 
	Complex cVcross_int(Fcomp cp,
                            Complex ceexp, double xoffs, Interval i) const;
	Complex cVcrossoverlap(Propdir hdir, const Mode& mv, Propdir vdir,
                               double zpos, Interval i) const;
	void cVcrossoverlap(const Mode& mv, double zpos, Interval i,
                            Complex& cff, Complex& cfb,
                            Complex& cbf, Complex& cbb) const;
 

private:
	// field evaluation, map: Fcomp -> cVfield
	Ivector cVzero;
	Ivector cVsfd;
	Cvector cVfac;
	Ivector cVdnq;
public:
	Ivector cVbrv;
	Propagation_type cVppt;	
};

/* --- mode overlaps ------------------------------------------------- */

/* integrate a fieldproduct of two modes over [i.x0, i.x1], real version */
double integrate(const Mode& sm1, Fcomp cp1,
                 const Mode& sm2, Fcomp cp2, Interval i);

/* ... integrate directional versions of a product of two mode profiles;
   cc1, cc2 == 1: take the complex conjugate of the component */
Complex integrate(const Mode& m1, Propdir d1, Fcomp c1, int cc1,
                  const Mode& m2, Propdir d2, Fcomp c2, int cc2, Interval i);

/* product 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   of two modes, alternatively forward or backward propagating */
Complex overlap(const Mode& m1, Propdir d1, const Mode& m2, Propdir d2); 
/* ... for vTE / vTM modes */
Complex cVoverlap(const Mode& m1, Propdir d1, const Mode& m2, Propdir d2); 

/* product 
   0.25\int(E_1x H_2y - E_1y H_2x + E_2x H_1y - E_2y H_1x)dx 
   (non complex conjugate fields)
   of two modes, alternatively forward or backward propagating */
Complex overlap_nc(const Mode& m1, Propdir d1, const Mode& m2, Propdir d2); 

/* product 
   0.5\int(E_1x H_2y - E_1y H_2x)dx 
   (non complex conjugate fields)
   of two modes, alternatively forward or backward propagating */
Complex overlap_ncns(const Mode& m1, Propdir d1, const Mode& m2, Propdir d2); 

/* for the BEP implementation: products 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   of two modes, combinations of forward and backward propagating versions */
void bidir_overlaps(const Mode& m1, const Mode& m2, 
                    Complex& off, Complex& ofb, Complex& obf, Complex& obb);

/* for the vBEP / vQUEP implementation: products 
   0.25\int(E_1x^* H_2y - E_1y^* H_2x + E_2x H_1y^* - E_2y H_1x^*)dx 
   of two modes, combinations of forward and backward propagating versions */
void cVbidir_overlaps(const Mode& m1, const Mode& m2,
                      Complex& off, Complex& ofb, 
                      Complex& obf, Complex& obb);


#endif  // SLAMODE_H