/* 
 * METRIC-Application:
 * Defect cavity in a short waveguide Bragg grating, QUEP simulation
 */

#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    2.00   // film refractive index
#define GPns    1.45   // substrate refractive index
#define GPt     0.2    // film thickness /mum
#define GPp     0.21   // reflector grating period /mum
#define GPs     0.11   // width of the holes /mum
#define GPd     0.6    // etching depth /mum (> GPt !)
#define GPN     4      // number of holes in each reflector 
#define GPw     0.2275 // width of the defect cavity /mum
#define Wavel   0.633  // vacuum wavelength /mum

#define CWbot   -3.0   // vertical computational window 
#define CWtop    3.0   // 
#define CWleft   3.0   // horizontal computational window,  
#define CWright  3.0   // distances from the outer dielectric interfaces 
#define NumMx    100   // number of spectral terms in the field discretization
#define NumMz    120   // 

#define DWx0    -2.0   // vertical display window
#define DWx1     2.0   //
#define DWz0     2.0   // horizontal display window, 
#define DWz1     2.0   // distances from the outer dielectric interfaces

/* short Bragg grating with central defect */
SegWgStruct defgrdef()
{
	Waveguide slab(1);
	slab.hx(1) = GPt;
	slab.hx(0) = 0.0;
	slab.n(2) = GPnc;
	slab.n(1) = GPnf;
	slab.n(0) = GPns;
	slab.lambda = Wavel;

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

	SegWgStruct gr(4*(GPN-1)+3);
	double z = 0.0; 
	int l = 0;
	gr(l) = slab; gr.hz(l) = z; ++l;
	for(int i=1; i<=GPN-1; ++i)
	{
		gr(l) = etched; z += GPs; gr.hz(l) = z; ++l;
		gr(l) = slab; z += GPp-GPs; gr.hz(l) = z; ++l;
	}
	gr(l) = etched; z += GPs; gr.hz(l) = z; ++l;
	gr(l) = slab; z += GPw; gr.hz(l) = z; ++l;
	for(int i=1; i<=GPN-1; ++i)
	{
		gr(l) = etched; z += GPs; gr.hz(l) = z; ++l;
		gr(l) = slab; z += GPp-GPs; gr.hz(l) = z; ++l;
	}
	gr(l) = etched; z += GPs; gr.hz(l) = z; ++l;
	gr(l) = slab; 
	return gr;
}

/* simulate the incidence of the fundamental guided mode */ 
int main()
{

	// the defect grating
	SegWgStruct gr = defgrdef();
	// a plot of the refractive index profile
	gr.plot('0', '0');

	// vertical computational window
        Interval vcw(CWbot, CWtop);
        // horizontal computational window
        Interval hcw(gr.hz(0)-CWleft, gr.hz(gr.nz)+CWright);

	// QUEP setup + spectral discretizaion
	QuepField qfld(gr, Pol, vcw, NumMx, hcw, NumMz);

        // input: the fundamental mode in the left channel
	qfld.input(LEFT, CC1, 0);

        // calculate the light propagation across the defect grating 
        qfld.quepsim();

	Complex a; 
	double p;
	fprintf(stderr, "\nScattered fundamental mode amplitudes:\n");
	a = qfld.Aout(LEFT, 0);
	p = qfld.Pout(LEFT, 0);
	fprintf(stderr, "R = %g + i %g,  |R|^2 = %g\n", a.re, a.im, p);
	a = qfld.Aout(RIGHT, 0);
	p = qfld.Pout(RIGHT, 0);
	fprintf(stderr, "T = %g + i %g,  |T|^2 = %g\n", a.re, a.im, p);

	double pb, pf, pd, pu;
        pb = qfld.Pout(LEFT);
        pf = qfld.Pout(RIGHT);
        pd = qfld.Pout(BOTTOM);
        pu = qfld.Pout(TOP);
        fprintf(stderr, "Total scattered power:\n");
        fprintf(stderr, "P_back      = %g\n", pb);
        fprintf(stderr, "P_forw      = %g\n", pf);
        fprintf(stderr, "P_down      = %g\n", pd);
        fprintf(stderr, "P_up        = %g\n", pu);
        fprintf(stderr, "Sum P_out   = %g\n", pb+pf+pd+pu);

        pb = qfld.Pgout(LEFT);
        pf = qfld.Pgout(RIGHT);
        fprintf(stderr, "Total scattered guided power:\n");
        fprintf(stderr, "P_back      = %g\n", pb);
        fprintf(stderr, "P_forw      = %g\n\n", pf);

	Fcomp fc = principalcomp(Pol);
	char pc = polchr(Pol);
	double dwl, dwr, dwb, dwt;
	dwb = DWx0;
	dwt = DWx1;
	dwl = gr.hz(0)-DWz0; 
	dwr = gr.hz(gr.nz)+DWz1;
	qfld.adjustphase(GPt/2.0, (gr.hz(gr.nz)-gr.hz(0))/2.0+GPw/4.0);
        qfld.viewer(dwb, dwt, dwl, dwr, 150, 150, 't', pc);
        qfld.plot(fc, REP, dwb, dwt, dwl, dwr, 150, 150, 't', pc);
        qfld.movie(dwb, dwt, dwl, dwr, 150, 150, 30, 't', pc);
        qfld.fplot(fc, dwb, dwt, dwl, dwr, 150, 150, 't', pc);
        qfld.fmovie(dwb, dwt, dwl, dwr, 150, 150, 30, 't', pc);

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