/* 
 * METRIC-Application:
 * A three layer waveguide: guided mode analysis, 
 * influence of defining parameters on propagation constants/effective indices, 
 * perturbational evaluation 
 */

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

#define Wgpns   1.5      // substrate refractive index
#define Wgpnf   2.0      // film refractive index
#define Wgpnc   1.0      // cover: air
#define Wgpt    1.0      // slab thickness 
#define Wavel   1.55     // vacuum wavelength 

                         // shifts of these parameters 
#define Deltans   0.01   // Wgpns  
#define Deltanf   0.01   // Wgpnf  
#define Deltanc   0.01   // Wgpnc  
#define Deltat    0.01   // Wgpt
#define Deltawl   0.01   // Wavel  

/* calculate the guided modes, apply perturbational expressions 
   for shifts of propagation constants / effective mode indices 
   due to alterations of the waveguide parameters */
int main()
{
	// define the waveguide
	Waveguide wg(1);
	wg.n(0) = Wgpns;
	wg.n(1) = Wgpnf;
	wg.n(2) = Wgpnc;
	wg.hx(0) = 0.0;
	wg.hx(1) = Wgpt;
	wg.lambda = Wavel;

		
	// display window
	Interval disp(wg.hx(0)-2.0, wg.hx(wg.nx)+2.0);
	// inspect the refractive index profile
	wg.plot(disp, '0', '0');

	// find all guided modes
	ModeArray tema, tmma, ma;
	modeanalysis(wg, TE, tema);
	modeanalysis(wg, TM, tmma);
	// merge the ModeArray into one
	ma = tema; ma.merge(tmma); tema.clear(); tmma.clear();

	// the unperturbed modes 
	for(int i=0; i<=ma.num-1; ++i)
	{	
		Mode m = ma(i);
		// propagation constant and effective index
		fprintf(stderr, "%s: beta = %g/mum, neff = %g, ", 
			m.ids, m.beta, m.neff); 

		// a profile plot of the principal component
		Fcomp cp = principalcomp(m.pol);
		m.plot(cp, ORG, disp, 500, dig10(i), dig1(i), 'L');
	}
	fprintf(stderr, "\n"); 

	// investigate the effect of altered parameters,
	// first order perturbation theory, 
	enum Parms {NS, NF, NC, TH, WL} par;
	for(par = NS; par <= WL; par = Parms(par+1))
	{
		Interval r;
		Perturbation p;
		
		// specify the Perturbation's,
		// consider the precise definitions in gengwed.c !
		switch(par)
		{
		case NS:
			// the substrate refractive index
			fprintf(stderr, "n_s: shift %g\n", Deltans);
			r = wg.layer(0);
			p = refindshift(Wgpns, Deltans, r);
			p.write(stderr);
			break;
		case NF:
			// the refractive index of the film layer
			fprintf(stderr, "n_f: shift %g\n", Deltanf);
			r = wg.layer(1);
			p = refindshift(Wgpnf, Deltanf, r);
			p.write(stderr);
			break;
		case NC:
			// the cover refractive index
			fprintf(stderr, "n_c: shift %g\n", Deltanc);
			r = wg.layer(2);
			p = refindshift(Wgpnc, Deltanc, r);
			p.write(stderr);
			break;
		case TH:
			// the film thickness 
			fprintf(stderr, "t: shift %g\n", Deltat);
			break;
		case WL:
			// the vacuum wavelength 
			fprintf(stderr, "lambda: shift %g\n", Deltawl);
			break;
		}

		// for each mode, evaluate ...
		for(int i=0; i<=ma.num-1; ++i)
		{	
			Mode m = ma(i);
			fprintf(stderr, " %s:  ", m.ids); 

			// ... the resulting shifts of the propagation constant
			// and of the effective mode index
			Complex deltabeta = CC0;
			double deltaneff = 0.0;
			switch(par)
			{
			case NS: case NF: case NC:
				deltabeta = m.phaseshift(p, FORW);
				deltaneff = deltabeta.re/m.k0;
				break;
			case TH:
				deltabeta = m.phaseshift(1, Deltat, FORW);
				deltaneff = deltabeta.re/m.k0;
				break;
			case WL:
				deltabeta = m.phaseshift(Deltawl, FORW);
				deltaneff = m.neffshift(Deltawl, FORW).re;
				break;
			}
			fprintf(stderr, "beta, shift: (%g+i%g)/mum  ", 
				deltabeta.re, deltabeta.im); 	
			fprintf(stderr, "neff, shift: %g\n", deltaneff); 	
		}

		// for comparison: rigorous calculation,
		// substract beta and neff of two modes related 
		// to the original and the modified waveguide  
		Waveguide wgp = wg;
		switch(par)
		{
		case NS:
			wgp.n(0) += Deltans;	
			break;
		case NF:
			wgp.n(1) += Deltanf;	
			break;
		case NC:
			wgp.n(2) += Deltanc;	
			break;
		case TH:
			wgp.hx(1) += Deltat;	
			break;
		case WL:
			wgp.lambda += Deltawl;	
			break;
		}
		ModeArray map;
		modeanalysis(wgp, TE, tema, 1); 
		modeanalysis(wgp, TM, tmma, 1);
		map = tema; map.merge(tmma); tema.clear(); tmma.clear();
		fprintf(stderr, " rigorous calculation:\n"); 
		for(int i=0; i<=ma.num-1; ++i)
		{	
			Mode m = ma(i);
			Mode mp = map(i);
			fprintf(stderr, " %s:  ", mp.ids); 

			// ... the resulting shifts of the propagation constant
			double deltabeta = mp.beta-m.beta;
			fprintf(stderr, "beta, diff: (%g+i0)/mum  ", 
				deltabeta); 	

			// ... and of the effective mode index
			double deltaneff = mp.neff-m.neff;
			fprintf(stderr, "neff, diff: %g\n", deltaneff); 	
		}

		// for comparison: "rigorous" tangent evaluation,
		// approximate the derivative of beta or neff with respect to 
		// the parameter as a "finite difference",
		// multiply by the actual parameter shift
		wgp = wg;
		double fd = 1.0e-6;
		double deltap = 0.0;
		switch(par)
		{
		case NS:
			wgp.n(0) += fd;	
			deltap = Deltans;
			break;
		case NF:
			wgp.n(1) += fd;	
			deltap = Deltanf;
			break;
		case NC:
			wgp.n(2) += fd;	
			deltap = Deltanc;
			break;
		case TH:
			wgp.hx(1) += fd;	
			deltap = Deltat;
			break;
		case WL:
			wgp.lambda += fd;	
			deltap = Deltawl;
			break;
		default: ;;
		}
		modeanalysis(wgp, TE, tema, 1); 
		modeanalysis(wgp, TM, tmma, 1);
		map = tema; map.merge(tmma); tema.clear(); tmma.clear();
		fprintf(stderr, " tangent evaluation, FD:\n"); 
		for(int i=0; i<=ma.num-1; ++i)
		{	
			Mode m = ma(i);
			Mode mp = map(i);
			fprintf(stderr, " %s:  ", mp.ids); 

			// ... the resulting shifts of the propagation constant
			double deltabeta = (mp.beta-m.beta)/fd*deltap;
			fprintf(stderr, "beta, shift: (%g+i0)/mum  ", 
				deltabeta); 	

			// ... and of the effective mode index
			double deltaneff = (mp.neff-m.neff)/fd*deltap;
			fprintf(stderr, "neff, shift: %g\n", deltaneff); 	
		}
		fprintf(stderr, "\n"); 
	}

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