#ifndef WMMCMT_H
#define WMMCMT_H

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

/*
 * wmmcmt.h
 * WMM Coupled Mode Theory 
 */

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

/* Coupled mode theory:
   setup coupling matrices, calculate supermode propagation constants
   and normalized amplitude vectors */
void CMTsetup(WMM_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 

/* compose CMT supermodes */
void CMTsupermodes(WMM_ModeArray& umode,  // modes to be coupled
		   Waveguide cpwg,        // the entire coupler structure
				          // as output by CMTsetup:
	           Dvector smpc,          //   supermode propagation constants
	           Dmatrix smav,          //   supermode amplitude vectors 
		   WMM_ModeArray& smode); // output: CMT supermodes 

/* Coupled mode theory:
   relative output amplitude of mode o, excitement in mode i */ 
Complex CMTamp(int i,              // number of input mode
               int o,              // number of output mode
	       double len,         // device length
	    		           // as output from CMTsetup:
	       Dmatrix& sigma,     //   power coupling matrix
	       Dvector& smpc,      //   supermode propagation constants
	       Dmatrix& smav);     //   supermode amplitude vectors 

/* Coupled mode theory:
   relative output power in mode o, excitement in mode i */ 
double CMTpower(int i,              // number of input mode
                int o,              // number of output mode
		double len,         // device length
				    // as output from CMTsetup:
	        Dmatrix& sigma,     //   power coupling matrix
	        Dvector& smpc,      //   supermode propagation constants
	        Dmatrix& smav);     //   supermode amplitude vectors 

/* Coupled mode theory:
   relative output amplitude of mode o, 
   two input modes i0, i1, with complex amplitudes c0, c1 */ 
double CMTpower2(int i0,            // number of first input mode
	         Complex c0,        //   its amplitude
                 int i1,            // number of second input mode
	         Complex c1,        //   its amplitude
                 int o,             // number of output mode
	         double len,        // device length
	   		            // as output from CMTsetup:
	         Dmatrix& sigma,    //   power coupling matrix
	         Dvector& smpc,     //   supermode propagation constants
	         Dmatrix& smav);    //   supermode amplitude vectors 

/* Coupled mode theory: evaluate field superposition, local intensity  
   umode's are assumed to be normalized !!! */
double CMTsz(double x,              // coordinates on the 
	     double y,              //   waveguide cross section
	     double z,              // propagation distance
	     WMM_ModeArray& umode,  // modes to be coupled
             int i,                 // number of input mode, power 1 at z=0
	  		            // as output from CMTsetup:
	     Dmatrix& sigma,        //   power coupling matrix
	     Dvector& smpc,         //   supermode propagation constants
	     Dmatrix& smav);        //   supermode amplitude vectors 

/* Coupled mode theory: evaluate field superposition, local intensity  
   on a mesh in the y-z-plane
   (npy+1)*(npz+1) points between y0<y<y1, 0<z<ldev, 
   s(yi, zi) = CMTsz(x, y0+(y1-y0)/npy*yi, z0+ldev/npz*zi),
   umode's are assumed to be normalized !!! */
void CMTprop(double x,              // x-level 
	     double y0,             // evaluate between 
	     double y1,             //   y0<y<y1, 
	     double ldev,           //   0<z<ldev 
	     int npy,               // number of
	     int npz,               //   mesh points
	     WMM_ModeArray& umode,  // modes to be coupled
             int i,                 // number of input mode, power 1 at z=0
	                            // as output from CMTsetup:
	     Dmatrix& sigma,        //   power coupling matrix
	     Dvector& smpc,         //   supermode propagation constants
	     Dmatrix& smav,         //   supermode amplitude vectors 
	     Dmatrix& s);           // output: power values

/* Coupled mode theory: evaluate field superposition, local intensity  
   on a mesh in the y-z-plane
   (npy+1)*(npz+1) points between y0<y<y1, 0<z<ldev, 
   s(yi, zi) = CMTsz(x, y0+(y1-y0)/npy*yi, z0+ldev/npz*zi),
   two input modes i0, i1, with complex amplitudes c0, c1, 
   umode's are assumed to be normalized !!!  */
void CMTprop2(double x,              // x-level 
	      double y0,             // evaluate between 
	      double y1,             //   y0<y<y1, 
	      double ldev,           //   0<z<ldev 
	      int npy,               // number of
	      int npz,               //   mesh points
	      WMM_ModeArray& umode,  // modes to be coupled
              int i0,                // number of first input mode
	      Complex c0,            //   its amplitude
              int i1,                // number of second input mode
	      Complex c1,            //   its amplitude
	                             // as output from CMTsetup:
	      Dmatrix& sigma,        //   power coupling matrix
	      Dvector& smpc,         //   supermode propagation constants
	      Dmatrix& smav,         //   supermode amplitude vectors 
	      Dmatrix& p);           // output: power values

#endif // WMMCMT_H