/* 
 * 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;
 * band structure analysis, 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 = 0.09;          // 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 Pmin   0.01            // Parameter scan: range and step size
#define Pmax   0.50            // 
#define Pstep  0.005           //  

/* 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(double pv = Pmin; pv <= Pmax; pv += Pstep)
{
	Wavel = GPp/pv;
	fprintf(stderr, "\n [lambda = %g, Lambda/lambda = %g] \n", Wavel, GPp/Wavel);

	// the grating, one period
	SegWgStruct gr = grdef();
	
	// 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.01;
	FBVOPLVL = 1.0e-4;
	FBModeArray ma;
	fbmsolve(gr, Pol, vcw, NumMx, 1, ma, 0);

	// data output: FB wavenumbers and directional power levels
	for(int j=0; j<=ma.num-1; ++j)
	{
		apptoxyf("bands", ma(j).beta*GPp/2.0/PI, GPp/Wavel);
		if(ma(j).beta > 0)
		{
			apptoxyf("pf", GPp/Wavel, ma(j).fld.Pout(SBRIGHT));
			apptoxyf("pb", GPp/Wavel, ma(j).fld.Pout(SBLEFT));
		}
	}

}
        fprintf(stderr, "\nOk.\n");
	return 0;
}