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

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

#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"waveg.h"
#include"maxweq.h"

/* Fcomp -> field character */
char fldchr(Fcomp cp)
{
	switch(cp)
	{
		case EX: case EY: case EZ: return 'E';
		case HX: case HY: case HZ: return 'H';
		case SZ: return 'S';
	}
	fprintf(stderr, "ERR: fldchr \n"); 
	return 'E';
}

/* Fcomp -> comp character */
char cpchr(Fcomp cp)
{
	switch(cp)
	{
		case EX: case HX: return 'x';
		case EY: case HY: return 'y';
		case EZ: case HZ: return 'z';
		case SZ: return 'z';
	}
	fprintf(stderr, "ERR: cpchr \n"); 
	return 'y';
}

/* 'EHS''xyz' -> Fcomp */
Fcomp fcomp(char eh, char cp)
{
	switch(eh)
	{
		case 'E': 
			switch(cp)
			{
			case 'x': return EX;
			case 'y': return EY;
			case 'z': return EZ;
			}
		case 'H': 
			switch(cp)
			{
			case 'x': return HX;
			case 'y': return HY;
			case 'z': return HZ;
			}
		case 'S': return SZ;
	}
	fprintf(stderr, "ERR: fcomp \n"); 
	return EY;
}

/* Polarization -> 't', 'v' */
char poltochar1(Polarization p)
{
	switch(p)
	{
	case QTE: return 't';
	case QTM: return 't';
	case VEC: return 'v';
	}
	return '*';
}
/* Polarization -> 'e', 'm', 'c' */
char poltochar2(Polarization p)
{
	switch(p)
	{
	case QTE: return 'e';
	case QTM: return 'm';
	case VEC: return 'c';
	}
	return '*';
}
/* Polarization -> principal field component */
Fcomp principalcomp(Polarization p)
{
	switch(p)
	{
	case QTE: return EY;
	case QTM: return HY;
	case VEC: 
		fprintf(stderr, 
			"... don't know principal field component ...\n"); 
		return EY;
	}
	fprintf(stderr, 
		"... don't know principal field component ...\n"); 
	return EY;
}

/* Symmetry -> 'n', 's', 'a' */
char symtochar(Symmetry s)
{
	switch(s)
	{
	case SYM: return 's';
	case ASY: return 'a';
	case NOS: return 'n';
	}
	return '*';
}

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

/* 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 */

/* evaluate variable factors */
double evalvfac(Vfac f, double b, double k0, double n)
{
	switch(f)
	{
	case SFID: return 1.0;
	case SFBETA: return b;
	case SFINVBETA: return 1.0/b;
	case SFEPS: return n*n;
	case SFIBMK0E: return 1.0/(b*b-k0*k0*n*n);
	case SFK0EPS: return k0*n*n;
	case SFK0: return k0;
	case SFINVOMEPSEP0: return val_invomep0(2.0*PI/k0)/(n*n);
	case SFINVOMMU0: return val_invommu0(2.0*PI/k0);
	}
	return 1.0;
}

/* the global representations ... */
Fieldrep qteEHrep[6];  // QTE
Fieldrep qtmEHrep[6];  // QTM
Fieldrep vefEHrep[6];  // vectorial EXEY formulation
Fieldrep vhfEHrep[6];  // vectorial HXHY formulation
Fieldrep vzfEHrep[6];  // vectorial EZHZ formulation

Fieldrep* EHrep;

/* Fcomp fc is mapped to EHrep[iFcomp(fc)] */
int iFcomp(Fcomp fc)
{
	switch(fc)
	{
		case EX: return 0;
		case EY: return 1;
		case EZ: return 2;
		case HX: return 3;
		case HY: return 4;
		case HZ: return 5;
		case SZ: fprintf(stderr, "ERR: iFcomp \n"); exit(1);
	}
	fprintf(stderr, "ERR: iFcomp \n"); 
	return 1;
}
/* Fcomp fc is mapped to a Fieldrep, depending on polarization an vec. form. */
Fieldrep getfieldrep(Fcomp cp, Polarization p, Vecform v)
{
	switch(p)
	{
	case QTE: return qteEHrep[iFcomp(cp)]; 
	case QTM: return qtmEHrep[iFcomp(cp)]; 
	case VEC: 
		switch(v)
		{
		case EXEY: return vefEHrep[iFcomp(cp)];
		case HXHY: return vhfEHrep[iFcomp(cp)];
		case EZHZ: return vhfEHrep[iFcomp(cp)];
		}
	}
	fprintf(stderr, "ERR: getfieldrep \n"); 
	return qteEHrep[iFcomp(cp)]; 
}

/* ... initialize ... */
void setEHrep(Fcomp fc,
	     int numc, Vfac glsfid, double glfac, 
	     Derlev fod1, int c1, Vfac sfid1, double fac1, 
	     Derlev fod2, int c2, Vfac sfid2, double fac2, 
	     Derlev fod3, int c3, Vfac sfid3, double fac3)
{
	int i;
	i = iFcomp(fc);
	EHrep[i].numc   = numc; 
	EHrep[i].glsfid = glsfid; 
	EHrep[i].glfac  = glfac;
	EHrep[i].ctr[0].c     = c1;
	EHrep[i].ctr[0].fod   = fod1; 
	EHrep[i].ctr[0].sfid  = sfid1; 
	EHrep[i].ctr[0].fac   = fac1;
	EHrep[i].ctr[1].c     = c2;
	EHrep[i].ctr[1].fod   = fod2; 
	EHrep[i].ctr[1].sfid  = sfid2; 
	EHrep[i].ctr[1].fac   = fac2;
	EHrep[i].ctr[2].c     = c3;
	EHrep[i].ctr[2].fod   = fod3; 
	EHrep[i].ctr[2].sfid  = sfid3; 
	EHrep[i].ctr[2].fac   = fac3;
	return;
}

/* initialize global field representation (vf,pol)EHrep */
void initfieldreppv(Polarization pol, Vecform vf)
{
	switch(pol)
	{
	case QTE:
		setEHrep(EX, 
			0, SFID, 0.0,
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(EY, 
			1, SFID, 1.0,
			FLD, 0, SFID, 1.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(EZ, 
			1, SFINVBETA, -1.0,
			DYF, 0, SFID, 1.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(HX, 
			2, SFINVOMMU0, 1.0,
			FLD, 0, SFBETA, -1.0, 
			DYY, 0, SFINVBETA, 1.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(HY, 
			1, SFINVOMMU0, 1.0,
			DXY, 0, SFINVBETA, -1.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(HZ, 
			1, SFINVOMMU0, 1.0,
			DXF, 0, SFID, 1.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		break;
	case QTM:
		setEHrep(EX, 
			2, SFINVOMEPSEP0, 1.0,
			FLD, 0, SFBETA, 1.0, 
			DYY, 0, SFINVBETA, -1.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(EY, 
			1, SFINVOMEPSEP0, 1.0,
			DXY, 0, SFINVBETA, 1.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(EZ, 
			1, SFINVOMEPSEP0, 1.0,
			DXF, 0, SFID, -1.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(HX, 
			0, SFID, 0.0,
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(HY, 
			1, SFID, 1.0,
			FLD, 0, SFID, 1.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		setEHrep(HZ, 
			1, SFINVBETA, -1.0,
			DYF, 0, SFID, 1.0, 
			FLD, 0, SFID, 0.0, 
			FLD, 0, SFID, 0.0); 
		break;
	case VEC:
		switch(vf)
		{
		case HXHY:
			setEHrep(EX, 
				3, SFINVOMEPSEP0, 1.0,
				FLD, 1, SFBETA, 1.0, 
				DXY, 0, SFINVBETA, -1.0, 
				DYY, 1, SFINVBETA, -1.0); 
			setEHrep(EY, 
				3, SFINVOMEPSEP0, 1.0,
				FLD, 0, SFBETA, -1.0, 
				DXX, 0, SFINVBETA, 1.0, 
				DXY, 1, SFINVBETA, 1.0); 
			setEHrep(EZ, 
				2, SFINVOMEPSEP0, 1.0,
				DYF, 0, SFID, 1.0, 
				DXF, 1, SFID, -1.0, 
				FLD, 0, SFID, 0.0); 
			setEHrep(HX, 
				1, SFID, 1.0,
				FLD, 0, SFID, 1.0, 
				FLD, 0, SFID, 0.0, 
				FLD, 0, SFID, 0.0); 
			setEHrep(HY, 
				1, SFID, 1.0,
				FLD, 1, SFID, 1.0, 
				FLD, 0, SFID, 0.0, 
				FLD, 0, SFID, 0.0); 
			setEHrep(HZ, 
				2, SFINVBETA, -1.0,
				DXF, 0, SFID, 1.0, 
				DYF, 1, SFID, 1.0, 
				FLD, 0, SFID, 0.0); 
			break;
		case EXEY:
			setEHrep(EX, 
				1, SFID, 1.0,
				FLD, 1, SFID, 1.0, 
				FLD, 1, SFID, 0.0, 
				FLD, 1, SFID, 0.0); 
			setEHrep(EY, 
				1, SFID, 1.0,
				FLD, 0, SFID, 1.0, 
				FLD, 1, SFID, 0.0, 
				FLD, 1, SFID, 0.0); 
			setEHrep(EZ, 
				2, SFINVBETA, -1.0,
				DXF, 1, SFID, 1.0, 
				DYF, 0, SFID, 1.0, 
				FLD, 1, SFID, 0.0); 
			setEHrep(HX, 
				3, SFINVOMMU0, 1.0,
				FLD, 0, SFBETA, -1.0, 
				DXY, 1, SFINVBETA, 1.0, 
				DYY, 0, SFINVBETA, 1.0); 
			setEHrep(HY, 
				3, SFINVOMMU0, 1.0,
				FLD, 1, SFBETA, 1.0, 
				DXX, 1, SFINVBETA, -1.0, 
				DXY, 0, SFINVBETA, -1.0); 
			setEHrep(HZ, 
				2, SFINVOMMU0, 1.0,
				DXF, 0, SFID, 1.0, 
				DYF, 1, SFID, -1.0, 
				FLD, 1, SFID, 0.0); 
			break;
		case EZHZ:
			setEHrep(EX, 
				2, SFIBMK0E, INVSQRTEP0,
				DYF, 0, SFK0, -1.0, 
				DXF, 1, SFBETA, -1.0, 
				FLD, 1, SFID, 0.0); 
			setEHrep(EY, 
				2, SFIBMK0E, INVSQRTEP0,
				DXF, 0, SFK0, 1.0, 
				DYF, 1, SFBETA, -1.0, 
				FLD, 1, SFID, 0.0); 
			setEHrep(EZ, 
				1, SFID, INVSQRTEP0,
				FLD, 1, SFID, 1.0, 
				FLD, 1, SFID, 0.0, 
				FLD, 1, SFID, 0.0); 
			setEHrep(HX, 
				2, SFIBMK0E, INVSQRTMU0,
				DYF, 1, SFK0EPS, 1.0, 
				DXF, 0, SFBETA, -1.0, 
				FLD, 1, SFID, 0.0); 
			setEHrep(HY, 
				2, SFIBMK0E, INVSQRTMU0,
				DXF, 1, SFK0EPS, -1.0, 
				DYF, 0, SFBETA, -1.0, 
				FLD, 1, SFID, 0.0); 
			setEHrep(HZ, 
				1, SFID, INVSQRTMU0,
				FLD, 0, SFID, 1.0, 
				FLD, 1, SFID, 0.0, 
				FLD, 1, SFID, 0.0); 
			break;
		}
		break;
	}
	return;
}

/* initialize global field representation (vf,pol)EHrep */
void initfieldrep()
{
	EHrep = qteEHrep; initfieldreppv(QTE, EXEY);
	EHrep = qtmEHrep; initfieldreppv(QTM, HXHY);
	EHrep = vefEHrep; initfieldreppv(VEC, EXEY);
	EHrep = vhfEHrep; initfieldreppv(VEC, HXHY);
	EHrep = vzfEHrep; initfieldreppv(VEC, EZHZ);
//	fprintf(stderr, "- Maxweq -\n");
	return;
}

/* ... invoke */
Maxweq::Maxweq()
{
	initfieldrep();
}

Maxweq MINIT;

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

/* typical initial field: linearly polarized Gaussian beam with circular 
   profile, propagating in a homogeneous medium in z-direction */ 

/* initialize */
Gaussianbeam::Gaussianbeam()
{
	l = 1.0;
	ng = 1.0;
	x0 = 0.0;
	y0 = 0.0;
	rg = 1.0;
	p0 = 1.0;
	alpha = 0.0;
}
Gaussianbeam::Gaussianbeam(double pl, double png, double px0, double py0, 
			   double prg, double pp0, double palpha)
{
	l = pl;
	ng = png;
	x0 = px0;
	y0 = py0;
	rg = prg;
	p0 = pp0;
	alpha = palpha;
}

/* real fields in A/um or V/um */ 
double Gaussianbeam::field(Fcomp cp, double x, double y)
{
	double f;
	if(cp == EZ) return 0.0;
	if(cp == HZ) return 0.0;
	if(cp == SZ) 
	{
		fprintf(stderr, "\n Gaussianbeam: SZ not implemented !\n");
		return 1.0;
	}
	f = sqrt(4.0/PI);
	f *= 1.0/rg;
	if(cp == HX || cp == HY) f *= 1.0/sqrt(val_k0(l)*val_invomep0(l)/ng);
	if(cp == EX || cp == EY) f *= sqrt(val_k0(l)*val_invomep0(l)/ng);
	f *= sqrt(p0);
	if(cp == HX) f *= -cos(alpha);
	if(cp == HY) f *=  sin(alpha);
	if(cp == EX) f *=  sin(alpha);
	if(cp == EY) f *=  cos(alpha);
	f *= exp(-((x-x0)*(x-x0)+(y-y0)*(y-y0))/rg/rg);
	return f;
}

/* 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)
{
	double ovl = 0.0; 
	int xi, yi;
	double x, y;
	double xs, ys;
	double fex, fhy, fey, fhx;
	xs = (r.x1-r.x0)/numpx;
	ys = (r.y1-r.y0)/numpy;
	for(yi=0; yi<=numpy; ++yi)
	{
		y = r.y0+((double) yi)*ys;
		for(xi=0; xi<=numpx; ++xi)
		{
			x = r.x0+((double) xi)*xs;
			fhy = fld2(HY, x, y);
			fex = 0; if(fhy != 0.0) fex = fld1(EX, x, y);
			fhx = fld2(HX, x, y);
			fey = 0; if(fhx != 0.0) fey = fld1(EY, x, y);
 			ovl += fex*fhy-fey*fhx; 
		}
	}
	return 0.5*ovl*xs*ys;
}
   	
/* ##################################################################### */

/* permittivity perturbation */

/* initialize */
Perturbation::Perturbation()
{
	r.x0 = r.y0 = r.x1 = r.y1 = 0.0;
	e = Cmatrix(3, 3);
	e.init(CC0);
}

/* initialize */
Perturbation::Perturbation(Rect rv)
{
	r = rv;
	e = Cmatrix(3, 3);
	e.init(CC0);
}

/* initialize */
Perturbation::Perturbation(Rect ri, double rxxv, double ryyv, double rzzv, 
			   double rxyv, double ixzv, double iyzv, double iabsv)
{
	r = ri;
	e = Cmatrix(3, 3);
	e(0, 0) = Complex(rxxv, iabsv);
	e(0, 1) = Complex(rxyv, 0.0);
	e(0, 2) = Complex(0.0,  ixzv);
	e(1, 0) = Complex(rxyv, 0.0);
	e(1, 1) = Complex(ryyv, iabsv);
	e(1, 2) = Complex(0.0,  iyzv);
	e(2, 0) = Complex(0.0, -ixzv);
	e(2, 1) = Complex(0.0, -iyzv);
	e(2, 2) = Complex(rzzv, iabsv);
}

/* initialize */
Perturbation::Perturbation(Rect ri, 
	                   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)	
{
	r = ri;
	e = Cmatrix(3, 3);
	e(0, 0) = p_xx;
	e(0, 1) = p_xy;
	e(0, 2) = p_xz;
	e(1, 0) = p_yx;
	e(1, 1) = p_yy;
	e(1, 2) = p_yz;
	e(2, 0) = p_zx;
	e(2, 1) = p_zy;
	e(2, 2) = p_zz;
}


/* output to file dat */
void Perturbation::write(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, "perturbed region: "); 
		r.write(dat);
		fprintf(dat, "\n"); 
		fprintf(dat, "permittivity perturbation:\n"); 
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(0,0).re,e(0,0).im,e(0,1).re,e(0,1).im,e(0,2).re,e(0,2).im);
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(1,0).re,e(1,0).im,e(1,1).re,e(1,1).im,e(1,2).re,e(1,2).im);
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(2,0).re,e(2,0).im,e(2,1).re,e(2,1).im,e(2,2).re,e(2,2).im);
	}
	else
	{
		comment("perturbation:", dat);
		r.write(dat);
		e.write(dat);
	}
	return;
}
		
/* input from file dat */
void Perturbation::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		r.x0 = r.y0 = r.x1 = r.y1 = 0.0;
		e = Cmatrix(3, 3);
		e.init(CC0);
	}
	else
	{
		r.read(dat);
		e.read(dat);
	}
	return;
}

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

/* 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)
{
	Perturbation p(rl); 
	double deps;
	deps = - n*alpha*lambda/2.0/PI;
	p.e.init(CC0);
	p.e(0, 0) = Complex(0.0, deps);
	p.e(1, 1) = Complex(0.0, deps);
	p.e(2, 2) = Complex(0.0, deps);
	return p;
}

/* 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)
{
	Perturbation p(rl); 
	double xi, w;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	w = sin(theta)*cos(phi);
	p.e(1, 2) = Complex(0.0, w*xi);
	p.e(2, 1) = Complex(0.0,-w*xi);
	w = sin(theta)*sin(phi);
	p.e(0, 2) = Complex(0.0,-w*xi);
	p.e(2, 0) = Complex(0.0, w*xi);
	w = cos(theta);
	p.e(0, 1) = Complex(0.0, w*xi);
	p.e(1, 0) = Complex(0.0,-w*xi);
	return p;
} 
/* as above, for polar orientation of the magnetization, vec(M) || x-axis, */
Perturbation magopt_polar(double sFrot, double n, double lambda, Rect rl)
{
	Perturbation p(rl); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(1, 2) = Complex(0.0, xi);
	p.e(2, 1) = Complex(0.0,-xi);
	return p;
}
/* for equatorial orientation of the magnetization, vec(M) || y-axis, */
Perturbation magopt_equat(double sFrot, double n, double lambda, Rect rl)
{
	Perturbation p(rl); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(0, 2) = Complex(0.0,-xi);
	p.e(2, 0) = Complex(0.0, xi);
	return p;
}
/* and for longitudinal orientation of the magnetization, vec(M) || z-axis */
Perturbation magopt_long(double sFrot, double n, double lambda, Rect rl)
{
	Perturbation p(rl); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(0, 1) = Complex(0.0, xi);
	p.e(1, 0) = Complex(0.0,-xi);
	return p;
}