/* 
 * METRIC-Application:
 * Oblique guided wave propagation towards a 90^o corner 
 * in a slab waveguide, vectorial 2-D problem, 
 * vQUEP solutions, incidence of semi-guided Gaussian wave bundles 
 */

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

#define MPng  3.4    // core refractive index
#define MPnb  1.45   // substrate/cladding refractive index
#define MPd   0.25   // core thickness /mum
#define Wavel 1.55   // vacuum wavelength  /mum
#define Pol   TE     // polarization of the incoming mode

#define MPth  41.    // incidence angle 
#define MPbw  13.22  // full field waist width of the beam /mum
#define MPy0  0.0    // focus coordinate y
#define MPz0  0.0    // focus coordinate z

                     // spectral discretization of the beam
#define MPNa   20    // number of angles
#define MPalim 0.01  // relative amplitudes: neglect below MPalim

#define CWxb  -3.1    // computational window: border positions
#define CWxt   3.1     
#define CWzl  -3.1     
#define CWzr   3.1     

                      // the discretized radiation field
#define NumM   150    // number of expansion terms for both directions

/* ------------------------------------------------------------------------ */

/* the corner configuration */
SegWgStruct stdef()
{
        Circuit cnr(1, 1);
        cnr.hx(0) = -MPd/2.0;
        cnr.hx(1) =  MPd/2.0;
        cnr.hz(0) = -MPd/2.0;
        cnr.hz(1) =  MPd/2.0;
        cnr.n(2, 0) = MPnb; cnr.n(2, 1) = MPng; cnr.n(2, 2) = MPnb;
        cnr.n(1, 0) = MPng; cnr.n(1, 1) = MPng; cnr.n(1, 2) = MPnb;
        cnr.n(0, 0) = MPnb; cnr.n(0, 1) = MPnb; cnr.n(0, 2) = MPnb;
        cnr.lambda = Wavel;

        SegWgStruct st(cnr);
        return st;
}

/* ------------------------------------------------------------------------ */

int main()
{
	// structure definition	
	SegWgStruct cnr = stdef();

	// computational window
        Interval cwx(CWxb, CWxt);
        Interval cwz(CWzl, CWzr);

	// container for a superposition of vQUEP solutions, 
	// specification of the incoming beam
	Bundle wave(cnr, LEFT, Pol, 0, 
	            MPth, MPbw, MPy0, MPz0, MPNa, MPalim, 
	            cwx, NumM, cwz, NumM);

	// containers for field plots, field discretization parameters
	double dwx0, dwx1, dwy0, dwy1, dwz0, dwz1;
	int ndpx, ndpy, ndpz;
	dwx0 =  -2.0; dwx1 =  15.0; ndpx = 200;
	dwy0 = -15.0; dwy1 =  15.0; ndpy = 300;
	dwz0 = -15.0; dwz1 =  2.0;  ndpz = 300;
	View vw;

	// longitudinal view, cross section
	vw = View(XZ, cnr, 
	         -1.7, 5.2, 200, 
	         -5.2, 1.7, 200, 0.0,
	         'l', 'g');
	wave.VA.add(vw);
	
	// bottom view, at the level of the horizontal core
	vw = View(YZ, cnr, 
	         dwy0, dwy1, ndpy, 
	         dwz0, dwz1, ndpz, 0.0,
	         'h', 'r');
	wave.VA.add(vw);

	// view of a plane parallel to the facet, 
	// at the center of the vertical core 
	vw = View(XY, cnr, 
	         dwx0, dwx1, ndpx, 
	         dwy0, dwy1, ndpy, 0.0, 
	         'c', 'r');
	wave.VA.add(vw);

	// containers for recording guided output powers
	Port pt;
	pt = Port(LEFT, TE, 0, 'R'); wave.PA.add(pt);
	pt = Port(LEFT, TM, 0, 'R'); wave.PA.add(pt);
	pt = Port(TOP,  TE, 0, 'T'); wave.PA.add(pt);
	pt = Port(TOP,  TM, 0, 'T'); wave.PA.add(pt);

	// carry out a series of vQUEP simulations for varying angles of incidence
        // accumulate fields as specified in the views
	wave.solve(5, 0.005, 0.01);
	
	// modal output
	fprintf(stderr, "\nRelative guided modal output power:\n");
	double sumP = 0.0;
	for(int j=0; j<=wave.PA.num-1; ++j)
	{
		fprintf(stderr, "%c_T%c%d = %g\n", 
			wave.PA(j).id, 
			polCHR(wave.PA(j).pol), 
			wave.PA(j).ord, 
			wave.PA(j).Pout/wave.Pin);
		sumP += wave.PA(j).Pout/wave.Pin;
	}
	fprintf(stderr, "  Sum = %g\n", sumP);

	// viewer and field plots (absolute values), output to .m-files 
	fprintf(stderr, "\n");
	for(int j=0; j<=wave.VA.num-1; ++j)
	{
		wave.VA(j).viewer();
		wave.VA(j).plot(mW);
		wave.VA(j).plot(mE);
		wave.VA(j).plot(mH);
	}
        fprintf(stderr, "\nOk.\n");
	return 0;
}