/* 
 * WMM-Application:
 * simple strip waveguide, perturbation theory:
 * - absorbing core 
 * - magnetooptic core 
 * - geometry variations
 */

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

/* waveguide parameters */
#define Wgpns 1.95  // substrate refractive index
#define Wgpnf 2.302 // film refractive index
#define Wgpnc 1.0   // cover: air
#define Wgpw  1.0   // strip width
#define Wgph  0.5   // strip height
#define Wgpl  1.3   // vacuum wavelength 

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

	g.hx(0) = 0.0;
	g.hx(1) = Wgph;
	g.hy(0) = -Wgpw/2.0;
	g.hy(1) =  Wgpw/2.0;

	g.n(0,0) = Wgpns;
	g.n(0,1) = Wgpns;
	g.n(0,2) = Wgpns;
	g.n(1,0) = Wgpnc;
	g.n(1,1) = Wgpnf;
	g.n(1,2) = Wgpnc;
	g.n(2,0) = Wgpnc;
	g.n(2,1) = Wgpnc;
	g.n(2,2) = Wgpnc;

	g.lambda = Wgpl;

	return g;
}

/* WMM analysis parameters */
WMM_Parameters pardef()
{
	WMM_Parameters p;

	p.vform = HXHY;
	p.vnorm = NRMMH;
	p.ccomp = CCALL;

	p.ini_d_alpha   = 0.05;
	p.ini_N_alpha   = 10;
	p.ini_alpha_max = 2.0;

	p.ini_steps = 10;
	p.ref_num   = 5;
	p.ref_exp   = 4.0;
	p.ref_sdf   = 0.5;

	p.fin_d_alpha   = 0.01;
	p.fin_N_alpha   = 30;
	p.fin_alpha_max = 3.0;

	p.btol   = 1.0e-7;
	p.mshift = 1.0e-10;

	return p;
}

/* perturbational analysis  based on fundamental QTE and QTM modes */
int main()
{
	WMM_Parameters par = pardef();
	WMM_ModeArray ma;
	WMM_Mode tem, tmm;
	Perturbation p, p1, p2;
	Complex db, db1, db2;
	double grad;
		
	// the basic waveguide 
	Waveguide wg = wgdef();

	// calculate and save the unperturbed semivectorial modes
	WMM_findfundmode(wg, QTE, SYM, 0.0, 0.0, par, '-', '-', ma);
	tem = ma(0);
	tem.write_def('-', '-');
	WMM_findfundmode(wg, QTM, SYM, 0.0, 0.0, par, '-', '-', ma);
	tmm = ma(0);
	tmm.write_def('-', '-');

	// second run should start here
//	tem.read_def(QTE, SYM, '-', '-');
//	tmm.read_def(QTM, SYM, '-', '-');

	// show the mode profiles
	Rect display(-Wgph, -Wgpw, Wgph*1.5, Wgpw);
	tem.mfile(EY, MOD, display, 75, 115, '0', '0', 'C');
	tmm.mfile(HY, MOD, display, 75, 115, '0', '0', 'C');

	// the core region
	Rect core(0.0, -Wgpw/2.0, Wgph, Wgpw/2.0);

fprintf(stderr, "\nAbsorbing core: (attenuation 0.001/mikrom, 10/cm)\n"); 
	// permittivity perturbation
	p = attenuation(0.001, Wgpnf, Wgpl, core);
	p.write(stderr);
	// TE:imaginary phaseshift
	db = tem.phaseshift(p);
	fprintf(stderr, "TE: phaseshift: %g+i%g /mikrom\n", db.re, db.im);
	fprintf(stderr, "    mode intensity attenuation: %g /cm\n", 
	                -2.0*db.im*10000.0);
	// make a plot ...
	Cvector pert(1);
	Cvector amp(1);
	ma.clear();
	ma.add(tem);
	pert(0) = db;
	amp(0) = CC1;
	ma.prop(amp, pert, Wgph/2.0, display.y0, display.y1, 75, 
	        0.0, 2000.0, 200, 't', 'e');
	// ... and a curve
	ma.writeiopower(tem, tem, pert, 500, 0.0, 2000.0, '-', '-');

	// the same procedure for the TM mode 
	db = tmm.phaseshift(p);
	fprintf(stderr, "TM: phaseshift: %g+i%g /mikrom\n", db.re, db.im);
	fprintf(stderr, "    mode intensity attenuation: %g /cm\n", 
	                -2.0*db.im*10000.0);
	ma.clear();
	ma.add(tmm);
	pert(0) = db;
	amp(0) = CC1;
	ma.prop(amp, pert, Wgph/2.0, display.y0, display.y1, 75, 
	        0.0, 2000.0, 200, 't', 'm');
	ma.writeiopower(tmm, tmm, pert, 500, 0.0, 2000.0, '-', '-');

fprintf(stderr, "\nMagnetooptic waveguide, double layer TM phase shifter:\n"); 
	// the two magnetooptic layers
        Rect bottoml(0.0, -Wgpw/2.0, 0.2, Wgpw/2.0);
        Rect topl(0.2, -Wgpw/2.0, Wgph, Wgpw/2.0);
	// permittivity perturbation
	p1 = magopt_equat(-3000.0, Wgpnf, Wgpl, bottoml);
	p1.write(stderr);
	p2 = magopt_equat(3000.0, Wgpnf, Wgpl, topl);
	p2.write(stderr);
	// TM phaseshift 
	db1 = tmm.phaseshift(p1);
	db2 = tmm.phaseshift(p2);
	db = db1+db2;
	fprintf(stderr, 
                "TM           phaseshift: %g+i%g /mikrom\n", 
                db.re, db.im);
	fprintf(stderr, 
                "   Nonreciprocal effect: %g /cm\n", 
                2.0*db.re*10000.0);

fprintf(stderr, "\nMagnetooptic waveguide, two domain TE phase shifter:\n"); 
	// the two magnetooptic domains
        Rect leftd(0.0, -Wgpw/2.0, Wgph, 0.0);
        Rect rightd(0.0, 0.0, Wgph, Wgpw/2.0);
	// permittivity perturbation
	p1 = magopt_polar(-3000.0, Wgpnf, Wgpl, leftd);
	p1.write(stderr);
	p2 = magopt_polar(3000.0, Wgpnf, Wgpl, rightd);
	p2.write(stderr);
	// TE phaseshift 
	db1 = tem.phaseshift(p1);
	db2 = tem.phaseshift(p2);
	db = db1+db2;
	fprintf(stderr, 
                "TE           phaseshift: %g+i%g /mikrom\n", 
                db.re, db.im);
	fprintf(stderr, 
                "   Nonreciprocal effect: %g /cm\n", 
                2.0*db.re*10000.0);

fprintf(stderr, "\nGeometry perturbations:\n"); 
	fprintf(stderr, "Gradients with respect to the strip height:\n"); 
	grad = tem.horgeovar(1, 1);
	fprintf(stderr, "TE: partial beta/partial h = %g mikrom^(-2)\n", grad);
	grad = tmm.horgeovar(1, 1);
	fprintf(stderr, "TM: partial beta/partial h = %g mikrom^(-2)\n", grad);

	fprintf(stderr, "Gradients with respect to the strip width:\n"); 
	grad = tem.vergeovar(1, 1);
	fprintf(stderr, "TE: partial beta/partial w = %g mikrom^(-2)\n", grad);
	grad = tmm.vergeovar(1, 1);
	fprintf(stderr, "TM: partial beta/partial w = %g mikrom^(-2)\n", grad);

	return 0;
}