/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * http://metric.computational-photonics.eu/ 
 */

/*
 * slaref.cpp
 * Reflection of plane waves from a dielectric multilayer stack
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"structure.h"
#include"gengwed.h"
#include"slamode.h"
#include"slamarr.h"
#include"slaref.h"

/* error message */
void slabreferror(const char *s)
{
	fprintf(stderr, "\nslaref: %s.\n", s);
	exit(1);
}

/* reflection of a plane wave from a dielectric multilayer stack:
   compute reflection + transmission coefficients, assemble the stationary
   electromagnetic field:
   wg: the multilayer structure, interface normal along the x-axis,
       needs at least 1 inner layer
   pol: input polarization, use TM for parallel polarization, 
                                TE for orthogonal polarization
   theta: incidence angle, in radiants
   in: intensity of the input wave, with respect to the surface (y-z-plane)
   ref (output): the reflected intensity, reflectivity for in=1
   trans (output): the transmitted intensity, transmittivity for in=1 
   fld (output): the electromagnetic field, basic wave profiles 
   amp (output): their amplitudes 
   rc, tc (output): complex relative amplitude coefficients 
                    of reflected, and transmitted waves */
void mlref(Waveguide wg, Polarization pol, double theta, 
           double in, 
           double& ref, Complex& rc,
           double& trans, Complex& tc, 
           ModeArray& fld, Cvector& amp)
{ 
	double beta, k0;
	int num;
	Complex df;
	Cvector dum(2);
	double err;

	num = wg.nx;
	k0 = val_k0(wg.lambda);
	beta = k0*wg.n(0)*sin(theta); 

	Cvector u(num+2);
	Cvector r(num+2);
	Dvector xu(num+2);
	Dvector xr(num+2);
	Dvector gamma(num+2);
	Cvector cg(num+2);
	Ivector s(num+2);
	Dvector eta(num+2);
	
	Cvector fca(num+2);
	Cvector fcb(num+2);
	Cvector fcc(num+2);
	Cvector fcd(num+2);
	Cvector fcA(num+2);
	Cvector fcB(num+2);
	Cvector fcC(num+2);
	Cvector fcD(num+2);
	Cvector fcR(num+2);

	for(int j=0; j<=num+1; ++j)
	{
		double g = k0*k0*wg.eps(j)-beta*beta;
		gamma(j)=sqrt(fabs(g));
		if(g >= 0.0)
		{
			s(j) = 1;
			cg(j) = Complex(gamma(j), 0.0);
			if(j==0) xu(j) = xr(j) = wg.hx(0);
			else
			{
				if(j==num+1) xu(j)=xr(j)=wg.hx(num);
				else xr(j)=xu(j)=0.5*(wg.hx(j-1)+wg.hx(j)); 
			}
		}
		else
		{
			s(j) = -1;
			cg(j) = Complex(0.0, -gamma(j));
			if(j >= 1  ) xu(j)=wg.hx(j-1); else xu(j)=wg.hx(0);
			if(j <= num) xr(j)=wg.hx(j);   else xr(j)=wg.hx(num);
		}
		if(pol == TE) eta(j)=1.0;
		else          eta(j)=1.0/wg.eps(j);
	}	
	for(int j=0; j<=num; ++j)
	{
		double dxu, dxr;
		Complex egu, egr;
		dxu = wg.hx(j)-xu(j);
		dxr = wg.hx(j)-xr(j);
		egu = exp(CCI*cg(j)*(-dxu));
		egr = exp(CCI*cg(j)*( dxr));
		fca(j) = egu; 
		fcb(j) = egr; 
		fcc(j) = CCI*cg(j)*egu*(-eta(j));
		fcd(j) = CCI*cg(j)*egr*( eta(j));

		dxu = wg.hx(j)-xu(j+1);
		dxr = wg.hx(j)-xr(j+1);
		egu = exp(CCI*cg(j+1)*(-dxu));
		egr = exp(CCI*cg(j+1)*( dxr));
		fcA(j) = egu; 
		fcB(j) = egr; 
		fcC(j) = CCI*cg(j+1)*egu*(-eta(j+1));
		fcD(j) = CCI*cg(j+1)*egr*( eta(j+1));
	}	

	fcR(num) = (fcc(num)/fcC(num)-fca(num)/fcA(num))/
	           (fcb(num)/fcA(num)-fcd(num)/fcC(num));
	for(int j=num-1; j>=0; --j)
	{	
		fcR(j) = ( fcc(j)/(fcC(j)+fcD(j)*fcR(j+1))
                          -fca(j)/(fcA(j)+fcB(j)*fcR(j+1)))/
		         ( fcb(j)/(fcA(j)+fcB(j)*fcR(j+1))
                          -fcd(j)/(fcC(j)+fcD(j)*fcR(j+1)));
	}

	if(pol == TE) u(0) = Complex(sqrt(2.0*fabs(in)
	                                  /gamma(0)/val_invommu0(wg.lambda)));
	else          u(0) = Complex(sqrt(2.0*fabs(in)*wg.eps(0)
	                                  /gamma(0)/val_invomep0(wg.lambda)));
	r(0) = fcR(0)*u(0);
	for(int j=1; j<=num; ++j)
	{
		u(j) = (fca(j-1)*u(j-1)+fcb(j-1)*r(j-1))/
		       (fcA(j-1)+fcB(j-1)*fcR(j));
		r(j) = fcR(j)*u(j);
	}
	u(num+1) = (fca(num)*u(num)+fcb(num)*r(num))/fcA(num);
	r(num+1) = CC0;

	err = 0.0;
	for(int j=0; j<=num; ++j)
	{
	df = 
	(u(j)*exp(CCI*cg(j)*(-(wg.hx(j)-xu(j))))
	+r(j)*exp(CCI*cg(j)*( (wg.hx(j)-xr(j)))))-
	(u(j+1)*exp(CCI*cg(j+1)*(-(wg.hx(j)-xu(j+1))))
	+r(j+1)*exp(CCI*cg(j+1)*( (wg.hx(j)-xr(j+1)))));
	err += df.abs();
	df = 
	(u(j)*exp(CCI*cg(j)*(-(wg.hx(j)-xu(j))))*CCI*cg(j)*(-eta(j))
	+r(j)*exp(CCI*cg(j)*( (wg.hx(j)-xr(j))))*CCI*cg(j)*( eta(j)))-
	(u(j+1)*exp(CCI*cg(j+1)*(-(wg.hx(j)-xu(j+1))))*CCI*cg(j+1)*(-eta(j+1))
	+r(j+1)*exp(CCI*cg(j+1)*( (wg.hx(j)-xr(j+1))))*CCI*cg(j+1)*( eta(j+1)));
	err += df.abs();
	}
//	if(err/(u(0).abs()) > 1.0e-8) slabreferror("Accuracy");
	if(err/(u(0).abs()) > 1.0e-8) fprintf(stderr, "\n(!) Acc: %g.\n", err/(u(0).abs()));

	Mode mfr, mfi;
	mfr.wg = wg;
	mfi.wg = wg;
	mfr.pol = pol;
	mfi.pol = pol;
	mfr.k0 = k0;
	mfi.k0 = k0;
	mfr.invommu0 = val_invommu0(wg.lambda);
	mfi.invommu0 = val_invommu0(wg.lambda);
	mfr.invomep0 = val_invomep0(wg.lambda);
	mfi.invomep0 = val_invomep0(wg.lambda);
	mfr.beta = beta;
	mfi.beta = beta;
	mfr.betaq = beta*beta;
	mfi.betaq = beta*beta;
	mfr.neff = beta/k0;
	mfi.neff = beta/k0;
	mfr.npcB = 0.0;
	mfi.npcB = 0.0;
	mfr.ppt = PROPAG;
	mfi.ppt = PROPAG;
	mfr.bdt_b = OPEN;
	mfi.bdt_b = OPEN;
	mfr.bp_b = -AWAY;
	mfi.bp_b = -AWAY;
	mfr.bdt_t = OPEN;
	mfi.bdt_t = OPEN;
	mfr.bp_t = AWAY;
	mfi.bp_t = AWAY;

	mfr.N = num+2;
	mfi.N = num+2;
	mfr.piece = new SmPiece[num+2];
	mfi.piece = new SmPiece[num+2];
	for(int l=0; l<=num+1; ++l)
	{
		mfr.piece[l].pol = pol;
		mfi.piece[l].pol = pol;
		Interval i = wg.layer(l);
		mfr.piece[l].bif = i.x0;
		mfi.piece[l].bif = i.x0;
		mfr.piece[l].tif = i.x1;
		mfi.piece[l].tif = i.x1;
		mfr.piece[l].k0 = k0;
		mfi.piece[l].k0 = k0;
		mfr.piece[l].k0nq = k0*k0*wg.eps(l);
		mfi.piece[l].k0nq = k0*k0*wg.eps(l);
		mfr.piece[l].inve = 1.0/wg.eps(l);
		mfi.piece[l].inve = 1.0/wg.eps(l);
		
		mfr.piece[l].g = gamma(l);
		mfi.piece[l].g = gamma(l);
		mfr.piece[l].xa = xr(l);
		mfi.piece[l].xa = xr(l);
		mfr.piece[l].xb = xu(l);
		mfi.piece[l].xb = xu(l);

		double gq = beta*beta-k0*k0*wg.eps(l);
		if(gq > 0.0) 
		{
			mfr.piece[l].bhv = HYP;
			mfi.piece[l].bhv = HYP;
		}
		else	
		{
			mfr.piece[l].bhv = OSC;
			mfi.piece[l].bhv = OSC;
		}
	}
	for(int j=0; j<=num+1; ++j)  
	{
		if(s(j) == 1)
		{
			mfr.piece[j].A = -(r(j)-u(j)).im;
			mfr.piece[j].B =  (r(j)+u(j)).re;
			mfi.piece[j].A =  (r(j)-u(j)).re;
			mfi.piece[j].B =  (r(j)+u(j)).im;
		}
		else
		{
			mfr.piece[j].A = r(j).re;
			mfr.piece[j].B = u(j).re;
			mfi.piece[j].A = r(j).im;
			mfi.piece[j].B = u(j).im;
		}
	}
	mfr.special = -1;
	mfi.special = -1;
	mfr.ids[0] = 'p'; mfr.ids[1] = 'w';
	mfr.ids[2] = 'R'; mfr.ids[3] = 'e';
	mfr.ids[4] = 0; mfr.ids[5] = 0; mfr.ids[6] = 0; mfr.ids[7] = 0;
	mfi.ids[0] = 'p'; mfi.ids[1] = 'w';
	mfi.ids[2] = 'I'; mfi.ids[3] = 'm';
	mfi.ids[4] = 0; mfi.ids[5] = 0; mfi.ids[6] = 0; mfi.ids[7] = 0;
	mfr.setfieldmap();
	mfi.setfieldmap();
	fld.clear();
	fld.add(mfr);
	fld.add(mfi);
	amp = Cvector(2);
	amp(0) = CC1;
	amp(1) = CCI;
	dum.init(CC0);

	df = fld.field(amp, dum, SX, wg.hx(0)-1.0, 0.0)
	    -fld.field(amp, dum, SX, wg.hx(num+1)+1.0, 0.0);      
	if(df.abs() > 1.0e-10) slabreferror("Field setup");

	if(pol == TE)
	{
		rc = r(0)/u(0);
		ref = 0.5*gamma(0)*val_invommu0(wg.lambda)*r(0).sqabs();
		if(s(num+1) > 0) 
		{
			tc = u(num+1)/u(0);
			trans = 0.5*gamma(num+1)*val_invommu0(wg.lambda)*u(num+1).sqabs();
		}
		else 
		{
			tc = CC0;
			trans = 0.0;
		}
	}
	else
	{
		rc = r(0)/u(0);
		ref = 0.5*gamma(0)*val_invomep0(wg.lambda)/wg.eps(0)*r(0).sqabs();
		if(s(num+1) > 0) 
		{
			tc = u(num+1)/u(0);
			trans = 0.5*gamma(num+1)/wg.eps(num+1)*val_invomep0(wg.lambda)*u(num+1).sqabs();
		}
		else 
		{
			tc = CC0;
			trans = 0.0;
		}
	}
	if(fabs((ref - rc.sqabs())/in) > 1.0e-10)
		slabreferror("Reflectance");

	if(s(num+1) > 0) 
	{
		double thetaT = asin(wg.n(0)/wg.n(num+1)*sin(theta));
		if(pol == TE)
		{
//			if(fabs((trans - wg.n(num+1)*cos(thetaT)/wg.n(0)/cos(theta)*tc.sqabs())/in) > 1.0e-10)
			if(fabs((trans - wg.n(num+1)*cos(thetaT)/wg.n(0)/cos(theta)*tc.sqabs())/in) > 1.0e-8)
				slabreferror("TE, transmittance");
		}
		else 
		{
			if(fabs((trans - wg.n(0)*cos(thetaT)/wg.n(num+1)/cos(theta)*tc.sqabs())/in) > 1.0e-10)
				slabreferror("TM, transmittance");
		}
	}

	if(fabs((in-ref-trans)/in) > 1.0e-10) 
	slabreferror("Energy conservation");
	return;
}
void mlref(Waveguide wg, Polarization pol, double theta, double in,
           double& ref, double& trans, ModeArray& fld, Cvector& amp)
{
	Complex rc, tc;
	mlref(wg, pol, theta, in, ref, rc, trans, tc, fld, amp);
	return;
}