/* 
 * METRIC-Application:
 * A circular step-index optical fiber, 
 * exact TE, TE, OAM modes and approximate LP modes
 */

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

#define Pnb    1.44   // cladding refractive index 
#define Png    1.45   // core refractive index 
#define PR     10     // core radius /mum 

#define Wavel  1.55   // vacuum wavelength

#define DW     21.0   // display window, x&y-extension [-DW, DW]

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

/* radial dielectric profile, interface positions relative to the core radius */
Waveguide fibdef()
{
	Waveguide fib(0);
	fib.hx(0) =  0.0; 
	fib.n(1) = Pnb;
	fib.n(0) = Png;
	fib.lambda = Wavel;
	return fib;
}

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

/* fiber mode analysis */
int main()
{
	FIModeArray fim;
	Waveguide fib = fibdef();
	Rect disp(-DW, -DW, DW, DW);

	// TE, TM, and OAM modes
 	fimsolve(fib, PR, 0.0, 0.0, 0, 5, fim, 0);
 	// for(int j=0; j<=fim.num-1; ++j) fim(j).reverse();
 	for(int j=0; j<=fim.num-1; ++j) fim(j).discretize(0.0, PR+DW, 500);
	for(int j=0; j<=fim.num-1; ++j)
	{
                fprintf(stderr, "\n[%d] %s: beta = %g, neff = %g, B = %g\n", j, fim(j).ids, fim(j).beta, fim(j).neff, fim(j).npcB);
		fim(j).plot(ER, MOD, disp, 350, 350, dig10(j), dig1(j));
		fim(j).plot(ET, MOD, disp, 350, 350, dig10(j), dig1(j));
		fim(j).plot(EZ, MOD, disp, 350, 350, dig10(j), dig1(j));
		fim(j).viewer(disp, 350, 350, dig10(j), dig1(j));
	}

/*
	// weak guidance approximation, LP modes (angular dependence: complex exponentials) 
 	fimsolveLP(fib, PR, 0.0, 0.0, 0, 5, fim, 0);
 	// for(int j=0; j<=fim.num-1; ++j) fim(j).reverse();
 	for(int j=0; j<=fim.num-1; ++j) fim(j).discretize(0.0, PR+DW, 500);
	for(int j=0; j<=fim.num-1; ++j)
	{
                fprintf(stderr, "\n[%d] %s: beta = %g, neff = %g, B = %g\n", j, fim(j).ids, fim(j).beta, fim(j).neff, fim(j).npcB);
		if(fim(j).pid == 'x') fim(j).plot(EX, MOD, disp, 350, 350, dig10(j), dig1(j));
		else                  fim(j).plot(EY, MOD, disp, 350, 350, dig10(j), dig1(j));
		fim(j).viewer(disp, 350, 350, dig10(j), dig1(j));
	}
*/	

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