#ifndef CMT_H
#define CMT_H

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

/*
 * cmt.h
 * Procedures related to the variants of coupled mode theory 
 */

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

/*
 * hybrid CMT:
 * analytical modal basis fields, 
 * combined with FEM-discretized amplitude functions
 */


/* types of CMT basis elements: guided mode, traveling horizontally or 
vertically, forward/backward or upward/downward, localized sources, 
Gaussian beams, rectangular resonances, the clockwise and anticlockwise
propagating bend modes associated with a circular cavity, a whispering gallery 
resonance of a circular cavity, and a resonance specified by a matrix of field 
values in an external file  */
enum HcmtBft {MHF, MHB, MVU, MVD, SRC, GBM, RES, CCC, CCA, WGM, MAT, FIM};
/* -> characters */
void hcmtbftchr(HcmtBft t, char& c0, char& c1, char &c2);

/* field component identifiers, for use by the HCMT solver */
enum FKcp {fcEy, fcHx, fcHz, dzEy, dxEy, dzHxmdxHz,
           fcHy, fcEx, fcEz, dzHy, dxHy, dzExmdxEz};
Fcomp fkcp2fcomp(FKcp  cp);
FKcp  fcomp2fkcp(Fcomp cp);
/* ... for the vectorial (3-D) variant */
enum C3cp {f3Ex, f3Ey, f3Ez, f3Hx, f3Hy, f3Hz, c3Ex, c3Ey, c3Ez, c3Hx, c3Hy, c3Hz};
Fcomp c3cp2fcomp(C3cp  cp);
C3cp  fcomp2c3cp(Fcomp cp);

/* for evaluation of Hcmt integrals & mode overlaps:
   parts of guided mode profiles will be neglected where the field strength 
   is below MCROP times the mode field profile maximum */
extern double MCROP;

/* for evaluation of Hcmt integrals & mode overlaps: 
   contributions of bend mode profiles will be neglected 
   for radii below BM_R_CROP_IN and radii beyond BM_R_CROP_OUT */
extern double BM_R_CROP_IN;
extern double BM_R_CROP_OUT;

/* for evaluation of Hcmt integrals & mode overlaps:
   contributions of WGM profiles will be neglected 
   for radii below  WM_R_CROP_IN and radii beyond WM_R_CROP_OUT */
extern double WM_R_CROP_IN;
extern double WM_R_CROP_OUT;

/* vHCMT, evanescent slab modes:
   amplitude discretization will be discontinued, if the absolute 
   slab mode level drops below EvSM_DISC_CROP */
extern double EvSM_DISC_CROP;


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

/* segments of the x-z-plane, for piecewise integration <-> FEM matrices */
#define MaxNumIsegs 100
enum IsegBt {CST, CRV, STL, LPC}; // different types of boundaries
class Isegment
{
public:
	// the integration interval along the z-axis
	double za; 
	double zb;
	// lower boundary: type & parametrization
	int lbt;
	double lbcx;
	double lbcz;
	double lbr;
	double lbsg;
	double ledta;
	
	// uper boundary: type & parametrization
	int ubt;
	double ubcx;
	double ubcz;
	double ubr;
	double ubsg;
	double uedta;

	// temporary data
	double sa, sb, x0, z0, alpha, cal, sal;
	double cwx0, cwx1, cwz0, cwz1;
	
	// part of the initialization
	void set_za_zb_al(Interval cwx, Interval cwz);
	int set_LPC_slopes(double &s1, double &s2);

	// the lower and upper boundary
	double lbxpos(double z) const;
	double ubxpos(double z) const;
		
	// for purposes of inspection:
	// save segment information to a file
	void write(FILE *dat);
	// MATLAB plot commands for the segment outline 
	void plot(FILE *dat, char *s);
};

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

/* basis element for the hybrid CMT, one of the interacting fields */
class HcmtElement 
{
public:
	// initialization, generic
	HcmtElement();
	// destroy
        ~HcmtElement();
        // copy constructor
        HcmtElement(const HcmtElement& he);
        // assignment
        HcmtElement& operator=(const HcmtElement& he);

	// initialization, guided modes
	// t: MHF or MHB for modes that travel forward / backward along z 
	//    MVU or MVD for modes that travel upward / downward along x 
	// m: the mode object
	// FEM discretization of the amplitude as a function of 
	// the x- or z-coordinate in [s_0, s_n], stepsize d_s, regular
	HcmtElement(HcmtBft t, Mode m, double s_0, double s_n, double d_s);

	// initialization, guided modes
	// t: MHF or MHB for modes that travel forward / backward along z 
	//    MVU or MVD for modes that travel upward / downward along x 
	// m: the mode object
	// FEM discretization of the amplitude as a function of 
	// the x- or z-coordinate in [s_0, s_n], irregular,
	// elements of maximum length d_s_disc symmetrically positioned at each 
	// respective interface in str, in-between a regular discretization
	// with maximum stepsize d_s_hom 
	HcmtElement(HcmtBft t, Mode m, Circuit str, 
	            double s_0, double s_n, double d_s_hom, double d_s_disc);
	HcmtElement(HcmtBft t, Mode m, SegWgStruct str, 
	            double s_0, double s_n, double d_s_hom, double d_s_disc);

	// initialization, localized source of polarization p, wavelength wl,  
	// outgoing Hankel functions H_2 of order o, localized at x_0, z_0,
	// embedded in a medium with refractive index n,
	// t: SRC,
	// FEM discretization of the amplitude as a function of 
	// the radial coordinate in [s_0, s_n], stepsize d_s, regular
	HcmtElement(HcmtBft t, Polarization p, double wl,
	            int o, double x_0, double z_0, double n,
	            double s_0, double s_n, double d_s);

	// initialization, Gaussian beam with polarization p, wavelength wl, 
	// width w at the waist, beam waist localized at x_0, z_0, a beam that 
	// propagates at an angle a with respect to the pos. z-axis in a 
	// medium with refractive index n,
	// t: GBM,
	// FEM discretization of the amplitude as a function of 
	// the coordinate along the beam axis in [s_0, s_n], stepsize d_s, reg.
	HcmtElement(HcmtBft t, Polarization p, double wl,
	            double w, double x_0, double z_0, double a, double n,
	            double s_0, double s_n, double d_s);

	// initialization, bimodal resonance of a dieletric square, 
	// established by the cores of modes m1, m2, 
	// center location at x_0, z_0, 
	// a single resonance amplitude is introduced
	HcmtElement(HcmtBft t, Mode m1, Mode m2, double x_0, double z_0);

	// initialization, bend mode traveling round a circular cavity 
	// t: CCC or CCA for modes that travel clockwise / anticlockwise
	// m: the bend mode object, origin at x_0, z_0
	// FEM discretization of the amplitude as a function of 
	// the angular coordinate in [0, 2 pi], stepsize d_p, regular
	// kR, optional: prescribed approxim. angular wavenumber
	HcmtElement(HcmtBft t, BDMode m, double x_0, double z_0, double d_p);
	HcmtElement(HcmtBft t, BDMode m, double x_0, double z_0, double d_p, int kR);

	// initialization, whispering gallery mode supported by a 
	// circular cavity 
	// t: WGM
	// wgm: the corresponding whispering gallery mode object
	HcmtElement(HcmtBft t, WGMode wgm);

	// initialization, guided mode supported by a circular step-index fiber 
	// t: FIM
	// fim: the corresponding fiber mode object
	HcmtElement(HcmtBft t, FIMode fim);
	// ... in an HcmtBundle: observe amplitude evolution 
	// at ny+1 equidistant positions ya <= y <= yb 
	HcmtElement(HcmtBft t, FIMode fim, double ya, double yb, int ny);

	// initialization, resonant field with complex frequency omega, 
	// discretized on a rectangular FD grid 
	// x_0 -> (numx) -> x_n, z_0 -> (numx) -> z_n, 
	// file input, t: MAT
	HcmtElement(HcmtBft t, Polarization p, Complex omega, 
	            const char *fname,
	            double x_0, double x_n, int numx, 
	            double z_0, double z_n, int numz);

	// Polarization
	Polarization pol = TE;
	// vacuum wavenumber
	double k0 = 1.0;
	// 1/(omega*mu0)
	double invommu0 = 1.0;
	// 1/(omega*epsilon0)
	double invomep0 = 1.0;

	// type of the CMT basis element
	HcmtBft type = MHF;

	// amplitude function, constant input / output levels 
	Complex amp(double s) const;
	// begin- and end-levels
	Complex ain = CC0;
	Complex aout = CC0;
	// circular amplitude function, end elements are identified 
	Complex camp(double s) const;

	// amplitude function -> .xyf data file 
   	//    foa: MOD, SQR, REP, IMP
   	//    disp: output window on the element coordinate axis
   	//    np: number of output points
   	//    ext0, ext1: filename id characters 
	void amp2file(Afo foa, Interval disp, int np, 
	              char ext0, char ext1) const;
	// as before, with default parameters, absolute value only 
	void amp2file(char ext0, char ext1) const;
	// as before, with default parameters, REP, IMP, MOD
	void amp2files(char ext0, char ext1) const;

	// FEM discretization parameters
	double s0 = 0.0;
	double sN = 1.0;
	double ds_hom = 1.0;
	double ds_disc = 1.0;

	// field values
	Complex field(Fcomp cp, double x, double z);

	// .. for use in the HCMT solver
	Complex el_field(FKcp cp, int k, double x, double z);
	Cvector el_field(int k, double x, double z);
	// .. for use in the vectorial HCMT solver
	Complex el_field_vec(C3cp cp, int k, double x, double z);
	Cvector el_field_vec(int k, double x, double z);

	// the coefficients of the FEM expansion, set by the HCMT solver
	Cvector u;

	// FEM representation of the amplitude function
	Dvector sbeg;	
	Dvector scen;	
	Dvector send;	
	Dvector sdsL;	
	Dvector sdsR;	
	Ivector given;
	int Ns = 0;
	void initFEMreg();
	void initFEMstr(Dvector p);
	double elf(int k, double s) const;	
	double d_elf(int k, double s) const;	
	double celf(int k, double s) const;	
	double d_celf(int k, double s) const;	
	Isegment **seg = NULL;
	Ivector nseg;
	int fillintsegs(int k, Isegment *lseg, Interval cwx, Interval cwz);

	// prepare for the HCMT solver
	void prepare(Interval cwx, Interval cwz);
	// polish, upon completed solver procedure
	void polish();

	// element information -> stderr
	void info() const;

	// index of the reciprocal element, 
	// waves traveling along the same channel in the reverse direction
	int recEli = 0;
	
	/* ---------------------------------------------------------- */
	// write (most) element data to file
	void write(FILE *dat);
	
	/* ---------------------------------------------------------- */

	// element specifications, only partially relevant for specific types
	Mode  mode;
	Mode rmode;
	BDMode bmode;
	WGMode wgmode;
	FIMode fimode;
private:
	Interval mext;
	Propdir dir = FORW;
	double d_beta = 0.0;
	double d_alpha = 0.0;
	double d_kappa = 0.0;
	Complex d_gamma = CC0;
	int ord = 1;
	double ordd = 1.0;
	int orda = 1;
	double x0 = 0.0;
	double xn = 1.0;
	int numx = 1;
	double dx = 1.0;
	double z0 = 0.0;
	double zn = 1.0;
	int numz = 1;
	double dz = 1.0;
	double baB = 1.0;
	double width = 1.0;
	double alpha = 1.0, sal = 1.0, cal = 1.0;
	double nb = 1.0;
	double kn = 1.0;
	double kqnq = 1.0;
	Cmatrix FDfEy;
	Cmatrix FDfHx;
	Cmatrix FDfHz;
	Cmatrix FDfdzEy;
	Cmatrix FDfdxEy;
	Cmatrix FDfdzHxmdxHz;

	/* --- vHCMT variant ----------------------------------------------- */

public:
	// problem type, 0: scalar TE/TM, 1: vectorial 
	int vM = 0;

	// for vM = 1: wavenumber parameter
	double ky = 1.0;

	// readjust wavenumber parameter
	void setky(double gky);

	// waveguide elements, normalized slab modes, 
	// potentially non-real kz: output power 
	double pout();

	// internal helpers
	Complex d_kz = CC1;
	double v_z0 = 1.0;

	// observe amplitude evolution of FIMs in a HcmtBundle
	Cvector pramp;  // projected modal amplitude 
	Dvector prpos;  // at position y;
};

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

/* container for the full hybrid CMT model, including the solver procedure */

enum HcmtStType {REC, OVL, MIC}; // different types of structure representation 
class HcmtField
{
public:
	// the structure in question
	HcmtStType sttype;

	// sttype = REC: rectangular waveguide circuit 
	SegWgStruct str; 

	// sttype = OVL: overlay of geometrical shapes
	OvlStruct ovl; 

	// sttype = MIC: symmetric circular microcavity 
	double st_w;  // access waveguides, width
	double st_g;  // gap between access waveguides and cavity
	double st_R;  // cavity radius, outer rim
	double st_d;  // bent core of the cavity, radial width 
	double st_x0; // cavity center position
	double st_z0; //
	double st_nb; // background refractive index 
	double st_ng; // core refractive index 
	double st_eb; // ... corresponding permittivity values
	double st_eg; //

	// map sttype -> local permittivity at x, z
 	double eps(double x, double z) const;

	// polarization: TE, TM
	Polarization pol;

	// vacuum wavelength and derived quantities
	double lambda;
	double k0;
	double invommu0;
	double invomep0;

	// the computational window, 
	Interval cwx;
	Interval cwz;
	// for the evaluation of the FEM integrals: 
        // number of integr. steps per wavelength (10pt Gaussian quadrature)
	int nispwl;
	
	// number of interacting fields
	int nel;

	// initialize, generic 
	HcmtField(); 
	// initialize, sttype = REC
	HcmtField(SegWgStruct s, Polarization p, double wl, 
	          Interval Cwx, Interval Cwz, int Nispwl);
	// initialize, sttype = OVL 
	HcmtField(OvlStruct s, Polarization p, double wl, 
	          Interval Cwx, Interval Cwz, int Nispwl);
	// initialize, sttype = MIC
	HcmtField(double w, double g, double R, double d,
	          double x0, double z0, double nb, double ng,
	          Polarization p, double wl, 
	          Interval Cwx, Interval Cwz, int Nispwl);
	// destroy
	~HcmtField();
	// copy constructor
	HcmtField(const HcmtField& hc);
	// assignment
	HcmtField& operator=(const HcmtField& hc);

	// free allocated memory
	void clear();

	// subscripting
	inline HcmtElement& operator() (int k) { return evec[k]; }
	inline HcmtElement operator() (int k) const { return evec[k]; }

	// add a constituting field; 
	// return value: the access index of fe in (*this)
	int addcf(HcmtElement fe);

	// add two reciprocal constituting fields, guided modes
	// t: MHF or MHB for modes that travel forward / backward along z 
	//    MVU or MVD for modes that travel upward / downward along x 
	// m: the mode object
	// FEM discretization of the amplitude as a function of 
	// the x- or z-coordinate in [s_0, s_n], stepsize d_s, regular
	// mif, mib (output): element indices of the of the forward/upwards 
	//                    and backward/downwards traveling waves
	void addchannel(HcmtBft t, Mode m, double s_0, double s_n, double d_s,
	                int& mif, int& mib);

	// HCMT FEM approach: integration over products of elements 
	// integral kernel 
	Cvector fkernel(int l, int m, double x, double z) const;
	// ... vectorial variant
	Cvector fkernel_vec(int l, int m, double x, double z) const;
	// integrate over x, one continuous piece 
	Cvector fpintegx(int l, int m, double xa, double xb, double z) const;
	// integrate over x, piecewise over parts of the structure 
	Cvector fintegx(int l, int m, double xa, double xb, double z) const;
	// integrate over z, one continuous piece 
	Cvector fpinteg(int l, int m, 
	                const Isegment& lseg, const Isegment& mseg, 
	                double za, double zb) const;
	// integrate over z, piecewise over parts of the structure 
	Cvector finteg(int l, int m, 
	               const Isegment& lseg, const Isegment& mseg) const;

	// HCMT system, quadratic boundary terms 
	Cvector btbeval(int l, int m) const;
	// HCMT system, linear boundary terms, inhomogeneities 
	Complex btreval(int l) const;

	// the HCMT solver, complex conjugated fields, Galerkin procedure
	void solve();

	// for the efficient evaluation of parameter scans: 
	// split solver procedure
	void init();
	void setupsystem(Cmatrix& M);
	void solve(Cmatrix M);

	// the HCMT solver, non complex conjugated fields, 
	//                  variational scheme incl. boundary terms;
	//              !!! bidirectional versions of all fields required, 
	//              <-> addchannel, field references  
	void varsolve();

	
	// HCMT supermode analysis, complex conjugate fields, Galerkin procedure;
	// returns a vector of supermode eigenfrequencies 
	Cvector supermodes();

	// filter the set of supermodes: returns indices of smodes 
	// with resonant wavelengths between minwl and maxwl, and
	// with attenuation below limal 
	Ivector sm_filter(double minwl, double maxwl, double limal);

	//  setup internal representation of supermode no. m, 
	//  return value: the respective: eigenfrequency
	Complex supermode(int m);

	// the resulting field 
	Complex field(Fcomp cp, double x, double z) const;
	// values on a rectangular mesh
	Cmatrix field(Fcomp cp, Interval wx, int npx, 
	                        Interval wz, int npz) const;

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

	// overlaps between the HCMT field and ...
	// ... mode m, traveling in direction d along the z-axis, at position z
	Complex overlap_x(const Mode& m, Propdir d, double z) const;
	// ... mode m, traveling in direction d along the x-axis, at position x
	Complex overlap_z(const Mode& m, Propdir d, double x) const;
	// coarse approximation for fields with pronounced discontinuities;
	// the routines assume continuous fields with continuous derivatives

	// ----------------------------------------------------------- 
	// preparation of physical field plots:
        // selection of the global phase / snapshot time 0 within one period

        // multiply the entire field by a phase factor exp(i phi),
        // the internal representation is modified accordingly
        void adjustphase(double phi);

        // adjust the global phase, such that a plot(cp, ORG) of the physical
        // basic field component (cp = EY for TE, cp = HY for TM) exhibits the
        // maximum field at position (x, z)
        void adjustphase(double x, double z);

	// adjust the global phase, such that a plot(cp, ORG) of the physical
	// basic field component (cp = EY for TE, cp = HY for TM) exhibits an
	// overall field maximum (-> time snapshot of an extremal state
	// in case of a configuration with partially standing waves)
	// ... determine the maximum on a regular mesh of numx * mumz points
	//     on the window given by xbeg, xend, zbeg, zend 
	void adjustphase(double xbeg, double xend, double zbeg, double zend,
	                 int numx, int numz);

	// ----------------------------------------------------------- 
	// visualization: output to MATLAB .m files  
	// evaluation of the mode interference pattern in the display window 
	// xbeg <= x <= xend, zbeg <= z <=zend 
	// npx, npz: number of points in output mesh
	// ext0, ext1: filename id characters for the .m file

	// image plot corresponding to field component cp
	// cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
	// foa: MOD, SQR, REP, IMP
	void plot(Fcomp cp, Afo foa,
                  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,
                  double xbeg, double xend, double zbeg, double zend, 
	          int npx, int npz, char ext0, char ext1) const;

	// animation of the light propagation
        // the frames show images of EY (TE) resp. HY (TM), at
        // ntfr equidistant times over one time period
        void movie(double xbeg, double xend, double zbeg, double zend,
                   int npx, int npz, int ntfr, char ext0, char ext1) const;

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

private:
	// the constituting fields
	HcmtElement *evec;
	// helper functions for overlap evaluations
	Complex ovlkernel_x(const Mode& m, Propdir d, double x, double z) const;
	Complex ovlkernel_z(const Mode& m, Propdir d, double x, double z) const;
	// for use by the HCMT solver:
	int NumU, u, g;
	Ivector U_bfi;
	Ivector U_eli;
	Cvector U_val;
	Ivector U_giv;
	Ivector U_r;
	int cckernel;
	// supermodes: internal representation
	Cvector om;
	Cmatrix smp;

	// ----------------------------------------------------------- 
	// vHCMT variant
public: 
	// problem type, 0: scalar TE/TM, 1: vectorial 
	int vM;

	// global wavenumber parameter 
	double ky;

	// change global wavenumber parameter
	void setky(double gky);

	// initialize, sttype = OVL 
	HcmtField(OvlStruct s, double gky, double wl, 
	          Interval Cwx, Interval Cwz, int Nispwl);
};

/* ------------------------------------------------------------------------   
   - Gaussian bundles of vectorial HCMT solutions                         - 
   ------------------------------------------------------------------------ */

/* container for field levels at specific points */
class cProbe
{
public:
	// location 
	double x; 
	double y;
	double z;

	// id characters
	char id; 

	// the local permittivity
	double eps;
	
	// the local e.m. field  
	Complex ex;
	Complex ey;
	Complex ez;
	Complex hx;
	Complex hy;
	Complex hz;

	// cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ
	Complex field(Fcomp cp) const;
	// foa: one of MOD, ORG, SQR, REP, IMP
	double field(Fcomp cp, Afo foa) const;
	// the field "strength", id: one of mE, mH, qE, qH, mW 
	double field(FSid id) const;

	// initialize 
	cProbe();
	cProbe(double cx, double cy, double cz, char ic);
};

/* an array of probes */
class cProbeArray
{
public:
	// number of probes included
	int num;
	// initialize
	cProbeArray();
	// destroy
	~cProbeArray();
	// copy constructor 
	cProbeArray(const cProbeArray& va);
	// assignment 
	cProbeArray& operator=(const cProbeArray& s);

	// free allocated memory
	void clear();

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

	// add a probe
	void add(cProbe m);

private:
	cProbe *pvec;
};

/* types of views: 
   cXZ: axis x (vertical) versus axis z (horizontal) at position y 
   cYZ: axis y (vertical) versus axis z (horizontal) at position x
   cXY: axis x (vertical) versus axis y (horizontal) at position z */
enum cVtype {cXZ, cYZ, cXY};

/* containers for discretized fields, plot routines */
class cView
{
public:
	// the underlying permittivity distribution,
	// adapted to the type of view in case of type YZ, XY
	OvlStruct st;
	
	// discretized field information on a regular rectangular grid 
	Cmatrix ex;
	Cmatrix ey;
	Cmatrix ez;
	Cmatrix hx;
	Cmatrix hy;
	Cmatrix hz;

	// axis orientation, x, y, z 
	// <-> plot axis horizontal, vertical, position
	cVtype type;
	// axis coordinates, characters 
	char va_chr;
	char ha_chr;

	// plot: vertical axis: va_N+1 points from va_beg to va_end
	double va_beg; 
	double va_end; 
	int    va_N;

	// plot: horizontal axis: ha_N+1 points from ha_beg to ha_end
	double ha_beg; 
	double ha_end; 
	int    ha_N;      

	// third axis: position
	double pos;

	// transform v, h, pos to x, y, z coordinates
	void trf(double v, double h, double& x, double& y, double& z) const; 

	// plot id characters
	char pid0; 
	char pid1; 
	
	// initialize: supply discretization information, setup matrices 
	cView();
	cView(cVtype t, OvlStruct s,
	     double v_beg, double v_end, int v_N, 
	     double h_beg, double h_end, int h_N, double p,
	     char ext0, char ext1);

	// ---------------------------------------------------------------
	// visualization: output to MATLAB .m files 
	// ext0, ext1: filename id characters for the m.file
	
	// image plot corresponding to field component cp
	// cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
	// foa: ORG, MOD, SQR, REP, IMP (ORG = REP)
	void plot(Fcomp cp, Afo foa) const; 

	// image plot of the field "strength", id: one of mE, mH, qE, qH, mW 
	void plot(FSid id) const;

	// electromagnetic energy density image
	void plot() const;

	// animation of the light propagation
	// the frames show images of component cp
	// at ntfr equidistant times over one time period
	void movie(Fcomp cp, int ntfr) const;

	// export full field data ("everything") into a viewer MATLAB file
	void viewer() const;
};

/* an array of containers for discretized fields */
class cViewArray 
{
public:
	// number of views included
	int num;
	// initialize
	cViewArray();
	// destroy
	~cViewArray();
	// copy constructor 
	cViewArray(const cViewArray& va);
	// assignment 
	cViewArray& operator=(const cViewArray& s);

	// free allocated memory
	void clear();

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

	// add a view
	void add(cView m);

private:
	cView *vvec;
};

/* registration of modal output power */
class cPort
{
public:
	// initialize	
	cPort();
	cPort(int e, char c1, char c2);
	// index of the HcmtElement of which the output is assembled 
	int eli;

	// identifier
	char id1;
	char id2;
	// upon completion of vHcmtBundle::solve():
	// Pout is the guided modal power carried across border bdr
	// by the mode of polarization pol, order ord
	double Pout;
};

/* an array of containers for registration of guided modal output power */
class cPortArray 
{
public:
	// number of ports included
	int num;

	// initialize
	cPortArray();
	// destroy
	~cPortArray();
	// copy constructor 
	cPortArray(const cPortArray& pa);
	// assignment 
	cPortArray& operator=(const cPortArray& s);

        // free allocated memory
        void clear();


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

	// add a port
	void add(cPort m);

private:
	cPort *pvec;

};

/* Bundles of vectorial HCMT solutions, container */

class vHcmtBundle 
{
public:
	// the structure under investigation
	OvlStruct str;
	
	// incoming mode, element index
	int inel;
	// effective index of incoming mode
	double Nin;
	// vacuum wavenumber
	double k0;
	
	// primary angle of incidence /^o
	double theta0;
	// beam half-width (1/e, field) at focus point /mum
	double bhw;
	// beam focus coordinates
	double y0;
	double z0;

	// lateral spectral discretization of the beam
	double ky0;  // primary lateral wavenumber
	int    Nky;  // number of discrete angles
	double Dky;  // angular window 
	double Wky;  // spectral half-width

	// resolve a resonance in the numerical integration
	double kyR0;  // in the lateral wavenumber interval
	double kyR1;
	int    NkyR;  // with the number of discrete angles

	// resolve output amplitude resoances 
	Ivector eliR; // in elements in eliR
	double  dkyR; // down to this wavenumber difference 

	// vHCMT computational window
	Interval cwx;
	Interval cwz;
	// number of vHCMT integration steps per wavelength
	int nispwl;

	// stepwise vHCMT: fill views with field data 
	void solve();
	// ... proceed from ka to kb in (roughly) nk steps 
	void solve(double ka, double kb, int nk);
	// ... resolve amplitude resonances 
	void solveR();

	// full setup:
	//    pre-filled field object f;
	// excitation:  
	//    a guided mode, HcmtElement eli, 
	//    coming in with unit amplitude, 
	//    at primary angle of incidence t (^o), 
	//    beam of full width w (1/e, field),
	//    represented by na angles between t-dt and t+dt,
	//    dt determined such that spectral contributions with relative 
	//    amplitude below dtlev are neglected; 
	//    positioned at focus coordinates fcy, fcz 
	vHcmtBundle(HcmtField f,
               int eli,
               double t, double w, double fcy, double fcz, int na, double dtlev);
	// ...
	// numerical integration needs to resolve a resonance 
	//    in angular range [tr-dtr, tr+dtr] 
	//    with nr equidisttant steps
	vHcmtBundle(HcmtField f,
               int eli,
               double t, double w, double fcy, double fcz, int na, double dtlev, 
	       double tr, double dtr, int nr);
	// ...
	// numerical integration needs to resolve resonances in 
	// the output amplitudes of elements with indices in reli, 
	// down to an angular stepsize dtr 
	vHcmtBundle(HcmtField f,
               int eli,
               double t, double w, double fcy, double fcz, int na, double dtlev, 
	       Ivector reli, double dtr);

	// containers for field inspection
	cViewArray VA;
	// containers for field levels at specific positions 
	cProbeArray FA;
	// containers for the registration of guided modal output power
	cPortArray PA;
	// registration of the guided modal input power
	double Pin;

	// a single lateral spectral component:
	double ky;
	HcmtField fld;
	Complex spfld(Fcomp cp, double x, double y, double z);
};

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

/* 
 * conventional coupled mode theory, 
 * codirectional guided wave propagation along parallel waveguide cores
 */

/* setup coupling matrices (called by CMTsetup) */
void CMTmats(ModeArray& umode,     // modes to be coupled
             Waveguide cpwg,       // the entire coupler structure
             Dmatrix& cmtS,        // output: power coupling matrix
             Dmatrix& cmtBK);      // output: coupling coefficients

/* setup coupling matrices, calculate supermode propagation constants
   and normalized amplitude vectors */
void CMTsetup(ModeArray& umode,     // modes to be coupled
	      Waveguide cpwg,       // the entire coupler structure
	      Dmatrix& sigma,       // output: power coupling matrix
	      Dvector& smpc,        // output: supermode propagation constants
	      Dmatrix& smav);       // output: supermode amplitude vectors 

/* coupled mode transfer matrix */
Cmatrix CMTtrafm(double len,     // length of the coupling segment
                                 // as output from CMTsetup:     
                 Dmatrix& sigma, // power coupling matrix
	         Dvector& smpc,  // supermode propagation constants
	         Dmatrix& smav); // supermode amplitude vectors 

#endif // CMT_H