#ifndef WGM_H
#define WGM_H

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

/*
 * wgm.h
 * whispering-gallery modes of circular cavities in 2-D 
 */

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

#define BF_CIRC // use routines for cylindrical functions adapted from the 
                // Circurs-package

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

/* a whispering-gallery mode */
class WGMode 
{
public:
	// Polarization
	Polarization pol;

	// cavity radius
	double R;
	// center position
	double x0;
	double z0;
	// refractive index profile; positions relative to radius R 
	Waveguide wg;
	// the cavity structure
	OvlStruct cav;	

	// angular order of the WGM
	int m;

	// complex eigenfrequency [rad/fs = rad/(10^(-15)s)] 
	inline Complex c_omega() const { return Complex(omega, alpha); };
	// related complex wavenumber  
	inline Complex c_k() const { return Complex(omega/SOLC0, alpha/SOLC0); }
	// eigenfrequency, real part [rad/fs = rad/(10^(-15)s)] 
	double omega;
	// eigenfrequency, imaginary part (attenuation) [rad/fs = rad/(10^(-15)s)] 
	double alpha;
	// 1/(omega_c * mu0) [A um/V]
	inline Complex invommu0() const 
		{ Complex o(omega, alpha); o *= MU0; o.inv(); return o; }
	// 1/(omega_c * ep0) [V um/A]
	inline Complex invomep0() const 
		{ Complex o(omega, alpha); o *= EP0; o.inv(); return o; }
	// Q factor
	double Q() const { return theQ(c_omega()); }
	// resonance wavelength
	double lambda_r() const { return val_lambda(omega); }
	// wavelength full-width-at-half-maximum of the resonance
	double fwhm_lambda() const { return FWHM_lambda(c_omega()); }

	// WGM field, time t=0 
	Complex field(Fcomp cp, double x, double z) const;
	Complex field(Fcomp cp, double r) const;
	void emfield(double x, double z,
		     Complex &ex, Complex &ey, Complex &ez,
		     Complex &hx, Complex &hy, Complex &hz) const;
	void emfield_r(double x, double z,
		       Complex &er, Complex &ey, Complex &ep,
		       Complex &hr, Complex &hy, Complex &hp) const;

	// ... values on a rectangular mesh
	Cmatrix field(Fcomp cp, Interval dwx, int npx,
				Interval dwz, int npz) const;

	// directional interference: mix with a*(reversed mode) 
	Complex field_ri(Fcomp cp, double x, double z, Complex a) const;
	void emfield_ri(double x, double z, Complex a,
		Complex &er, Complex &ey, Complex &ep,
		Complex &hr, Complex &hy, Complex &hp) const;

	// radial order: 
	// number of interior nodes in the real part of the radial profile  
        int ord;  
	// token of the mode 
	char ids[16];

	// determine the radial order: rough lookup; set mode id  
	void classify();

	// translate the WGM: shift the origin by dx, dz
	void translate(double dx, double dz);

	// direction of propagation: FORW (clockwise, default), BACK (anticlockwise) 
	Propdir dir;
	// reverse the WGM: clockwise propagation -> anticlockwise propagation
	void reverse();

	// ----------------------------------------------------------- 
	// visualization: output to MATLAB .m files  
	// evaluation of the field 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

	// profile plot corresponding to component cp (static, t=0)
	// cp: EX, EY, EZ, HX, HY, HZ, ER, EP, HR, HP 
	// 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;

	// export full field data ("everything") into a viewer MATLAB file
	// animation ignoring the field decay
        void viewer(double xbeg, double xend, double zbeg, double zend,
                    int npx, int npz, char ext0, char ext1) const;
	
	// piecewise internal representation by Bessel functions,  
	// coefficients
	Cvector cA;
	Cvector cB;

	// evaluate transverse resonance condition, set cA, cB accordingly
	Complex travres();
	Complex travres(Complex o);
	double tr_error;
        Complex tr_Ainit;
        void tr_init(Complex o);
	// converge transverse resonance condition towards root;
	// return 1: success; 0: failure 
	int converge(Complex z1, Complex z2, int quiet);

	// radial profile, principal component 
	Complex profile(double r) const;
	// ... and its derivative
	Complex d_profile(double r) const;

	// typical interior and exterior crop radii
        double cRi;
        double cRe;

	// use discretized versions (1) of the radial profile
	int disc_prof;
	// on radial interval [ir < wg.hx(.) < re] with nr steps 
        void discretize(double ri, double re, int nr);

private: 
	// storage space for discretized profiles 
	Ivector n_rad;
	Dmatrix v_rad;
	Cmatrix v_prf;
	Cmatrix v_dpr;
	Cmatrix d_prf;
	Cmatrix d_dpr;
	double d_ri;
        double d_re;
};

/* default profile crop radii: */
/* contributions for radii below 
   (radius of the innermost interface)-WGMp_R_CROP_IN*(vacuum wavelength) 
   will be neglected */
extern double WGMp_R_CROP_IN;
/* contributions for radii above  
   (radius of the outermost interface)+WGMp_R_CROP_OUT*(vacuum wavelength)
   will be neglected */
extern double WGMp_R_CROP_OUT;
/* relative field levels below WGMp_R_CROP_LVL will be neglected */
extern double WGMp_R_CROP_LVL;
/* mode profile evaluation: neglect fields at radii below 
   min(R, vac.wavelength)*WGMp_R_tiny */
extern double WGMp_R_TINY;
/* mode profile classification: relative node detection level */
extern double WGMp_NodeLVL;


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

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

	// free allocated memory
	void clear();

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

	// add a mode
	void add(WGMode m);
	// delete a mode entry
	void remove(int i);
	// add an entire WGModeArray nma 
	void merge(WGModeArray& nma);

	// sort the array by WGM attenuation: lowest first
	void sort();

	// prune the array:
	// remove modes with identical order m, with closeby frequencies, 
	// or with identical radial order
	void prune();
	// among modes with identical radial order 
	// select the one closest to resonance wavelength twl
	void prune(double twl);
	// discard all modes with resonant wavelengths outside interval twi
	void prune(Interval twi);
	
private:
	WGMode *mvec;
};

/* ------------------------------------------------------------------------
   WGM mode analysis procedures
   ------------------------------------------------------------------------ */

/* Heuristics for mode finding: 
   mode order identification by solving an equivalent straight problem 
   with Dirichlet boundary conditions */

// staircase approximation of the equiv. st. problem, number of slices
extern int WGMsP_NumSlices;
// radial distance of interior boundary to pos. R, in vacuum wavelengths
extern double WGMsP_IntBd;
// radial distance of exterior boundary to pos. R, in vacuum wavelengths
extern double WGMsP_ExtBd;
// number of spectral terms examined
extern int WGMsP_NSpT;
// minimum effective index (real part) considered
extern double WGMsP_MinNeff;
// relative extension of effective index range towards lower values
extern double WGMsP_NRangeExt;
// acceptable relative attenuation for prospective confined modes 
extern double WGMsP_RelAttLim;
// range of angular mode indices around each trial number  
extern int WGMsP_mRange;


/* complex root finding */

// initial wavelength interval, half-width in mum */
extern double WGMsP_IniDLambda;
// initial attenuation, upper level */
extern double WGMsP_IniMaxAtt;
// initial attenuation, lower level */
extern double WGMsP_IniMinAtt;
// secant method, terminal distance of root prototypes
extern double WGMsP_RootAcc;
// number of trials for convergence with modified initial conditions
extern int WGMsP_LimNumIts;
// maximum growth of trial value distance permitted per iteration step
extern double WGMsP_ConvGrowth;
// mimimum distance in abs. frequency for resonances to be considered different
extern double WGMsP_MinOmDiff;
// maximum acceptable level of absolute TRC violation 
extern double WGMsP_MaxTRC;


/* solver procedures */ 

/* for the cavity with radial layering w, radius r, centered at (xc, zc),  
   try to find pol-polarized WGMs of low radial order, 
   with eigenfrequencies next to val_omega(twl), for the target wavelength twl, 
   consider ntre * ntim trial values for complex eigenfrequencies,
   lorm=1: consider modes of lowest radial order only,  
   quiet = 0, 1, 2: all, less, no  log output  */
int wgmsolve(Waveguide w, double r, double xc, double zc, Polarization pol,
             double twl, int lorm, int ntre, int ntim, WGModeArray& wgma, int quiet);

/* for the cavity with radial layering w, radius r, centered at (xc, zc),  
   try to find pol-polarized WGMs of low radial order, 
   with eigenfrequencies corresponding to wavelengths in the interval twi, 
   consider ntre * ntim trial values at central resonance wavelength,  
   lorm=1: consider modes of lowest radial order only,  
   quiet = 0, 1, 2: all, less, no  log output  */
int wgmsolve(Waveguide w, double r, double xc, double zc, Polarization pol,
             Interval twi, int lorm, int ntre, int ntim, WGModeArray& wgma, int quiet);

/* ... generic: 
   stype = 't': search resonances at target wavelength twl,
   stype = 'i': find resonances within wavelength interval twi */
int wgmsolve(char stype, Waveguide w, double r, double xc, double zc, Polarization pol,
             double twl, Interval twi, int lorm, int ntre, int ntim, WGModeArray& wgma, int quiet);

/* Bessel functions, integer order, complex argument, wrapper functions */
Complex BfJ(int o, Complex a);
Complex BfY(int o, Complex a);
Complex dBfJ(int o, Complex a);
Complex dBfY(int o, Complex a);

#endif // WGM_H