/* 
 * METRIC-Application:
 * Oblique guided wave propagation towards a slab facet,
 * vectorial 2-D problem including out-of plane losses, 
 * vQUEP solutions, incidence of semi-guided Gaussian wave bundles 
 */

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

#define MPnf  2.0    // core refractive index
#define MPns  1.5    // substrate refractive index
#define MPnc  1.0    // cladding refractive index
#define MPd   0.5    // core thickness /mum
#define Wavel 1.55   // vacuum wavelength  /mum
#define Pol   TE     // polarization of the incoming mode

#define MPth  30.0   // primary incidence angle 
#define MPbw  10.0   // 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  -5.1    // computational window: border positions
#define CWxt   5.1     
#define CWzl  -8.1     
#define CWzr   8.1     

                      // the discretized radiation field
#define NumMx   100   // number of expansion terms, x-direction
#define NumMz   150   // number of expansion terms, z-direction

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

/* the facet configuration */
SegWgStruct stdef()
{
	Circuit fac(1, 1);
	fac.hx(0) = 0.0;
	fac.hx(1) = MPd;
	fac.hz(0) = -MPd;
	fac.hz(1) = 0.0;
	fac.n(2, 0) = MPnc; fac.n(2, 1) = MPnc; fac.n(2, 2) = MPnc;
	fac.n(1, 0) = MPnf; fac.n(1, 1) = MPnf; fac.n(1, 2) = MPnc;
	fac.n(0, 0) = MPns; fac.n(0, 1) = MPns; fac.n(0, 2) = MPns;
	fac.lambda = Wavel;

	SegWgStruct st(fac);
	return st; 
}

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

int main()
{
	
	// structure definition
	SegWgStruct fac = 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(fac, LEFT, Pol, 0, 
	            MPth, MPbw, MPy0, MPz0, MPNa, MPalim, 
	            cwx, NumMx, cwz, NumMz);

	// containers for field plots, field discretization parameters
	double dwx0, dwx1, dwy0, dwy1, dwz0, dwz1;
	int ndpx, ndpy, ndpz;
	double plds = 0.1;
	dwx0 =  -1.8; dwx1 =   1.5; ndpx = ((int) floor((dwx1-dwx0)/plds))*5;
	dwy0 = -20.0; dwy1 =  20.0; ndpy = ((int) floor((dwy1-dwy0)/plds));
	dwz0 = -20.0; dwz1 =  20.0; ndpz = ((int) floor((dwz1-dwz0)/plds));
	View vw;

	// longitudinal view, cross section
	vw = View(XZ, fac, 
	         dwx0, dwx1, ndpx, 
	         dwz0, dwz1, ndpz, 0.0,
	         'l', 'g');
	wave.VA.add(vw);
	
	// bottom view, at a level in the substrate 
	vw = View(YZ, fac, 
	         dwy0, dwy1, ndpy, 
	         dwz0, dwz1, ndpz, -MPd/2.0,
	         'l', 'o');
	wave.VA.add(vw);

	// bottom view, at the core center 
	vw = View(YZ, fac, 
	         dwy0, dwy1, ndpy, 
	         dwz0, dwz1, ndpz, MPd/2.0,
	         'h', 'r');
	wave.VA.add(vw);

	// bottom view, at a level in the cladding
	vw = View(YZ, fac, 
	         dwy0, dwy1, ndpy, 
	         dwz0, dwz1, ndpz, MPd+MPd/2.0,
	         't', 'p');
	wave.VA.add(vw);

	// view of a plane parallel to the facet, at a position in the core
	vw = View(XY, fac, 
	         dwx0, dwx1, ndpx, 
	         dwy0, dwy1, ndpy, -MPd/2.0, 
	         'c', '0');
	wave.VA.add(vw);

	// view of a plane parallel to the facet at a position beyond the core
	vw = View(XY, fac, 
	         dwx0, dwx1, ndpx, 
	         dwy0, dwy1, ndpy,  MPd/2.0, 
	         'c', '1');
	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);
	
	// carry out a series of vQUEP simulations for varying angles of incidence
	// accumulate fields as specified in the views
	wave.solve(3, 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;
}