/* 
 * METRIC-Application:
 * Propagation of Gaussian beams in a homogeneous dielectric medium
 * - crossing & interference of wide beams originating from the 
 *   left and lower border
 */

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

#define Pol    TE     // light polarization 
#define FPn    1.45   // refractive index of the medium
#define Wavel  1.55   // vacuum wavelength /mum

#define CW     20.0   // computational window, x- and z-extension 
#define NumM   50     // number of spectral terms in the field discretizations

#define DW     20.0   // display window, x- and z-extension

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

/* specification of the incoming Gaussian beams */
#define BPwidth   3.0   // half with 
#define BPx0     -1.0   // center position 
#define BPtheta  10.0   // tilt angle, in degrees 
Gaussianbeam Beam1(Pol, BPwidth,  BPx0,  BPtheta, Wavel, FPn);
Gaussianbeam Beam2(Pol, BPwidth, -BPx0, -BPtheta, Wavel, FPn);

/* the computational domain, homogeneous space */
Circuit circuitdef()
{
	Circuit cr;
	cr = Circuit(1, 1);
	cr.hx(0) = -CW/6.0;
	cr.hx(1) =  CW/6.0;
	cr.hz(0) = -CW/6.0;
	cr.hz(1) =  CW/6.0;
	cr.n(0, 0) = FPn; cr.n(0, 1) = FPn; cr.n(0, 2) = FPn;
	cr.n(1, 0) = FPn; cr.n(1, 1) = FPn; cr.n(1, 2) = FPn;
	cr.n(2, 0) = FPn; cr.n(2, 1) = FPn; cr.n(2, 2) = FPn;
	cr.lambda = Wavel;
	return cr;
}

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

/* QUEP simulation */
int main()
{
	fprintf(stderr, "\nPropagation of Gaussian beams:\n");	
	fprintf(stderr, "------------------------------\n");	
	fprintf(stderr, "lambda = %g mum\n", Wavel);
	fprintf(stderr, "n      = %g\n", FPn);
	fprintf(stderr, "w      = %g mum\n", BPwidth);
	fprintf(stderr, "x0     = %g mum\n", BPx0);
	fprintf(stderr, "theta  = %g^o\n", BPtheta);

	// the structure under investigation 
	Circuit fsp = circuitdef();
	// vertical and horizontal computational window
	Interval cw(-CW/2.0, CW/2.0);
	// QUEP setup 
	QuepField qfld(fsp, Pol, cw, NumM, cw, NumM);

	// input: the beams as defined above, at the left and lower borders  
	qfld.input(LEFT,   Beam1, CC1);
	qfld.input(BOTTOM, Beam2, CC1);
	// QUEP solution of the Helmholtz problem
	qfld.quepsim(); 

	double pl, pr, pb, pt;
	fprintf(stderr, "\nTotal incoming power:\n");
	pl = qfld.Pin(LEFT);
	pr = qfld.Pin(RIGHT);
	pb = qfld.Pin(BOTTOM);
	pt = qfld.Pin(TOP);
	fprintf(stderr, "P_left    = %g\n", pl);
	fprintf(stderr, "P_right   = %g\n", pr);
	fprintf(stderr, "P_bottom  = %g\n", pb);
	fprintf(stderr, "P_top     = %g\n", pt);
	fprintf(stderr, "Sum P_in  = %g\n", pl+pr+pb+pt);

	fprintf(stderr, "\nTotal outgoing power:\n");
	pl = qfld.Pout(LEFT);
	pr = qfld.Pout(RIGHT);
	pb = qfld.Pout(BOTTOM);
	pt = qfld.Pout(TOP);
	fprintf(stderr, "P_left    = %g\n", pl);
	fprintf(stderr, "P_right   = %g\n", pr);
	fprintf(stderr, "P_bottom  = %g\n", pb);
	fprintf(stderr, "P_top     = %g\n", pt);
	fprintf(stderr, "Sum P_out = %g\n", pl+pr+pb+pt);
	
	char pc = polchr(Pol);
	Fcomp fc = principalcomp(Pol);
	Mlop_Print = NO;
	Mlop_Colour = YES;

	qfld.viewer(-0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 't', pc);
	qfld.plot(fc, MOD, -0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 't', pc);
	qfld.movie(-0.5*DW,0.5*DW,-0.5*DW,0.5*DW, 250, 250, 30, 't', pc); 

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