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

/* 
 * waveg.cpp
 * waveguide definition 
 */

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

/* error message */
void wavegerror(const char *m) 
{
	fprintf(stderr, "\n ERR: %s\n", m);
	exit(1);
}

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

/* Rect: rectangular section from the waveguide cross section */
/* initialize */
Rect::Rect()
{
	x0 = 0.0;
	y0 = 0.0;
	x1 = 1.0;
	y1 = 1.0;
}
Rect::Rect(double xa, double ya, double xb, double yb)
{
	if(xa <= xb) { x0 = xa; x1 = xb; }
	else         { x0 = xb; x1 = xa; }
	if(ya <= yb) { y0 = ya; y1 = yb; }
	else         { y0 = yb; y1 = ya; }
}

/* output to FILE dat */
void Rect::write(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, "[%g,%g|%g,%g]", x0, y0, x1, y1);
	}
	else     
	{
		fputdouble(x0, dat);
		fputdouble(y0, dat);
		fputdouble(x1, dat);
		fputdouble(y1, dat);
	}
	return;
}

/* input from FILE dat */
void Rect::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) 
	{
		wavegerror("Rect: read");
	}
	else 
	{
		x0 = fgetdouble(dat);
		y0 = fgetdouble(dat);
		x1 = fgetdouble(dat);
		y1 = fgetdouble(dat);
	}
	return;
}
 
/* ------------------------------------------------------------------- */

/* Waveguide: waveguide geometry + permittivity profile; vacuum wavelength */
#define AWAY 1.0e+10

/* initialize */
Waveguide::Waveguide()
{
	nx = 0;
	ny = 0;	
	hx = Dvector();
	hy = Dvector();
	lambda = 1.0;
	n = Dmatrix();
	xyplane = Rect(-AWAY,-AWAY,AWAY,AWAY);
}

/* initialize */
Waveguide::Waveguide(int vnx, int vny) 
{
	nx = vnx;
	ny = vny;
	hx = Dvector(nx+1);
	hy = Dvector(ny+1);
	lambda = 1.0;
	n = Dmatrix(nx+2, ny+2);
	xyplane = Rect(-AWAY,-AWAY,AWAY,AWAY);
}

/* initialize */
Waveguide::Waveguide(int vnx, Dvector vhx,
		     int vny, Dvector vhy, 
		     double l, Dmatrix vn)
{
	nx = vnx;
	ny = vny;
	hx = vhx;
	hy = vhy;
	lambda = l;
	n = vn;
	xyplane = Rect(-AWAY,-AWAY,AWAY,AWAY);
}

/* free allocated memory */
void Waveguide::strip()
{
	nx = 0;
	ny = 0;
	hx.strip();
	hy.strip();
	n.strip();
}

/* set Rectangle index (l, m) corresponding to position (x,y) */ 
void Waveguide::rectidx(double x, double y, int& l, int& m) const
{
	l = 0;
	while(l<=nx && x>hx(l)) ++l;
	m = 0;
	while(m<=ny && y>hy(m)) ++m;
	return;
}

/* get rectangle boundaries corresponding to index (l,m) */  
Rect Waveguide::rectbounds(int l, int m) const
{
	Rect r;
	if(l <= 0) r.x0 = xyplane.x0;
	else r.x0 = hx(l-1);
	if(l >= nx+1) r.x1 = xyplane.x1;
	else r.x1 = hx(l);
	if(m <= 0) r.y0 = xyplane.y0;
	else r.y0 = hy(m-1);
	if(m >= ny+1) r.y1 = xyplane.y1;
	else r.y1 = hy(m);
	return r;
}

/* get rectangle boundaries corresponding to position (x,y) */ 
Rect Waveguide::rectbounds(double x, double y) const
{
	int l, m;
	rectidx(x, y, l, m);
	return rectbounds(l, m);
}

/* test neighbourhood of two rectangles */
int Waveguide::testconnect(int l0, int m0, int l1, int m1) const
{
	int ld, md;
	ld = l1-l0; if(ld<0) ld = -ld;	
	md = m1-m0; if(md<0) md = -md;	
	if(ld==1 && md==0) return 1;
	if(ld==0 && md==1) return 1;
	if(ld==0 && md==0) return 1;
	return 0;
}

/* test matching of boundary positions with another waveguide */
#define HDIST 1.0e-8
int Waveguide::bdmatch(Waveguide w) const
{
	int j, k, e;
	for(j=0; j<=w.nx; ++j)
	{
		e = 0;
		for(k=0; k<=nx; ++k)
		{
			if(fabs(hx(k)-w.hx(j))<=HDIST) e = 1; 
		}
		if(e != 1) return 1;	
	}
	for(j=0; j<=w.ny; ++j)
	{
		e = 0;
		for(k=0; k<=ny; ++k)
		{
			if(fabs(hy(k)-w.hy(j))<=HDIST) e = 1; 
		}
		if(e != 1) return 1;	
	}
	return 0;
}

/* permittivity on rectangle l, m */
double Waveguide::eps(int l, int m) const
{
	double e;
	e = n(l,m);
	return e*e;
}

/* permittivity at position x, y */
double Waveguide::eps(double x, double y) const
{
	double e;
	int l, m;
	rectidx(x, y, l, m);
	e = n(l,m);
	return e*e;
}

/* upper bound for effective mode indices, default value */
double Waveguide::defaultneffmax() const
{
	int l, m;
	double neffmax = 0.0;
	for(l=0; l<=nx+1; ++l)
	{
		for(m=0; m<=ny+1; ++m)
		{
			if(neffmax<n(l,m)) neffmax = n(l,m);
		}
	}
	return neffmax;
}

/* lower bound for effective mode indices, default value */
double Waveguide::defaultneffmin() const
{
	int l, m;
	double neffmin = 0.0;
	double nmax;
	nmax = defaultneffmax();
	for(l=0; l<=nx+1; ++l)
	{
		m=0;    if(neffmin<n(l,m) && n(l,m)<nmax) neffmin = n(l,m);
		m=ny+1; if(neffmin<n(l,m) && n(l,m)<nmax) neffmin = n(l,m);
	}
	for(m=1; m<=ny; ++m)
	{
		l=0;    if(neffmin<n(l,m) && n(l,m)<nmax) neffmin = n(l,m);
		l=nx+1; if(neffmin<n(l,m) && n(l,m)<nmax) neffmin = n(l,m);
	}
	return neffmin;
}

/* translate: hx -> hx+dx, hy -> hy+dy */
void Waveguide::translate(double dx, double dy)
{
	int l, m;
	for(l=0; l<=nx; ++l) hx(l) += dx;
	for(m=0; m<=ny; ++m) hy(m) += dy;
	xyplane.x0 += dx;
	xyplane.x1 += dx;
	xyplane.y0 += dy;
	xyplane.y1 += dy;
	return;
}

/* check mirror symmetry of the waveguide w.r.t. central x-z-plane, return value m
   m == 0: no symmetry
   m >= 1: symmetric structure with respect to slice m  */
int Waveguide::checksymmetry_pmy() const
{
	if(ny % 2 == 0) return 0;
	double c = (hy(ny)+hy(0))/2.0;
	int cm, dl;
        rectidx(0.0, c, dl, cm);
	for(int m=0; m<=ny; ++m)
	{
		if(fabs((c-hy(m)) - (hy(ny-m)-c))> 1.0e-7) return 0;
	}
	for(int l=0; l<=nx+1; ++l)
	{
		for(int m=0; m<=ny+1; ++m)
		{
			if(fabs(n(l,m)-n(l,ny-m+1)) > 1.0e-7) return 0;
		}
	}
	return cm;
}

/* check mirror symmetry of the waveguide w.r.t. central y-z-plane, return value l
   l == 0: no symmetry
   l >= 1: symmetric structure with respect to layer l  */
int Waveguide::checksymmetry_pmx() const
{
	if(nx % 2 == 0) return 0;
	double c = (hx(nx)+hx(0))/2.0;
	int cl, dm;
        rectidx(c, 0.0, cl, dm);
	for(int l=0; l<=nx; ++l)
	{
		if(fabs((c-hx(l)) - (hx(nx-l)-c))> 1.0e-7) return 0;
	}
	for(int l=0; l<=nx+1; ++l)
	{
		for(int m=0; m<=ny+1; ++m)
		{
			if(fabs(n(l,m)-n(nx-l+1,m)) > 1.0e-7) return 0;
		}
	}
	return cl;
}

/* output to FILE dat */
void Waveguide::write(FILE *dat)
{
	int i, j;
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, "Nx: %d, Hx: ", nx);
		for(i=0; i<=nx; ++i) fprintf(dat, "%g ", hx(i));
		fprintf(dat, "\nNy: %d, Hy: ", ny);
		for(i=0; i<=ny; ++i) fprintf(dat, "%g ", hy(i));
		fprintf(dat, "\n");
		fprintf(dat, "lambda: %g\n", lambda);
		fprintf(dat, "n:\n");
		for(i=nx+1; i>=0; --i) 
		{
			for(j=0; j<=ny+1; ++j) 
				fprintf(dat, "%g  ", n(i,j));
			fprintf(dat, "\n");
		}
		fprintf(dat, "eps:\n");
		for(i=nx+1; i>=0; --i) 
		{
			for(j=0; j<=ny+1; ++j) 
				fprintf(dat, "%g  ", eps(i,j));
			fprintf(dat, "\n");
		}
	}
	else
	{
		comment("nx", dat);
		fputint(nx, dat);
		comment("ny", dat);
		fputint(ny, dat);
		comment("hx[0..nx]", dat);
		hx.write(dat);
		comment("hy[0..ny]", dat);
		hy.write(dat);
		comment("lambda", dat);
		fputdouble(lambda, dat);
		comment("n[0..nx+1,0..ny+1]", dat);
		n.write(dat);
	}
	return;
}

/* input from FILE dat */
void Waveguide::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) 
	{
		wavegerror("Waveguide: read");
	}
	nx = fgetint(dat);
	ny = fgetint(dat);
	hx.read(dat);
	hy.read(dat);
	lambda = fgetdouble(dat);
	n.read(dat);
	return;
}