/* 
 * METRIC-Application:
 * Plane wave reflection from a dielectric multilayer stack,
 * defect grating, configuration with a transmission resonance
 */

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

#define Pol        TE    // light polarization
#define SPnb     2.25    // backgound refractive index
#define SPng     4.50    // refractive index of high index regions 
double  SPg   = 0.125;   // gap thickness 
double  SPt   = 0.125;   // thickness of the high index regions 
int     SPN   = 19;      // number of periods in the reflection grating
double  Wavel = 2*PI/3.8599162;     // vacuum wavelength 
double  Theta =  0.0;    // incidence angle, in degrees

#define DWx    2.0    // display window x, distances from the outer boundaries 
#define DWz    1.0    // display window z, half width 

/* layer stack definition */
Waveguide stdef()
{
	Waveguide g(2*SPN+1);
	int l = 0;
	g.n(l)  = SPnb;
	g.hx(l) =  0.0;
	++l;
	for(int p=1; p<=SPN; ++p)
	{
		g.n(l)  = SPng;
		g.hx(l) = g.hx(l-1)+SPt;
		++l;
		g.n(l)  = SPnb;
		g.hx(l) = g.hx(l-1)+SPg;
		++l;
	}
	g.n(l)  = SPng;
	g.hx(l) = g.hx(l-1)+SPt;
	++l;
	g.n(l)  = SPnb;
	g.n(SPN+1) = SPng;
	g.lambda = Wavel;
	return g;
}

/* calculate the plane wave reflection & transmission */
int main()
{
	ModeArray fld;
	Cvector amp(2);
	double in, ref, trans; 
	Fcomp fc = principalcomp(Pol);
		
	// define the layer stack
	Waveguide st = stdef();

	// solve the problem
	in = 1.0;
	mlref(st, Pol, Theta*PI/180.0, in, ref, trans, fld, amp);

	fprintf(stderr, "I = %g\n", in);
	fprintf(stderr, "R = %g\n", ref);
	fprintf(stderr, "T = %g\n", trans);

	// field plots 
	Interval disp(st.hx(0)-DWx, st.hx(2*SPN+1)+DWx);
	fld(0).plot(fc, ORG, disp, 500, 'r', 'e', 'L');
	fld(1).plot(fc, ORG, disp, 500, 'i', 'm', 'V');

	Cvector bamp(amp.nel);
	bamp.init(CC0);
	fld.plot(amp, bamp, 
	         fc, REP, disp.x0, disp.x1, -DWz, DWz, 150, 150, '0', '0');
	fld.plot(amp, bamp, 
	         fc, IMP, disp.x0, disp.x1, -DWz, DWz, 150, 150, '0', '0');
	fld.plot(amp, bamp, 
	         SX, REP, disp.x0, disp.x1, -DWz, DWz, 150, 150, '0', '0');
	fld.plot(amp, bamp,
	         SZ, REP, disp.x0, disp.x1, -DWz, DWz, 150, 150, '0', '0');
	fld.movie(amp, bamp, 
	          disp.x0, disp.x1, -DWz, DWz, 300, 50, 30, '0', '0');
	fld.viewer(amp, bamp, 
	           disp.x0, disp.x1, -DWz, DWz, 300, 50, '0', '0');
	
	fprintf(stderr, "\nOk.\n");
	return 0;
}