/* 
 * METRIC-Application:
 * Crossing between perpendicular waveguides,
 * modeling by hybrid analytical /  numerical coupled mode theory
 */

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

#define Pol    TE
#define CPnb   1.45   // background refractive index 
#define CPng   3.40   // waveguide cores: refractive index 
#define CPw    0.20   // width of the horizontal waveguide /mum 
#define CPv    0.45   // width of the vertical waveguide /mum 
#define Wavel  1.55   // vacuum wavelength /mum

#define DWz    3.0    // display window
#define DWx    3.0    // (symmetrical in x, z, distance from the center) 

// HCMT basis fields, parameters of the FEM discretization of amplitudes
#define FEMstep 0.025
#define Z0     (-CPv/2.0-1.5)
#define ZN     ( CPv/2.0+1.5)
#define X0     (-CPw/2.0-1.5)
#define XN     ( CPw/2.0+1.5)

// evaluation of HCMT integrals: computational window
#define CWz0   -4.0   
#define CWz1    4.0 
#define CWx0   -4.0 
#define CWx1    4.0 

#define Nispwl 4   // number of integration steps (10 pt Gaussian quadrature)
                   // per wavelength

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

/* a three layer slab waveguide of thickness t */
Waveguide wgdef(double t)
{
	Waveguide g(1);
	g.hx(0) = -t/2.0;
	g.hx(1) =  t/2.0;
	g.n(0) = CPnb;
	g.n(1) = CPng;
	g.n(2) = CPnb;
	g.lambda = Wavel;
	return g;
}

/* the waveguide crossing */
Circuit circuitdef()
{
	Circuit cr;
	cr = Circuit(1, 1);
	cr.hx(0) = -CPw/2.0;
	cr.hx(1) =  CPw/2.0;
	cr.hz(0) = -CPv/2.0;
	cr.hz(1) =  CPv/2.0;
	cr.n(0, 0) = CPnb; cr.n(0, 1) = CPng; cr.n(0, 2) = CPnb;
	cr.n(1, 0) = CPng; cr.n(1, 1) = CPng; cr.n(1, 2) = CPng;
	cr.n(2, 0) = CPnb; cr.n(2, 1) = CPng; cr.n(2, 2) = CPnb;
	cr.lambda = Wavel;
	return cr;
}

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

/* crossing simulation */
int main()
{

	// MCROP = 1.0e-6;

	// the full straucture
	SegWgStruct cross = SegWgStruct(circuitdef());
	// HCMT field container, initialization
	HcmtField fld(cross, Pol, Wavel, 
	            Interval(CWx0, CWx1), Interval(CWz0, CWz1), Nispwl);

	// computaion of the HCMT basis fields
	Waveguide wg;
	HcmtElement el;
	ModeArray mh, mv;

	// modes of the horizontal core
	wg = wgdef(CPw);
	modeanalysis(wg, Pol, mh);
	// remember the indices in the array of basis elements
	Ivector mhf_eli(mh.num);
	Ivector mhb_eli(mh.num);

	for(int m=0; m<=mh.num-1; ++m) 
	{
		// add basis fields: 
 		// unidirectional modes one by one,  
		// el = HcmtElement(MHF, mh(m), Z0, ZN, FEMstep);
		// el.ain = CC0; if(m==0) el.ain = CC1;
		// mhf_eli(m) = fld.addcf(el);
		// el = HcmtElement(MHB, mh(m), Z0, ZN, FEMstep);
		// el.ain = CC0; ;
		// mhb_eli(m) = fld.addcf(el);

		// or fields for bidirectional propagation in one step 
		fld.addchannel(MHF, mh(m), Z0, ZN, FEMstep, 
		               mhf_eli(m), mhb_eli(m));
		fld(mhf_eli(m)).ain = CC0;
		if(m==0) fld(mhf_eli(m)).ain = CC1;
		fld(mhb_eli(m)).ain = CC0;
	}

	// modes of the vertical core
	wg = wgdef(CPv);
	modeanalysis(wg, Pol, mv);
	// remember the indices in the array of basis elements
	Ivector mvu_eli(mv.num);
	Ivector mvd_eli(mv.num);
	for(int m=0; m<=mv.num-1; ++m) 
	{
		// add basis fields: 
 		// unidirectional modes one by one,  
		// el = HcmtElement(MVU, mh(m), Z0, ZN, FEMstep);
		// el.ain = CC0; 
		// mvu_eli(m) = fld.addcf(el);
		// el = HcmtElement(MVD, mh(m), Z0, ZN, FEMstep);
		// el.ain = CC0; ;
		// mvd_eli(m) = fld.addcf(el);

		// or fields for bidirectional propagation in one step 
		fld.addchannel(MVU, mv(m), X0, XN, FEMstep, 
		               mvu_eli(m), mvd_eli(m));
		fld(mvu_eli(m)).ain = CC0;
		fld(mvd_eli(m)).ain = CC0;
	}

	// solve the problem according to the HCMT scheme
	fld.solve();

	
	// maximum amplitude in time snapshot plots of the principal field
	fld.adjustphase(0.0, 0.0);
	Fcomp pc = principalcomp(Pol);

	fprintf(stderr, "\n");
        Complex f = fld.field(pc, 0.0, 0.0);
        fprintf(stderr, "field(%c%c, 0, 0) = %g + i %g\n\n", 
		fldchr(pc), cpchr(pc), f.re, f.im);

	fprintf(stderr, "\n");
	double tp = 0.0, p;
	fprintf(stderr, "Squared output amplitudes: \n");
	for(int j=0; j<=mh.num-1; ++j)
	{
		p = (fld(mhb_eli(j)).aout).sqabs(); tp += p;
		fprintf(stderr, " B%d: %g\n", j, p);
		p = (fld(mhf_eli(j)).aout).sqabs(); tp += p;
		fprintf(stderr, " F%d: %g\n", j, p);
	}
	for(int j=0; j<=mv.num-1; ++j)
	{
		p = (fld(mvd_eli(j)).aout).sqabs(); tp += p;
		fprintf(stderr, " D%d: %g\n", j, p);
		p = (fld(mvu_eli(j)).aout).sqabs(); tp += p;
		fprintf(stderr, " U%d: %g\n", j, p);
	}
	fprintf(stderr, "Sum: %g\n", tp);

	// amplitude functions -> data files
	fprintf(stderr, "\n");
	for(int j=0; j<=mh.num-1; ++j)
	{
		fld(mhf_eli(j)).amp2files(dig10(j), dig1(j));
		fld(mhb_eli(j)).amp2files(dig10(j), dig1(j));
	}
	for(int j=0; j<=mv.num-1; ++j)
	{
		fld(mvu_eli(j)).amp2files(dig10(j), dig1(j));
		fld(mvd_eli(j)).amp2files(dig10(j), dig1(j));
	}

	fprintf(stderr, "\n");
	fld.plot(pc, REP, -DWx, DWx, -DWz, DWz, 200, 200, '0', '0');
	fld.viewer(-DWx, DWx, -DWz, DWz, 200, 200, '0', '0');

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