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

/*
 * gengwed.cpp
 * Electrodynamic constants, basic relations, general definitions
 * concerning 2D optical wave propagation and eigenmode analysis
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"structure.h"
#include"gengwed.h"
#include"cylfunc.h"

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

/* passive cavities: translation between characterizing quantities */

/* Q factor, given the complex eigenfrequency */
double theQ(Complex omega) { return 0.5*omega.re/omega.im; };

/* cavity lifetime [fs], given the imag. part of the compl. freq. [rad/fs] */
double tauC(double alpha) { return 0.5/alpha; };

/* outgoing radiation, power vs. frequency, FWHM [rad/fs] */
double FWHM_omega(double alpha) { return 2.0*alpha; }

/* outgoing radiation, power vs. wavelength, FWHM [mum] */
double FWHM_lambda(Complex omega) { return val_lambda(omega.re)*2.0*omega.im/omega.re; }

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

/* Fcomp -> field character 'E', 'H', 'S' */
char fldchr(Fcomp cp)
{
	switch(cp)
	{
		case EX: case EY: case ER: case EZ: case EP: case ET: return 'E';
		case HX: case HY: case HR: case HZ: case HP: case HT: return 'H';
		case SX: case SY: case SR: case SZ: case SP: case ST: return 'S';
	}
	fprintf(stderr, "ERR: fldchr \n"); 
	return '?';
}

/* Fcomp -> comp character 'x', 'y', 'z', 'r', 'p' */
char cpchr(Fcomp cp)
{
	switch(cp)
	{
		case EX: case HX: case SX: return 'x';
		case EY: case HY: case SY: return 'y';
		case ER: case HR: case SR: return 'r';
		case EZ: case HZ: case SZ: return 'z';
		case EP: case HP: case SP: return 'p';
		case ET: case HT: case ST: return 't';
	}
	fprintf(stderr, "ERR: cpchr \n"); 
	return '?';
}

/* 'EHS''xyrzp' -> Fcomp */
Fcomp fcomp(char eh, char cp)
{
	switch(eh)
	{
		case 'E': 
			switch(cp)
			{
			case 'x': return EX;
			case 'y': return EY;
			case 'r': return ER;
			case 'z': return EZ;
			case 'p': return EP;
			case 't': return ET;
			}
		case 'H': 
			switch(cp)
			{
			case 'x': return HX;
			case 'y': return HY;
			case 'r': return HR;
			case 'z': return HZ;
			case 'p': return HP;
			case 't': return HT;
			}
		case 'S':
			switch(cp)
			{
			case 'x': return SX;
			case 'y': return SY;
			case 'r': return SR;
			case 'z': return SZ;
			case 'p': return SP;
			case 't': return ST;
			}
	}
	fprintf(stderr, "ERR: fcomp \n"); 
	return EY;
}

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

/* FSid ->  character 'E', 'H', 'W' */
char idchr(FSid id)
{
	switch(id)
	{
		case mE: case qE: return 'E';
		case mH: case qH: return 'H';
		case mW: return 'W';
	}
	fprintf(stderr, "ERR: idchr \n"); 
	return 'W';
}

/* FSid ->  character 'm', 'q' */
char mqchr(FSid id)
{
	switch(id)
	{
		case mE: case mH: case mW: return 'm';
		case qE: case qH:          return 'q';
	}
	fprintf(stderr, "ERR: idchr \n"); 
	return 'm';
}

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

/* Polarization -> 'e', 'm' */
char polchr(Polarization p)
{
	if(p == TM) return 'm';
	else return 'e';
}
/* Polarization -> 'E', 'M' */
char polCHR(Polarization p)
{
	if(p == TM) return 'M';
	else return 'E';
}

/* Polarization -> principal field component */
Fcomp principalcomp(Polarization p)
{
	if(p == TM) return HY;
	else return EY;
}

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

/* convert complex field values f to real numbers, 
   depending on the attribute foa */
double realize(Complex f, Afo foa)
{
 	switch(foa)
	{
	case MOD:
		return f.abs();
		break;
	case SQR:
		return f.sqabs();
		break;
	case REP: case ORG:
		return f.re;
		break;
	case IMP:
		return f.im;
		break;       
	}
	return 1.0;
}

Dvector realize(Cvector f, Afo foa)
{
 	switch(foa)
	{
	case MOD:
		return f.abs();
		break;
	case SQR:
		return f.sqabs();
		break;
	case REP: case ORG:
		return f.re();
		break;
	case IMP:
		return f.im();
		break;       
	}
	return f.abs();
}

Dmatrix realize(Cmatrix f, Afo foa)
{
 	switch(foa)
	{
	case MOD:
		return f.abs();
		break;
	case SQR:
		return f.sqabs();
		break;
	case REP: case ORG:
		return f.re();
		break;
	case IMP:
		return f.im();
		break;       
	}
	return f.abs();
}

/* Afo -> 'm', 'o', 's', 'r', 'i' */
char afochr(Afo foa)
{
	switch(foa)
        {
                case MOD: return 'm'; break;
                case ORG: return 'f'; break;
                case SQR: return 'q'; break;
                case REP: return 'r'; break;
                case IMP: return 'i'; break;
        }                          
	return '*';
}

/* Afo + Fcomp -> '|Hx|', 'Re Ey', 'Im Hz' */
int afocpdesc(Afo foa, Fcomp cp, char *desc)
{
	int i=0;
	switch(foa)
        {
                case MOD: case SQR: 
			desc[i++] = '|'; 
			break;
                case REP: 
			desc[i++] = 'R'; desc[i++] = 'e'; desc[i++] = ' '; 
			break;
                case IMP: 
			desc[i++] = 'I'; desc[i++] = 'm'; desc[i++] = ' '; 
			break;
                case ORG: 	
			break;
        }                          
	desc[i++] = fldchr(cp);
	desc[i++] = cpchr(cp);
	switch(foa)
        {
                case MOD: 
			desc[i++] = '|'; 
			break;
                case SQR: 
			desc[i++] = '|'; desc[i++] = '^'; desc[i++] = '2'; 
			break;
                case REP: case IMP: case ORG: 	
			break;
        }                          
	desc[i] = 0;
	return i;
}

/* FSid -> '|E|', '|E|^2', 'W', etc.; filled into desc */
int iddesc(FSid id, char *desc)
{
	int i=0;
	switch(id)
        {
                case mE: 
			desc[i++] = '|'; desc[i++] = 'E'; desc[i++] = '|'; 
			break;
                case mH: 
			desc[i++] = '|'; desc[i++] = 'H'; desc[i++] = '|'; 
			break;
                case qE: 
			desc[i++] = '|'; desc[i++] = 'E'; desc[i++] = '|'; 
			desc[i++] = '^'; desc[i++] = '2'; 
			break;
                case qH: 
			desc[i++] = '|'; desc[i++] = 'H'; desc[i++] = '|'; 
			desc[i++] = '^'; desc[i++] = '2'; 
			break;
                case mW: 
			desc[i++] = 'W';
			break;
        }                          
	desc[i] = 0;
	return i;
}

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

/* permittivity perturbations, planar waveguides, 1D cross sections */

/* initialize */
Perturbation::Perturbation()
{
	pr.x0 = pr.x1 = 0.0;
	e = Cmatrix(3, 3);
	e.init(CC0);
}

/* initialize */
Perturbation::Perturbation(Interval iv)
{
	pr = iv;
	e = Cmatrix(3, 3);
	e.init(CC0);
}

/* initialize */
Perturbation::Perturbation(Interval iv, double rxxv, double ryyv, double rzzv, 
			   double rxyv, double ixzv, double iyzv, double iabsv)
{
	pr = iv;
	e = Cmatrix(3, 3);
	e(0, 0) = Complex(rxxv, iabsv);
	e(0, 1) = Complex(rxyv, 0.0);
	e(0, 2) = Complex(0.0,  ixzv);
	e(1, 0) = Complex(rxyv, 0.0);
	e(1, 1) = Complex(ryyv, iabsv);
	e(1, 2) = Complex(0.0,  iyzv);
	e(2, 0) = Complex(0.0, -ixzv);
	e(2, 1) = Complex(0.0, -iyzv);
	e(2, 2) = Complex(rzzv, iabsv);
}

/* initialize */
Perturbation::Perturbation(Interval iv, 
	                   Complex p_xx, Complex p_xy, Complex p_xz,	
	                   Complex p_yx, Complex p_yy, Complex p_yz,	
	                   Complex p_zx, Complex p_zy, Complex p_zz)	
{
	pr = iv;
	e = Cmatrix(3, 3);
	e(0, 0) = p_xx;
	e(0, 1) = p_xy;
	e(0, 2) = p_xz;
	e(1, 0) = p_yx;
	e(1, 1) = p_yy;
	e(1, 2) = p_yz;
	e(2, 0) = p_zx;
	e(2, 1) = p_zy;
	e(2, 2) = p_zz;
}

/* output to file dat */
void Perturbation::write(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, " perturbed region: "); 
		pr.write(dat);
		fprintf(dat, " permittivity perturbation:\n"); 
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(0,0).re,e(0,0).im,e(0,1).re,e(0,1).im,e(0,2).re,e(0,2).im);
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(1,0).re,e(1,0).im,e(1,1).re,e(1,1).im,e(1,2).re,e(1,2).im);
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(2,0).re,e(2,0).im,e(2,1).re,e(2,1).im,e(2,2).re,e(2,2).im);
	}
	else
	{
		comment("perturbation:", dat);
		pr.write(dat);
		e.write(dat);
	}
	return;
}
		
/* input from file dat */
void Perturbation::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) return;
	else
	{
		pr.read(dat);
		e.read(dat);
	}
	return;
}

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

/* special perturbations */

/* a change of the local refractive index:  
   r: the perturbed region on the waveguide cross section 
   n: uniform refractive index of the material in r
   deltan: the shift of the refractive index 
   The perturbation leads to an isotropic real first order permittivity shift 
   of 2*n*deltan */
Perturbation refindshift(double n, double deltan, Interval r)
{
	Perturbation p(r); 
	double deps;
	deps = 2.0*n*deltan;
	p.e.init(CC0);
	p.e(0, 0) = Complex(deps, 0.0);
	p.e(1, 1) = Complex(deps, 0.0);
	p.e(2, 2) = Complex(deps, 0.0);
	return p;
}

/* isotropic absorbtion: 
   r: the lossy region on the waveguide cross section 
   alpha: intensity attenuation [1/micrometer] of the material in r,
          the material damps the intensity of a plane wave according to
          exp(-alpha z), where z is the propagation distance
   n: refractive index of the material in r
   lambda: vacuum wavelength [micrometer]
   The attenuation leads to a complex refractive index n - i alpha lambda/4/pi,
   or a permittivity perturbation of - i n alpha lambda / (2 pi) */
Perturbation attenuation(double alpha, double n, double lambda, Interval r)
{
	Perturbation p(r); 
	double deps;
	deps = - n*alpha*lambda/2.0/PI;
	p.e.init(CC0);
	p.e(0, 0) = Complex(0.0, deps);
	p.e(1, 1) = Complex(0.0, deps);
	p.e(2, 2) = Complex(0.0, deps);
	return p;
}

/* magnetooptic permittivity contribution:
   r: the magnetooptic region on the waveguide cross section 
   sFrot: specific Faraday rotation of the material [degrees/cm] (!)
   theta, phi: orientation of the magnetization in spherical coordinates,
               vec(M) = M(sin(theta)cos(phi), sin(theta)sin(phi), cos(theta)) 
   n: refractive index of the material in r
   lambda: vacuum wavelength [micrometer] */
Perturbation magopt(double sFrot, double theta, double phi, 
                    double n, double lambda, Interval r)
{
	Perturbation p(r); 
	double xi, w;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	w = sin(theta)*cos(phi);
	p.e(1, 2) = Complex(0.0, w*xi);
	p.e(2, 1) = Complex(0.0,-w*xi);
	w = sin(theta)*sin(phi);
	p.e(0, 2) = Complex(0.0,-w*xi);
	p.e(2, 0) = Complex(0.0, w*xi);
	w = cos(theta);
	p.e(0, 1) = Complex(0.0, w*xi);
	p.e(1, 0) = Complex(0.0,-w*xi);
	return p;
} 
/* as above, for polar orientation of the magnetization, vec(M) || x-axis, */
Perturbation magopt_polar(double sFrot, double n, double lambda, Interval r)
{
	Perturbation p(r); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(1, 2) = Complex(0.0, xi);
	p.e(2, 1) = Complex(0.0,-xi);
	return p;
}
/* for equatorial orientation of the magnetization, vec(M) || y-axis, */
Perturbation magopt_equat(double sFrot, double n, double lambda, Interval r)
{
	Perturbation p(r); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(0, 2) = Complex(0.0,-xi);
	p.e(2, 0) = Complex(0.0, xi);
	return p;
}
/* and for longitudinal orientation of the magnetization, vec(M) || z-axis */
Perturbation magopt_long(double sFrot, double n, double lambda, Interval r)
{
	Perturbation p(r); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(0, 1) = Complex(0.0, xi);
	p.e(1, 0) = Complex(0.0,-xi);
	return p;
}

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

/* permittivity perturbations, channel waveguides, 2D cross sections */

/* initialize */
ChWgPert::ChWgPert()
{
	r.x0 = r.y0 = r.x1 = r.y1 = 0.0;
	e = Cmatrix(3, 3);
	e.init(CC0);
}

/* initialize */
ChWgPert::ChWgPert(Rect rv)
{
	r = rv;
	e = Cmatrix(3, 3);
	e.init(CC0);
}

/* initialize */
ChWgPert::ChWgPert(Rect ri, double rxxv, double ryyv, double rzzv, 
	 		    double rxyv, double ixzv, double iyzv, double iabsv)
{
	r = ri;
	e = Cmatrix(3, 3);
	e(0, 0) = Complex(rxxv, iabsv);
	e(0, 1) = Complex(rxyv, 0.0);
	e(0, 2) = Complex(0.0,  ixzv);
	e(1, 0) = Complex(rxyv, 0.0);
	e(1, 1) = Complex(ryyv, iabsv);
	e(1, 2) = Complex(0.0,  iyzv);
	e(2, 0) = Complex(0.0, -ixzv);
	e(2, 1) = Complex(0.0, -iyzv);
	e(2, 2) = Complex(rzzv, iabsv);
}

/* initialize */
ChWgPert::ChWgPert(Rect ri, 
	           Complex p_xx, Complex p_xy, Complex p_xz,	
	           Complex p_yx, Complex p_yy, Complex p_yz,	
	           Complex p_zx, Complex p_zy, Complex p_zz)	
{
	r = ri;
	e = Cmatrix(3, 3);
	e(0, 0) = p_xx;
	e(0, 1) = p_xy;
	e(0, 2) = p_xz;
	e(1, 0) = p_yx;
	e(1, 1) = p_yy;
	e(1, 2) = p_yz;
	e(2, 0) = p_zx;
	e(2, 1) = p_zy;
	e(2, 2) = p_zz;
}


/* output to file dat */
void ChWgPert::write(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		fprintf(dat, "perturbed region: "); 
		r.write(dat);
		fprintf(dat, "\n"); 
		fprintf(dat, "permittivity perturbation:\n"); 
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(0,0).re,e(0,0).im,e(0,1).re,e(0,1).im,e(0,2).re,e(0,2).im);
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(1,0).re,e(1,0).im,e(1,1).re,e(1,1).im,e(1,2).re,e(1,2).im);
		fprintf(dat, " %g+i%g   %g+i%g   %g+i%g\n",
		  e(2,0).re,e(2,0).im,e(2,1).re,e(2,1).im,e(2,2).re,e(2,2).im);
	}
	else
	{
		comment("perturbation:", dat);
		r.write(dat);
		e.write(dat);
	}
	return;
}
		
/* input from file dat */
void ChWgPert::read(FILE *dat)
{
	if(dat == stderr || dat == stdout)
	{
		r.x0 = r.y0 = r.x1 = r.y1 = 0.0;
		e = Cmatrix(3, 3);
		e.init(CC0);
	}
	else
	{
		r.read(dat);
		e.read(dat);
	}
	return;
}

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

/* special perturbations */

/* isotropic absorbtion: 
   rl: the lossy area on the waveguide cross section 
   alpha: intensity attenuation of the material in rl, [1/micrometer]
          the material damps the intensity of a plane wave according to
          exp(-alpha z)
   n: refractive index of the material in rl
   lambda: vacuum wavelength [micrometer]
   The attenuation leads to a complex refractive index n - i alpha lambda/4/pi,
   or a permittivity perturbation of - i n alpha lambda / (2 pi) */
ChWgPert attenuation(double alpha, double n, double lambda, Rect rl)
{
	ChWgPert p(rl); 
	double deps;
	deps = - n*alpha*lambda/2.0/PI;
	p.e.init(CC0);
	p.e(0, 0) = Complex(0.0, deps);
	p.e(1, 1) = Complex(0.0, deps);
	p.e(2, 2) = Complex(0.0, deps);
	return p;
}

/* magnetooptic permittivity contribution:
   rl: the magnetooptic area on the waveguide cross section 
   sFrot: specific Faraday rotation of the material [degrees/cm] (!)
   theta, phi: orientation of the magnetization in spherical coordinates,
               vec(M) = M(sin(theta)cos(phi), sin(theta)sin(phi), cos(theta)) 
   n: refractive index of the material in rl
   lambda: vacuum wavelength [micrometer] */
ChWgPert magopt(double sFrot, double theta, double phi, 
                double n, double lambda, Rect rl)
{
	ChWgPert p(rl); 
	double xi, w;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	w = sin(theta)*cos(phi);
	p.e(1, 2) = Complex(0.0, w*xi);
	p.e(2, 1) = Complex(0.0,-w*xi);
	w = sin(theta)*sin(phi);
	p.e(0, 2) = Complex(0.0,-w*xi);
	p.e(2, 0) = Complex(0.0, w*xi);
	w = cos(theta);
	p.e(0, 1) = Complex(0.0, w*xi);
	p.e(1, 0) = Complex(0.0,-w*xi);
	return p;
} 
/* as above, for polar orientation of the magnetization, vec(M) || x-axis, */
ChWgPert magopt_polar(double sFrot, double n, double lambda, Rect rl)
{
	ChWgPert p(rl); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(1, 2) = Complex(0.0, xi);
	p.e(2, 1) = Complex(0.0,-xi);
	return p;
}
/* for equatorial orientation of the magnetization, vec(M) || y-axis, */
ChWgPert magopt_equat(double sFrot, double n, double lambda, Rect rl)
{
	ChWgPert p(rl); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(0, 2) = Complex(0.0,-xi);
	p.e(2, 0) = Complex(0.0, xi);
	return p;
}
/* and for longitudinal orientation of the magnetization, vec(M) || z-axis */
ChWgPert magopt_long(double sFrot, double n, double lambda, Rect rl)
{
	ChWgPert p(rl); 
	double xi;
	xi = n*lambda/PI*sFrot*PI/180.0/10000.0;
	p.e.init(CC0);
	p.e(0, 1) = Complex(0.0, xi);
	p.e(1, 0) = Complex(0.0,-xi);
	return p;
}

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

/* convert a complex permittivity eps into a complex refractive index n
   such that eps = n*n, eps != 0 required; 
   eps = eps'-i eps'', n = n'-i n'', sign convention: n' > 0 */
Complex cepston(Complex e)
{
	double ed = e.re;;
	double edd = -e.im;
	double ea = sqrt(ed*ed+edd*edd);
	double nd, ndd;
	nd  = sqrt(0.5*(ea+ed));
	ndd = sqrt(0.5*(ea-ed));
	if(edd < 0.0) ndd = -ndd;
	return Complex(nd, -ndd);
}
 
/* ---------------------------------------------------------------------- */

/* Gaussian beams, input specification for the Helmholtz solvers,
   paraxial representation of moderately tilted beams, Gaussian profile
   prescribed on the x-axis, propagation in the positive z-direction
   with specific wavelength in a homogeneous dielectric medium */

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

/* initialize */
Gaussianbeam::Gaussianbeam(Polarization tp, double tw, double tx0, 
                           double ttheta, double tlambda, double tng)
{
	pol = tp;
	w = tw;
	x0 = tx0;
	theta = ttheta;
	lambda = tlambda;
	ng = tng;
	
	double k0 = val_k0(lambda);
	double f = 1.0;
	if(pol == TE) f = val_invommu0(lambda);
	else          f = val_invomep0(lambda)/(ng*ng);
	tef = CCI*val_invommu0(lambda);
	tmf = CCI*val_invomep0(lambda)/(ng*ng);
	if(theta<=-180.0 || theta>=180.0) gberror("initialize: theta");
	kncth = k0*ng*cos(theta*PI/180.0);
	knsth = k0*ng*sin(theta*PI/180.0);
	amp = 1.0/sqrt(f*sqrt(2.0*PI)/4.0*w*kncth);
}

/* field profile at the initial z-position */
Complex Gaussianbeam::field(Fcomp cp, double x) const
{
	Complex f;
	switch(pol) 
	{
	case TE:
		switch(cp)
		{
		case EY: return phi(x);
		case HX: f = tef*dzphi(x); f.neg(); return f;
		case HZ: f = tef*dxphi(x); return f;
		case EX: case EZ: case HY: return CC0;
		case SZ: f = tef*phi(x)*(dzphi(x)).conj();
		         return Complex(-0.5*f.re, 0.0);
		default: gberror("field: cp not implemented");
		}
		break;
	case TM:
		switch(cp)
		{
		case HY: return phi(x);
		case EX: f = tmf*dzphi(x); return f;
		case EZ: f = tmf*dxphi(x); f.neg(); return f;
		case HX: case HZ: case EY: return CC0;
		case SZ: f = tmf*(phi(x)).conj()*dzphi(x);
		         return Complex(0.5*f.re, 0.0);
		default: gberror("field: cp not implemented");
		}
		break;
	default: gberror("field: pol");
		break;
	}
	return CC0;
}

/* helper function: princical field component */
Complex Gaussianbeam::phi(double x) const
{
	double dx = x-x0;
	double dxdwq = dx/w;
	dxdwq *= dxdwq;
	return expi(-knsth*dx)*exp(-dxdwq)*amp;
}

/* helper function: princical field component, z-derivative */
Complex Gaussianbeam::dzphi(double x) const
{
	Complex f = CCmI*kncth+2.0/w/w*(x-x0)*knsth/kncth;
	return f*phi(x);
}

/* helper function: princical field component, x-derivative */
Complex Gaussianbeam::dxphi(double x) const
{
	Complex f = CCmI*knsth-2.0/w/w*(x-x0);
	return f*phi(x);
}

/*
 * helper functions describing a Gaussian beam propagating along the z-axis
 */

Complex gbm_phi(double kn, double w, double x, double z)
{
	Complex f;
	double zr = w*w*kn/2.0;
	double wi = w*sqrt(1.0+z*z/zr/zr);
	double rr = z/(z*z+zr*zr);
	f = Complex(sqrt(w/wi)*exp(-x*x/wi/wi), 0.0);
	f *= expi(-(kn*(z+0.5*x*x*rr)-0.5*atan(z/zr)));
	return f;
}
Complex gbm_dxphi(double kn, double w, double x, double z)
{
	double t1 = w*w;
	double t5 = kn*kn;
	double t8 = t1*t1;
	double t12 = z*z;
	double t16 = x*x;
	double t24 = atan(2.0*z/t1/kn);
	double t26 = t8*t5;
	double t35 = t26+4.0*t12;
	double t45 = sqrt(sqrt(t35));
	Complex t39 = exp(-((2.0*CCI)*t5*kn*z*t8+(8.0*CCI)*kn*t12*z+(4.0*CCI)*kn*z*t16+(CCmI)*t24*t26+(4.0*CCmI)*t24*t12+2.0*t16*t1*t5)/t35/2.0);
	double t42 = sqrt(kn);
	Complex f = -2.0*(t1*kn+(2.0*CCI)*z)*t39*x*t42*kn*w/t45/t35;
	return f;
}
Complex gbm_dzphi(double kn, double w, double x, double z)
{
	double t1 = kn*kn;
	double t2 = t1*t1;
	double t4 = w*w;
	double t5 = t4*t4;
	double t6 = t5*t5;
	double t8 = t1*kn;
	Complex t12 = (2.0*CCI)*t8;
	double t13 = x*x;
	double t17 = z*z;
	double t20 = z*t5;
	double t23 = t4*t13;
	double t31 = t17*t17;
	double t36 = t17*z;
	double t50 = atan(2.0*z/t4/kn);
	double t52 = t5*t1;
	double t60 = t52+4.0*t17;
	Complex t64 = exp(-(t12*t20+8.0*CCI*kn*t36+4.0*CCI*kn*z*t13+CCmI*t50*t52+4.0*CCmI*t50*t17+2.0*t23*t1)/t60/2.0);
	double t66 = sqrt(kn);
	double t69 = t60*t60;
	double t70 = sqrt(sqrt(t60));
	Complex f = -(CCI*t2*kn*t6+CCmI*t8*t5*t4+t12*t13*t5+8.0*CCI*t8*t5*t17+2.0*t20*t1-8.0*t23*z*t1+4.0*CCmI*kn*t4*t17+16.0*CCI*kn*t31+8.0*CCmI*kn*t13*t17+8.0*t36)*t64*t66*w/t70/t69;
	return f;
}
Complex gbm_dxxphi(double kn, double w, double x, double z)
{
	double t1 = w*w;
	double t2 = t1*t1;
	double t4 = kn*kn;
	double t5 = t4*kn;
	double t7 = x*x;
	double t12 = t2*t4;
	double t15 = t7*z;
	double t19 = z*z;
	double t25 = t19*z;
	double t40 = atan(2.0*z/t1/kn);
	double t50 = t12+4.0*t19;
	Complex t54 = exp(-((2.0*CCI)*t5*z*t2+(8.0*CCI)*kn*t25+(4.0*CCI)*kn*t15+(-CCI)*t40*t12+(-4.0*CCI)*t40*t19+2.0*t7*t1*t4)/t50/2.0);
	double t56 = sqrt(kn);
	double t60 = t50*t50;
	double t61 = sqrt(sqrt(t50));
	Complex f = -2.0*(t2*t1*t5-2.0*t5*t7*t2+(2.0*CCI)*z*t12+(-8.0*CCI)*t1*t15*t4+4.0*t1*kn*t19+8.0*kn*t7*t19+(8.0*CCI)*t25)*t54*t56*kn*w/t61/t60;
	return f;
}
Complex gbm_dzzphi(double kn, double w, double x, double z)
{
	double t1 = kn*kn;
	double t2 = t1*t1;
	double t3 = t2*t2;
	double t5 = w*w;
	double t6 = t5*t5;
	double t7 = t6*t6;
	double t8 = t7*t7;
	double t11 = t1*kn;
	double t12 = x*x;
	double t13 = t11*t12;
	double t14 = z*z;
	double t15 = t14*t14;
	double t16 = t15*z;
	double t19 = t7*t6;
	double t20 = t2*t1;
	double t23 = t6*t5;
	double t28 = t2*kn;
	double t29 = t5*t7;
	double t32 = t15*t15;
	double t38 = t6*t1;
	double t45 = t6*t2;
	double t46 = t12*t12;
	double t50 = t1*t12;
	double t54 = t12*t20;
	double t58 = t23*t2;
	double t66 = t2*t11;
	double t69 = t14*z;
	Complex t82 = t3*t1*t8+96.0*CCI*t69*t11*t23+3.0*t19*t20-2.0*t3*t7*t23+4.0*CCmI*t19*z*t66+256.0*t1*t32+12.0*t7*t2*t14-48.0*t38*t15-64.0*t2*t12*t6*t15-96.0*t45*t46*t14+576.0*t50*t15*t5+16.0*t54*t7*t14+96.0*t58*t12*t14+192.0*CCI*t16*kn*t5+192.0*CCmI*t69*t13*t6+12.0*CCI*z*t28*t29+384.0*CCI*t16*kn*t12+256.0*CCI*t5*t13*t16+48.0*CCmI*t7*t69*t28;
	double t87 = t28*t12;
	double t121 = t15*t14;
	Complex t139 = 16.0*CCI*t29*t66*t12*z+128.0*CCI*t23*t87*t69+128.0*CCmI*t5*t11*t46*t69+72.0*CCmI*z*t87*t7+192.0*CCmI*t6*t16*t11+256.0*CCmI*t15*t69*kn+4.0*t3*t12*t19-96.0*t58*t15+4.0*t20*t46*t7-12.0*t54*t29+16.0*t3*t19*t14+96.0*t20*t7*t15+256.0*t45*t121-24.0*t20*t29*t14-256.0*t50*t121-128.0*t1*t5*t121+64.0*t1*t46*t15+32.0*CCI*t23*t28*t46*z-192.0*t121;
	double t154 = atan(2.0*z/t5/kn);
	double t164 = t38+4.0*t14;
	Complex t168 = exp(-(2.0*CCI*t11*z*t6+8.0*CCI*kn*t69+4.0*CCI*kn*z*t12+CCmI*t154*t38+4.0*CCmI*t154*t14+2.0*t12*t5*t1)/t164/2.0);
	double t170 = sqrt(kn);
	double t173 = t164*t164;
	double t174 = t173*t173;
	double t175 = sqrt(sqrt(t164));
	Complex f = -(t82+t139)*t168*t170*w/t175/t174;
	return f;
}

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

/* a signed square root */
double sgn_sqrt(double x)
{
	if(x >= 0) return sqrt(x);
	return -sqrt(-x);
}

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

/*
 * for computing parameter scans of resonant structures, narrow peaks in an
 * otherwise regular background:
 * given a range r of arguments and an initial stepsize h0, determine the
 * next value of a series or arguments x (not ordered!, x is being extended
 * by one value), such that initially appearing minima and maxima in the curve
 * x, y are resolved with a minimum stepsize of hm upon completion
 * return value != 0: finished; 
 * ext_type: -1: resolve minima in y(x),
 * ext_type:  1: resolve maxima in y(x), default: resolve both;
 * minydev: minimum absolute deviation in y for an extremum to be identified
 * rough, inefficient, ...
 */
int next(Interval r, double h0, double hm, Dvector &x, Dvector y, int ext_type,
         double minydev)
{
	if(r.x1-r.x0 < h0) return 1;
	if(x.nel == 0) 
	{
		x.append(r.x0);
		return 0;
	}
	double xm;
	xm = x.max();
	if(xm <= r.x1-h0) 
	{
		x.append(xm+h0);
		return 0;
	}
	if(x.nel != y.nel) return 1;

	Dvector xc, yc;
	xc = x; yc = y; sort(xc, yc);

	if(xc.nel <= 2) return 1;
	int p=0;
	while(p+2 <= xc.nel-1)
	{
		double x0, x1, x2;
		double y0, y1, y2;
		x0 = xc(p);
		x1 = xc(p+1);
		x2 = xc(p+2);
		if(x1-x0>hm || x2-x1>hm)
		{
			y0 = yc(p);
			y1 = yc(p+1);
			y2 = yc(p+2);
			if(y0<y1-minydev && y1-minydev>y2 && ext_type != -1) 
			{
				if(x1-x0 >= x2-x1)	
				{
					x.append(0.5*(x0+x1));
					return 0;
				}
				else
				{
					x.append(0.5*(x1+x2));
					return 0;
				}
			}
			if(y0>y1+minydev && y1+minydev<y2 && ext_type != 1) 
			{
				if(x1-x0 >= x2-x1)	
				{
					x.append(0.5*(x0+x1));
					return 0;
				}
				else
				{
					x.append(0.5*(x1+x2));
					return 0;
				}
			}
		}
		++p;
	}
	return 1;
}
int next(Interval r, double h0, double hm, Dvector &x, Dvector y, int ext_type)
{
	return next(r, h0, hm, x, y, ext_type, 0.0);
} 
int next(Interval r, double h0, double hm, Dvector &x, Dvector y)
{
	return next(r, h0, hm, x, y, 0);
}