#ifndef QUEPFLD_H
#define QUEPFLD_H

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

/*
 * quepfld.h
 * Quadridirectional mode expansion, crossed sequences of arrays of modes
 */

class QuepField 
{
public:
	// two crossed sequences of waveguide segments
	BepField h;    // horizontal 
	BepField v;    // vertical 
	// computational window
	Interval cwx;  // vertical
	Interval cwz;  // horizontal

	// initialize
	QuepField();
	QuepField(SegWgStruct hst, Interval wx, Interval wz);
	QuepField(Circuit rc, Interval wx, Interval wz);

	// full QUEP initialization: setup the container for the 
	// quadridirectional mode expansion & spectral discretization
	//    rc: the structure under investigation,
	//    p:  fields of polarization p,
	//    wx: vertical computational window,
	//    Mx: number of expansion terms per vertical slice,
	//    wz: horizontal computational window,
	//    Mx: number of expansion terms per horizontal layer,
	//    suppress log output, if quiet == 1 */
        QuepField(SegWgStruct rc, Polarization p,
                  Interval wx, int Mx, Interval wz, int Mz, int quiet);
	QuepField(SegWgStruct rc, Polarization p,
	          Interval wx, int Mx, Interval wz, int Mz);
        QuepField(Circuit rc, Polarization p,
                  Interval wx, int Mx, Interval wz, int Mz, int quiet);
	QuepField(Circuit rc, Polarization p,
	          Interval wx, int Mx, Interval wz, int Mz);

	// destroy
	~QuepField();
	// copy constructor 
	QuepField(const QuepField& quef);
	// assignment 
	QuepField& operator=(const QuepField& quef);

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

	// ---------------------------------------------------------------
	// the waveguides associated with the borders of the circuit
	Waveguide port(CBorder bdr) const;
	// the arrays of in- and outgoing modes crossing the circuit borders
	ModeArray modes(CBorder bdr) const;

	// ---------------------------------------------------------------
	// specification of input fields, additive:

	// mode m with amplitude A at border bdr
	void input(CBorder bdr, Complex A, const Mode& m); 
	// mode of order o with amplitude A at border bdr
	void input(CBorder bdr, Complex A, int o); 
	// tilted Gaussian beam gb, multiplied by amplitude A, at border bdr
	void input(CBorder bdr, const Gaussianbeam& gb, Complex A); 
	// ... a normalized beam, unspecified phase 
	void input(CBorder bdr, const Gaussianbeam& gb); 

	// vM: amplitude of incoming mode of order o passing border bdr, 
	// with polarization pol, evaluated at propagation coordinate p   
	Complex Ain(CBorder bdr, Polarization pol, int o, double p) const;

	// ---------------------------------------------------------------
	// output amplitudes & powers:

	// amplitude of outgoing mode m at border bdr 
	Complex Aout(CBorder bdr, const Mode& m) const;
	// amplitude of outgoing mode of order o (number o) at border bdr 
	Complex Aout(CBorder bdr, int o) const;
	// vM: amplitude of outgoing mode of order o at border bdr, 
	// with polarization pol  
	Complex Aout(CBorder bdr, Polarization pol, int o) const;
	// vM: amplitude of outgoing mode of order o passing border bdr, 
	// with polarization pol, evaluated at propagation coordinate p   
	Complex Aout(CBorder bdr, Polarization pol, int o, double p) const;

	// output power associated with mode m at border bdr
	double Pout(CBorder bdr, const Mode& m) const;
	// output power associated with mode of order o (number o) at border bdr
	double Pout(CBorder bdr, int o) const;
	// output power associated with mode of order o at border bdr 
	// with polarization pol (vM) 
	double Pout(CBorder bdr, Polarization pol, int o) const;

	// guided output power crossing border bdr
	double Pgout(CBorder bdr) const;
	// guided output power crossing border bdr, 
	// associated with polarization pol (vM) 
	double Pgout(CBorder bdr, Polarization pol) const;

	// total output power crossing border bdr (guided and nonguided fields)
	double Pout(CBorder bdr) const;
	// total output power crossing border bdr (guided and nonguided fields)
	// associated with polarization pol (vM) 
	double Pout(CBorder bdr, Polarization pol) const;


	// ---------------------------------------------------------------
	// input power (readout only):

	// total input power crossing border bdr (guided and nonguided fields)
	double Pin(CBorder bdr) const;

	// ---------------------------------------------------------------
	// inspect modal expansions & amplitudes at the border regions 
	void inspect(CBorder bdr) const;
	void inspect() const;

	// ---------------------------------------------------------------
	// compute the discretized quadridirectional mode spectra: 
	// modes of polarization p, 
        //    Mx expansion terms per slice, 
	//    Mz terms per layer  
	//    suppress log output, if quiet == 1
	void discretize(Polarization p, int Mx, int Mz, int quiet);

	// ---------------------------------------------------------------
	// check the QuepField for consistency ... 
	void check();	
	// ... & set global values for
	Polarization QF_pol;     // polarization (obsolete, if QF_vM) 
	double       QF_lambda;  // vacuum wavelength
	int          QF_vM;      // vTE and vTM modes 
	Complex      QF_ky;      // with uniform wavenumber ky

	// ---------------------------------------------------------------
	// quadridirectional mode interference evaluation 

	// field superposition at point x, z 
	// cp: EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
	Complex field(Fcomp cp, double x, double z) const;

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

	// local electromagnetic energy density at x, z
	double endens(double x, double z) const;

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

	// ---------------------------------------------------------------

	// quadridirectional mode expansion solver:
	// given the structure with the accompanying mode expansions,
	// where input mode amplitudes are specified by 
        // h.setf(0, f), h.setb(h.s.nz+1, b), 
	// v.setf(0, u), v.setb(v.s.nz+1, d), or the input methods,
	// find the remaining mode coefficients.
	void quepsim();

	// as before, with a test whether the total relative power balance is 
	// violated by not more than pviol,
	// return value 1: ok, 0: no success. 
	int quepsim(double pviol);

	// as before, but with a movable computational window:
	// try at most numtr times to obtain a result with a violation of the
	// total relative power balance of not more than pviol,
	// by repeating the calculation with computational window boundary
	// positions shifted by small integer multiples of +/- bshift times 
	// the window width, return value 1: ok, 0: no success.
	int quepsim_s(int numtr, double pviol, double bshift);
	
	// as before, with default boundary shift parameters
	int quepsim_s() { return quepsim_s(5, 0.005, 0.01); };

	// ---------------------------------------------------------------
	// vQUEP implementation

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

	// compute the discretized quadridirectional mode spectra:
	// vTE and vTM modes, uniform wavenumber ky,
	// 2*Mx expansion terms per vertical slice,
	// 2*Mz terms per horizontal layer 
	void discretize(Complex ky, int Mx, int Mz, int quiet);

	// full vQUEP initialization: setup the container for the 
	// quadridirectional mode expansion & spectral discretization
	//    rc: the structure under investigation,
	//    ky: fields of vTE and vTM polarizations, with y-wavenumber ky 
	//    wx: vertical computational window,
	//    Mx: number of expansion terms per vertical slice,
	//    wz: horizontal computational window,
	//    Mx: number of expansion terms per horizontal layer,
	//    suppress log output, if quiet == 1 */
        QuepField(SegWgStruct rc, Complex ky,
                  Interval wx, int Mx, Interval wz, int Mz, int quiet);
	QuepField(SegWgStruct rc, Complex ky,
	          Interval wx, int Mx, Interval wz, int Mz);
        QuepField(Circuit rc, Complex ky,
                  Interval wx, int Mx, Interval wz, int Mz, int quiet);
	QuepField(Circuit rc, Complex ky,
	          Interval wx, int Mx, Interval wz, int Mz);
	// ... with ky specified by the input field: 
	// guided mode of polarization pol, order ord, 
	// coming in at border bdr, with amplitude amp, 
	// at angle of incidence theta (^o) 
        QuepField(SegWgStruct rc, CBorder bdr, 
	          Polarization pol, int ord, Complex amp, double theta,
                  Interval wx, int Mx, Interval wz, int Mz, int quiet);
        QuepField(SegWgStruct rc, CBorder bdr, 
	          Polarization pol, int ord, Complex amp, double theta,
                  Interval wx, int Mx, Interval wz, int Mz);
        QuepField(Circuit rc, CBorder bdr, 
	          Polarization pol, int ord, Complex amp, double theta,
                  Interval wx, int Mx, Interval wz, int Mz, int quiet);
        QuepField(Circuit rc, CBorder bdr, 
	          Polarization pol, int ord, Complex amp, double theta,
                  Interval wx, int Mx, Interval wz, int Mz);

	// ---------------------------------------------------------------
	// 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 computational window
	void adjustphase(int numx, int numz);
	// ... as before, with (coarse) defaults for numx, numz
	void adjustphase();

	// adjust the global phase, such that a plot(cp, ORG) of the physical
	// component cp exhibits the maximum field at position (x, z) 
	void adjustphase(Fcomp cp, double x, double z);

	// adjust the global phase, such that a plot(cp, ORG) of the physical
	// component cp 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 computational window
	void adjustphase(Fcomp cp, int numx, int numz);
	// ... as before, with (coarse) defaults for numx, numz
	void adjustphase(Fcomp cp);

	// ---------------------------------------------------------------
	// 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(Fcomp cp, Interval dwx, int npx, 
	                        Interval dwz, int npz) const;
	Dmatrix field(Fcomp cp, Afo foa, Interval dwx, int npx, 
	                                 Interval dwz, int npz) const;

	// ... np equidistant field values on the straight line 
	// between (x0, z0) and (x1, z1) 
	Cvector field(Fcomp cp, double x0, double z0, 
	                        double x1, double z1, int np) const;
	Dvector field(Fcomp cp, Afo foa, double x0, double z0, 
	                                 double x1, double z1, int np) const;

	// ---------------------------------------------------------------
	// 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 plot points in the output mesh
	// ext0, ext1: filename id characters for the m.file
        // assumption: uniform wavelength of all included modes 
	
	// 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, 
                  double xbeg, double xend, double zbeg, double zend,
                  int npx, int npz, char ext0, char ext1) const; 

	// image plot of the relative phase
	// cp: one of EX, EY, EZ, HX, HY, HZ
	void phasemap(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, 
	          double xbeg, double xend, double zbeg, double zend,
                  int npx, int npz, char ext0, char ext1) const;

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

	// fancy style surface plot corresponding to field component cp
	// cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ
	void fplot(Fcomp cp,
	           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 component cp, or EY (TE) | HY (TM), resp.,
	// at ntfr equidistant times over one time period
	void movie(Fcomp cp, 
	           double xbeg, double xend, double zbeg, double zend,
                   int npx, int npz, int ntfr, char ext0, char ext1) const;
	void movie(double xbeg, double xend, double zbeg, double zend,
                   int npx, int npz, int ntfr, char ext0, char ext1) const;

	// animation of the light propagation, fancy style ;-)
	// the frames show images of component cp, or EY (TE) | HY (TM), resp.,
	// at ntfr equidistant times over one time period
	void fmovie(Fcomp cp, 
	            double xbeg, double xend, double zbeg, double zend,
	            int npx, int npz, int ntfr, char ext0, char ext1) const;
	void fmovie(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;
};

#endif // QUEPFLD_H