/* 
 * METRIC-Application:
 * Oblique guided wave propagation towards a 90^o corner 
 * in a slab waveguide, vQUEP solution 
 */

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

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

                      // the discretized radiation field
#define NumM   300    // 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, control output to .m-file
	SegWgStruct cnr = stdef();
	cnr.plot('0', '0');

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

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

	// and solution
        qfld.quepsim_s(5, 0.003, 0.01);

	// modal output
	Complex a;
	double r, phi;
	double xr = -MPd/2.0;
	double zr = -MPd/2.0;
	fprintf(stderr, "\nScattered fundamental mode amplitudes:\n");
	a = qfld.Aout(LEFT, TE, 0, zr)/qfld.Ain(LEFT, Pol, 0, zr);
	r = a.sqabs();
	phi = a.arg();
	fprintf(stderr, "Left: TE: R = %g, phi = %g\n", r, phi);
	a = qfld.Aout(LEFT, TM, 0, zr)/qfld.Ain(LEFT, Pol, 0, zr);
	r = a.sqabs();
	phi = a.arg();
	fprintf(stderr, "      TM: R = %g, phi = %g\n\n", r, phi);
	a = qfld.Aout(TOP, TE, 0, xr)/qfld.Ain(LEFT, Pol, 0, zr);
        r = a.sqabs();
	phi = a.arg();
	fprintf(stderr, "Top:  TE: T = %g, phi = %g\n", r, phi);
	a = qfld.Aout(TOP, TM, 0, xr)/qfld.Ain(LEFT, Pol, 0, zr);
        r = a.sqabs();
	phi = a.arg();
	fprintf(stderr, "      TM: T = %g, phi = %g\n\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 = -2.7;
	dwr =  1.7;
	dwb = -1.7;
	dwt =  2.7;

	// field viewer and static field plots, .m-files, 
	// absolute values of electric and magnetic fields, energy density 
	int ndpx, ndpz;
	ndpx = 180;
	ndpz = 180;
	char id1, id2;
	id1 = 't';
	id2 = polchr(Pol);
        qfld.viewer(dwb, dwt, dwl, dwr, ndpx, ndpz, id1, id2);
        qfld.plot(mW, dwb, dwt, dwl, dwr, ndpx, ndpz, id1, id2);
        qfld.plot(mE, dwb, dwt, dwl, dwr, ndpx, ndpz, id1, id2);
        qfld.plot(mH, dwb, dwt, dwl, dwr, ndpx, ndpz, id1, id2);

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