#ifndef CPLXWG_H
#define CPLXWG_H

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

/* 
 * cplxwg.h
 * Modes of waveguides with complex permittivity
 */

#include <stdio.h>

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

/* waveguide geometry + permittivity profile; vacuum wavelength ..., 
   complex media, including losses or gain */
class cWaveguide
{
public:
	// number of inner layers
	int nx;
	// positions of the dielectric interfaces
	Dvector hx;

	// get interval index corresponding to position x  
	int layeridx(double x) const;
	// get layer boundaries corresponding to index l  
	Interval layer(int l) const;
	// get layer boundaries corresponding to position x  
	Interval layer(double x) const;

	// vacuum wavelength 
	double lambda;
	// refractive index profile
	Cvector n;
	// permittivity on layer l
	Complex eps(int l) const;
	// permittivity at position x
	Complex eps(double x) const;

	// translate: hx -> hx+dx
	void translate(double dx);

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

	// initialize
	cWaveguide();
	cWaveguide(int vnx); 
	cWaveguide(int vnx, Dvector vhx, double l, Cvector n);
	cWaveguide(Waveguide wg); 

	// free allocated memory
	void strip();

	// compare the present object to wg, return 1 if equal, 0 otherwise
	int equal(const cWaveguide& wg) const;

	// remove unnecessary boundaries from waveguide
	void optimize();

	// check consistency of the waveguide object: error&exit, if not ok
	void consistency() const;

	// check the symmetry of the waveguide, return value l 
        // l == 0: no symmetry 
        // l >= 1: symmetric structure with respect to layer l
        int checksymmetry() const;

        // split the present waveguide into two, at layer sl 
        // the pieces are stored in wg1 and wg2
        void split(int sl, cWaveguide& wg1, cWaveguide& wg2) const;

        // expand the waveguide towards a structure 
        // that is symmetric with respect to xs; xs > hx(nx) 
        cWaveguide expand(double xs);

        // mirror the waveguide at xs; xs > hx(nx) 
        cWaveguide mirror(double xs);

	// distance in permittivity between this and twg, 0: waveguides are equal
	double distance(const cWaveguide twg) const;

	// an intermediate between this and twg, at ratio r
	cWaveguide intermediate(const cWaveguide twg, double r) const;
	
	// refractive index profile -> MATLAB .m file  
	// disp: output region on the x axis; ext0, ext1: filename id chars 
	void plot(Interval disp, char ext0, char ext1) const;
	void plot(char ext0, char ext1) const;

	// the real waveguide, with permittivity Re(eps), special (!)
	Waveguide re() const;

	// check, whether this is a (real) Waveguide
	bool isreal() const;
	
	// bounds for effective mode permittivities, default values
	double defaultepseffmin() const;
	double defaultepseffmax() const;

	// modify the outer regions for guidance, if possible
	void unleak(double epsmin);
};

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

/* internal mode representation: field on a homogeneous layer */
class cSmPiece
{
public:
	// Polarization
	Polarization pol;
	// the layer interfaces 
	double bif; // bottom
	double tif; // top
	// vacuum wavenumber 
	double k0;
	// two elementary amplitudes
	Complex A;
	Complex B;
	// "transverse" wave number
	Complex kappa;
	// local coordinate offsets
	double xa;
	double xb;
	// the local permittivity
	Complex eps;
	
	// initiate: set polarization, interfaces, permittivity
	void init(Polarization p, Interval i, Complex e, double lambda);

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

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

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

	// maximum absolute level of the principal field component within a layer
	double maxf(Interval i, double &ph) const;

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

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

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

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

	// 1: reverse sign of root of kappaq, 0: don't
	int krootrev;
	// translate longitudinal to lateral wavenumbers
	void gqtokappa(Complex gq);
	// on boundary pieces, boundary conditions realized:
	// 0: bound, 1: leaky
	int bdktype;
};

/* a mode of a complex waveguide */
class cMode
{
	friend class cModeArray;
public:
	// initialize
	cMode();
	cMode(const cWaveguide& g,
             Boundary_type bdtb, double bpb, Boundary_type bdtt, double bpt,
             Polarization p);
	cMode(const Mode m);
	// destroy 
	~cMode();
	// copy constructor
	cMode(const cMode& m);
	// assignment
	cMode& operator=(const cMode& m);

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

	// squared propagation constant
	Complex gammaq;
	// propagation constant, gamma = beta - i alpha
	Complex gamma;
	// effective mode index 
	Complex neff;
	// phase propagation constant
	double beta;
	// attenuation constant
	double alpha;
	// propagation length, Lp = 1/(2*alpha)
	double Lp();

	// mode order
	int ord;

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

	// 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];
	// exterior field type identifier  
	char btids[8];
	// set these strings
	void setids();
	// ... for given order 
	void setids(int o);

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

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

	// the mode profile, 
	// values of component cp: EX - HZ at transverse coordinate x
	// ... for the forward traveling mode
	Complex field(Fcomp cp, double x) const;
	// ... for forward or backward traveling modes
	Complex field(Fcomp cp, Propdir d, double x) 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;

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

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

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

	// integrate a fieldproduct cp1^* cp2 over the interval [i.x0, i.x1]
	Complex integrate(Fcomp cp1, Fcomp cp2, Interval i) const;
	// integrate a fieldproduct of two modes over Interval i,
	// use the complex conjugate of *this' contribution 
	Complex integrate(Fcomp cp, const cMode& m, Fcomp cpm, Interval i) const;

	// determine, whether the absolute squared field can be integrated over region r
	// r = 'a': the entire x-axis,
	// r = 'b': lowermost segment of the waveguide,
	// r = 't': uppermost segment of the waveguide
	bool isintegrable(char r) 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) const;

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

	// output to file dat
	void write(FILE *dat);
	// 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 cMode& m) const;

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

	// transverse resonance evaluation 
        Complex travres(Complex tgq); 

	// reverse transverse resonance evaluation 
        Complex revtravres(Complex tgq); 

	// set the relevant variables for an actual modal solution
	void polish(Complex gq);

	// follow propagation constants into the complex plane 
	int traceto(Complex egq, int kr_b, int kr_t, int quiet, double &ftrE);

	// 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;	
	cSmPiece *piece;

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

	// comptational steps that led to this mode
	// -1: undetermined, unfinished mode object
	// 0: real mode, direct real mode analysis
	// 1: complexificaton of a real mode
	// 2: complexification and leakification of a real mode
	int cst;
	
	// computation by complexification of the real mode
	Mode orgM;
};

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

/* arrays of Modes, mode interference evaluation */
class cModeArray 
{
public:
	// number of modes included
	int num;
	// initialize
	cModeArray();
	// destroy
	~cModeArray();
	// copy constructor 
	cModeArray(const cModeArray& ma);
	// assignment 
	cModeArray& operator=(const cModeArray& s);

	// free allocated memory
	void clear();

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

	// subscripting
	inline cMode& operator() (int i) { return mvec[i]; }
	inline cMode operator() (int i) const { return mvec[i]; }

	// add a mode
	void add(cMode m);
	// delete a mode entry
	void remove(int i);
	// add an entire cModeArray nma, nma is cleared ! 
	void merge(cModeArray& nma);

	// sort the array by squared propagation constants, highest first
	void sort();

	// remove doublets from the array
	void prune();

	// compare the present object to ma, return 1 if equal, 0 otherwise
	int equal(cModeArray& ma);

	// -----------------------------------------------------------
	// mode interference evaluation, bidirectional mode superpositions,
	// all modes are assumed to belong to the same waveguide ! 
	//
	// field and energy density at a point (x, z), 
	// amplitudes of forward and backward propgating mode variants must 
	// be supplied, the interpretation depends on whether a z-interval 
	// is present as an optional argument: 
	// argument list ends with x, z:
	//    famp: complex amplitudes of the forward traveling modes at z=0,
	//    bamp: complex amplitudes of the backward traveling modes at z=0,
	//    (note that z can be negative);
	// argument list ends with x, z, z0, z1:
	//    famp: complex amplitudes of the forward traveling modes at z=z0,
	//    bamp: complex amplitudes of the backward traveling modes at z=z1.
	    	
	// the local field, 
	// cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, or SZ 
	Complex field(const Cvector& famp, const Cvector& bamp, 
	              Fcomp cp, double x, double z) const;
	Complex field(const Cvector& famp, const Cvector& bamp, 
	              Fcomp cp, double x, double z, double z0, double z1) const;

	// local field "strength",  id: one of mE, mH, qE, qH, mW 
	double lfs(const Cvector& famp, const Cvector& bamp,
		   FSid id, double x, double z) const;
	double lfs(const Cvector& famp, const Cvector& bamp,
		   FSid id, double x, double z, double z0, double z1) const;

	// local energy density 
	double endens(const Cvector& famp, const Cvector& bamp, 
	              double x, double z) const;
	double endens(const Cvector& famp, const Cvector& bamp, 
	              double x, double z, double z0, double z1) const;

	// the six electric and magnetic components at x, z
        void emfield(const Cvector& famp, const Cvector& bamp,
	             double x, double z,
                     Complex& ex, Complex& ey, Complex& ez,
                     Complex& hx, Complex& hy, Complex& hz) const;
        void emfield(const Cvector& famp, const Cvector& bamp,
	             double x, double z, double z0, double z1,
                     Complex& ex, Complex& ey, Complex& ez,
                     Complex& hx, Complex& hy, Complex& hz) const;

	
	// the optical field, discretized on a regular rectangular mesh,  
	// npx x npz points on a display window dwx x dwz
	// cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
	// foa: ORG, MOD, SQR, REP, IMP (ORG = REP) 
	Cmatrix field(const Cvector& famp, const Cvector& bamp, 
		      Fcomp cp, Interval dwx, int npx, 
				Interval dwz, int npz) const;
	Dmatrix field(const Cvector& famp, const Cvector& bamp, 
		      Fcomp cp, Afo foa, Interval dwx, int npx, 
					 Interval dwz, int npz) const;

	// ----------------------------------------------------------- 
	// visualization: output to MATLAB .m files  
	// evaluation of the mode interference pattern in the display window 
	// xbeg <= x <= xend, zbeg <= z <=zend 
	// amp_f: amplitudes of the forward propagating modes at z=0 
	// amp_b: amplitudes of the backward propagating modes at z=0
	// npx, npz: number of points in output mesh
	// ext0, ext1: filename id characters for the .m file
        // assumption: uniform polarization and wavelength of all included modes

	// image plot corresponding to field component cp
	// cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
	// foa: MOD, SQR, REP, IMP
	void plot(const Cvector& amp_f, const Cvector& amp_b, 
                  Fcomp cp, Afo foa,
                  double xbeg, double xend, double zbeg, double zend, 
	          int npx, int npz, char ext0, char ext1) const;

	// relative phase, image plot 
	void phasemap(const Cvector& amp_f, const Cvector& amp_b, 
                      Fcomp cp, 
                      double xbeg, double xend, double zbeg, double zend, 
	              int npx, int npz, char ext0, char ext1) const;

	// image plot of the field "strength", id: one of mE, mH, qE, qH, mW  
	void plot(FSid id, const Cvector& amp_f, const Cvector& amp_b, 
                  double xbeg, double xend, double zbeg, double zend, 
	          int npx, int npz, char ext0, char ext1) const;

	// electromagnetic energy density image 
	void plot(const Cvector& amp_f, const Cvector& amp_b, 
                  double xbeg, double xend, double zbeg, double zend, 
	          int npx, int npz, char ext0, char ext1) const;

	// export full field data ("everything") into a viewer MATLAB file
        void viewer(const Cvector& amp_f, const Cvector& amp_b,
                    double xbeg, double xend, double zbeg, double zend,
                    int npx, int npz, char ext0, char ext1) const;

	// ----------------------------------------------------------- 
private:
	cMode *mvec;
};

/* - Mode analysis of complex waveguides -------------------------------- */

/* complexification & leakyfication: 
   maximum number of steps in waveguide change */
extern double CPLXs_WgSteps;

/* bound mode analysis: follow propagation constants into the complex plane */
cModeArray complexify(const cMode o, const cWaveguide twg, int nst, int quiet);
cModeArray complexify(const cMode o, const cWaveguide twg, int quiet);

/* leaky mode analysis: follow propagation constants into the complex plane */
cModeArray leakify(const cMode o, const cWaveguide twg, int nst, int quiet);
cModeArray leakify(const cMode o, const cWaveguide twg, int quiet);

/* tolerance for fixing the propagation constants */
extern double CPLXs_GqTol;
double gqTol(Complex gq);

/* minimum distance at which propagation constants are considered different */
extern double CPLXs_GqSep;
double gqSep(Complex gq);

/* tolerance for the maximum relative error in the continuity of the 
   basic field component for accepting a field as a valid mode */
extern double CPLXs_DiscErrTol;

/* potential SPsolutions: increase in search range for eps_eff */
extern double CPLXs_EeInc;

/* --- mode analysis, helper functions, realized problem --------------- */

/* trap propagation constants */
int trapCR(ModeArray& ma, Mode& m,
          double bqmin, double bqmax, double sf, int sm, int lfsc,
          int verb);

/* --- mode analysis --------------------------------------------------- */

/* determine modes of the complex waveguide, 
   starting with real squared propagation constants between qqmin and qqmax, 
   for the problem specified by wg, pol, 
   boundaries of type bdtb at bpb, of type bdtt at bpt; 
   verb: 0: no control output, 1: progress, 2: all */ 
void findmodesC(cWaveguide wg, Polarization pol,
                Boundary_type bdtb, double bpb,
                Boundary_type bdtt, double bpt,
                double bqmin, double bqmax, cModeArray& ma, int verb);

/* --- mode analysis, simpler driver routines -------------------------- */

/* guided mode analysis, returns the number of found modes 
   wg: the structure under consideration
   pol: polarization type, TE or TM
   ma: the modes found (output)
   quiet == 1: suppress log output */
int modeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma, int quiet); 
int modeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma); 
/* leaky mode analysis 
   epsout: permittivity applied to outer layers for initial bound mode computation */
int leakymodeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma, double epsout, int quiet);
int leakymodeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma, int quiet);
int leakymodeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma, double epsout);
int leakymodeanalysis(cWaveguide wg, Polarization pol, cModeArray& ma);

#endif  // CPLXWG_H