#ifndef WMMARR_H
#define WMMARR_H

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

/*
 * wmmarr.h
 * arrays of WMM modes, mode interference evaluation
 */

/* mode pointer */
typedef  WMM_Mode*  WMM_ModeP;

/* array of WMM_Modes */
#define MAXNUMMODES 100   // maximum number of modes in array

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

	// free allocated memory
	void clear();

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

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

	// add a mode
	void add(WMM_Mode m);
	// delete a mode entry
	void remove(int i);
	// add an entire WMM_ModeArray nma, nma is cleared ! 
	void merge(WMM_ModeArray& 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 fieldmat(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 fieldvec(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 WMM_Mode& imode, int m, const WMM_Mode& omode);
   
	// single relative power level for device length l    
	double iopower(const WMM_Mode& imode, const WMM_Mode& omode, 
                       Cvector pert, double l);

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

	// ... write to file     
	void writeiopower(const WMM_Mode& imode, const WMM_Mode& 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 mfile(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 fancymfile(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 WMM_ModeArray& imodes, Waveguide iwg,
                    const WMM_ModeArray& omodes, Waveguide 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 WMM_Mode& imode, 
                    const WMM_ModeArray& omodes, Waveguide 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 WMM_Mode& imode, 
                    const WMM_Mode& 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:
	WMM_ModeP mvec[MAXNUMMODES];
};

#endif // WMMARR_H