/* 
 * METRIC-Application:
 * A short asymmetric waveguide Bragg grating (parameter scans), 
 * vEIM, conventional EIM, and QUEP simulations 
 */

#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   // grating period /mum
#define GPs     0.11   // width of the holes /mum
#define GPd     0.6    // etching depth /mum (> GPt !)
#define GPN     6      // number of holes
double  Wavel = 0.633; // vacuum wavelength /mum

/* scan options (adjust main loop and output commands in the main program): */

// a single wavelength
#define NumP     1     
double Pval[NumP] = {0.633};

// a series of predefined wavelengths
// #define NumP     3     
// double Pval[NumP] = {0.40, 0.47, 0.60};

// scan over a wavelength interval
// #define NumP 300
// #define Pbeg 0.3 
// #define Pend 0.8

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

/* QUEP simulations, parameters */
#define CWbot   -1.5   // vertical computational window,
#define CWtop    1.5   //
#define CWleft   2.0   // horizontal computational window,
#define CWright  2.0   // distances from the outer dielectric interfaces
#define NumMx    100   // number of spectral terms in the field discretization
#define NumMz    130   //

/* EIM, effective index for hole regions */
#define H_Neff   1.0

/* short Bragg grating */
SegWgStruct grdef()
{
	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(2*GPN-1);
	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; 
	return gr;
}

/* simulate the incidence of the fundamental guided mode */ 
int main()
{
for(int pi=0; pi<=NumP-1; ++pi)
{
	
	double pv = Pval[pi];
//	double pv = Pbeg + ((double)pi)/((double)NumP)*(Pend-Pbeg);
	
	Wavel = pv;
	fprintf(stderr, "\n[%d]:  Wavel = %g\n", pi, Wavel); 

	// the grating
	SegWgStruct gr = grdef();

	// vEIM 

	EimField fld(gr, Pol, Wavel);
	fld.solve(0);

	apptoxyf("vtr", pv, fld.trans);
	apptoxyf("vrf", pv, fld.ref);
 
	fld.adjustphase(GPt/2.0, -1.0);
        double zbeg = gr.hz(0)-DWz0;
        double zend = gr.hz(gr.nz)+DWz1;
	Fcomp cp = principalcomp(Pol);
	fld.plot(cp, REP, DWx0, DWx1, zbeg, zend, 100, 300, 'v', 'a'+pi);
	fld.viewer(DWx0, DWx1, zbeg, zend, 100, 300, 'v', 'a'+pi);

	// EIM, conventional 

        fprintf(stderr, "\nEIM:\n");

	Waveguide rs = gr(0);
	ModeArray ma;
	modeanalysis(rs, Pol, ma);
	if(ma.num < 1) exit(1);
	Mode xp = ma(0);
	// xp.plot('x', 'p');

	Waveguide hp(gr.nz);
	hp.hx = gr.hz;
	hp.lambda = gr(0).lambda;
	hp.special = 1;
	Circuit grc = gr.circuit();
	double s_nef = xp.neff;
	double h_nef = H_Neff;

	for(int s=0; s<=gr.nz+1; ++s)
	{
		if(gr(s).equal(gr(0))) hp.n(s) = s_nef*s_nef;
		else                   hp.n(s) = h_nef*h_nef;
	}
	fprintf(stderr, "eps_eff_s: %g\n", hp.n(0));
	// hp.plot('e', 'f');

	double in, ref, trans;
	ModeArray efld;
	Cvector amp(2);
	in = 1.0;
	mlref(hp, Pol, 0.0, in, ref, trans, efld, amp);

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

	apptoxyf("etr", pv, trans);
	apptoxyf("erf", pv, ref);

	// QUEP 

        Interval vcw(CWbot, CWtop);
        Interval hcw(gr.hz(0)-CWleft, gr.hz(gr.nz)+CWright);
	QuepField qfld(gr, Pol, vcw, NumMx, hcw, NumMz);
	qfld.input(LEFT, CC1, 0);
        int r;
	r = qfld.quepsim_s(5, 0.01, 0.0097);
	if(r != 0)
	{
		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 = |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 = |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", pf);
		fprintf(stderr, "Loss        = %g\n", 1.0-pb-pf);

		apptoxyf("qtr", pv, pf);
		apptoxyf("qrf", pv, pb);
		apptoxyf("qlo", pv, 1.0-pf-pb);
	
		qfld.adjustphase(GPt/2.0, -1.0);
		qfld.plot(cp, REP, DWx0, DWx1, zbeg, zend, 100, 300, 
		          'q', 'a'+pi);
		qfld.viewer(DWx0, DWx1, zbeg, zend, 100, 300, 'q', 'a'+pi);
	}
}
        fprintf(stderr, "\nOk.\n");
	return 0;
}