#ifndef MAXWEQ_H
#define MAXWEQ_H

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

/* 
 * maxweq.h
 * relations concerning the electromagnetic fields for waveguide modes  
 */

/* compiler xlC on IBM SP needs this statement ?!? */
#undef HZ

#define PI	3.14159265358979323846
#define PI_2	1.57079632679489661923
#define PI_4	0.78539816339744830962
#define PIM2    6.28318530717958647693

/* vacuum speed of light [um/fs] */
#define	SOLC0       2.99792458E-1
/* 1/sqrt(mu0) [sqrt(A*um/(V*fs))] */
#define INVSQRTMU0  2.82094792E-2 
/* 1/sqrt(ep0) [sqrt(V*um/(A*fs))] */
#define	INVSQRTEP0 10.62736593 

/* vacuum permittivity [10^(-15) A s/(V um) = A fs/(V um) ] */
#define EP0 8.854187817E-3
/* vacuum permeability [10^(-15) V s/(A um) = V fs/(A um) ] */
#define MU0 1.2566370614E+3 

/* vacuum wavenumber [1/um], function of the vacuum wavelength [um] */
inline double val_k0(double lambda) { return 6.283185307179586/lambda; }; 
/* 1/omega/epsilon_0 [V um / A], function of the vacuum wavelength [um] */
inline double val_invomep0(double lambda) { return 59.95849160*lambda; };
/* 1/omega/mu_0 [A um / V], function of the vacuum wavelength [um] */
inline double val_invommu0(double lambda) { return (4.224638618e-4)*lambda; }; 
/* angular frequency [radiants/fs], function of the vacuum wavelength [um] */
inline double val_omega(double lambda) { return 1.883651567/lambda; };
/* corresponding period in time [fs], function of the vacuum wavelength [um] */
inline double val_period_T(double lambda) { return lambda*3.335640952; }; 
/* vacuum wavelength [um], function of the angular frequency [radiants/fs] */
inline double val_lambda(double omega) { return 1.883651567/omega; };

/* field components */ 
enum Fcomp {EX, EY, EZ, HX, HY, HZ, SZ};

/* Fcomp -> field character */
char fldchr(Fcomp cp);

/* Fcomp -> comp character */
char cpchr(Fcomp cp);

/* 'EHS''xyz' -> Fcomp */
Fcomp fcomp(char eh, char cp);

/* polarization parameter */ 
enum Polarization {QTE, QTM, VEC};

/* Polarization -> 't', 'v' */
char poltochar1(Polarization p);
/* Polarization -> 'e', 'm', 'c' */
char poltochar2(Polarization p);
/* Polarization -> principal field component */
Fcomp principalcomp(Polarization p);

/* symmetry parameter */
enum Symmetry {SYM, ASY, NOS};

/* Symmetry -> 'n', 's', 'a' */
char symtochar(Symmetry s);

/* Attributes for field output */ 
// MOD: absolute value 
// ORG: the original field 
// SQR: the squared field
// REP: real part 
// IMP: imaginary part 
enum Afo {MOD, ORG, SQR, REP, IMP};

/* derivation levels */ 
enum Derlev {FLD=0, DXF=1, DYF=2, DXY=3, DXX=4, DYY=5};

/* vectorial mode formulation */
enum Vecform {EXEY, HXHY, EZHZ};

/* ##################################################################### */

/* representation of a field component as a sum of derivations of one or
   two basic fields, depending on the polarization and formulation of the 
   vectorial mode problem */

/* variable factors */
enum Vfac {SFID, SFBETA, SFINVBETA, SFEPS, SFIBMK0E,
	   SFK0EPS, SFK0, SFINVOMEPSEP0, SFINVOMMU0};
/* ... evaluate */
double evalvfac(Vfac f, double b, double k0, double n);

/* the global representations ... */
typedef struct
{
	int      numc;         /* Number of contributions; 1, 2, or 3 */
	Vfac   glsfid;         /* global factors: variabel */
	double  glfac;         /*                 fixed    */
	struct
	{
		int        c;         /* basic component 0/1 */ 
		Derlev   fod;         /* its derivation level */
		Vfac    sfid;         /* variable factor */
		double   fac;         /* fixed factor */
	} ctr[3];
} Fieldrep;

extern Fieldrep qteEHrep[6];  // QTE
extern Fieldrep qtmEHrep[6];  // QTM
extern Fieldrep vefEHrep[6];  // vectorial EXEY formulation
extern Fieldrep vhfEHrep[6];  // vectorial HXHY formulation
extern Fieldrep vzfEHrep[6];  // vectorial EZHZ formulation

/* Fcomp fc is mapped to EHrep[iFcomp(fc)] */
int iFcomp(Fcomp fc);
/* Fcomp fc is mapped to a Fieldrep, depending on polarization an vec. form. */
Fieldrep getfieldrep(Fcomp cp, Polarization p, Vecform v);

/* initialize global field representations */
void initfieldreps();

class Maxweq
{
public:
	// initialize
	Maxweq();
};

extern Maxweq MINIT;

/* ##################################################################### */

/* typical initial field: linearly polarized Gaussian beam with circular 
   profile, propagating in a homogeneous medium in z-direction */ 
class Gaussianbeam
{
public:
	double l;     // vacuum wavelength, given in um
	double ng;    // refractive index 
	double x0;    // location of intensity maximum
	double y0;    
	double rg;    // beam radius
	double p0;    // total power in W
	double alpha; // polarization angle versus TE 

	// initialize 
	Gaussianbeam();
	Gaussianbeam(double l, double ng, double x0, double y0, double rg,
		     double p0, double alpha);

	// real fields in A/um  or V/um  
	double field(Fcomp cp, double x, double y);
};

/* overlap 0.5*intint fld1_EX*fld2_HX - fld1_EY*fld2_HY dxdy of 
   real fields, integrals are evaluated as sums over step-functions on an 
   equidistant mesh with (numx+1)*(numy+1) points on rectangle r */
double overlap(double (*fld1)(Fcomp cp, double x, double y),
               double (*fld2)(Fcomp cp, double x, double y),
	       Rect r, int numpx, int numpy);

/* ##################################################################### */

/* permittivity perturbation */
class Perturbation
{
public:
	// part of x-y-plane with nonzero perturbation 	
	Rect r;          

	// entries of the permittivity perturbation tensor, complex 3x3 matrix 
	Cmatrix e;

	// initialize 
	Perturbation();
	Perturbation(Rect r);
	Perturbation(Rect r, 
	             Complex p_xx, Complex p_xy, Complex p_xz,	
	             Complex p_yx, Complex p_yy, Complex p_yz,	
	             Complex p_zx, Complex p_zy, Complex p_zz);	
	Perturbation(Cmatrix p);
	// relevant quantities for isotropic absorbing part 
        // and Hermitian anisotropic part  
	Perturbation(Rect r, double rxx, double ryy, double rzz, 
			     double rxy, double ixz, double iyz, double iabs);

	// input
	void read(FILE *dat);
	//output
	void write(FILE *dat);
};

/* ##################################################################### */

/* special perturbations */

/* isotropic absorbtion: 
   rl: the lossy area on the waveguide cross section 
   alpha: intensity attenuation of the material in rl, [1/micrometer]
          the material damps the intensity of a plane wave according to
          exp(-alpha z)
   n: refractive index of the material in rl
   lambda: vacuum wavelength [micrometer]
   The attenuation leads to a complex refractive index n - i alpha lambda/4/pi,
   or a permittivity perturbation of - i n alpha lambda /(2 pi) */
Perturbation attenuation(double alpha, double n, double lambda, Rect rl);

/* magnetooptic permittivity contribution:
   rl: the magnetooptic area on the waveguide cross section 
   sFrot: specific Faraday rotation of the material [degrees/cm] (!)
   theta, phi: orientation of the magnetization in spherical coordinates,
               vec(M) = M(sin(theta)cos(phi), sin(theta)sin(phi), cos(theta)) 
   n: refractive index of the material in rl
   lambda: vacuum wavelength [micrometer] */
Perturbation magopt(double sFrot, double theta, double phi, 
                    double n, double lambda, Rect rl);

/* as above, for polar orientation of the magnetization, vec(M) || x-axis, */
Perturbation magopt_polar(double sFrot, double n, double lambda, Rect rl);

/* for equatorial orientation of the magnetization, vec(M) || y-axis, */
Perturbation magopt_equat(double sFrot, double n, double lambda, Rect rl);

/* and for longitudinal orientation of the magnetization, vec(M) || z-axis */
Perturbation magopt_long(double sFrot, double n, double lambda, Rect rl);

/* ##################################################################### */

#endif // MAXWEQ_H