/* 
 * METRIC-Application:
 * A short waveguide Bragg grating, parameter scans 
 */

#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.200  // film thickness /mum
double  GPp;           // grating period /mum
double  GPs;           // width of the holes /mum
#define GPd     0.6    // etching depth /mum (> GPt !)
#define GPN     4      // number of holes 
#define Wavel   0.633  // vacuum wavelength /mum

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

#define Pmin     0.1   // parameter scan: minimum,
#define Pmax     0.3   //                 maximum,
#define NumP     20    //                 number of parameter values

/* 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; ++pi)
{
	double pv = Pmin+((double) pi)/((double) NumP)*(Pmax-Pmin);
	GPp = pv;
	GPs = pv/2.0;
	fprintf(stderr, "\n[%d]:  p = %g,  s = %g\n", pi, GPp, GPs); 

	// the grating
	SegWgStruct gr = grdef();

	// 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 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", pf);
	fprintf(stderr, "Loss        = %g\n", 1.0-pb-pf);

	apptoxyf("tr", pv, pf);
	apptoxyf("rf", pv, pb);
	apptoxyf("lo", pv, 1.0-pf-pb);
}
        fprintf(stderr, "\nOk.\n");
	return 0;
}