/* 
 * METRIC-Application:
 * A three layer waveguide, magnetooptic core, equatorial configuration,
 * magnetooptic phase shift (TM), dispersion curves
 */

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

#define Wgpns   1.95    // substrate refractive index
#define Wgpnf   2.302   // film refractive index
#define Wgpnc   1.0     // cover: air
double  Wgpt  = 0.5;    // slab thickness /mum 
#define Wavel   1.3     // vacuum wavelength /mum
#define WgpFrot 3000.0  // specific Faraday rotation of the core material
                        //                                     (in ^o/cm)

#define MaxOrd 100      // highest order of a recorded mode

#define Pbeg 0.01       // Parameter variations: minimum value,
#define Pend 2.0        //                       maximum value,
#define NumP 200        //                       number of points

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

/* waveguide definition */
Waveguide wgdef()
{
	Waveguide g(1);

	g.hx(0) = 0.0;
	g.hx(1) = Wgpt;
	g.n(0) = Wgpns;
	g.n(1) = Wgpnf;
	g.n(2) = Wgpnc;
	g.lambda = Wavel;
	return g;
}

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

/* calculate effective indices and the nonreciprocal effect 
   (the difference between propagation constants of modes traveling in the 
    forward and backward directions through the magnetooptic waveguide)
   versus the layer thickness parameter */
int main()
{
	ModeArray ma;
	char name[20] = "tm__neff";
	char psname[20] = "tm__nrps";

	Dmatrix neff(MaxOrd+1, NumP+1); 
	Dmatrix nrps(MaxOrd+1, NumP+1); 
	Dvector pval(NumP+1); 
	int omax = -1;
	neff.init(0.0);
	
	fprintf(stderr, "\nDispersion curves:\n");
	for(int pi=0; pi<=NumP; ++pi)
	{
		pval(pi) = Pbeg+((double)pi)/((double)NumP)*(Pend-Pbeg);
		Wgpt = pval(pi);	

		// define the waveguide
		Waveguide wg = wgdef();

		// the magnetooptic permittivity perturbation 
		// in the core layer
		Perturbation pert = magopt_equat(WgpFrot, wg.n(1), wg.lambda, 
		                                 wg.layer(1));

		// find the guided modes
		modeanalysis(wg, TM, ma, 1);
		
		// store effective mode indices, evaluate the perturbation
		for(int j=0; j<=ma.num-1; ++j)
		{
			int o = ma(j).nodes(); 
			if(o <= MaxOrd)
			{
				neff(o, pi) = ma(j).neff;	
				// nonreciprocal phase shift in 1/cm
				nrps(o, pi) = ma(j).phaseshift(pert, FORW).re
				              *2.0*10000.0; 
				if(o > omax) omax = o;
			}
		}
		int s = NumP/60;
                if(pi % s == 0) fprintf(stderr, ".");
	}
	fprintf(stderr, "\n");
	// write the data to files, individually per mode order, for ...
	Dvector oarg(NumP+1); 
	Dvector oval(NumP+1); 
	Dvector nval(NumP+1); 
	for(int o=0; o<=omax; ++o)
	{
		int num = 0;
		name[2] = dig10(o);
		name[3] =  dig1(o);
		psname[2] = dig10(o);
		psname[3] =  dig1(o);
		for(int pi=0; pi<=NumP; ++pi)
		{
			if(neff(o, pi) > 0.0)	
			{
				oarg(num) = pval(pi);
				oval(num) = neff(o, pi);
				nval(num) = nrps(o, pi);
				++num;
			}
		}
		if(num >= 1)
		{
		// ... effective mode indices
		toxyf(name, oarg.subvector(num, 0), oval.subvector(num, 0));
		// ... and magnetooptic phase shifts
		toxyf(psname, oarg.subvector(num, 0), nval.subvector(num, 0));
		}
	}
	fprintf(stderr, "Ok.\n");
	return 0;
}