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

#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 MPth  30.0   // angle of incidence
#define Wavel 1.55   // vacuum wavelength  /mum
#define Pol   TE     // polarization of the incoming mode

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

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

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

/* 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, control output	to .m-file
	SegWgStruct fac = stdef();
	fac.plot('0', '0');

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

	// vQUEP setup
	QuepField qfld(fac,  
	               LEFT, Pol, 0, CC1, MPth,
	               cwx, NumM, cwz, NumM);

	// ... and solution
        qfld.quepsim_s(3, 0.005, 0.01);

	// modal output 
	Complex a;
	double r, phi;
	fprintf(stderr, "\nScattered fundamental mode amplitudes:\n");
	a = qfld.Aout(LEFT, TE, 0, 0.0)/qfld.Ain(LEFT, Pol, 0, 0.0);
        r = a.sqabs();
	phi = a.arg();
        fprintf(stderr, "TE: R = %g, phi = %g\n", r, phi);
	a = qfld.Aout(LEFT, TM, 0, 0.0)/qfld.Ain(LEFT, Pol, 0, 0.0);
        r = a.sqabs();
	phi = a.arg();
        fprintf(stderr, "TM: R = %g, phi = %g\n", r, phi);

	// power balance
	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);
	double pbTE, pfTE, pdTE, puTE;
	double pbTM, pfTM, pdTM, puTM;
        pbTE = qfld.Pout(LEFT, TE);
        pfTE = qfld.Pout(RIGHT, TE);
        pdTE = qfld.Pout(BOTTOM, TE);
        puTE = qfld.Pout(TOP, TE);
        pbTM = qfld.Pout(LEFT, TM);
        pfTM = qfld.Pout(RIGHT, TM);
        pdTM = qfld.Pout(BOTTOM, TM);
        puTM = qfld.Pout(TOP, TM);
        fprintf(stderr, "... per polarization:\n");
        fprintf(stderr, "TE: P_back = %g\n", pbTE);
        fprintf(stderr, "    P_forw = %g\n", pfTE);
        fprintf(stderr, "    P_down = %g\n", pdTE);
        fprintf(stderr, "    P_up   = %g\n", puTE);
        fprintf(stderr, "    Sum    = %g\n", pbTE+pfTE+pdTE+puTE);
        fprintf(stderr, "TM: P_back = %g\n", pbTM);
        fprintf(stderr, "    P_forw = %g\n", pfTM);
        fprintf(stderr, "    P_down = %g\n", pdTM);
        fprintf(stderr, "    P_up   = %g\n", puTM);
        fprintf(stderr, "    Sum    = %g\n", pbTM+pfTM+pdTM+puTM);
        fprintf(stderr, "Sum        = %g\n\n", pbTE+pfTE+pdTE+puTE+pbTM+pfTM+pdTM+puTM);

	// display window
	double dwl, dwr, dwb, dwt;
	dwl = -4.2;
	dwr =  3.2;
	dwb = -1.8;
	dwt =  1.5;

	// field viewer and field component plots (absolute values), .m-files 
	qfld.adjustphase(HZ);
        qfld.viewer(dwb, dwt, dwl, dwr, 150, 260, 't', 'e');
        qfld.plot(EX, MOD, dwb, dwt, dwl, dwr, 150, 260, 't', 'e');
        qfld.plot(EY, MOD, dwb, dwt, dwl, dwr, 150, 260, 't', 'e');
        qfld.plot(EZ, MOD, dwb, dwt, dwl, dwr, 150, 260, 't', 'e');
        qfld.plot(HX, MOD, dwb, dwt, dwl, dwr, 150, 260, 't', 'e');
        qfld.plot(HY, MOD, dwb, dwt, dwl, dwr, 150, 260, 't', 'e');
        qfld.plot(HZ, MOD, dwb, dwt, dwl, dwr, 150, 260, 't', 'e');

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