/* 
 * METRIC-Application:
 * Grating in a symmetric three layer waveguide, 
 * (infinite sequence of equally spaced identical dielectric squares),
 * parameters: J.D. Joannopoulos et al., "Photonic Crystals: 
 *             Molding the flow of light", Princeton University Press, 2008;
 * computation of Floquet-Bloch modes, BEP.
 */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"metric.h"

#define Pol     TE             // light polarization
#define GPnc    1.0            // cover refractive index
#define GPnf    (sqrt(12.0))   // film refractive index
#define GPp     1.0            // grating period /mum
#define GPt     0.4            // film thickness /mum
#define GPs     0.6            // width of the holes /mum
double  Wavel = 5.0;           // vacuum wavelength /mum

double  CWbot = -5.123;        // vertical computational window 
double  CWtop =  5.123;        // 
int     NumMx =  70;           // number of spectral terms in the field discretization,
#define NumMxMin 20            // adapted to a minimum to ensure convergence of the eigensolver 
#define NSDTPWL  20            // number of initial (spectral) discretiz. terms per wavelength

#define DWx0    -1.8           // vertical display window
#define DWx1     1.8           //
#define DWz0    -1.8           // horizontal display window, 
#define DWz1    11.8           //
#define NPx       90           // field plots: number of discretization points
#define NPz      300           //

                               // differnt configurations to be analyzed
#define NumP 5
double  Pval[NumP] = {0.1, 0.2, 0.2433, 0.4086, 0.4641} ;

/* periodic dielectric waveguide: one Period */
SegWgStruct grdef()
{
	Waveguide slab(1);
	slab.hx(1) =  GPt/2.0;
	slab.hx(0) = -GPt/2.0;
	slab.n(2) = GPnc;
	slab.n(1) = GPnf;
	slab.n(0) = GPnc;
	slab.lambda = Wavel;

	Waveguide etched(1);
	etched.hx(1) =  GPt/2.0;
	etched.hx(0) = -GPt/2.0;
	etched.n(2) = GPnc;
	etched.n(1) = GPnc;
	etched.n(0) = GPnc;
	etched.lambda = Wavel;

	SegWgStruct gr(3);
	// the period extends from gr.hz(0) to gr.hz(3)
	double z = 0.0; 
	int l = 0;
	gr(l) = slab; gr.hz(l) = z; ++l;
	gr(l) = slab; z += 0.5*(GPp-GPs); gr.hz(l) = z; ++l;
	gr(l) = etched; z += GPs; gr.hz(l) = z; ++l;
	gr(l) = slab; z += 0.5*(GPp-GPs); gr.hz(l) = z; ++l;
	gr(l) = slab; 
	return gr;
}

int main()
{
for(int v=0; v<=NumP-1; ++v)
{
	Wavel = GPp/Pval[v];
	fprintf(stderr, "\n\n --- lambda = %g,  Lambda/lambda = %g ---\n",
		Wavel, GPp/Wavel);

	// the grating, one period
	SegWgStruct gr = grdef();
	// prepare a plot 
	// SegWgStruct pgr = gr.unwrap(DWz0-1.0, DWz1+1.0);
	// pgr.write(stderr);
	// pgr.plot(DWx0, DWx1, DWz0, DWz1, '0', '0');
	// exit(1);

	
	// vertical computational window, adapt to vacuum wavelength
	double cwt = CWtop;
	double cwb = CWbot;
	if((cwt-cwb) < 10.0*Wavel)
	{	
		double cw0 = 0.5*(cwt+cwb);
		cwt = cw0+5.0123*Wavel;
		cwb = cw0-5.0123*Wavel;
	}	 
	NumMx = (int) ceil((cwt-cwb)/Wavel*((double)NSDTPWL));
	Interval vcw(cwb, cwt);

	// invoke the BEP-Floquet-Bloch solver
	FBMAXATTLVL = 0.1;
	FBVOPLVL = 0.0001;
	FBModeArray ma;
	fbmsolve(gr, Pol, vcw, NumMx, 1, ma, 0);

	// field plots and animations for one directional mode per configuration
	int done = 0;
	for(int j=0; j<=ma.num-1; ++j)
	{
		if(done==0 && ma(j).power() > 0.0)
		{
			fprintf(stderr, "[%d] P = %g  beta*Lambda/2/pi = %g Lambda/lambda = %g\n",
				j, ma(j).power(), ma(j).beta*GPp/2.0/PI, GPp/Wavel);
			ma(j).adjustphase();
			Fcomp fc = principalcomp(Pol);
			ma(j).viewer(DWx0, DWx1, DWz0, DWz1, NPx, NPz, dig10(v), dig1(v)); 
			ma(j).plot(fc, MOD, DWx0, DWx1, DWz0, DWz1, NPx, NPz, dig10(v), dig1(v));
			ma(j).plot(fc, REP, DWx0, DWx1, DWz0, DWz1, NPx, NPz, dig10(v), dig1(v));
			ma(j).movie(DWx0, DWx1, DWz0, DWz1, NPx, NPz, 30, dig10(v), dig1(v));
			done = 1;
		}
	}
}
        fprintf(stderr, "\nOk.\n");
	return 0;
}